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:

               Aeq'  x  = a_eq
               Aneq  x >= b_neq

   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');

SEE ALSO
   PDL, PDL::Opt::NonLinear

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

AUTHOR
   Mark Grimes, <[email protected]>

COPYRIGHT AND LICENSE
   This software is copyright (c) 2014 by Mark Grimes, <[email protected]>.

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

You are viewing proxied material from ftp.icm.edu.pl. The copyright of proxied material belongs to its original authors. Any comments or complaints in relation to proxied material should be directed to the original authors of the content concerned. Please see the disclaimer for more details.