NAME
PDL::PP::Include::Fit::ExpRate  buildtime includes to expose ExpRate C functions
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);
...
},
DESCRIPTION
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 PDL::PPbased files (i.e. .pd files) so that your own C code can access the same C functions used by PDL::Fit::ExpRate.
working here  mention that these are function pointers, not functions
This module accomplishes all of this by calling a couple of ppfunctions when you say use PDL::PP::Include::Fit::ExpRate
. These ppfunctions 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 pp_addxs
and 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 preallocated 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:
returntype funcname (
/* Input */ type var1, type var2, ...
/* Output */ type * var1, type * var2, ...
);
 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 backsubsitituion. This effectively achieves matrix inversion, solving for x the equation
A * x = y
 exprate_fit_quadratic

void exprate_fit_quadratic ( /* Input */ double * xs, double * ys, int N, /* Output */ double * coefs );
Given N xy data points, performs a leastsquares quadratic fit to the data, storing the results in the coefs array. The coefs array should have enough room for three elements.
 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 xy 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] ) \ /
 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 xy 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);
 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 xy 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 santychecked against the original values of A, B, and lambda. If the adjustment for A is larger than
trust_radius * A
, then the whole adjustment needs to be scaled down until it is no larger thantrust_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.