view release on metacpan or  search on metacpan

lib/Math/Brent.pm  view on Meta::CPAN

    print "Minimum is at sinc($x) = $y\n";

In either case the output will be C<Minimum is at sinc(4.4934094397196) = -.217233628211222>

This module has implementations of Brent's method for one-dimensional
minimisation of a function without using derivatives. This algorithm
cleverly uses both the Golden Section Search and parabolic
interpolation.

Anonymous subroutines may also be used as the function reference:

    my $cubic_ref = sub {my($x) = @_; return 6.25 + $x*$x*(-24 + $x*8));};

    my($x, $y) = Minimise1D(3, 1, $cubic_ref);
    print "Minimum of the cubic at $x = $y\n";

In addition to finding the minimum, there is also an implementation of the
Van Wijngaarden-Dekker-Brent Method, used to find a function's root without
using derivatives.

    use Math::Brent qw(Brentzero);

    my $tolerance = 1e-7;
    my $itmax = 80;

    sub wobble
    {
        my($t) = @_;
        return $t - cos($t);
    }

    #
    # Find the zero somewhere between .5 and 1.
    #
    $r = Brentzero(0.5, 1.0, \&wobble, $tolerance, $itmax);

=head1 EXPORT

Each function can be exported by name, or all may be exported by using the tag 'all'.

=head2 FUNCTIONS

The functions may be imported by name, or by using the export
tag "all".

=cut

=head3 Minimise1D()

Provides a simple interface to the L</BracketMinimum()> and L</Brent()>
routines.

Given a function, an initial guess for the function's
minimum, and its scaling, this routine converges
to the function's minimum using Brent's method.

    ($x, $y) = Minimise1D($guess, $scale, \&func);

The minimum is reached within a certain tolerance (defaulting 1e-7), and
attempts to do so within a maximum number of iterations (defaulting to 100).
You may override them by providing alternate values:

    ($x, $y) = Minimise1D($guess, $scale, \&func, 1.5e-8, 120);

=cut

sub Minimise1D
{
    my ($guess, $scale, $func, $tol, $itmax) = @_;
    my ($a, $b, $c) = BracketMinimum($guess - $scale, $guess + $scale, $func);

    return Brent($a, $b, $c, $func, $tol, $itmax);
}

#
# BracketMinimum
#
# BracketMinimum is MNBRAK minimum bracketing routine from section 10.1
# of Numerical Recipies
#

my $GOLD = 0.5 + sqrt(1.25); # Default magnification ratio for intervals is phi.
my $GLIMIT = 100.0; # Max magnification for parabolic fit step
my $TINY = 1E-20;

=head3 BracketMinimum()

    ($ax, $bx, $cx, $fa, $fb, $fc) = BracketMinimum($ax, $bx);

Given a function reference B<\&func> and distinct initial points B<$ax>
and B<$bx>, this routine searches in the downhill direction and returns
a list of the three points B<$ax>, B<$bx>, B<$cx> which bracket the
minimum of the function, along with the function values at those three
points, $fa, $fb, $fc.

The points B<$ax>, B<$bx>, B<$cx> may then be used in the function Brent().

=cut

sub BracketMinimum
{
    my ($ax, $bx, $func) = @_;
    my ($fa, $fb) = (&$func($ax), &$func($bx));

    #
    # Swap the a and b values if we weren't going in
    # a downhill direction.
    #
    if ($fb > $fa)
    {
	my $t = $ax; $ax = $bx; $bx = $t;
	$t = $fa; $fa = $fb; $fb = $t;
    }

    my $cx = $bx + $GOLD * ($bx - $ax);
    my $fc = &$func($cx);

    #
    # Loop here until we bracket
    #

lib/Math/Brent.pm  view on Meta::CPAN

	    }

	    $u = $cx + $GOLD * ($cx - $bx);
	    $fu = &$func($u);
	}
	elsif (($cx - $u) * ($u - $ulim) > 0)
	{
	    # parabolic  fit between C and limit
	    $fu = &$func($u);

	    if ($fu < $fc)
	    {
		$bx = $cx; $cx = $u;
		$u = $cx + $GOLD * ($cx - $bx);
		$fb = $fc; $fc = $fu;
		$fu = &$func($u);
	    }
	}
	elsif (($u - $ulim) * ($ulim - $cx) >= 0)
	{
	    # Limit parabolic U to maximum
	    $u = $ulim;
	    $fu = &$func($u);
	}
	else
	{
	    # eject parabolic U, use default magnification
	    $u = $cx + $GOLD * ($cx - $bx);
	    $fu = &$func($u);
	}

	# Eliminate oldest point and continue
	$ax = $bx; $bx = $cx; $cx = $u;
	$fa = $fb; $fb = $fc; $fc = $fu;
    }

    return ($ax, $bx, $cx, $fa, $fb, $fc);
}

#
# The complementary step is (3 - sqrt(5))/2, which resolves to 2 - phi.
#
my $CGOLD = 2 - $GOLD;
my $ZEPS = 1e-10;

=head3 Brent()

Given a function and a triplet of abcissas B<$ax>, B<$bx>, B<$cx>, such that

=over 4

=item 1. B<$bx> is between B<$ax> and B<$cx>, and

=item 2. B<func($bx)> is less than both B<func($ax)> and B<func($cx)>),

=back

Brent() isolates the minimum to a fractional precision of about B<$tol>
using Brent's method.

A maximum number of iterations B<$itmax> may be specified for this search - it
defaults to 100. Returned is a list consisting of the abcissa of the minum
and the function value there.

=cut

sub Brent
{
    my ($ax, $bx, $cx, $func, $tol, $ITMAX) = @_;
    my ($d, $u, $x, $w, $v); # ordinates
    my ($fu, $fx, $fw, $fv); # function evaluations

    $ITMAX //= 100;
    $tol //= 1e-8;

    my $a = min($ax, $cx);
    my $b = max($ax, $cx);

    $x = $w = $v = $bx;
    $fx = $fw = $fv = &$func($x);
    my $e = 0.0; # will be distance moved on the step before last
    my $iter = 0;

    while ($iter < $ITMAX)
    {
	my $xm = 0.5 * ($a + $b);
	my $tol1 = $tol * abs($x) + $ZEPS;
	my $tol2 = 2.0 * $tol1;

	last if (abs($x - $xm) <= ($tol2 - 0.5 * ($b - $a)));

	if (abs($e) > $tol1)
	{
	    my $r = ($x-$w) * ($fx-$fv);
	    my $q = ($x-$v) * ($fx-$fw);
	    my $p = ($x-$v) * $q-($x-$w)*$r;

	    $p = -$p if (($q = 2 * ($q - $r)) > 0.0);

	    $q = abs($q);
	    my $etemp = $e;
	    $e = $d;

	    unless ( (abs($p) >= abs(0.5 * $q * $etemp)) or
		($p <= $q * ($a - $x)) or ($p >= $q * ($b - $x)) )
	    {
                #
	        # Parabolic step OK here - take it.
                #
	        $d = $p/$q;
	        $u = $x + $d;

	        if ( (($u - $a) < $tol2) or (($b - $u) < $tol2) )
	        {
		    $d = copysign($tol1, $xm - $x);
	        }
	        goto dcomp; # Skip the golden section step.
	    }
	}

        #