NAME

PDL::PP::Include::Fit::ExpRate - build-time 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::PP-based 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 pp-functions when you say 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 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 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, ...
);

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

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.

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]  )
                    \                                  /

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);

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 trust_radius * A, then the whole adjustment needs to be scaled down until it is no larger than 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.