File Coverage

blib/lib/PDL/Opt/Simplex.pm

Criterion	Covered	Total	%
statement	64	67	95.5
branch	19	24	79.1
condition	3	5	60.0
subroutine	6	6	100.0
pod	0	2	0.0
total	92	104	88.4

package PDL::Opt::Simplex;
use strict;
use warnings;
use PDL;
use PDL::Exporter;

our @ISA = qw/PDL::Exporter/;
our @EXPORT_OK = qw/simplex make_simplex/;
our %EXPORT_TAGS = ( Func => [@EXPORT_OK] );
our $VERSION = '2.097';

*simplex = \&PDL::simplex;

sub make_simplex {
  my ($init, $initsize) = @_;
  my ( $nd, $nd2, @otherdims ) = ( $init->dims, 1 );
  pop @otherdims if @otherdims; # drop the spurious "1"
  return unless $nd2 == $nd + 1 or $nd2 == 1;
  return $init if $nd2 == $nd + 1;
  my $simp = PDL->zeroes( $nd, $nd + 1, @otherdims );
  # Constructing a tetrahedron:
  # At step n (starting from zero)
  # take vertices 0..n and move them 1/(n+1) to negative dir on axis n.
  # Take vertex n+1 and move it n/(n+1) to positive dir on axis n
  my $uppertri = $simp->copy;
  ones($nd * ($nd+1) / 2)->tritosquare($uppertri->slice(",:-2")->t);
  my $iseq = sequence($nd);
  my $pjseq = $iseq / ( $iseq + 1 );
  $simp .= $init - $initsize * $pjseq * $uppertri;
  $simp->slice(",1:")->diagonal(0,1) += $initsize * (1-$pjseq);
  $simp;
}

sub PDL::simplex {
  my ( $init, $initsize, $minsize, $maxiter, $sub, $logsub, $t ) = @_;
  my $simp = make_simplex($init, $initsize) // return;
  my $nd = $simp->dim(0);
  my $vals = $sub->($simp);
  my $ssize;
  while ($maxiter--) {
    $ssize = ( $simp - $simp->slice(":,0") )->magnover->maxover;
    $logsub->( $simp, $vals, $ssize ) if $logsub;
    last unless $ssize > $minsize;
    my $valsn = !$t ? $vals : $vals - $t * log( $vals->random + 0.00001 );
    my $minind = $valsn->minimum_ind;
    my $maxind = $valsn->maxover_n_ind(2);
    my $maxind0 = $maxind->at(0);
    my @worstvals = $valsn->index($maxind)->list;
    my $bestval = $valsn->at($minind);
    my $ssum = ($simp->t->sumover - (my $worst = $simp->slice(":,($maxind0)"))) / $nd;
    my $new = 2 * $ssum - $worst;
    my $val = $sub->($new)->at(0);
    $val += $t * log( rand() + 0.00001 ) if $t;
    my $removetop = 0;
    if ( $val < $bestval ) {
      my $newnew = $new + $ssum - $worst;
      my $val2   = $sub->($newnew);
      if ( $val2->at(0) < $val ) { # CASE1 Reflection and Expansion
        $new = $newnew;
        $val = $val2;
      } # else CASE2 Reflection
      $removetop = 1;
    } elsif ( $val < $worstvals[1] ) { # CASE3 Reflection
      $removetop = 1;
    } else {
      my $newnew = 0.5 * ($ssum + $worst);
      my $val2   = $sub->($newnew);
      if ( $val2->at(0) < $worstvals[0] ) { # CASE4 Contraction
        $new = $newnew;
        $val = $val2;
        $removetop = 1;
      }
    }
    if ($removetop) {
      $worst .= $new;
      $vals->slice( "($maxind0)" ) .= $val;
    } else { # CASE5 Multiple Contraction
      $simp .= 0.5 * ($simp->slice(":,$minind") + $simp);
      my $idx = which( sequence($nd+1) != $minind );
      $vals->index($idx) .= $sub->($simp->dice_axis(1,$idx));
    }
  }
  my $mmind = $vals->minimum_ind;
  return ( $simp->slice(":,$mmind"), $ssize, $vals->index($mmind) );
}

1;

__END__