PRO scale_pattern, data, weights, patt_arr, fsol, fsol_err, model, chisq, status, danidl_cppcode, CONSTANT = constant, ERRBAR = errbar, $
; Description: This module fits a linear combination of a set of patterns, stored in the pattern array
; "patt_arr", to a set of data "data". Mathematically, the data "data" are modelled as:
;
; M[i] = f_1*P_1[i] + f2*P_2[i] + ... + f_N*P_N[i]
;
; where M[i] is the model for the ith data point D[i], P_j represents the jth pattern,
; f_j is the jth scale factor, and N is the number of patterns. The model can also be
; adjusted to include an additive constant C by setting the keyword CONSTANT. If required,
; the fit can be weighted via the parameter "weights". The scale factors (and constant,
; if requested) are returned to the user via the solution vector "fsol" along with the
; model M via the parameter "model".
; In the case that the weights of the fit "weights" represent the inverse variance
; weights for the dependent data "data", then the uncertainties on the fitted parameters
; are represented by "fsol_err", and the chi squared of the fit is represented by "chisq".
; The module also provides the option of treating the weights as representing the
; 1-sigma uncertainties on the data values "data", in which case the actual weights that
; are used in the fit are the inverse variance weights, and again the uncertainties on
; the fitted parameters are represented by "fsol_err", and the chi squared of the fit is
; represented by "chisq".
; The problem described above is that of a general linear least-squares problem. The
; method used to solve for the scale factors "fsol" is that of the construction and
; solution of the normal equations, which are represented by the following matrix equation:
;
; A x = v
;
; where A is a square symmetric positive-definite matrix (the Hessian matrix), x is the
; vector of scale factors to be determined, and v is a vector. The solution proceeds
; via Cholesky factorisation of A = G transpose(G), and then solving G y = v and
; transpose(G) x = y in turn (i.e. solution via forwards and back substitution). This
; method is more efficient and numerically stable than simply inverting the matrix A
; and calculating x = A^(-1) v. The module returns the least-squares matrix A via the
; keyword MATRIX and the least-squares vector v via the keyword VECTOR.
; The inverse of A happens to be the covariance matrix. Specifically, in the case that
; the fit is weighted by the inverse variances, the square root of the diagonal elements
; of the covariance matrix represent the uncertainties on the fitted parameters. Hence
; the module will also calculate the inverse of A, independently of the factorisation
; method used to solve the equation A x = v, in order to obtain the covariance matrix and
; the uncertainties on the fitted parameters "fsol_err". The covariance matrix is returned
; via the keyword COV_MATRIX. The calculation of the covariance matrix and the parameter
; uncertainties is optional and may be controlled via the use of the keyword NO_INVERSION.
; For more details on the theory of the general linear least-squares problem, please
; see Chapter 15.4 of Numerical Recipes 3rd Edition (Press et al. 2007).
; The module also implements an important optimisation in the calculation of the least
; squares matrix for N > 1. This optimisation is applicable if the set of patterns have
; the functional form P_j = P^r_j where P is an underlying pattern. In this case, there
; is the possibility that some of the least squares matrix elements are repeated and
; therefore less calculations are required in order to construct the full matrix. An
; example of such a linear least-squares problem is that of fitting a polynomial in one
; dimension. This feature is accessible through the keyword POWERS1 which may be used to
; supply the set of powers r_j.
; It may also be the case that the set of patterns have the functional form
; P_j = P^r_j Q^s_j where P and Q are underlying patterns, in which case a similar
; optimisation is applicable in the calculation of the least squares matrix for N > 1.
; An example of such a linear least-squares problem is that of fitting a polynomial in
; two dimensions. This feature is accessible through the use of both of the keywords
; POWERS1 and POWERS2 used to supply the sets of powers r_j and s_j respectively.
; For further flexibility, the module provides the option of iterating the fitting
; process. By setting the keyword ITERATION_PAR to an IDL structure that contains the
; six tags "method", "sigma", "maxrej", "permanent", "tol" and "maxiter" defined below
; (see under "Keywords"), the module will do the following.
; If "sigma" is positive, then the sigma-clip threshold S is taken to be the value
; of "sigma". If "sigma" is negative or zero, then the sigma-clip threshold S is
; calculated as sqrt(ln(Ngood)), where Ngood is the number of data values in "data" with
; positive weights. This expression for the sigma-clip threshold S may be derived from
; considerations of the Bayesian information criterion (BIC).
; If the tag "method" is set to "abs_resid", then, during each iteration of the fitting
; process, all data values with absolute residuals greater than S times the standard
; deviation of the residuals are flagged as candidate outliers. If the tag "method" is
; set to "norm_abs_resid", then, during each iteration of the fitting process, all data
; values with normalised absolute residuals greater than S, where the normalising factors
; are the square roots of the corresponding weights, are flagged as candidate outliers.
; If the tag "maxrej" is set to an INTEGER/LONG type positive number, then this number
; is considered to be the maximum number of data values that can be rejected in any one
; iteration, and consequently, only this number of candidate outliers with the largest
; residuals are rejected. If the tag "maxrej" is set to a FLOAT/DOUBLE type positive
; number that is less than or equal to "1.0", then this number is considered to be the
; maximum fraction of data values that can be rejected in any one iteration, and
; consequently, only the corresponding number of candidate outliers with the largest
; residuals are rejected. If the tag "maxrej" is set to anything else, then there is no
; limit on the number of data values that can be rejected in any one iteration.
; If the tag "permanent" is set to "1", then data values that are rejected will remain
; rejected for all subsequent iterations, meaning that the set of data values used at each
; iteration will only diminish. If the tag "permanent" is set to any other value, then
; data values will not be permanently rejected so that if a data value is initially
; rejected, then it may be subsequently included in the fitting process during the later
; iterations. This protects against non-outlier data values being permanently excluded
; from the calculation of the fit parameters when real outlier data values skew the
; parameter estimates during the early iterations.
; The iterations stop when the maximum fractional change in the fit parameters is less
; than or equal to "tol", or when there are not enough good data values to be used in the
; fit. To protect against the possibility of the fitting process not converging (or
; entering into an oscillatory regime), a limit of "maxiter" iterations is placed on the
; number of iterations that are performed.
; Useful extra information about the fitting process is returned via the keyword
; EXTRA_INFO, and the subscripts of the data values that were used in the (final iteration
; of the) fit are returned via the keyword DATA_USED_SUBS.
; This module will use C++ code if possible to achieve a faster execution than the
; default IDL code. The C++ code can be up to a factor of ~1.20 times faster than the IDL
; code in this case.
;
; N.B: In this module, the definition of the standard deviation of the residuals that is
; used is given by:
;
; Stddev(R) = sqrt( [ SUM_i ( R_i )^2 ] / (N_d - N_p) )
;
; where R is the set of residuals R_i = x_i - Model(X), X is the set of data with
; N_d elements, x_i is the ith element of X, Model(X) is the set of model values, and
; N_p is the number of fit parameters. The above definition of the standard deviation
; is employed because the parameters are estimated from the data itself.
;
; Input Parameters:
;
; data - FLOAT/DOUBLE SCALAR/VECTOR/ARRAY - A scalar/vector/array of input data to which the linear
; combination of patterns are to be fitted.
; weights - FLOAT/DOUBLE SCALAR/VECTOR/ARRAY - A scalar/vector/array with the same number of elements
; as "data", that contains the weights to be used when
; fitting the linear combination of patterns. If this
; parameter does not have the correct number type or the
; same number of elements as "data", then all the data
; values are weighted equally. Data points to be ignored
; can be pre-flagged by setting the corresponding weights
; in this scalar/vector/array to zero (or negative values).
; If the keyword ERRBAR is set, then the module will
; assume that the values in this parameter represent the
; 1-sigma uncertainties on the data values, and it will
; calculate the actual weights to be used in the fit as
; the inverse variance weights.
; patt_arr - FLOAT/DOUBLE SCALAR/VECTOR/ARRAY - An array, (usually) with at least two dimensions,
; containing the set of patterns P_1, P_2, ..., P_N to be
; fitted as a linear combination to the input data "data".
; The first dimension of "patt_arr" should always be of
; size N, which is the number of patterns to be fitted.
; The remaining dimensions of "patt_arr" may be of any
; configuration with the condition that each sub-array
; "patt_arr[j,*]" (or pattern) contains the same number
; of elements as "data". In the case that N = 1 and there
; is only one data value in "data", then "patt_arr" may
; be supplied as a scalar value or a single element
; vector. To ensure a stable solution, none of the
; patterns in "patt_arr" should be identical, or a
; linear combination of any other set of patterns. The
; module will fail if any one pattern consists solely of
; zeros for the data points with positive weights.
; danidl_cppcode - STRING - The full directory path indicating where the DanIDL C++ code is installed.
; The shared library of DanIDL C++ routines should exist as
; "/dist/libDanIDL.so" within the installation.
;
; Output Parameters:
;
; fsol - DOUBLE VECTOR - A vector of scale factors of length N elements representing the solutions f_1,
; f_2, ..., f_N. If the keyword CONSTANT is set, then this vector will actually
; be of length N+1, and the last element represents the fitted additive constant
; C.
; fsol_err - DOUBLE VECTOR - A vector of the same length as "fsol" with values representing the
; uncertainties on the values stored in "fsol". The values in this vector
; only represent uncertainties in the case where the weights "weights" are
; the inverse variances for the dependent data "data", or if the keyword
; ERRBAR is set (see below) and the weights represent the 1-sigma
; uncertainties on the dependent data "data". If the keyword NO_INVERSION
; is set and more than one pattern is being fitted, then the covariance
; matrix will not be calculated and therefore the uncertainties on the
; fitted parameters will also not be calculated. In this case "fsol_err" is
; returned as a vector of the same length as "fsol" where are all elements
; are set to the non-sensical value of "-1.0".
; model - DOUBLE SCALAR/VECTOR/ARRAY - A scalar/vector/array (of the same dimensions as the input
; parameter "data") representing the fitted model M.
; chisq - DOUBLE - The chi squared of the fit ignoring bad data values flagged in "weights" and rejected
; data values identified during the iterative fitting process. This quantity only
; represents the chi squared in the case where the input weights "weights" represent
; the inverse variances for the dependent data "data", or if the keyword ERRBAR is set
; (see below) and the weights represent the 1-sigma uncertainties on the dependent data
; "data".
; status - INTEGER - If the fit parameters were calculated successfully and no iterations were performed,
; then "status" is returned with a value of "4". If the fit parameters were calculated
; successfully and the iterations were stopped because the algorithm tolerance was
; achieved, then "status" is returned with a value of "3". If the fit parameters were
; calculated successfully but the iterations were stopped because the maximum number
; of iterations "maxiter" had been performed, then "status" is returned with a value
; of "2". If the fit parameters were calculated successfully but the iterations were
; stopped because there were not enough good data values to be used in the next
; iteration, then "status" is returned with a value of "1". If the fit parameters were
; not calculated successfully, then "status" is returned with a value of "0".
;
; Keywords:
;
; If the keyword CONSTANT is set (as "/CONSTANT"), then the module will include an additive constant C
; in the model such that:
;
; M[i] = f_1*P_1[i] + f2*P_2[i] + ... + f_N*P_N[i] + C
;
; The fitted value of C and its uncertainty will be returned as the last elements of the vectors "fsol"
; and "fsol_err" respectively.
;
; If the keyword ERRBAR is set (as "/ERRBAR"), then the module will assume that the weights "weights"
; represent the 1-sigma uncertainties on the data values "data", and it will calculate the actual
; weights to be used in the fit as the inverse variance weights.
;
; If the keyword NO_INVERSION is set (as "/NO_INVERSION") and more than one pattern is being fitted,
; then the module will not calculate the covariance matrix (which is the inverse of the least squares
; matrix A). Consequently the covariance matrix is not returned via the keyword COV_MATRIX and the
; parameter uncertainties (which are given by the square root of the diagonal elements of the covariance
; matrix) are not calculated. This option is useful for avoiding the relatively costly calculation of
; the inverse of A when the user is not interested in the covariance matrix or the parameter
; uncertainties, and it can therefore be used to build more efficient code.
;
; If the keyword POWERS1 is set to a one-dimensional vector r_j of INTEGER/LONG type numbers with N
; elements, and the keyword POWERS2 is not set to a one-dimensional vector of INTEGER/LONG type
; numbers with N elements, then the module will assume that the set of patterns have the functional
; form P_j = P^r_j where P is an underlying pattern. This assumption will then be used to speed up
; the calculation of the least squares matrix by avoiding redundant calculations. Note that if N = 1
; then this keyword is ignored.
;
; If both of the keywords POWERS1 and POWERS2 are each set to one-dimensional vectors r_j and s_j of
; INTEGER/LONG type numbers each with N elements, then the module will assume that the set of patterns
; have the functional form P_j = P^r_j Q^s_j where P and Q are underlying patterns. This assumption
; will then be used to speed up the calculation of the least squares matrix by avoiding redundant
; calculations. Note that if N = 1 then this keyword is ignored.
;
; If the keyword MATRIX is set to a named variable, then, on return, this variable will contain a
; two-dimensional array of size NxN elements that represents the matrix A in the normal equations (see
; description above). If the keyword CONSTANT is set, then this matrix will actually be of size
; (N+1)x(N+1) elements. ;;;;;;;;;WHEN N = 1????
;
; If the keyword VECTOR is set to a named variable, then, on return, this variable will contain a
; one-dimensional vector of length N elements that represents the vector v in the normal equations (see
; description above). If the keyword CONSTANT is set, then this vector will actually be of length
; (N+1) elements.
;
; If the keyword COV_MATRIX is set to a named variable, then, on return, this variable will contain a
; two-dimensional array of size NxN elements that represents the covariance matrix of the solution. If
; the keyword CONSTANT is set, then this matrix will actually be of size (N+1)x(N+1) elements. If the
; keyword NO_INVERSION is set and more than one pattern is being fitted, then the covariance matrix will
; not be calculated.
;
; If the keyword ITERATION_PAR is set to an IDL structure with the six tags described below, then the
; module will iterate the fitting process (as decribed under "Description"). The tags that must be
; present are:
;
; (i) method - STRING - This tag has two acceptable values. If "method" is set to "abs_resid", then,
; during each iteration of the fitting process, all data values with absolute
; residuals greater than S times (see "sigma") the standard deviation of the
; residuals are flagged as candidate outliers. If "method" is set to
; "norm_abs_resid", then, during each iteration of the fitting process, all
; data values with normalised absolute residuals greater than S, where the
; normalising factors are the square roots of the corresponding weights, are
; flagged as candidate outliers.
;
; (ii) sigma - FLOAT/DOUBLE - This tag specifies the sigma-clip rejection threshold to be used. If
; "sigma" is positive, then the sigma-clip threshold S is taken to be
; the value of "sigma". If "sigma" is negative or zero, then the
; sigma-clip threshold S is calculated as sqrt(ln(Ngood)), where Ngood
; is the number of data values in "data" with positive weights. This
; expression for the sigma-clip threshold S may be derived from
; considerations of the Bayesian information criterion (BIC).
;
; (iii) maxrej - INTEGER/LONG or FLOAT/DOUBLE - If this tag is set to an INTEGER/LONG type positive
; number, then the value of this tag is considered to
; be the maximum number of data values that can be
; rejected in any one iteration. If this tag is set
; to a FLOAT/DOUBLE type positive number that is less
; than or equal to "1.0", then the value of this tag
; is considered to be the maximum fraction of data
; values that can be rejected in any one iteration.
; If this tag is set to anything else, then there is
; no limit on the number of data values that can be
; rejected in any one iteration.
;
; (iv) permanent - INTEGER/LONG - If this tag is set to "1", then data values that are rejected will
; remain rejected for all subsequent iterations, meaning that the
; set of data values used at each iteration will only diminish. If
; this tag is set to any other value, then data values will not be
; permanently rejected so that if a data value is initially excluded
; then it may be subsequently included in the fitting process during
; the later iterations.
;
; (v) tol - FLOAT/DOUBLE - This tag specifies the tolerance of the iterative fitting process. The
; iterations will stop if the maximum fractional change in the fit parameters
; between subsequent iterations is less than or equal to "tol". This parameter
; must have a value greater than or equal to "0.0", and less than "1.0". A
; value for "tol" of "0.0" means that the iterations will continue until none
; of the fit parameters change from one iteration to the next.
;
; (vi) maxiter - INTEGER/LONG - This tag specifies the maximum number of iterations to be performed
; before returning the fit parameters. If this parameter is set to a
; negative value or zero, then there will be no constraint on the
; maximum number of iterations that may be performed.
;
; If the keyword EXTRA_INFO is set to a named variable, then, on return, this variable will contain an
; IDL structure with the following tags:
;
; (i) ndata - LONG - The number of data values in the input parameter "data".
; (ii) ndata_good - LONG - The initial number of data values with positive weights.
; (iii) ndata_used - LONG - The number of data values used in the calculation of the final fit
; parameters.
; (iv) niter - LONG - The number of iterations performed in the fitting process.
;
; If the keyword DATA_USED_SUBS is set to a named variable, then, on return, this variable will contain
; a VECTOR of LONG type numbers representing the subscripts for the data values that were used in the
; calculation of the final fit parameters.
;
; Author: Dan Bramich (dan.bramich@hotmail.co.uk)
;
; History:
;
; 28/10/2011 - Added the keyword NO_INVERSION (dmb).
; 27/10/2011 - Implemented Cholesky factorisation with forwards and back subsititution in order to solve
; the normal equations (dmb).
; 06/04/2011 - Keywords POWERS1 and POWERS2 added (dmb).
; 31/01/2011 - Module expanded to allow iteration of the fitting process (dmb).
; 02/12/2010 - Module expanded to cope with input data of any dimensions (dmb).
; 25/09/2010 - Module fully rewritten to take advantage of the matrix multiplication operations within
; IDL. A great improvement in speed and efficiency of the code was achieved (dmb).
; 16/07/2010 - Module created (dmb).
;Set the default output parameter values
;Check that "data" is of the correct number type
;If "weights" does not have the same number of elements as "data" or the correct number type, then weight all
;data points equally, otherwise set any negative weights equal to zero while checking that not all data points
;are flagged as bad. Also, if the keyword ERRBAR is set, then convert the 1-sigma uncertainties on the data
;values to inverse variances.
;Check that "patt_arr" is of the correct number type
;If there is only one data value
;Check that the keyword CONSTANT is not set
;Check that "patt_arr" has only one element
;Calculate the least squares matrix and the vector on the right hand side
;Check that the least squares matrix element is not zero
;Compute the inverse of the least squares matrix
;Calculate the scale factor and its uncertainty
;Construct the model for the input data
;Determine the subscripts of the data values used in the fit
;Set "status" to "4"
;If there is more than one data value, then check that "patt_arr" has at least two dimensions. Also, if the
;keyword CONSTANT is set, then include the required pattern to enable the fit of the constant parameter.
;Check that the number of data points to be fitted that have positive weights is at least as great as the
;number of parameters to be determined
;Determine if the user has indicated that the set of patterns are formed from powers of one or two underlying
;patterns, and if the keyword CONSTANT is set, then include the zeroth power in the relevant power vector(s)
;Determine if the fitting process should be iterated
;Determine the subscripts of the diagonal elements of the least squares matrix
;If the keyword NO_INVERSION is not set, then prepare an identity matrix
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;Iterative fitting loop
;Start the iterative fitting loop
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;Calculate the least squares matrix and vector
;If this is the first iteration
;Calculate the least squares matrix and the vector on the right hand side
;If the iteration schema requires it, then save the original least squares matrix and the vector on
;the right hand side
;If this is not the first iteration
;Record the previous set of fit parameters
;If data values are not being permanently rejected, then reset the least squares matrix and the vector
;on the right hand side
;In this case, it is faster to subtract the contributions of the rejected data values from the least
;squares matrix and the vector on the right hand side, rather than recalculating everything from scratch
;Check that there are no zero diagonal elements in the least squares matrix
;;;;;;;;;;;;;;;;;;;;;;;Don't generate subscripts here. Count instead using long(total())
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;Solve the normal equations
;If there is only one pattern to fit
;Compute the inverse of the least squares matrix
;Solve the system of linear equations
;Determine the uncertainties on the fitted parameters
;If there are multiple patterns to fit
;Compute the Cholesky factorisation of the least squares matrix
;Solve the system of linear equations
;If the keyword NO_INVERSION is not set
;Solve an identity system of linear equations in order to compute the inverse of the least squares matrix
;Determine the uncertainties on the fitted parameters
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;Calculate the model and residuals
;Construct the current model for the input data
;Calculate the squared residuals of the current fit
;Calculate the chi squared of the current fit
;Update the number of iterations performed
;If no iterations are to be performed, then exit from the iterative fitting loop
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;Iterative schema
;Save a copy of the current weights
;If this is not the first iteration
;Calculate the maximum fractional change in the fit parameters
;If the maximum fractional change in the fit parameters is less than or equal to "tol", then
;exit the iterative fitting loop
;If the maximum number of iterations has been reached, then exit the iterative fitting loop
;If iteration is being performed based on the absolute residuals of the fit
;Only perform data point rejection if there are more data points than fit parameters
;Estimate the variance of the residuals
;If data values are not to be rejected permanently, then reset the current weights to the original weights
;Determine a set of candidate outlier data values
;If iteration is being performed based on the normalised absolute residuals of the fit
;If data values are not to be rejected permanently, then reset the current weights to the
;original weights
;Determine a set of candidate outlier data values
;If no (further) data values are to be rejected, then exit the iterative fitting loop
;If a limit has been set on the number of data points that may be rejected in any one iteration,
;then determine what this limit is for the current iteration
;If not all candidate outlier data values are to be rejected in this iteration
;Determine the required set of candidate outlier data values with the largest residuals
;Reject the outlier data values
;Determine the new number of good data values
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;Set the output parameters
;Set the output parameter values
;If the fitting process was iterated
;Set the output parameter values in the IDL structure "extra_info"
;Determine the subscripts of the data values used in the final fit
;If the iterative fitting loop was exited because the fit tolerance was achieved, then set "status"
;to "3" and finish
;If the iterative fitting loop was exited because the maximum number of iterations was reached, then
;set "status" to "2" and finish
;If the iterative fitting loop was exited because there were not enough good data values to be used
;in the next iteration, then set "status" to "1"
;If the fitting process was not iterated
;Set the output parameter values in the IDL structure "extra_info"
;Determine the subscripts of the data values used in the fit
;Set "status" to "4"