package PDL::PP::Include::Fit::ExpRate;
use strict;
use warnings;

=head1 NAME

PDL::PP::Include::Fit::ExpRate - build-time includes to expose ExpRate C functions

=head1 SYNOPSIS

 # In your PDL::PP file:
 use PDL::PP::Include::Fit::ExpRate 'PDL';
 
 # Later in your PP code
 pp_def ('your_pp_function', Pars => 'xs(n); ys(n); ...',
     Code => q{
         ...
         double A, B, lambda, sum_sq_error;
         
         /* Get the initial parameter estimate for an exponential fit */
         exprate_estimate_parameters($P(xs), $P(ys), $SIZE(n), 
             &A, &B, &lambda, &sum_sq_error);
         ...
     },

=head1 DESCRIPTION

L<PDL::Fit::ExpRate> provides a method for fitting a noisy time series
to an exponential decay. The method is built upon a number of simpler
C functions which themselves have PDL methods. This module is meant to
be included in your L<PDL::PP>-based files (i.e. .pd files) so that your
own C code can access the same C functions used by L<PDL::Fit::ExpRate>.



working here - mention that these are function pointers, not functions



This module accomplishes all of this by calling a couple of pp-functions
when you say C<use PDL::PP::Include::Fit::ExpRate>. These pp-functions
add code to the generated XS files generated by PDL::PP. Having "use"d this
module, you can safely call the following C functions in the code for your
C<pp_addxs> and C<pp_def> declarations.

The documentation for each function lists the arguments to each function.
Since C only allows a single output, and because C does not perform garbage
collection, a number of these functions take pre-allocated containers as
their input, writing the results of their computations to these containers.
The argument lists always begin with the inputs and end with the outputs.
The basic format looks like this:

  return-type func-name (
      /* Input  */ type var1, type var2, ...
      /* Output */ type * var1, type * var2, ...
  );

=over

=item exprate_three_by_three_Householder

  void exprate_three_by_three_Householder (
      /* Input  */ double A[3][3], double y[3],
      /* Output */ double x[3]
  );

This function performs Householder diagonalization followed by
back-subsitituion. This effectively achieves matrix inversion, solving for
x the equation

 A * x = y

=item exprate_fit_quadratic

  void exprate_fit_quadratic (
      /* Input  */ double * xs, double * ys, int N,
      /* Output */ double * coefs
  );

Given N x-y data points, performs a least-squares quadratic fit to the data,
storing the results in the coefs array. The coefs array should have enough
room for three elements.

=item exprate_accum_sum_sq_err

  double exprate_accum_sum_sq_err (
      /* Input  */ double * xs, double * ys, int N,
                   double A, double B, double lambda
      /* Output  --  none, returns the result */
  );

Given N x-y data points and the parameters for the exponential curve A, B,
and lambda, this calculates the sum of the squared error, which can be
written as

                     /           x[i] * lambda          \
 sum_sq_error = sum (  A + B * e                 - y[i]  )
                     \                                  /

=item exprate_estimate_parameters

  void exprate_estimate_parameters (
      /* Input  */ double * xs, double * ys, int N,
      /* Output */ double * A, double * B, double * lambda,
                   double * sum_sq_error
  );

Given N x-y data points, this function uses a quadratic fit to estimate the
parameters A, B, and lambda for an exponential fit. It also calculates the
sum of the squared errors. To call this function, you should pass it the
address of the variables for A, B, lambda, and sum_sq_error where you want
the stored results. For example:

 double A, B, lambda, sum_sq_error;
 exprate_estimate_parameters(xs, ys, N, &A, &B, &lambda, &sum_sq_error);

=item exprate_newton_method_step

  int exprate_newton_method_step (
      /* Input  */ double * xs, double * ys, int N, double trust_radius,
      /* InOut  */ double * A, double * B, double * lambda,
      /* Output */ double * sum_sq_error
  );

Performs a single guarded Newton update step on the parameters A, B, and
lambda given N x-y data points and a specified trust radius. Note that,
unlike many of the previous functions, the values for A, B, and lambda are
both used as input and modified as output. The return value is a boolean
value that indicates whether or not the step had to be shortened to fit
within the trust radius.

This function works as follows. First, we assume that A, B, and lambda are
estimates of your data's actual coefficients. How should we modify these
values to obtain better estimates? Given the data and the current value of
the coefficients, this function analyticaly computes the gradient and
Hessian matrix in the coefficient space, and then inverts the Hessian to
obtain the adjustments to A, B, and lambda:

 Hessian * adjustments = -gradient

This would give a perfect set of adjustments if the function were a quadratic.
However, the function of interest is an exponential, which is not perfectly
approximated by a quadratic fit. Thus the adjustments need to be santy-checked
against the original values of A, B, and lambda. If the adjustment for A is
larger than C<trust_radius * A>, then the whole adjustment needs to be scaled
down until it is no larger than C<trust_radius * A>. Similar checks are
performed on B and lambda. Once we have obtained an adjustment that leads to
corrections that all fall within the trust radius, we apply that adjustment
to the coefficients and compute the sum_sq_error.

If it happens that the adjustment had to be corrected, the function returns 
a value of 1. If the adjustment fell within the trust radius, the function
returns a value of 0.

=back

=cut

#########################
# Function delcarations #
#########################

our $function_declarations = q{
	void exprate_three_by_three_Householder (double A[3][3], double y[3], double x[3]);
	void exprate_fit_quadratic (double *, double *, int, double *);
	double exprate_accum_sum_sq_err (double *, double *, int, double, double, double);
	void exprate_estimate_parameters (double *, double *, int, double *, double *, double *, double *);
	int exprate_newton_method_step (double *, double *, int, double, double *, double *, double *, double *);
};

#########################################
# Convert those to pointer declarations #
#########################################

our $pointer_declarations = $function_declarations;
$pointer_declarations =~ s/(\t\w+ )(\w+)/$1(*$2)/g;

##############################################
# XS code that assigns the function pointers #
##############################################

our $boot_code = q{

	/* Ensure that PDL::Fit::ExpRate is loaded */
	load_module(PERL_LOADMOD_NOIMPORT, newSVpv("PDL::Fit::ExpRate", 0), 0);
	
	/* Assign the pointers */
	SV * the_sv = get_sv("PDL::Fit::ExpRate::__householder_func_addr", 0);
	exprate_three_by_three_Householder = INT2PTR(void*, SvIV(the_sv));
	the_sv = get_sv("PDL::Fit::ExpRate::__quadratic_func_addr", 0);
	exprate_fit_quadratic = INT2PTR(void*, SvIV(the_sv));
	the_sv = get_sv("PDL::Fit::ExpRate::__accum_func_addr", 0);
	exprate_accum_sum_sq_err = INT2PTR(void*, SvIV(the_sv));
	the_sv = get_sv("PDL::Fit::ExpRate::__estimate_func_addr", 0);
	exprate_estimate_parameters = INT2PTR(void*, SvIV(the_sv));
	the_sv = get_sv("PDL::Fit::ExpRate::__newton_func_addr", 0);
	exprate_newton_method_step = INT2PTR(void*, SvIV(the_sv));

};

sub import {
	if (grep /PDL/, @_) {
		# Make sure they already included PDL::PP
		croak("You must include PDL::PP before PDL::PP::Include::Fit::Exprate")
			unless defined $PDL::PP::VERSION;
		
		PDL::PP::pp_addhdr $pointer_declarations;
		PDL::PP::pp_add_boot $boot_code;
	}
}

1;