PRO calculate_least_squares_matrix_and_vector, patt_arr, weights, data, matrix, vector, status, danidl_cppcode, POWERS1 = powers1, POWERS2 = powers2, NO_VECTOR = no_vector
; Description: This module calculates the least squares matrix and vector for a general linear
; least-squares problem, which is the problem of fitting a linear combination of a set
; of patterns to a set of data. Mathematically, the 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 corresponding
; normal equations are represented by the following matrix equation:
;
; A x = v
;
; where A is a square symmetric matrix (the least squares or Hessian matrix), x is the
; vector of scale factors to be determined, and v is a vector. The solution of this
; system is achieved by inverting the matrix A:
;
; x = A^(-1) v
;
; 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).
; This module calculates the least squares matrix A from the set of patterns "patt_arr"
; and the set of data weights "weights", and returns the matrix A via the parameter
; "matrix". The module also calculates the least squares vector v from the set of patterns
; "patt_arr", the set of data weights "weights" and the set of data values "data", and
; returns the vector v via the parameter "vector". Note that the values of the elements
; of the least squares matrix do not depend on the actual data values themselves.
; 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.
; The module provides the option of not calculating the least squares vector v by use
; of the keyword NO_VECTOR. If this keyword is set appropriately, then the input parameter
; "data" will be ignored, and the output parameter "vector" will not be calculated. This
; feature is useful if the user only requires the calculation of the least squares matrix.
; 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.5 times faster than the IDL
; code in this case.
;
; Input Parameters:
;
; patt_arr - FLOAT/DOUBLE VECTOR/ARRAY - An array, with at least two dimensions, containing the
; set of patterns P_1, P_2, ..., P_N to be used to calculate
; the least squares matrix and vector. 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 such that each
; sub-array "patt_arr[j,*]" has a total of M elements. Note
; that since IDL treats an N by 1 array as an N-element
; one-dimensional vector, this module accepts such an
; N-element one-dimensional vector to represent the case when
; M is equal to "1".
; weights - FLOAT/DOUBLE SCALAR/VECTOR/ARRAY - A scalar/vector/array of data weights with M elements.
; If this parameter does not have the correct number
; type or M elements, then all the weights are assumed
; to be "1.0".
; data - FLOAT/DOUBLE SCALAR/VECTOR/ARRAY - A scalar/vector/array of data values with M elements.
; If the keyword NO_VECTOR is set (as "/NO_VECTOR"), then
; this parameter is ignored.
; 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:
;
; matrix - DOUBLE VECTOR/ARRAY - A two-dimensional array of size NxN elements that represents the least
; squares matrix A. ;;;;;;;;Can come back as vector if N = 1 right? CHECK
; vector - DOUBLE VECTOR - A one-dimensional vector of length N elements that represents the least
; squares vector v. If the keyword NO_VECTOR is set (as "/NO_VECTOR"), then
; this vector is not calculated, and this parameter is set to "[0.0D]".
; status - INTEGER - If the module successfully calculated the least squares matrix and vector, then
; "status" is returned with a value of "1", otherwise it is returned with a value
; of "0".
;
; Keywords:
;
; 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 NO_VECTOR is set (as "/NO_VECTOR"), then the input parameter "data" will be ignored,
; and the least squares vector "vector" will not be calculated. This option is useful if the user
; only requires the calculation of the least squares matrix.
;
; Author: Dan Bramich (dan.bramich@hotmail.co.uk)
;
; History:
;
; 01/04/2011 - Module created (dmb).
;Set the default output parameter values
;Check that "patt_arr" is of the correct number type, and determine the number of patterns
;If "weights" does not have M elements or the correct number type, then assume that all the data weights have
;the value "1.0"
;If the least squares vector is to be calculated, then check that "data" is of the correct number type and
;that it has the correct number of elements
;If there is only one pattern
;Calculate the least squares matrix
;Calculate the least squares vector if not instructed otherwise
;Set "status" to "1" and return
;Determine if the shared library of DanIDL C++ routines is installed
;Reform the input arrays
;Determine if the user has indicated that the set of patterns are formed from powers of one or two underlying
;patterns
;If the set of patterns are not formed from powers of one or two underlying patterns
;Set up a cross reference matrix
;If the set of patterns are formed from powers of one underlying pattern
;Set up a cross reference matrix
;If the set of patterns are formed from powers of two underlying patterns
;Set up a cross reference matrix
;Set up the least squares matrix and vector
;Define the subscripts of the first pattern in the pattern array
;For each matrix column
;If the diagonal element in the current matrix column is still to be calculated
;Extract the pattern corresponding to the current matrix column
;Calculate a temporary vector
;Calculate the diagonal matrix element
;Fill in all the matrix elements that have the same value as the current matrix element
;Calculate the current element of the least squares vector if not instructed otherwise
;Determine how many elements in the current matrix column are still to be calculated
;If the diagonal element in the current matrix column has already been calculated
;Determine how many elements in the current matrix column are still to be calculated
;If there are no more elements in the current matrix column still to be calculated
;Calculate the current element of the least squares vector if not instructed otherwise
;Move on to the next matrix column
;If there is at least one element in the current matrix column still to be calculated
;Calculate a temporary vector
;Calculate the current element of the least squares vector if not instructed otherwise
;If the shared library of DanIDL C++ routines is installed, then use the C++ code to speed up the calculation
;of the remaining matrix elements in the current matrix column
;For each matrix element in the current matrix column that is still to be calculated
;Check that the current matrix element is still to be calculated
;Calculate the current matrix element
;Fill in all the matrix elements that have the same value as the current matrix element
;If the shared library of DanIDL C++ routines is not installed, then use the default IDL code
;For each matrix element in the current matrix column that is still to be calculated
;Check that the current matrix element is still to be calculated
;Calculate the current matrix element
;Fill in all the matrix elements that have the same value as the current matrix element
;Set "status" to "1"