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

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};

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:

                A' x  = a
                B  x >= b

    This routine solves the quadratic programming optimization problem

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

        optionally constrained by:

                A' x  = a
                B  x >= b

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

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

    sub qp_orig { my ($Dmat, $dvec, $Amat, $bvec, $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( ! defined $bvec ){    # if b is undef, create it
        $bvec = zeros($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 bvec are incompatible!")
        if $q != $bvec->nelem;                # Check A is [_ x q]
      croak("Value of meq is invalid!")
        if ($meq > $q) || ($meq < 0 );

      my $iact = zeros($q);              # Store which constraints are active
      my $nact = pdl(0);                 # Store number of active constraints
      my $r    = $n < $q ? $n : $q;      # Used to size work space
      my $sol  = zeros($n->sclr);              # Store the solution [n]
      my $lagr = zeros($q->sclr);              # Store the Lagranges for the constraints
      my $crval= pdl(0);                 # Value at min
      my $work = zeros((2*$n+$r*($r+5)/2+2*$q+1)->sclr);  # Work space
      my $iter = zeros(2);               # Store info about interations
      my $ierr = pdl(0);                 # Input: 1=Factorized; Output: error flag

      my $res = qpgen2(
                       $Dmat->copy, $dvec->copy,
                       $n, $n,
                       $sol, $lagr,
                       $crval,
                       $Amat->transpose->copy,
                       $bvec->copy, $n,
                       $q, $meq,
                       $iact, $nact,
                       $iter, $work,
                       $ierr
            );

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

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

      # TODO: process/return the results
      #
      # From R implementation:
      #
      # list(solution=res1$sol,
      #      value=res1$crval,
      #      unconstrained.solution=res1$dvec,
      #      iterations=res1$iter,
      #      Lagrangian = res1$lagr,
      #      iact=res1$iact[1:res1$nact])   
    }
    EOD

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

    sub qp { my ($Dmat, $dvec, $Amat, $avec, $Bmat, $bvec) = @_;

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

      if( ! defined $bvec ){    # if b is undef, create it
        $bvec = zeros($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("Bmat and dvec are incompatible!")
        if $n != $Bmat->dim(1);               # Check B is [n x _]
      croak("Amat and avec are incompatible!")
        if $Amat->dim(0) != $avec->nelem;                # Check A is [_ x q]
      croak("Bmat and bvec are incompatible!")
        if $Bmat->dim(0) != $bvec->nelem;                # Check B is [_ x q]

      my $iact = zeros($q);              # Store which constraints are active
      my $nact = pdl(0);                 # Store number of active constraints
      my $r    = $n < $q ? $n : $q;      # Used to size work space
      my $sol  = zeros($n->sclr);              # Store the solution [n]
      my $lagr = zeros($q->sclr);              # Store the Lagranges for the constraints
      my $crval= pdl(0);                 # Value at min
      my $work = zeros((2*$n+$r*($r+5)/2+2*$q+1)->sclr);  # Work space
      my $iter = zeros(2);               # Store info about interations
      my $ierr = pdl(0);                 # Input: 1=Factorized; Output: error flag

      my $A = $Amat->glue( 0, $Bmat )->copy;
      my $a = $avec->glue( 0, $bvec )->copy;  # ->flat?
      my $meq = $Amat->dim(0);

      print "\$A' = ", $A->transpose;
      print "\$a = ", $a;
    print "n = $n";
    print "q = $q";
    print "meq = $meq";
    print "\n";

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

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

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

      # TODO: process/return the results
      #
      # From R implementation:
      #
      # list(solution=res1$sol,
      #      value=res1$crval,
      #      unconstrained.solution=res1$dvec,
      #      iterations=res1$iter,
      #      Lagrangian = res1$lagr,
      #      iact=res1$iact[1:res1$nact])   
    }
    EOD

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

SEE ALSO
    PDL, PDL::Opt::NonLinear

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

AUTHOR
    Mark Grimes, <mgrimes@cpan.org>

COPYRIGHT AND LICENSE
    This software is copyright (c) 2013 by Mark Grimes, <mgrimes@cpan.org>.

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