Creating Non-linear Fit Plugins

Kst provides header files designed to simplify the creation of non-linear least-squares fit plugins. The following sections detail the use of the header files.

Header Files and Definitions

The non-linear fit header files are located in kst/plugins/fits_nonlinear/ of the Kst source tarball. The files are named non_linear.h and non_linear_weighted.h for unweighted and weighted fits, respectively. To use these files, include only one of them in the source code for your plugin:

#include <../non_linear.h>

#include <../non_linear_weighted.h>

(by convention, we will place the source code for the plugin one directory below where the header files are).

As non-linear fitting is an iterative process, you must also define the maximum number of iterations that should be performed. The non-linear fitting algorithm will stop when at least one of the following conditions is true:

The maximum number of iterations has been reached.
A precision of 10^-4 has been reached.

In addition, you need to define the number of parameters in the model, as it is not passed to the fitting function explicitly. To define these two values, include the following at the top of your source code:

#define NUM_PARAMS [num1]
#define MAX_NUM_ITERATIONS [num2]

replacing [num1] with the number of parameters in the model, and [num2] with the maximum number of iterations to perform.

Implementing Required Functions

To use the header files for non-linear fits, you must provide the function to use as the model, the partial derivatives of the function with respect to each parameter, and initial estimates of the best-fit parameters. To do this, three functions must be implemented. These functions are described below.

double function_calculate( double dX, double* pdParameters)

This function calculates the y value of the fitting model for a given x value dX, using the supplied array of parameters pdParameters. The order of parameters in pdParameters is arbitrary, but should be consistent with the other two implemented functions. For example, for an exponential model, function_calculate could be implemented as follows:

double function_calculate( double dX, double* pdParameters ) {
  double dScale	 = pdParameters[0];
  double dLambda = pdParameters[1];
  double dOffset = pdParameters[2];
  double dY;

  dY  = ( dScale * exp( -dLambda * dX ) ) + dOffset;

  return dY;
}

void function_derivative( double dX, double* pdParameters, double* pdDerivatives )

This function calculates the partial derivatives of the model function for a give value of x dX. The partial derivatives should be returned in pdDerivatives. The order of the partial derivatives in the array pdDerivatives should correspond to the order of the parameters in pdParameters (i.e. if pdParameters[0] contains the parameter lambda for an exponential model, pdDerivatives[0] should contain the derivative of the model with respect to lambda).

void function_initial_estimate( const double* pdX, const double* pdY, int iLength, double* pdParameterEstimates )

This function provides an initial estimate of the best-fit parameters to the fitting function. The array of x values and y values of the data points are provided in pdX and pdY respectively, and the number of data points is provided by iLength. You can use any or none of these parameters at your discretion to calculate the initial estimate. The function should put the calculated initial estimates in pdParameterEstimates, with the order of the estimates corresponding to the order of the parameters in pdParameters of function_calculate and function_derivative. Keep in mind that the initial estimate is important in determining whether or not the fitting function converges to a solution.

Calling the Fitting Functions

Once all the required functions have been implemented, the fitting function from the included header file can be called:

kstfit_nonlinear( inArrays, inArrayLens,
		  inScalars, outArrays,
		  outArrayLens, outScalars );

kstfit_nonlinear_weighted( inArrays, inArrayLens,
                           inScalars, outArrays,
                           outArrayLens, outScalars );

depending on whether you are implementing a non-weighted fit or a weighted fit.

The function will return 0 on success, or -1 on error, so it is simplest to set the return value of the exported C function to be equal to the return value of the fitting function. To maintain simplicity, the code for the plugin can simply pass the arguments given to the exported C function to the fitting function. Note, however, that inArrays must be structured as follows:

inArrays[0] must contain the array of x coordinates of the data points
inArrays[1] must contain the array of y coordinates of the data points
inArrays[2] only exists if kstfit_linear_weighted is being called, and must contain the array of weights to use for the fit.

The easiest way to ensure that inArrays is structured correctly is to specify the correct order of input vectors in the XML file for the plugin.

After kstfit_linear_unweighted or kstfit_linear_weighted is called, outArrays and outScalars will be set as follows:

outArrays[0] will contain the array of fitted y values.
outArrays[1] will contain the array of residuals.
outArrays[2] will contain the array of best-fit parameters that were estimated.
outArrays[3] will contain the covariance matrix, returned row after row in an array.
outScalars[0] will contain chi^2/nu, where chi^2 is the weighted sum of squares of the residuals, and nu is the degrees of freedom.

outArrayLens will be correctly set to indicate the length of each output array.

Ensure that the specified outputs in the XML file match those that the exported C function returns (which in most cases will simply be the outputs returned by the fitting function).

Example

The following is an example of a non-linear fit plugin that performs a fit to an exponential model.

/*
 *  Exponential fit plugin for KST.
 *  Copyright 2004, The University of British Columbia
 *  Released under the terms of the GPL.
 */

#define NUM_PARAMS 3
#define MAX_NUM_ITERATIONS 500

#include "../non_linear.h"

void function_initial_estimate( const double* pdX, const double* pdY,
                                int iLength, double* pdParameterEstimates ) {
  KST_UNUSED( pdX )
  KST_UNUSED( pdY )
  KST_UNUSED( iLength )

  pdParameterEstimates[0] =  1.0;
  pdParameterEstimates[1] =  0.0;
  pdParameterEstimates[2] =  0.0;
}

double function_calculate( double dX, double* pdParameters ) {
  double dScale	 = pdParameters[0];
  double dLambda = pdParameters[1];
  double dOffset = pdParameters[2];
  double dY;

  dY  = ( dScale * exp( -dLambda * dX ) ) + dOffset;

  return dY;
}

void function_derivative( double dX, double* pdParameters, double* pdDerivatives ) {
  double dScale	 = pdParameters[0];
  double dLambda = pdParameters[1];
  double dExp;
  double ddScale;
  double ddLambda;
  double ddOffset;

  dExp     = exp( -dLambda * dX );
  ddScale  = dExp;
  ddLambda = -dX * dScale * dExp;
  ddOffset = 1.0;

  pdDerivatives[0] = ddScale;
  pdDerivatives[1] = ddLambda;
  pdDerivatives[2] = ddOffset;
}

extern "C" int kstfit_exponential(const double *const inArrays[], const int inArrayLens[],
		const double inScalars[],
		double *outArrays[], int outArrayLens[],
		double outScalars[]);

int kstfit_exponential(const double *const inArrays[], const int inArrayLens[],
		const double inScalars[],
		double *outArrays[], int outArrayLens[],
		double outScalars[])
{
  return kstfit_nonlinear( inArrays, inArrayLens,
                           inScalars, outArrays,
			   outArrayLens, outScalars );
}