File Coverage

blib/lib/Algorithm/CurveFit.pm

Criterion	Covered	Total	%
statement	123	135	91.1
branch	27	48	56.2
condition	13	36	36.1
subroutine	11	11	100.0
pod	1	1	100.0
total	175	231	75.7

package Algorithm::CurveFit;

use 5.006;
use strict;
use warnings;

our $VERSION = '1.06';

require Exporter;

our @ISA = qw(Exporter);

our %EXPORT_TAGS = (
    'all' => [ qw( curve_fit ) ]
);

our @EXPORT_OK = ( @{ $EXPORT_TAGS{'all'} } );

our @EXPORT = qw();

use Carp qw/confess/;
use Math::Symbolic qw/parse_from_string/;
use Math::MatrixReal;
use Data::Dumper;

# machine epsilon
use constant EPS => 2.2e-16;
use constant SQRT_EPS => sqrt(EPS);

sub curve_fit {
    shift @_ if not ref $_[0] and defined $_[0] and $_[0] eq 'Algorithm::CurveFit';

    confess('Uneven number of arguments to Algorithm::CurveFit::curve_fit.')
      if @_ % 2;

    my %args = @_;

    # Formula
    confess("Missing 'formula' parameter.") if not defined $args{formula};
    my $formula;
    if (ref($args{formula}) =~ /^Math::Symbolic/) {
        $formula = $args{formula};
    }
    else {
        eval { $formula = parse_from_string( $args{formula} ); };
        confess( "Cannot parse formula '" . $args{formula} . "'. ($@)" )
          if not defined $formula or $@;
    }

    # Variable (optional)
    my $variable = $args{variable};
    $variable = 'x' if not defined $variable;
    confess("Formula '"
          . $args{formula}
          . "' not explicitly dependent on "
          . "variable '$variable'." )
      if not grep { $_ eq $variable } $formula->explicit_signature();

    # Parameters
    my $params = $args{params};

    confess("Parameter 'params' has to be an array reference.")
      if not defined $params
      or not ref($params) eq 'ARRAY';
    my @parameters = @$params;
    confess('No parameters specified.') if not @parameters;
    confess('Individual parameters need to be array references.')
      if grep { not defined $_ or not ref($_) eq 'ARRAY' } @parameters;
    foreach my $p (@parameters) {
        confess("Weird parameter\n'"
              . Dumper($p)
              . "' Should have the format\n"
              . "[ NAME_STRING, GUESSED_VALUE, ACCURACY ]\n"
              . "With the accuracy being optional. See docs." )
          if @$p > 3
          or @$p < 2
          or grep { not defined $_ } @$p;

        confess("Formula '"
              . $args{formula}
              . "' not explicitly dependent on "
              . "parameter '"
              . $p->[0]
              . "'." )
          if not grep { $_ eq $p->[0] } $formula->explicit_signature();
    }

    # XData
    my $xdata = $args{xdata};
    confess('X-Data missing.')
      if not defined $xdata
      or not ref($xdata) eq 'ARRAY'
      or not @$xdata;
    my @xdata = @$xdata;

    # YData
    my $ydata = $args{ydata};
    confess('Y-Data missing.')
      if not defined $ydata
      or not ref($ydata) eq 'ARRAY'
      or not @$ydata;
    confess('Y-Data and X-Data need to have the same number of elements.')
      if not @$ydata == @xdata;
    my @ydata = @$ydata;

    # Max_Iter (optional)
    my $max_iter = $args{maximum_iterations};
    $max_iter = 0 if not defined $max_iter;

    # Add third element (dlamda) to parameter arrays in case they're missing.
    foreach my $param (@parameters) {
        push @$param, 0 if @$param < 3;
    }

    # Array holding all first order partial derivatives of the function in respect
    # to the parameters in order.
    my @derivatives;
    my @param_names = ($variable, map {$_->[0]} @parameters);
    foreach my $param (@parameters) {
        my $deriv =
          Math::Symbolic::Operator->new( 'partial_derivative', $formula,
            $param->[0] );
        $deriv = $deriv->simplify()->apply_derivatives()->simplify();
        my ($sub, $trees) = Math::Symbolic::Compiler->compile_to_sub($deriv, \@param_names);
        if ($trees) {
            push @derivatives, $deriv; # residual trees, need to evaluate
        } else {
            push @derivatives, $sub;
        }
    }

    # if not compilable, close over a ->value call for convenience later on
    my $formula_sub = do {
        my ($sub, $trees) = Math::Symbolic::Compiler->compile_to_sub($formula, \@param_names);
        $trees
        ? sub {
              $formula->value(
                map { ($param_names[$_] => $_[$_]) } 0..$#param_names
              )
          }
        : $sub
    };

    my $dbeta;

    # Iterative approximation of the parameters
    my $iteration = 0;

    # As long as we're under max_iter or maxiter==0
    while ( !$max_iter || ++$iteration < $max_iter ) {
        # Generate Matrix A
        my @cols;
        my $pno = 0;
        my @par_values = map {$_->[1]} @parameters;
        foreach my $param (@parameters) {
            my $deriv = $derivatives[ $pno++ ];

            my @ary;
            if (ref $deriv eq 'CODE') {
                foreach my $x ( 0 .. $#xdata ) {
                    my $xv = $xdata[$x];
                    my $value = $deriv->($xv, @par_values);
                    if (not defined $value) { # fall back to numeric five-point stencil
                        my $h = SQRT_EPS*$xv; my $t = $xv + $h; $h = $t-$xv; # numerics. Cf. NR
                        $value = $formula_sub->($xv, map { $_->[1] } @parameters)
                    }
                    push @ary, $value;
                }
            }
            else {
                $deriv = $deriv->new; # better safe than sorry
                foreach my $x ( 0 .. $#xdata ) {
                    my $xv = $xdata[$x];
                    my $value = $deriv->value(
                      $variable => $xv,
                      map { ( @{$_}[ 0, 1 ] ) } @parameters # a, guess
                    );
                    if (not defined $value) { # fall back to numeric five-point stencil
                        my $h = SQRT_EPS*$xv; my $t = $xv + $h; $h = $t-$xv; # numerics. Cf. NR
                        $value = $formula_sub->($xv, map { $_->[1] } @parameters)
                    }
                    push @ary, $value;
                }
            }
            push @cols, \@ary;
        }

        # Prepare matrix of datapoints X parameters
        my $A = Math::MatrixReal->new_from_cols( \@cols );

        # transpose
        my $AT = ~$A;
        my $M  = $AT * $A;

        # residuals
        my @beta =
          map {
            $ydata[$_] - $formula_sub->(
                $xdata[$_],
                map { $_->[1] } @parameters
              )
          } 0 .. $#xdata;
        $dbeta = Math::MatrixReal->new_from_cols( [ \@beta ] );

        my $N = $AT * $dbeta;

        # Normalize before solving => better accuracy.
        my ( $matrix, $vector ) = $M->normalize($N);

        # solve
        my $LR = $matrix->decompose_LR();
        my ( $dim, $x, $B ) = $LR->solve_LR($vector);

        # extract parameter modifications and test for convergence
        my $last = 1;
        foreach my $pno ( 1 .. @parameters ) {
            my $dlambda = $x->element( $pno, 1 );
            $last = 0 if abs($dlambda) > $parameters[ $pno - 1 ][2];
            $parameters[ $pno - 1 ][1] += $dlambda;
        }
        last if $last;
    }

    # Recalculate dbeta for the squared residuals.
    my @beta =
      map {
        $ydata[$_] - $formula_sub->(
            $xdata[$_],
            map { $_->[1] } @parameters
          )
      } 0 .. $#xdata;
    $dbeta = Math::MatrixReal->new_from_cols( [ \@beta ] );

    my $square_residual = $dbeta->scalar_product($dbeta);
    return $square_residual;
}

1;
__END__