package GoldenSection;

# Craig Kelley -- November 20, 2001

# Performs a golden section search to discover the minimum value of a
# bracketed linear function

use strict;
use Carp qw(cluck confess);

# export constants to the user of this class, if they want them
require Exporter;
my @ISA = qw(Exporter);
my @EXPORT_OK = qw ($GOLDEN_RATIO $GOLDEN_INVERSE);

# Some constants
my $GOLDEN_RATIO = (3 - sqrt(5)) / 2;
my $GOLDEN_INVERSE = 1- $GOLDEN_RATIO;

#
# Public Methods
#

# CONSTRUCTORS

# new_with_point
#
# Creates a new object; a golden section search attempts to minimize a linear parabolic function
#
# The following are mandatory arguments:
#
#  code reference    : A subroutine that performs the linear function in question
#  float             : A starting point on the function
#  float             : A positive number that describes an acceptable variance from the
#                      actual minumum; if this is set to '0' an absolute minimum will be
#                      found, but the algorithm may run forever due to float limitations. A
#                      good heruistic to find an absolute minimum is the square root of
#                      your machine's floating point precision[* See Credits]
sub new_with_point {

  my ($that, $function, $starting_point, $delta) = @_;
  my $class = ref($that) || $that;
  my $self = {};

  get_set_function($self, $function);
  get_set_delta($self, $delta);
  find_bounds($self);

  bless $self, $class;
}

# new_with_bracket
#
# Creates a new object; a golden section search attempts to minimize a linear parabolic function
#
# The following are mandatory arguments:
#
#  code reference    : A subroutine that performs the linear function in question
#  float             : A positive number that describes an acceptable variance from the
#                      actual minumum; if this is set to '0' an absolute minimum will be
#                      found, but the algorithm may run forever due to float limitations. A
#                      good heruistic to find an absolute minimum is the square root of
#                      your machine's floating point precision[* See Credits]
#  float             : The left bound such that the minimum lies in the positive direction
#                      of this point
#  float             : The right bound such that the minimum lies in the negative direction
#                      of this point
sub new_with_bracket {

  my ($that, $function, $delta, $left, $right) = @_;
  my $class = ref($that) || $that;
  my $self = {};

  get_set_function($self, $function);
  get_set_bounds($self, $left, $right);
  get_set_delta($self, $delta);

  bless $self, $class;
}

# ACCESSORS / MUTATORS

# get_set_function
#
# returns the current function (CODE ref) and sets it to the passed-in argument, if present
sub get_set_function {

  my ($self, $new_function) = @_;

  if (defined $new_function) {
    # sanity checks
    confess "Passed in function must be a CODE reference"
      unless (ref($new_function) eq "CODE");
  }
  my $retval = $self->{FUNCTION};
  if (defined $new_function) {
    $self->{FUNCTION} = $new_function;
  }

  return $retval;
}

# get_set_delta
#
# returns the current delta value and sets it to the passed-in argument, if present
sub get_set_delta {

  my ($self, $new_delta) = @_;

  if (defined $new_delta) {
    # sanity checks
    confess "Delta needs to be a positive number"
      if ($new_delta < 0);
    cluck "Delta is set to zero, it is possible that the search will never end"
      if ($new_delta == 0);
  }
  my $retval = $self->{DELTA};
  if (defined $new_delta) {
    $self->{DELTA} = $new_delta;
  }

  return $retval;
}

# get_set_bounds
#
# returns the current bounds (left, right) and sets either to the passed-in arguments, if present
sub get_set_bounds {

  my ($self, $new_left, $new_right) = @_;

  my @retval = ($self->{LEFT_BOUND}, $self->{RIGHT_BOUND});
  if (defined $new_left) {
    $self->{LEFT_BOUND} = $new_left;
  }
  if (defined $new_right) {
    $self->{RIGHT_BOUND} = $new_right;
  }

  if ($self->{LEFT_BOUND} >= $self->{RIGHT_BOUND}) {
    ($self->{LEFT_BOUND}, $self->{RIGHT_BOUND}) = @retval;
    cluck "get_set_delta ($new_left, $new_right) would leave object in a bad state; left must" . 
      "be less than right : reseting values";
  }

  return @retval;
}

# minimize
#
# Minimizes the function; returns an acceptable minimization and takes no arguments
sub minimize {

  my $self = shift;

  my $min = undef;		# return value

  my $a = $self->{LEFT_BOUND};
  my $c = $self->{RIGHT_BOUND};
  my $b = $a + (($c - $a) / 2);

  my @X = ($a,
	   $b - $GOLDEN_RATIO * ($b - $a),
	   $b,
	   $c);

  my @Fvalues = ( $self->{FUNCTION}->($X[1]),
		  $self->{FUNCTION}->($X[2]));

  while (abs ($X[3] - $X[0]) > ($self->{DELTA} *
					    (abs($X[1]) +
					     abs($X[2])))) {
    if ($Fvalues[1] < $Fvalues[0]) {
      $X[0] = $X[1];
      $X[1] = $X[2];
      $X[2] = $GOLDEN_INVERSE * $X[1] + $GOLDEN_RATIO * $X[3];
      shift @Fvalues;
      push @Fvalues,  $self->{FUNCTION}->($X[2]);
    }
    else {
      $X[3] = $X[2];
      $X[2] = $X[1];
      $X[1] = $GOLDEN_INVERSE * $X[2] + $GOLDEN_RATIO * $X[0];
      pop @Fvalues;
      unshift @Fvalues, $self->{FUNCTION}->($X[1]);
    }

  }
  if ($Fvalues[0] < $Fvalues[1]) {
    return ($Fvalues[0], $X[1]);
  }
  else {
    return ($Fvalues[1], $X[2]);
  }

  return $min;

}

#
# PRIVATE METHODS
#

# find_bounds
#
# blindly bracket the minimum point in a linear function using the smiley
# heuristic in 3 points ^_^
sub find_bounds {

  my ($self, $point, $step_size) = @_;

  my $step_size = 2
    unless (defined $step_size);
  my $MAX_STEP_SIZE = 100;

  my ($left, $middle, $right);
  # initial points to try
  if (defined $point) {
    $left = $middle - $step_size;
    $right = $middle + $step_size;
  }
  else {
    $left = -1;
    $middle = 0;
    $right = 1;
  }

  my $f_left = $self->{FUNCTION}->($left);
  my $f_middle = $self->{FUNCTION}->($middle);
  my $f_right = $self->{FUNCTION}->($right);

  my $min_found = undef;
  while (!defined $min_found) {

    # do we have a smiley face?
    if ( ( $f_left > $f_middle) &&
	 ( $f_right > $f_middle) ) {
      $min_found = 1;
    }
    # downhill to the right?
    elsif ($f_left > $f_right) {
      $left = $middle;
      $f_left = $f_middle;
      $middle = $right;
      $f_middle = $f_right;
      $right += $step_size;
      $f_right = $self->{FUNCTION}->($right);
    }
    # downhill to the left?
    elsif ($f_left < $f_right) {
      $right = $middle;
      $f_right = $f_middle;
      $middle = $left;
      $f_middle = $f_left;
      $left -= $step_size;
      $f_left = $self->{FUNCTION}->($left);
    }
    else {
      # equal (very very unlikely to ever happen)
      $right -= ($step_size / (rand(5) + 1));
      $f_left = $self->{FUNCTION}->($left);
      $right += ($step_size / (rand(5) + 1));
      $f_right = $self->{FUNCTION}->($right);
    }
    $step_size *= 1.5;
    if ($step_size > $MAX_STEP_SIZE) {
      $step_size = $MAX_STEP_SIZE;
    }
  }

  get_set_bounds($self, $left, $right);
  return 1;

}

# CREDITS
#
# "Golden Section Search in One Direction"
# Numerical Recipies in C: The Art of Scientific Computing (ISBN 0-521-43108-5)
# Copyright 1988-1992, Cambridge University Press Programs
# Copyright 1988-1992, Numerical Recipies Software

# make perl happy
1;