package PDL::Opt::Powell;

use strict;
use warnings;
use Math::Brent;
use PDL;

# ($optimum,$ssize,$optval) = powell( $init, $initsize, $minsize, $maxiter, $sub, ( $log_sub, $U ) );
#

sub powell {
    my $init     = shift;
    my $initsize = shift;
    my $minsize  = shift;
    my $maxiter  = shift;
    my $sub      = shift;
    my $log_sub  = shift;

    my $nd        = $init->dim(0);
    my $U    = shift(@_) || identity($nd); # columns of directions

    my $maxiter1d = $maxiter;
    my $minsize1d = $minsize;

    if ( !ref $initsize ) {
        $initsize = ones($nd) * $initsize; # initial step size
    }

    # ($nd-1) indexes of the columns of $U 
    my $minidx = zeroes($nd-1);

    my $P0   = $init;
    my $f0   = $sub->($P0);
    my $Pcur = $P0;
    my $fcur = $f0;
    my $ssize; # RMS distance
    while ( $maxiter-- ) {
        my @deltaf; # decrements at each direction
        for my $i ( 0 .. $nd - 1 ) {
            my $Ui    = $U->slice("($i)"); # currently moving direction
            # 1D sub for 1D minimization
            my $sub1d = sub {
                my $res = $sub->( $Pcur + $_[0] * $Ui );
                return ( ( ref $res ) ? $res->sclr : $res );
            };
            my $initsize1d = inner( $initsize, $Ui )->sclr;
            my ( $Ucoeff, $fnext ) =
              Math::Brent::Minimise1D( 0, $initsize1d, $sub1d, $minsize1d,
                $maxiter1d );

            $Pcur = $Pcur + $Ucoeff * $Ui;
            push @deltaf, $fcur - $fnext;
            $fcur = $fnext;
        }
        $ssize = sqrt( sum( ( $P0 - $Pcur )**2 ) ); #dRMS

        if ($log_sub) {
            # returns variables, value, dRMS, vector_set
            # for display
            $log_sub->( $Pcur, $fcur, $ssize, $U );
        }
        last if ( $ssize < $minsize ); # termination criterion

        my $deltaf = pdl(@deltaf);
        my $dfmax  = $deltaf->max;

        my $new  = 2 * $Pcur - $P0;
        my $fnew = $sub->($new);

        if (
            # This acceptance rule is from Numerical Recipes 3rd Ed.
            ( $fnew < $f0 )
            &&
            (
                2 * ( $f0 - 2*$fcur + $fnew ) * ( $f0 - $fcur - $dfmax )**2 <
                ( $f0 - $fnew )**2 * $dfmax
            )
          )
        {
            my $Unew = norm( $Pcur - $P0 );
#            print "New direction is accepted $Unew\n";
            minimum_n_ind( $deltaf, $minidx );
            $U = append( $Unew->transpose, $U->dice_axis( 0, $minidx ) );
        }

        $P0 = $Pcur;
        $f0 = $fcur;
    }

    return ( $P0, $ssize, $fcur );
}

1;