################################################################
package optimize;
################################################################
=pod

=head1 NAME

  optimize

=head1 SUBROUTINES

  goldensearch()

=head1 SYNOPSIS

I<goldensearch> is a unidimensional minimization function.  It
is not very fast, but it is very solid.

I<goldensearch> was essentially adapted from Numerical Recipees.
It adds some checking, plus pass-through of fixed parameters (as
provided by readdata) into the function.

=head1 USAGE

  goldensearch( very-negative-x, mid-x, very-positive-x, tolerance, function-reference, passthrough-to-function )

where

  both very-negative-x and very-positive-x should have a value above mid-x, function-reference
  is something like &functionname, and passthrough can be anything.

=cut

################################################################

$VERSION = 1.0;
use constant DEBUG=>0;		# set to 1 if you want to see more about the search

use strict;

sub goldensearch {
  my $ax= shift(@_);  my $bx= shift(@_);  my $cx= shift(@_);
  my $tol= shift(@_);
  my $function= shift(@_);

  # the remainder are fixed parameters that are passed through

  use constant R => 0.6180339;
  use constant C => 0.3819661;

  if ($function->($ax, @_) < $function->($bx,  @_)) {
    die "1: sorry, but f($ax)=".($function->($ax, @_))." is less than f($bx)=".($function->($bx, @_))."!\n"; }
  if ($function->($cx,  @_) < $function->($bx,  @_)) { 
    die "2: sorry, but f($cx)=".($function->($cx,  @_))." is less than f($bx)=".($function->($bx,  @_))."!\n"; }

  my ($f0,$f3);

  my $x0= $ax;
  my $x3= $cx;
  my ($x1, $x2);
  if (abs($cx-$bx) > abs($bx-$ax)) {
    $x1= $bx;
    $x2=$bx+C*($cx-$bx);
  }
  else {
    $x2= $bx;
    $x1= $bx-(C*($bx-$ax));
  }
  my $f1= $function->($x1, @_);
  my $f2= $function->($x2, @_);
  while (abs($x3-$x0) > $tol*(abs($x1)+abs($x2))) {
    if ($f2<$f1) {
      $x0=$x1; $x1=$x2; $x2=R*$x1+C*$x3;
      $f0=$f1; $f1=$f2; $f2= $function->($x2, @_);
      if (DEBUG) { print STDERR "[golden: trying $x2 -> $f2]\n"; }
    }
    else {
      $x3=$x2; $x2=$x1; $x1=R*$x2+C*$x0;
      $f3=$f2; $f2=$f1; $f1= $function->($x1, @_);
      if (DEBUG) { print STDERR "[golden: trying $x1 -> $f1]\n"; }
    }
  }
  if ($f1<$f2) { return ($f1, $x1); }
  else { return ($f2, $x2); }
}
  
1;