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