The Perl Toolchain Summit needs more sponsors. If your company depends on Perl, please support this very important event.
our $VERSION = '0.27';
use PDL::Exporter;

pp_setversion($VERSION);

# Cargo culted from PDL::Opt::NonLinear
pp_addhdr('
pdl    *pdl1, *pdl2, *pdl3, *pdl4, *pdl5;
SV    *sv_pdl1, *sv_pdl2, *sv_pdl3, *sv_pdl4, *sv_pdl5;
#include <math.h>
#include <stdio.h>
#include <string.h>
/* Change names when fixing glibc-2.1 bug */
#ifdef MY_FIXY0
#define y0(a) fixy0(a)
extern double fixy0(double a);
#endif
#ifdef MY_FIXYN
#define yn(a,b) fixyn(a,b)
extern double fixyn(int a, double b);
#endif
');

## handle various cases of 'finite'
#
if ($^O =~ /MSWin/) {
# _finite in VC++ 4.0
pp_addhdr('
#define finite _finite
#include <float.h>
/* avoid annoying warnings */
typedef long int logical;
typedef long int integer;
typedef long int ftnlen;
#ifdef __cplusplus
typedef float (*paramf)(...);
typedef double (*paramd)(...);
typedef void (*paramv)(...);
#else
typedef float (*paramf)();
typedef double (*paramd)();
typedef void (*paramv)();
#endif
');
} else {
pp_addhdr('
/* avoid annoying warnings */
typedef int logical;
typedef int integer;
typedef int ftnlen;
#ifdef __cplusplus
typedef float (*paramf)(...);
typedef double (*paramd)(...);
typedef void (*paramv)(...);
#else
typedef float (*paramf)();
typedef double (*paramd)();
typedef void (*paramv)();
#endif
')
}

pp_addpm({At=>'Top'},<<'EOD');
use strict;
use warnings;
use PDL::Ufunc;
use PDL::Ops;
use PDL::NiceSlice;
use Carp;

# ABSTRACT: Quadratic programming solver for PDL

=head1 NAME

PDL::Opt::QP - Quadratic programming solver for PDL

=head1 SYNOPSIS

    use PDL;
    use PDL::NiceSlice;
    use PDL::Opt::QP;

    my $mu   = pdl(q[ 0.0427 0.0015 0.0285 ])->transpose; # [ n x 1 ]
    my $mu_0 = 0.0427;
    my $dmat = pdl q[ 0.0100 0.0018 0.0011 ;
                      0.0018 0.0109 0.0026 ;
                      0.0011 0.0026 0.0199 ];
    my $dvec = zeros(3);
    my $amat = $mu->glue( 0, ones( 1, 3 ) )->copy;
    my $bvec = pdl($mu_0)->glue( 1, ones(1) )->flat;
    my $meq  = pdl(2);

    my $sol = qp( $dmat, $dvec, $amat, $bvec, $meq );
    say "Solution: ", $sol->{x};

=head1 DESCRIPTION

This routine uses Goldfarb/Idnani algorithm to solve the
following minimization problem:

           minimize  f(x) = 0.5 * x' D x  -  d' x
              x

    optionally constrained by:

            Aeq'  x  = a_eq
            Aneq  x >= b_neq

=cut

EOD

# Fortran function signature:
#  subroutine qpgen2( ->dmat, ->dvec, fddmat, n, <-sol, <-lagr,
#      <-crval, ->amat, ->bvec, fdamat, q, ->meq, <-iact, <-nact,
#      <-iter, work, <->ierr)
#  integer n, i, j, l, l1,
# *     info, q, iact(*), iter(*), it1,
# *     ierr, nact, iwzv, iwrv, iwrm, iwsv, iwuv, nvl,
# *     r, fdamat, iwnbv, meq, fddmat
#  double precision dmat(fddmat,*), dvec(*), lagr(*), sol(*), bvec(*)
# $     ,work(*), temp, sum, t1, tt, gc, gs, crval,nu, amat(fdamat,*)
# $     , vsmall, tmpa, tmpb

pp_def("qpgen2",
    HandleBad => 0,
    Pars => 'dmat(m,m); dvec(m);
        [o]sol(m); [o]lagr(q); [o]crval();
        amat(m,q); bvec(q); int meq();
        int [o]iact(q); int [o]nact();
        int [o]iter(s=2); [t]work(z); int [o]ierr();
    ',
    RedoDimsCode => pp_line_numbers(__LINE__, q{
        /* Calculate z */
        int m_size = $PDL(dmat)->dims[0];
        int q_size = $PDL(amat)->dims[1];

        /* r = max(n, q) */
        int r_size = m_size < q_size ? q_size : m_size;

        $SIZE(z) = 2*m_size + r_size*(r_size+5)/2 + 2*q_size + 1;
    }),
    GenericTypes => [D],
    Code => pp_line_numbers(__LINE__, q{
        extern int qpgen2_(
            double *dmat, double *dvec, integer *fddmat, integer *n,

            double *sol, double *lagr, double *crval, double *amat,
            double *bvec, integer *fdamat, integer *q, integer *meq,
            integer *iact, integer *nact, integer *iter, double *work,
            integer *ierr
        );

        int m_size = $SIZE(m);
        int q_size = $SIZE(q);
        int factor_and_output_err;

        threadloop %{
            /* set the factor input to zero to indicate that the intput isn't
             * factored. */
            factor_and_output_err = 0;
            qpgen2_(
                $P(dmat),
                $P(dvec),
                &m_size,
                &m_size,
                $P(sol),
                $P(lagr),
                $P(crval),
                $P(amat),
                $P(bvec),
                &m_size,
                &q_size,
                $P(meq),
                $P(iact),
                $P(nact),
                $P(iter),
                $P(work),
                &factor_and_output_err
            );
            /* Store the error results */
            $ierr() = factor_and_output_err;
        %}
    }),
    Doc => q{

=for ref

This routine solves the quadratic programming optimization problem

           minimize  f(x) = 0.5 x' D x  -  d' x
              x

    optionally constrained by:

            Aeq'  x  = a_eq
            Aneq  x >= b_neq


.... more docs to come ....
});

pp_add_exported('', 'qp_orig');
pp_addpm({At=>'Bot'},<<'EOD');

sub qp_orig {
  my ($Dmat, $dvec, $Amat, $avec, $meq) = @_;

  my $n = pdl $Dmat->dim(1);    # D is an [n x n] matrix
  my $q = pdl $Amat->dim(0);    # A is an [n x q] matrix

  if( $avec->isnull ){ $avec = zeros(1,$q); }

  croak("Dmat is not square!")
    if $n != $Dmat->dim(0);               # Check D is [n x n]
  croak("Dmat and dvec are incompatible!")
    if $n != $dvec->nelem;                # Check d is [n]
  croak("Amat and dvec are incompatible!")
    if $n != $Amat->dim(1);               # Check A is [n x _]
  croak("Amat and avec are incompatible!")
    if $q != $avec->nelem;                # Check A is [_ x q]
  croak("Value of meq is invalid!")
    if ($meq > $q) || ($meq < 0 );

  #  Pars => 'dmat(m,m); dvec(m);
  #      [o]sol(m); [o]lagr(q); [o]crval();
  #      amat(m,q); bvec(q); int meq();
  #      int [o]iact(q); int [o]nact();
  #      int [o]iter(s=2); [t] work(z); int [o]ierr();

  my ( $sol, $lagr, $crval, $iact, $nact, $iter, $ierr ) = qpgen2(
                   $Dmat->copy, $dvec->copy,
                   $Amat->transpose->copy,
                   $avec->copy,
                   $meq,
        );

  croak "qp: constraints are inconsistent, no solution!"
      if any($ierr == 1);
  croak "qp: matrix D in quadratic function is not positive definite!"
      if any($ierr == 2);
  croak "qp: some problem with mininization" if any($ierr);

  return {
    x     => $sol,
    lagr  => $lagr,
    crval => $crval,
    iact  => $iact,
    nact  => $nact,
    iter  => $iter,
    ierr  => $ierr,
  };

}
EOD

pp_add_exported('', 'qp');
pp_addpm pp_line_numbers(__LINE__, q{

sub qp {
  my ($Dmat, $dvec, %args) = @_;

  my $col   = 0;
  my $row   = 1;
  my $stack = 2;

  my $n = pdl $Dmat->dim($row);    # D is an [n x n] matrix
  my $s = pdl $Dmat->dim($stack);  # ... of $s stacked problems

  # Default handling for A_eq and A_neq
  my $A_eq  = exists $args{A_eq}  ? $args{A_eq}  : zeroes(0, $n);
  my $A_neq = exists $args{A_neq} ? $args{A_neq} : zeroes(0, $n);

  my $m = pdl $A_eq->dim($col);    # A is an [n x m] matrix
  my $p = pdl $A_neq->dim($col);   # A is an [n x p] matrix

  # Default handling for a_eq and a_neq
  # These have to be [?x1] matrixes to ensure threading works if glue,
  # otherwise they are relicated not just attached.
  my $a_eq  = exists $args{a_eq}  ? $args{a_eq}  : zeroes(1, $m);
  my $a_neq = exists $args{a_neq} ? $args{a_neq} : zeroes(1, $p);

  check_dimmensions( 'Dmat', $Dmat, $n, $n, $s );
  check_dimmensions( 'dvec', $dvec, $n, $s );
  check_dimmensions( 'A_eq', $A_eq, $m, $n, $s );
  check_dimmensions( 'a_eq', $a_eq, $m, $s );
  check_dimmensions( 'A_neq', $A_neq, $p, $n, $s );

  if( $p == 0 ){ # If a_neq has zero elems (per stack), then it will be [0x1]
      check_dimmensions( 'a_neq', $a_neq, 1, $p, $s );
  } else {       # otherwise it will be a [p] vector
      check_dimmensions( 'a_neq', $a_neq, $p, $s );
  }

  # Combine _eq and _neq, use meq to specifiy number of equality constratins
  my $A = $A_eq->glue( 0, $A_neq );
  my $a = $a_eq->glue( 0, $a_neq );
  my $meq = $A_eq->dim($col);

  #  Pars => 'dmat(m,m); dvec(m);
  #      [o]sol(m); [o]lagr(q); [o]crval();
  #      amat(m,q); bvec(q); int meq();
  #      int [o]iact(q); int [o]nact();
  #      int [o]iter(s=2); [t] work(z); int [io]ierr();

  my ( $sol, $lagr, $crval, $iact, $nact, $iter, $ierr ) = qpgen2(
                   $Dmat->copy, $dvec->copy,
                   $A->transpose->copy,
                   $a->copy,
                   $meq,
        );

  croak "qp: constraints are inconsistent, no solution!"
      if any($ierr == 1);
  croak "qp: matrix D in quadratic function is not positive definite!"
      if any($ierr == 2);
  croak "qp: some problem with mininization" if any($ierr);

  return {
    x     => $sol,
    lagr  => $lagr,
    crval => $crval,
    iact  => $iact,
    nact  => $nact,
    iter  => $iter,
    ierr  => $ierr,
  };

}

sub check_dimmensions {
    my ($name, $pdl, @dims) = @_;

    pop @dims if $dims[-1] == 1;    # remove last dim if a dimension of 1
    my $got      = pdl $pdl->dims;
    my $expected = pdl @dims;

    croak( sprintf( "dimmension check failed for %s (is %s)\nexpected [%s]\n%s",
            $name, $pdl->info, join(',',@dims), qp_usage() )
     ) unless all $got == $expected;
}

sub qp_usage {
    return q{
     usage: qp( Dmat [n x n], dvec [n], A_eq  => [n x m], a_eq  => [m],
                                        a_neq => [n x p], a_neq => [p] )
        where A_eq and a_eq are the equality constraints
          and A_neq and a_neq are the inequality constraints (>=)
     All inputs can add one addition dimmension to "stack" multiple
     Optimization problems. The returned solution (x) will have stack
     dimmensions.
    };
}
});

pp_addpm({At=>'Bot'},<<'EOD');

=head1 SEE ALSO

L<PDL>, L<PDL::Opt::NonLinear>

=head1 BUGS

Please report any bugs or suggestions at L<http://rt.cpan.org/>

=head1 AUTHOR

Mark Grimes, E<lt>mgrimes@cpan.orgE<gt>

=head1 COPYRIGHT AND LICENSE

This software is copyright (c) 2014 by Mark Grimes, E<lt>mgrimes@cpan.orgE<gt>.

This is free software; you can redistribute it and/or modify it under
the same terms as the Perl 5 programming language system itself.

=cut
EOD

pp_done();  # you will need this to finish pp processing