lib/PDL/Opt/QP.pd - metacpan.org

our $VERSION = '0.21';
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); int fddmat(); int n();
        [o]sol(m); [o]lagr(q); [o]crval();
        amat(m,q); bvec(q); int fdamat(); int q(); int meq();
        int [o]iact(q); int [o]nact();
        int [o]iter(s=2); [t]work(z); int [io]ierr();
    ',
    GenericTypes => [D],
    Code => '
        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
        );

        qpgen2_(
            $P(dmat),
            $P(dvec),
            $P(fddmat),
            $P(n),
            $P(sol),
            $P(lagr),
            $P(crval),
            $P(amat),
            $P(bvec),
            $P(fdamat),
            $P(q),
            $P(meq),
            $P(iact),
            $P(nact),
            $P(iter),
            $P(work),
            $P(ierr)
        );

        // Not sure if we will need to process the solutions here
        // for (i = 0; i < $SIZE(n); i++)
        //  $x(n=>i) = xtmp[i];
        // $maxit()=it;
',
    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 = pdl()->reshape(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 );

  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

  # print "\$A = ", $Amat;
  # print "\$a = $avec\n";
  # print "n = $n\n";
  # print "q = $q\n";
  # print "meq = $meq\n";

  my $res = qpgen2(
                   $Dmat->copy, $dvec->copy,
                   $n, $n,
                   $sol, $lagr,
                   $crval,
                   $Amat->transpose->copy,
                   $avec->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, $A_eq, $a_eq, $A_neq, $a_neq) = @_;

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

  my $n = pdl $Dmat->dim($row);    # D is an [n x n] matrix
  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
  my $q = $m + $p;    # A is an [n x q] matrix

  if( $A_eq->isnull ){ $A_eq = pdl()->reshape(0,$n); }
  if( $a_eq->isnull ){ $a_eq = pdl()->reshape(1,$m); }
  if( $A_neq->isnull ){ $A_neq = pdl()->reshape(0,$n); }
  if( $a_neq->isnull ){ $a_neq = pdl()->reshape(1,$p); }

  croak("dimmension check failed: Dmat [n x n*] is not square")
    if $Dmat->dim($col) != $n;
  croak("dimmension check failed: Dmat [n x n] and dvec [n* x 1]")
    if $dvec->nelem != $n;
  croak("dimmension check failed: A_eq [n* x m] and a_eq [n x 1]")
    if $A_eq->dim($row) != $n;
  croak("dimmension check failed: A_eq [n x m] and a_eq [m* x 1]")
    if $a_eq->nelem != $m;
  croak("dimmension check failed: A_neq [n* x p] and a_neq [p x 1]")
    if $A_neq->dim($row) != $n;
  croak("dimmension check failed: A_neq [n x p] and a_neq [p* x 1]")
    if $a_neq->nelem != $p;

  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 = $A_eq->glue( 0, $A_neq )->copy;
  my $a = $a_eq->glue( 0, $a_neq )->copy;  # ->flat?
  my $meq = $A_eq->dim(0);

  # print "\$A = ", $A;
  # print "\$a = $a\n";
  # print "n = $n\n";
  # print "q = $q\n";
  # print "meq = $meq\n";

  my $res = qpgen2(
                   $Dmat->copy, $dvec->copy,
                   $n, $n,
                   $sol, $lagr,
                   $crval,
                   $A->transpose->copy,
                   $a->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_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) 2013 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
	Global
`s`	Focus search bar
`?`	Bring up this help dialog
	GitHub
`g` `p`	Go to pull requests
`g` `i`	Go to GitHub issues (only if GitHub is preferred repository)
	POD
`g` `a`	Go to author
`g` `c`	Go to changes
`g` `i`	Go to issues
`g` `d`	Go to dist
`g` `r`	Go to repository/SCM
`g` `s`	Go to source
`g` `b`	Go to file browse
Search terms
module: (e.g. module:Plugin)
distribution: (e.g. distribution:Dancer auth)
author: (e.g. author:SONGMU Redis)
version: (e.g. version:1.00)
Keyboard Shortcuts
