view release on metacpan or  search on metacpan

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

I<L<Numerical Recipes Third Edition|http://numerical.recipes/aboutNR3book.html>>
book, section 9.1.

The function supports the following parameters:

=over

=item C<max_iterations>

How many iterations of our algorithm will be applied at most while trying to
bracket the given function. This gives an order of magnitude of the number of
times that C<$func> will be evaluated. Defaults to I<100>.

=item C<do_outward>

Whether the function will try to bracket a root in an interval larger than the
one given by C<[$x1, $x2]>. Defaults to I<1>.

=item C<do_inward>

Whether the function will try to bracket a root in an interval smaller than the
one given by C<[$x1, $x2]>. Defaults to I<1>.

One of C<do_inward> or C<do_outward> at least must be a  true value.

=item C<inward_split>

Tuning parameter describing the starting number of intervals into which the
starting interval is split when looking inward for a bracket. Defaults to I<3>.

Note that the algorithm may change and this parameter may stop working or may
take a different meaning in the future.

=item C<inward_factor>

Tuning parameter describing a factor by which the inwards interval are split
at each iteration. Defaults to I<3>.

Note that the algorithm may change and this parameter may stop working or may
take a different meaning in the future.

=item C<outward_factor>

Tuning parameter describing how much the starting interval is grown at each
iteration when looking outward for a bracket. Defaults to I<1.6>.

Note that the algorithm may change and this parameter may stop working or may
take a different meaning in the future.

=back

=cut

Readonly my $DEFAULT_INWARD_SPLIT => 3;
Readonly my $DEFAULT_INWARD_FACTOR => 3;
Readonly my $DEFAULT_OUTWARD_FACTOR => 1.6;

sub _create_bracket_inward_state ($x1, $x2, $f1, %params) {
  my $s = {ret => undef};
  $s->{split} = $params{inward_split} // $DEFAULT_INWARD_SPLIT;
  croak 'inward_split must be at least 2' if $s->{split} < 2;
  $s->{factor} = $params{inward_factor} // $DEFAULT_INWARD_FACTOR;
  croak 'inward_factor must be at least 2' if $s->{factor} < 2;
  @{$s}{'x1', 'x2'} = ($x1, $x2);
  $s->{f1} = $f1;
  lock_keys(%{$s});
  return $s;
}

sub _do_bracket_inward ($f, $s) {
  my $dx = ($s->{x2} - $s->{x1}) / $s->{split};
  my $xa = $s->{x1};
  my ($fa, $fb) = ($s->{f1});
  for my $j (1 .. $s->{split}) {
    my $xb = $xa + $dx;
    $fb = $f->($xb);
    if ($fa * $fb < 0) {
      $s->{ret} = [$xa, $xb, $fa, $fb];
      return 1;
    }
    ($xa, $fa) = ($xb, $fb);
  }
  $s->{split} *= $s->{factor};
  return 0;
}

sub _create_bracket_outward_state ($f, $x1, $x2, $f1, %params) {
  my $s = {ret => undef};
  $s->{factor} = $params{outward_factor} // $DEFAULT_OUTWARD_FACTOR;
  croak 'outward_factor must be larger than 1' if $s->{factor} <= 1;
  @{$s}{'x1', 'x2'} = ($x1, $x2);
  @{$s}{'f1', 'f2'} = ($f1, $f->($x2));
  lock_keys(%{$s});
  return $s;
}

sub _do_bracket_outward ($f, $s) {
  if ($s->{f1} * $s->{f2} < 0) {
    $s->{ret} = [@{$s}{'x1', 'x2', 'f1', 'f2'}];
    return 1;
  }
  if (abs($s->{f1}) < abs($s->{f2})) {
    $s->{x1} += $s->{factor} * ($s->{x1} - $s->{x2});
    $s->{f1} = $f->($s->{x1});
  } else {
    $s->{x2} += $s->{factor} * ($s->{x2} - $s->{x1});
    $s->{f2} = $f->($s->{x2});
  }
  return 0;
}

sub bracket ($func, $x1, $x2 = undef, %params) {
  if (!defined $x2 || $x1 == $x2) {
    Readonly my $LARGISH_FACTOR => 1000;
    $x2 += $LARGISH_FACTOR * $EPS;
  }
  my $max_iter = $params{max_iterations} // $DEFAULT_MAX_ITERATIONS;
  croak 'max_iterations must be positive' if $max_iter <= 0;

  my $f = _wrap_func($func);
  my $f1 = $f->($x1);