
MatrixLibContents1. Introduction1.1 Matrix Data Types 1.2 Prefixes of MatrixLib Functions 1.3 C/C++ Specifics 1.3.1 MatObj, the objectoriented interface for matrix functions 1.3.2 Directpointer form of matrix function calls 1.4 Pascal/Delphi Specifics 2. Management of Dynamically Generated Matrices 3. Initialization of Matrices 4. DataType Conversions 5. Symmetry Operations (Transposing, Rotating, Reflecting), Interpolation, Augmenting, Deleting, Extracting, and Filling Parts of a Matrix 6. Arithmetic Operations Performed on a Single Row, Column, or the Diagonal 7. Operations Performed Along All Rows or All Columns Simultaneously, or the Diagonal; Center of Gravity of a Matrix 8. Operations Involving Two Rows or Two Colums 9. WholeMatrix Arithmetics: Addition, Multiplication 10. Linear Algebra 11. Eigenvalues and Eigenvectors, Matrix SquareRoot 12. TwoDimensional FourierTransform Methods 13. Data Fitting 13.1 Polynomials 13.2 General Linear Model Functions 13.3 NonLinear Models 13.4 Fitting Multiple Data Sets 13.5 Helper Functions for Nonlinear Fits 13.6 Summary of Fitting Functions 14. Matrix Input and Output 15. Graphical Representation of Matrices 1. IntroductionThe MatrixLib part of OptiVec builds upon the principles established in the VectorLib part. You may wish to refer to the introductory chapters of the VectorLib documentation, before reading on for MatrixLib.Back to MatrixLib Table of Contents OptiVec home 1.1 Matrix Data TypesAs VectorLib defines vector data types, here are the matrix data types, defined in MatrixLib:
The ordering of elements is the same as in the twodimensional arrays provided by the respective target compilers. This means that the matrices are stored rowwise in MatrixLib versions for C and C++ compilers, but columnwise in the Pascal/Delphi versions. While we recommend to exclusively use these dynamically allocated matrix types, static matrices defined, e.g., as
Back to MatrixLib Table of Contents OptiVec home 1.2 Prefixes of MatrixLib FunctionsEach MatrixLib function has a prefix defining the data type on which it acts:
In a few cases, the prefix is augmented by a threeletter code denoting special matrix properties:
Back to MatrixLib Table of Contents OptiVec home 1.3 C/C++ Version Specifics:1.3.1 MatObj, the objectoriented interface for matrix functionsSimilarly to the vector objects of VecObj, the objectoriented interface for the matrix functions is contained in the includefiles <fMatObj.h>, <dMatObj.h> etc., with one includefile for each of the datatypes supported in OptiVec. Remember that, in order to get the whole interface (for all data types at once),#include <OptiVec.h>. For access to any of the vector graphics functions, always include <OptiVec.h>. By analogy with the alias names of the vector objects, the matrix objects get the alias names fMatObj, dMatObj, and so on, with the datatype signalled by the first one or two letters of the class name. The matrix constructors are:
For the matrix classes, the arithmetic operators
If you ever need to process a MatObj matrix in a "classic" plainC OptiVec function, you may use the member functions The syntax of all MatObj functions is described below together with the basic MatrixLib functions for which tMatObj serves as a wrapper. Please note the following restrictions, coming from the tight control of matrix dimensions for compatibility:
1.3.2 Directpointer form of matrix function callsIn the C/C++ version of MatrixLib, all functions taking matrix arguments exist in a second form. In this form, all matrix arguments are replaced by pointers to the first elements of the respective matrices. You can explicitly use the second form by leaving away the underbar of the functionname prefix. CallingMF_mulM( MC, MA, MB, htA, lenA, lenB ); is equivalent to calling MFmulM(&(MC[0][0]),&(MA[0][0]),&(MB[0][0]),htA,lenA,lenB); or MFmulM( MC[0], MA[0], MB[0], htA, lenA, lenB ); Actually, the runtime routines are in the second form, and the macros defined in <MFstd.h> etc. convert calls to the first form into calls to the second form. Back to MatrixLib Table of Contents OptiVec home 1.4 Pascal/Delphi Version Specifics:In the Pascal/Delphi version of MatrixLib, matrices of all data types are actually defined as pointers to the element M[0][0]. This means you can pass static matrices (like MX: array[0..5][0..5] of Single; ) to MatrixLib functions with the address operator:MF_equ1( @(MX[0][0]), 6 ); The situation is somewhat different for dynamic arrays in Delphi. While, in the onedimensional case, they they can be used with all VectorLib functions by simply typecasting fArrays into fVectors, twodimensional arrays require a different approach: Dynamic twodimensional arrays can be created in Delphi by declaring them as "array of fArray", "array of dArray", etc. Shortcut definitions for these types are also given in the VecLib unit, as f2DArray, d2DArray, etc. As a consequence of the described way of defining 2DArrays in Delphi, each row of a matrix is stored in a separate vector. This means that the matrix elements do no longer occupy a single, contiguous memory space. Therefore, 2DArrays cannot be used directly with MatrixLib functions. Instead, they have to be copied into matrices by calling MF_2DArrayToMatrix etc., before MatrixLib functions can be used on them. The reverse conversion is available as MF_MatrixTo2DArray. In the following chapters, a brief summary of all MatrixLib function is given, ordered into groups according to their functionality. Click on the links to obtain more information about the individual functions mentioned. Back to MatrixLib Table of Contents OptiVec home 2. Management of Dynamically Generated Matrices
OptiVec's dynamically allocated matrices are aligned on 32byte boundaries, which allows for optimum cacheline matching for the Pentium processor and its currently available successors and competitors. C/C++ version only:
Both C/C++ and Pascal/Delphi versions:
To assign a value to a specific matrix element, you should use the syntax *(MF_Pelement( MX, ht, len, 3, 4 )) = 3.0; /* for C/C++*/ MF_Pelement( MX, ht, len, 3, 4 )^ := 3.0; (* for Pascal/Delphi *) Back to MatrixLib Table of Contents OptiVec home 3. Initialization of MatricesIn order to initialize all matrix elements with the same value, or to perform the same arithmetic operation on all matrix elements simultaneously, please use the respective VectorLib function. You can do so, because all matrices are in fact stored as vectors, occupying a contiguous space in memory. All you have to do is to pass the first row of the matrix (rather than the matrix itself) to the vector function. Fore example, you can initialize all matrix elements with a constant C by callingVF_equC( MA[0], ht*len, C ); /* C/C++ */ VF_equC( MA, ht*len, C ); (* Pascal/Delphi *) As you shall see, some of the most common operations of this kind are also explicitly defined in MatrixLib, like initialization with zero, MF_equ0, or multiplication by a constant, available as MF_mulC (see chapter 9). Here are the typical matrix initialization functions:
Twodimensional windows for spectral analysis are provided by:
Back to MatrixLib Table of Contents OptiVec home 4. DataType ConversionsMatrices of every data type can be converted into every other. Only a few examples are given; the rest should be obvious.
Back to MatrixLib Table of Contents OptiVec home 5. Symmetry Operations (Transposing, Rotating, Reflecting),

MF_transpose  transpose a matrix: MB = MA^{T} 
MCF_hermconj  Hermitian conjugate: MB = MA^{T*} 
MF_rotate90  clockwise rotation by 90° 
MF_rotate180  rotation by 180° 
MF_rotate270  clockwise rotation by 270° (or counterclockwise rotation by 90 °) 
MF_Rows_rev  reverse the element order along rows; this corresponds to a reflection of the matrix at the Y axis 
MF_Cols_rev  reverse the element order along columns; this corresponds to a reflection of the matrix at the X axis 
MF_Rows_reflect  set the upper halves of all rows equal to their reversed lower halves 
MF_Cols_reflect  set the upper halves of all columns equal to their reversed lower halves 
MF_UequL  copy lowerdiagonal elements into upperdiagonal by indexreflection, so as to get a symmetric matrix 
MF_LequU  copy upperdiagonal elements into lowerdiagonal 
MF_polyinterpol  polynomial interpolation 
MF_ratinterpol  rational interpolation 
MF_natCubSplineInterpol  natural cubic spline interpolation 
MF_equMblock  extract a block, i.e. a submatrix for which the sampling interval both along rows and along columns is 1 
MF_equMblockT  extract the transpose of a block 
MF_submatrix  extract a submatrix with sampling intervals along rows and/or columns possibly different from 1 
MF_block_equM  copy a matrix back into a block of another (normally larger) matrix 
MF_block_equMT  copy the transpose of a matrix into a block of another (normally larger) matrix 
MF_submatrix_equM  copy a submatrix back into another (normally larger) matrix 
MF_Row_extract  extract a single row and copy it into a vector 
MF_Col_extract  extract a single column and copy it into a vector 
MF_Dia_extract  extract the diagonal and copy it into a vector 
MF_Trd_extract  extract a tridiagonal matrix from a general matrix 
MF_Row_insert  augment a matrix by insertion of one row 
MF_Col_insert  augment a matrix by insertion of one column 
MF_Row_delete  delete one row of a matrix 
MF_Col_delete  delete one column of a matrix 
Back to MatrixLib Table of Contents OptiVec home
MF_Row_neg  multiply all elements of a specific row by 1 
MF_Col_neg  multiply all elements of a specific column by 1 
MF_Row_addC  add a constant to all elements of a specific row 
MF_Col_addC  add a constant to all elements of a specific column 
MF_Dia_addC  add a constant to all diagonal elements 
MF_Row_addV  add corresponding vector elements to all elements of a specific row 
MF_Col_addV  add corresponding vector elements to all elements of a specific column 
MF_Dia_addV  add corresponding vector elements to the diagonal elements 
MF_Row_subC  subtract a constant from all elements of a specific row 
MF_Col_subrC  reverse substraction: difference between column elements and a constant 
MF_Dia_mulV  multiply the diagonal elements with corresponding vector elements 
MF_Row_divV  divide all elements of a specific row by corresponding vector elements 
MF_Col_divrC  reverse division: division of a constant by the individual column elements 
Back to MatrixLib Table of Contents OptiVec home
MF_Rows_max  store the maxima of all rows in a column vector 
MF_Cols_max  store the maxima of all colums in a row vector 
MF_Dia_max  return the maximum of the diagonal as a scalar 
MF_Rows_absmax  store the absolute maxima of all rows in a column vector 
MF_Cols_absmax  store the absolute maxima of all colums in a row vector 
MF_Dia_absmax  return the absolute maximum of the diagonal as a scalar 
MF_Rows_absmin  store the absolute minima of all rows in a column vector 
MF_Cols_absmin  store the absolute minima of all colums in a row vector 
MF_Dia_absmin  return the absolute minimum of the diagonal as a scalar 
MF_Rows_sum  sum, taken along rows and stored in a column vector 
MF_Cols_sum  sum, taken along colums and stored in a row vector 
MF_Dia_sum  sum of the diagonal elements 
MF_Rows_prod  product, taken along rows and stored in a column vector 
MF_Cols_prod  product, taken along colums and stored in a row vector 
MF_Dia_prod  product of the diagonal elements 
MF_Rows_runsum  running sum along rows 
MF_Cols_runsum  running sum along columns 
MF_Rows_runprod  running product along rows 
MF_Cols_runprod  running product along columns 
MF_Rows_rotate  rotate all rows by a specified number of positions 
MF_Cols_rotate  rotate all columns by a specified number of positions 
MF_Rows_reflect  set the upper halves of all rows equal to their reversed lower halves 
MF_Cols_reflect  set the upper halves of all columns equal to their reversed lower halves 
Please note that multiplying all rows or all columns by one and the same vector is equivalent to a multiplication by a diagonal matrix, which is provided by MF_mulMdia and MF_TmulMdia.
For all of the above functions involving maxima (...max, ...absmax), the corresponding minima are found by the functions named ...min, ...absmin.
For complex numbers, maxima can be defined by various criteria. The following table summarizes the functions finding them (again, of course, the corresponding functions for the minima exist as well):
MCF_Rows_absmax  store the maxima of the absolute values of all rows in a (realvalued) column vector 
MCF_Cols_absmax  store the maxima of the absolute values of all colums in a (realvalued) row vector 
MCF_Dia_absmax  return the maximum of the absolute values of the diagonal elements as a (real) scalar 
MCF_Rows_absmaxReIm  find the maximum absolute values of all real and of all imaginary parts separately along the rows of a matrix; merge the maxima into complex numbers, and store them in a column vector 
MCF_Cols_absmaxReIm  find the maximum absolute values of all real and of all imaginary parts separately along the columns of a matrix; merge the maxima into complex numbers, and store them in a row vector 
MCF_Dia_absmaxReIm  find the maximum absolute values of all real and of all imaginary parts separately along the diagonal of a square matrix; merge these two maxima into one complex number and return it as a scalar 
MCF_Rows_cabsmax  find the complex numbers of largest magnitude along rows and store these row maxima in a (complex) column vector 
MCF_Cols_cabsmax  find the complex numbers of largest magnitude along columns and store these column maxima in a (complex) row vector 
MCF_Dia_cabsmax  find the complex number of largest magnitude along the diagonal of a square matrix and return it as a (complex) scalar 
MCF_Rows_maxReIm  find the maxima of all real and of all imaginary parts separately along the rows of a matrix; merge the maxima into complex numbers, and store them in a column vector 
MCF_Cols_maxReIm  find the maxima of all real and of all imaginary parts separately along the columns of a matrix; merge the maxima into complex numbers, and store them in a row vector 
MCF_Dia_maxReIm  find the maximum of all real and of all imaginary parts separately along the diagonal of a square matrix; merge these two maxima into one complex number and return it as a scalar 
MCF_Rows_sabsmax  find the complex numbers of largest sum Re+Im along rows and store these row maxima in a (complex) column vector 
MCF_Cols_sabsmax  find the complex numbers of largest sum Re+Im along columns and store these column maxima in a (complex) row vector 
MCF_Dia_sabsmax  find the complex number of largest sum Re+Im along the diagonal of a square matrix and return it as a (complex) scalar 
To determine the center of gravity of a matrix, you have the choice between the following two functions:
MF_centerOfGravityInd  center of gravity, returned as interpolated element indices 
MF_centerOfGravityV  center of gravity of an MZ matrix with explicitly given X and Y axes 
Back to MatrixLib Table of Contents OptiVec home
MF_Rows_exchange  exchange two rows 
MF_Cols_exchange  exchange two columns 
MF_Rows_add  add one row to another (destination += source) 
MF_Cols_add  add one column to another 
MF_Rows_sub  subtract one row from another (destination = source) 
MF_Cols_sub  subtract one column from another 
MF_Rows_Cadd  add scaled row to another (destination += C * source) 
MF_Cols_Cadd  add scaled column to another 
MF_Rows_lincomb  linear combination of two rows 
MF_Cols_lincomb  linear combination of two columns 
Back to MatrixLib Table of Contents OptiVec home
MF_addM  add two matrices 
MF_addMT  add one matrix and the transpose of another matrix MC = MA + MB^{T} 
MF_subM  subtract one matrix from another 
MF_subMT  subtract a transposed matrix MC = MA  MB^{T} 
MF_subrMT  subtract a matrix from another, transposed, matrix MC = MB^{T}  MA 
MF_mulC  multiply all matrix elements by a constant 
MCF_mulReC  multiply all elements of a complex matrix by a real number 
MF_divC  divide all matrix elements by a constant 
MCF_divReC  divide all elements of a complex matrix by a real number 
MFs_addM  scaled addition of two matrices: MC = c * (MA + MB) 
MFs_addMT  scaled addition of one matrix and the transpose of another: MC = c * (MA + MB^{T}) 
MFs_subM  scaled subtraction of two matrices: MC = c * (MA  MB) 
MFs_subMT  scaled subtraction of one matrix and the transpose of another: MC = c * (MA  MB^{T}) 
MFs_subrMT  scaled reverse subtraction of one matrix and the transpose of another: MC = c * (MB^{T}  MA) 
MF_lincomb  linear combination: MC = ca * MA + cb * MB 
b)Matrix multiplication:
MF_mulV  multiply a matrix by a column vector: Y = MA * X 
MF_TmulV  multiply the transpose of a matrix by a column vector: Y = MA^{T} * X 
VF_mulM  multiply a row vector by a matrix: Y = X * MA 
VF_mulMT  multiply a row vector by the transpose of a matrix: Y = X * MA^{T} 
MF_mulM  multiply two matrices: MC = MA * MB 
MF_mulMT  multiply one matrix by the transpose of another matrix: MC = MA * MB^{T} 
MF_TmulM  multiply the transpose of a matrix by another matrix: MC = MA^{T} * MB 
MF_TmulMT  multiply the transpose of one matrix by the transpose of another matrix: MC = MA^{T} * MB^{T} 
MCF_mulMH  multiply one matrix by the hermitian conjugate of another matrix: MC = MA * MB^{T *} 
MCF_HmulM  multiply the hermitian conjugate of a matrix by another matrix: MC = MA^{T *} * MB 
c) Multiplications of general matrices by diagonal matrices or vice versa
The diagonal matrix is passed to the respective function as a vector.
MF_mulMdia  multiply a general matrix by a diagonal matrix: MC = MA * MBDia 
MF_TmulMdia  multiply the transpose of a matrix by a diagonal matrix: MC = MA^{T} * MBDia 
MFdia_mulM  multiply a diagonal matrix by a general matrix: MC = MADia * MB 
MFdia_mulMT  multiply a diagonal matrix by the transpose of a general matrix: MC = MADia * MB^{T} 
Back to MatrixLib Table of Contents OptiVec home
MF_solve  solve simultaneous linear equations (using LU decomposition) 
MF_inv  invert a matrix 
MF_det  determinant of a matrix 
MF_solveBySVD  solve simultaneous linear equations, using Singular Value Decomposition 
MF_safeSolve  tries first solution by LUD; if that fails, SVD is done 
MF_LUdecompose  decompose into LU form 
MF_LUimprove  improve accuracy of LU decomposition by iteration 
MF_LUDresult  check if MF_LUdecompose was successful 
MF_LUDsetEdit  set editing threshold for MF_LUdecompose; may be used to work around singularities 
MF_LUDgetEdit  retrieve currently set threshold 
MF_LUsolve  solve simultaneous linear equations, given the matrix in LU form 
MF_LUinv  invert matrix already composed into LU form 
MF_LUdet  determinant of matrix already composed into LU form 
MF_SVdecompose  Singular Value Decomposition 
MF_SVsolve  solve SVdecomposed set of linear equations 
MF_SVDsetEdit  set threshold for Singular Value editing 
MF_SVDgetEdit  retrieve current threshold for Singular Value editing 
Back to MatrixLib Table of Contents OptiVec home
MFsym_eigenvalues  eigenvalues with or without eigenvectors of a symmetric real matrix 
MFsym_sqrt  squareroot of a symmetric, positive definite matrix 
Back to MatrixLib Table of Contents OptiVec home
MF_FFTtoC  Forward Fast Fourier Transform (FFT) of a real matrix; the result is a cartesian complex matrix 
MF_FFT  Forward and backward FFT of real matrices; using symmetry relations, the complex result is packed into a real matrix of the same size as the input matrix 
MCF_FFT  Forward and backward FFT of complex matrices 
MF_convolve  Convolution with a spatial response function 
MF_deconvolve  Deconvolution 
MF_filter  Spatial filtering 
MF_autocorr  Spatial autocorrelation 
MF_xcorr  Spatial crosscorrelation 
MF_spectrum  Spatial frequency spectrum 
MF_Rows_FFT  Fourier Transform along rows 
MF_Cols_FFT  Fourier Transform along columns 
Back to MatrixLib Table of Contents OptiVec home
VF_linregress  equallyweighted linear regression on XY data 
VF_linregresswW  the same with nonequal weighting 
VF_polyfit  fitting of one XY data set to a polynomial 
VF_polyfitwW  the same for nonequal datapoint weighting 
VF_linfit  fitting of one XY data set to an arbitrary function linear in its parameters 
VF_linfitwW  the same for nonequal datapoint weighting 
VF_nonlinfit  fitting of one XY data set to an arbitrary, possibly nonlinear function 
VF_nonlinfitwW  the same for nonequal datapoint weighting 
VF_multiLinfit  fitting of multiple XY data sets to one common linear function 
VF_multiLinfitwW  the same for nonequal datapoint weighting 
VF_multiNonlinfit  fitting of multiple XY data sets to one common nonlinear function 
VF_multiNonlinfitwW  the same for nonequal datapoint weighting 
MF_linfit  fit XYZ data to a function linear in its parameters 
MF_linfitwW  the same with nonequal weighting of individual data points 
MF_nonlinfit  fit XYZ data to an arbitrary, possibly nonlinear function 
MF_nonlinfitwW  the same for nonequal datapoint weighting 
MF_multiLinfit  fit multiple XYZ data sets to one common linear function 
MF_multiLinfitwW  the same for nonequal datapoint weighting 
MF_multiNonlinfit  fit multiple XYZ data sets to one common nonlinear function 
MF_multiNonlinfitwW  the same for nonequal datapoint weighting 
void LinModel( fVector BasFuncs, float x, unsigned nfuncs )
{
BasFuncs[0] = sin(x);
BasFuncs[1] = cos(x);
}
Note that the coefficients a_{i} are not used in the model function itself. The argument nfuncs (which is neglected in the above example) allows to use a variable number of basis functions. It is also possible to switch fitting parameters "on" and "off", i.e., "free" or "fixed". To this end, the routine VF_linfit takes a "status" array as an argument. For all members of the parameter vector that are to be treated as free, the corresponding "status" entry must be set to 1. If status_{i} is set to 0, the corresponding parameter, a_{i}, is treated as fixed at its value upon input.
Internally, VF_linfit employs a Singular Value Decomposition algorithm to obtain a solution even for (near)singular linear systems. Thereby, coefficients a_{i} whose significance is lower than a threshold, Thresh, are
set to 0 instead of infinity. This threshold can be modified by calling VF_setLinfitNeglect. The current threshold can be retrieved by VF_getLinfitNeglect.
Having noted above that functions like
y = sin(a_{0}x)+ cos(a_{1}x)
require their own treatment, we arrive at the last and most general class of model functions:
void NonlinModel( fVector Y, fVector X, ui size, fVector A)
{
for( ui i=0; i<size; i++ )
Y[i] = sin( A[0]*X[i] ) + cos( A[1]*X[i] );
}
VF_nonlinfit will call this model function again and again, hundreds and thousands of times, with different values in the parameter vector A. This means that a nonlinear fitting routine will spend most of its time in your model function (and its derivatives, see below), and you might wish to optimize this model function as far as possible, using vector functions. Suppose you have declared and allocated the fVector YHelp somewhere within your main program. Then you can write:
void OptNonlinModel( fVector Y, fVector X, ui size, fVector A)
{
VFx_sin( YHelp, X, size, A[0], 0.0, 1.0 );
VFx_cos( Y, X, size, A[1], 0.0, 1.0 );
VF_addV( Y, Y, YHelp, size );
}
The parameter array A you just saw being used by VF_nonlinfit to call your model function, is the one you have to pass as the first argument to VF_nonlinfit. Upon input, "A" must (!) contain your initial "best guess" of the parameters. The better your guess, the faster VF_nonlinfit will converge. If your initial guess is too far off, convergence might be very slow or even never attained. Upon output, "A" will contain the best set of parameters that could be found.
All nonlinfit functions return the figureofmerit of the best parameter array "A" found. To perform the fit, these functions need not only the fitting function itself, but also its partial derivatives with respect to the individual coefficients. Therefore, an argument "Derivatives" is required in calling the nonlinfit functions. If you happen to know the partial derivatives analytically, Derivatives should point to a function calculating them. If you do not know them, call with Derivatives = NULL (nil in Pascal/Delphi). In the latter case, all derivatives will be calculated numerically inside the functions. In cases where you know some, but not all of the partial derivatives of Y with respect to the coefficients a_{i}, you could also calculate dY / da_{i} wherever you have an analytic formula, and call VF_nonlinfit_autoDeriv for the remaining coefficients. To demonstrate this possibility, here is a function coding the derivatives for our nonlinear sine and cosine example:
void DerivModel( fVector dYdAi, fVector X, ui size, unsigned iPar, fVector A, VF_NONLINFITWORKSPACE *ws )
{
switch( iPar )
{
case 0:
for( ui i=0; i<size; i++ )
dYdAi[i] = X[i] * cos( A[0]*X[i] );
break;
case 1: /* say we don't know this derivative: */
VF_nonlinfit_autoDeriv( dYdAi, X, size, ipar, A, &ws );
}
}
As you might have guessed, the only purpose of the input argument "ws" for DerivModel is to be passed on to VF_nonlinfit_autoDeriv, in case this becomes necessary. Think twice, however, before resorting to VF_nonlinfit_autoDeriv: The function calls needed to determine the derivatives numerically may easily slow down your fit by a factor of 3 to 10! Anyway, as described above for the model function itself, also its derivatives should be implemented using vector functions. This leads us finally to the optimized version of the derivatives function:
void OptDerivModel( fVector dYdAi, fVector X, ui size, unsigned iPar, fVector A, VF_NONLINFITWORKSPACE *ws )
{
switch( iPar )
{
case 0:
VFx_cos( dYdAi, X, size, A[0], 0.0, 1.0 );
VF_mulV( dYdAi, dYdAi, X, size ); break;
case 1:
VFx_sin( dYdAi, X, size, A[1], 0.0, 1.0 );
VF_mulV( dYdAi, dYdAi, X, size );
}
}
Actual working examples for the polynomial, general linear, and general nonlinear cases are contained in the OptiVec shareware versions you can download.
The nonlinear fitting routines are highly sophisticated and offer the user a lot of different options. These options may be set by the function V_setNonlinfitOptions. To retrieve current settings, use
V_getNonlinfitOptions.
All options are packed into a structure (struct in C/C++, record in Pascal/Delphi) named VF_NONLINFITOPTIONS.
It has the following fields:
int FigureOfMerit;  0: least squares fitting
1: robust fitting, optimizing for minimum absolute deviation (default: 0) 
float AbsTolChi;  absolute change of c^{2} (default: EPSILON) 
float FracTolChi;  fractional change of c^{2} (default: SQRT_EPSILON) 
float AbsTolPar;  absolute change of all parameters (default: SQRT_MIN) 
float FracTolPar;  fractional change of all parameters (default: SQRT_EPSILON)
The four xxxTolChi and xxxTolPar parameters describe the convergence conditions: if the changes achieved in successive iterations are smaller than demanded by these criteria, this signals convergence. Criteria which are not applicable should be set to 0.0. 
unsigned HowOftenFulfill;  how often fulfill one of the above conditions before convergence is considered achieved (default: 3) 
unsigned LevelOfMethod;  1: LevenbergMarquardt method,
2: Downhill Simplex (Nelder and Mead) method, 3: both methods alternating; add 4 to this in order to try breaking out of local minima; 0: no fit, calculate only c^{2} (and Covar) (default: 1) 
unsigned LevMarIterations;  max.number of successful iterations of LevMar during one run (default: 100) 
unsigned LevMarStarts;  number of LevMar restarts per run (default: 2) 
float LambdaStart,
LambdaMin, LambdaMax, LambdaDiv, LambdaMul;  treatment of LevMar parameter lambda (don't touch, unless you are an expert!) 
unsigned DownhillIterations;  max. number of successful iterations in Downhill Simplex method (default: 200) 
float DownhillReflection,
DownhillContraction, DownhillExpansion;  treatment of reshaping of the simplex in Downhill Simplex method (don't touch unless you are an expert!) 
unsigned TotalStarts;  max. number of LevMar/Downhill pairs (default: 16) 
fVector UpperLimits, LowerLimits;  impose upper and/or lower limits on parameters. Default for both: NULL or nil 
void (*Restrictions)(void);  pointer to a userdefined function, implementing restrictions on the parameters which cannot be formulated as simple upper/lower limits. The function must check the whole parameter vector and edit the parameters as needed. Default: NULL or nil.
A typical application of this function would be a model, in which the parameters A_{i} have to be in a specific order (ascending or descending). So your Restrictions might call VF_sort to enforce that order. Apart from such very special cases, it is generally better not to use Restrictions at all, but let the nonlinfit routine run its course without too much interfering in it. 
As described above, all nonlinear fitting functions need a set of variables for internal use which, for the VF_nonlinfit family of functions, is contained in a struct (Delphi: record) VF_NONLINFITWORKSPACE. It is passed by its address to the respective function. This set of variables need not be initialized. It does not contain userretrievable information. In the MF_nonlinfit family of functions, a similar set of internal variables is needed as MF_NONLINFITWORKSPACE.
A typical call of VF_nonlinfit would look like this in C/C++:
VF_NONLINFITWORKSPACE ws;
VF_NONLINFITOPTIONS fopt;
V_getNonlinfitOptions( &fopt );
// at this point, modify fopt as desired...
VF_nonlinfit( ParValues, AStatus, nParameters, X, Y, sz, ModelFunc, DerivativeFuncs, &ws, &fopt );
or in Delphi:
ws: VF_NONLINFITWORKSPACE;
fopt: VF_NONLINFITOPTIONS;
V_getNonlinfitOptions( @fopt );
// at this point, modify fopt as desired...
VF_nonlinfit( ParValues, AStatus, nParameters, X, Y, sz, @ModelFunc, @DerivativeFuncs, @ws, @fopt );
If you replace the parameter &ws with NULL / nil, VF_nonlinfit will generate its own workspace.
fVector X, Y;  X and Y vectors (independent variables) 
fMatrix MZ;  measured data z=f(x,y) 
fMatrix MInvVar;  inverse variances of the individual matrix elements (needed only for the "with weights" functions) 
unsigned htZ, lenZ;  matrix dimensions 
float WeightOfExperiment;  individual weight to be assigned to the whole experiment (again needed only for the weighted variants) 
fVector X, Y;  X and Y vectors 
fVector InvVar;  inverse variances of the individual vector elements (needed only for the "with weights" functions) 
ui size;  the number of vector elements 
float WeightOfExperiment;  individual weight to be assigned to the whole experiment (again needed only for the weighted variants) 
type VF_EXPERIMENT = record
X, Y, InvVar: fVector;
size: UInt;
WeightOfExperiment: Single;
end;
type PVF_EXPERIMENT = ^VF_EXPERIMENT;
type MF_EXPERIMENT = record
X, Y: fVector;
MZ, MInvVar: fMatrix;
htZ, lenZ: UInt;
WeightOfExperiment: Single;
end;
type PMF_EXPERIMENT = ^MF_EXPERIMENT;
Both in VF_EXPERIMENT and MF_EXPERIMENT, InvVar and WeightOfExperiment are needed only for the weighted variants of the multifit functions.
VF_nonlinfit_getBestValues  best set of parameters found so far 
VF_nonlinfit_getFOM  the best figureofmerit (c^{2}, chisquare, or c, chiabs, depending on the VF_NONLINFITOPTIONS member "FigureOfMerit") obtained so far 
VF_nonlinfit_getTestDir  returns the test direction (+1 for upwards, 1 for downwards) during "breakout" attempts (levelofmethod greater than 3, see the description of VF_NONLINFITOPTIONS below) 
VF_nonlinfit_getTestPar  index of the parameter currently under "breakout" investigation 
VF_nonlinfit_getTestRun  index of the current "breakout" test run (for each fitted parameter, one test run is performed) 
VF_nonlinfit_stop  stops fitting after completion of the current LevenbergMarquardt or Downhill cycle 
There are two possibilities how these functions can be called. The first way is to call them from within your model function. For example, you might install a counter into your model function and, upon reaching a certain number of calls, retrieve the current chi2, parameter set, etc. The second way is open only for multithreaded programs: a thread, different from the one performing the fit, may call any of the above functions to control the progress of the fit.
Both of these ways require the VF_NONLINFITWORKSPACE to be globally accessible. When a nonlinfit function is called from multiple threads simultaneously, this poses an obvious problem, especially if you wish to call the helper function from your ModelFunc: How does the model function "know" in real time, from which thread it was called? The solution is: Store a thread index in your parameter vector A. For example, in your main program, you might have something like:
float A1[11], A2[11], ...;
int AStatus1[11], AStatus2[11], ...;
VF_NONLINFITWORKSPACE ws1, ws2, ...;
....
VF_random( A1, 10, 1, 1.0, +1.0 ); // any better guess for starting values would be welcome
VI_equC( AStatus1, 10, 1 ); // fist 10 parameters are to be fitted
A1[10] = 1; // thread index = 1
AStatus1[10] = 0; // of course, this parameter is frozen
... throw thread #1 with VF_nonlinfit taking A1, AStatus1, and ws1 as parameters
VF_random( A2, 10, 1, 1.0, +1.0 ); // starting values for second fit
VI_equC( AStatus2, 10, 1 );
A2[10] = 2; // thread index = 2
AStatus2[10] = 0; // again, this parameter is frozen
... throw thread #2 with VF_nonlinfit taking A2, AStatus2, and ws2 as parameters
In your model function then, you would have something along the lines of:
void ModelFunc( fVector Y, fVector X, ui size, fVector A )
{
static unsigned counter1=0, counter2=0, ...;
float CurrentFOM;
....
switch( (unsigned)(A[10]) )
{
case 1: if( (++counter1 % 1024) == 0 )
{ CurrentFOM = VF_nonlinfit_getFOM( ws1 );
printf( "\nCount 1: %u, FOM: %f", counter1, CurrentFOM );
} break;
case 1: if( (++counter2 % 1024) == 0 )
{ CurrentFOM = VF_nonlinfit_getFOM( ws2 );
printf( "\nCount 2: %u, FOM: %f", counter2, CurrentFOM );
} break;
....
}
....
}
Back to MatrixLib Table of Contents OptiVec home
MF_cprint  print a matrix to the screen. If necessary, rows are cut off at the screen boundaries. If there are more rows than screen lines, proper paging is applied. (DOS and Delphi only) 
MF_print  print a matrix to the screen (without paging or row cutoff); (DOS, EasyWin and Delphi only) 
MF_fprint  print a matrix in ASCII format to a stream 
MF_store  store in binary format 
MF_recall  retrieve in binary format 
MF_write  write in ASCII format in a stream 
MF_read  read from an ASCII file 
Back to MatrixLib Table of Contents OptiVec home
MF_xyzAutoDensityMap  Color density map for z=f(x,y) with automatic scaling of the x and y axes and of the color density scale between mincolor and maxcolor 
MF_xyzDataDensityMap  z=f(x,y) color density map, plotted into an existing axis frame, and using the color density scale set by the last call to an AutoDensityMap function 
MF_zAutoDensityMap  Color density map for z=f(i,j) with automatic scaling of the x and y axes and of the color density scale between mincolor and maxcolor. i and j are the indices in x and y direction, respectively 
MF_zDataDensityMap  Color density map for z=f(i,j), plotted into an existing axis frame, and using the color density scale set by the last call to an AutoDensityMap function 
M_setDensityBounds  Set a color scale for matrix colordensity plots. 
M_setDensityMapBounds  Set a color scale and draw an XY coordinate system for matrix colordensity plots. 
M_findDensityMapBounds  Calculate a color scale and draw an XY coordinate system for matrix colordensity plots, ensuring that the grid lines of the coordinate system correspond to exact (rather than only rounded) values. 
Back to MatrixLib Table of Contents OptiVec home
Copyright © 19962020 OptiCode  Dr. Martin Sander Software Development