| |
| MPFIT: A MINPACK-1 Least Squares Fitting Library in C |
| |
| Original public domain version by B. Garbow, K. Hillstrom, J. More' |
| (Argonne National Laboratory, MINPACK project, March 1980) |
| |
| Tranlation to C Language by S. Moshier (moshier.net) |
| |
| Enhancements, documentation and packaging by C. Markwardt |
| (comparable to IDL fitting routine MPFIT |
| see http://cow.physics.wisc.edu/~craigm/idl/idl.html) |
| Copyright (C) 2003, 2004, 2006, 2007, 2009, 2010, 2013 Craig B. Markwardt |
| |
| |
| SUMMARY of CHANGES |
| 16 Feb 2009 - version 1.0 - initial release |
| 18 Feb 2009 - version 1.1 - add 'version' field to 'results' struct |
| 22 Nov 2009 - - allow to compile with C++ compiler |
| - change to memset() instead of nonstandard bzero() |
| - for Microsoft, proprietary equivalent of finite() |
| 04 Oct 2010 - - add static declarations, remove some compiler warnings |
| (reported by Lars Kr. Lundin) |
| 13 Nov 2010 - version 1.2 - additional documentation, cleanup of mpfit.h |
| 23 Apr 2013 - version 1.3 - add MP_NO_ITER; bug fix mpside=2 when debugging |
| (thanks to M. Wojdyr) |
| |
| $Id: README,v 1.18 2016/06/02 19:14:16 craigm Exp $ |
| |
| INTRODUCTION |
| |
| MPFIT uses the Levenberg-Marquardt technique to solve the |
| least-squares problem. In its typical use, MPFIT will be used to |
| fit a user-supplied function (the "model") to user-supplied data |
| points (the "data") by adjusting a set of parameters. MPFIT is |
| based upon MINPACK-1 (LMDIF.F) by More' and collaborators. |
| |
| For example, a researcher may think that a set of observed data |
| points is best modelled with a Gaussian curve. A Gaussian curve is |
| parameterized by its mean, standard deviation and normalization. |
| MPFIT will, within certain constraints, find the set of parameters |
| which best fits the data. The fit is "best" in the least-squares |
| sense; that is, the sum of the weighted squared differences between |
| the model and data is minimized. |
| |
| The Levenberg-Marquardt technique is a particular strategy for |
| iteratively searching for the best fit. This particular |
| implementation is drawn from a robust routine called MINPACK-1 (see |
| NETLIB). This version allows upper and lower bounding constraints |
| to be placed on each parameter, or the parameter can be held fixed. |
| |
| The user-supplied function should compute an array of weighted |
| deviations between model and data. In a typical scientific problem |
| the residuals should be weighted so that each deviate has a |
| gaussian sigma of 1.0. If X represents values of the independent |
| variable, Y represents a measurement for each value of X, and ERR |
| represents the error in the measurements, then the deviates could |
| be calculated as follows: |
| |
| for (i=0; i<m; i++) { |
| deviates[i] = (y[i] - f(x[i])) / err[i]; |
| } |
| |
| where m is the number of data points, and where F is the |
| function representing the model evaluated at X. If ERR are the |
| 1-sigma uncertainties in Y, then the sum of deviates[i] squared will |
| be the total chi-squared value, which MPFIT will seek to minimize. |
| |
| Simple constraints can be placed on parameter values by using the |
| "pars" parameter to MPFIT, and other parameter-specific options can be |
| set. See below for a description of this optional parameter |
| descriptor. |
| |
| MPFIT does not perform more general optimization tasks. MPFIT is |
| customized, based on MINPACK-1, to the least-squares minimization |
| problem. |
| |
| BASIC CALLING INTERFACE of MPFIT |
| |
| At the very least, the user must supply a 'user function,' which |
| calculates the residuals as described above, and one or more |
| parameters. The calling interface to mpfit is: |
| |
| int mpfit(mp_func funct, int m, int npar, |
| double *xall, mp_par *pars, mp_config *config, |
| void *private_data, |
| mp_result *result); |
| |
| The user function <b>funct</b> is a C function which computes the |
| residuals using any user data that is required. The number of |
| residuals is specified by the integer <b>m</b>. The nature and |
| properties of user function are described in greater detail below. |
| |
| The user function parameters are passed to mpfit as the array |
| double xall[npar]; |
| where <b>npar</b> is the number of parameters of the user function. |
| It must be the case that m > npar, i.e. there are more data points |
| than free parameters. The user must pass an initial "best guess" of |
| the user parameters to mpfit(), and upon successful return, the xall[] |
| array will contain the "best fit" parameters of the fit (and the original |
| values are overwritten). |
| |
| The user function is responsible to compute an array of generic |
| residuals. Usually these residuals will depend on user-dependent |
| quantities such as measured data or independent variables such "x" |
| when fitting a function of the form y(x). The user should pass these |
| quantities using the optional <b>private_data</b> pointer. If |
| private_data is not used then a null pointer should be passed. |
| Otherwise, it can be any C pointer, typically a pointer to a structure |
| such as this: |
| struct example_private_data { /* EXAMPLE: fitting y(x) */ |
| double *x; /* x - independent variable of model */ |
| double *y; /* y - measured "y" values */ |
| double *y_error; /* y_error - measurement uncertainty in y */ |
| }; |
| The user would be responsible to declare such a structure, and to fill |
| it with pointers to the relevant data arrays before calling mpfit(). |
| mpfit() itself does not inspect or modify the private_data contents, |
| but merely passes it along to the user function |
| |
| The structure <b>pars</b> is optional, and allows the user to specify |
| additional information and constraints about each user function |
| parameter. For example, the user can specify simple bounding |
| constraints. If passed, then pars[] must be dimensioned as, |
| mp_par pars[npar]; |
| where mp_par is a structure defined in mpfit.h. If no special |
| parameter information is necessary, then the user can pass |
| pars == 0. |
| |
| The optional structure <b>config</b> configures how mpfit() behaves, |
| and is described in greater detail below. By passing a null pointer, |
| the user obtains the default behavior. |
| |
| The structure <b>result</b> contains results of the fit, returned by |
| mpfit(). The user should pass a pointer to a structure of type |
| 'mp_result' (which should be zeroed), and upon return, the structure |
| is filled with various diagnostic quantities about the fitting run. |
| These quantities are described in greater detail below. If these |
| diagnostics are not required, then the user may pass a null pointer. |
| |
| |
| USER FUNCTION |
| |
| The user must define a function which computes the appropriate values |
| as specified above. The function must compute the weighted deviations |
| between the model and the data. The user function may also optionally |
| compute explicit derivatives (see below). The user function should |
| be of the following form: |
| |
| int myfunct(int m, int n, double *p, double *deviates, |
| double **derivs, void *private) |
| { |
| int i; |
| /* Compute function deviates */ |
| for (i=0; i<m; i++) { |
| deviates[i] = {function of p[] and private data}; |
| } |
| |
| return 0; |
| } |
| |
| The user function parameters are defined as follows: |
| int m - number of data points |
| int n - number of parameters |
| double *p - array of n parameters |
| double *deviates - array of m deviates to be returned by myfunct() |
| double **derivs - used for user-computed derivatives (see below) |
| (= 0 when automatic finite differences are computed) |
| |
| User functions may also indicate a fatal error condition by returning |
| a negative error code. Error codes from -15 to -1 are reserved for |
| the user function. |
| |
| |
| EXAMPLE 1 - USER FUNCTION y(x) |
| |
| Here is a sample user function which computes the residuals for a simple |
| model fit of the form y = f(x). The residuals are defined as |
| (y[i] - f(x[i]))/y_error[i]. We will use the example_private_data |
| structure defined above. The user function would appear like this: |
| |
| int myfunct_y_x(int m, int n, double *p, double *deviates, |
| double **derivs, struct example_private_data *private) |
| { |
| int i; |
| double *x, *y, *y_error; |
| |
| /* Retrieve values of x, y and y_error from private structure */ |
| x = private.x; |
| y = private.y; |
| y_error = private.y_error; |
| |
| /* Compute function deviates */ |
| for (i=0; i<m; i++) { |
| deviates[i] = (y[i] - f(x[i], p)) / y_error[i]; |
| } |
| |
| return 0; |
| } |
| |
| In this example, the user would be responsible to define the function |
| f(x,p) where x is the independent variable and p is the array of user |
| parameters. Although this example shows y = f(x,p), it is trivial to |
| extend it to additional independent variables, such as u = f(x,y,p), |
| by adding additional elements to the private data structure. |
| |
| The setup and call to mpfit() could look something like this: |
| int m = 10; /* 10 data points */ |
| double x[10] = {1,2,3,4,5,6,7,8,9,10}; /* Independent variable */ |
| double y[10] = {-0.42,1.59,0.8,3.2,2.09,6.39,6.29,9.01,7.47,7.58};/* Y measurements*/ |
| double y_error[10] = {1,1,1,1,1,1,1,1,1,1}; /* Y uncertainties */ |
| struct example_private_data private; |
| |
| double xall[2] = { 1.0, 1.0 }; /* User parameters - initial guess */ |
| /* Fill private data structure to pointers with user data */ |
| private.x = x; |
| private.y = y; |
| private.y_error = y_error; |
| |
| /* Call mpfit() */ |
| status = mpfit(&myfunct_y_x, m, npar, xall, 0, 0, &private, 0); |
| |
| USER-COMPUTED DERIVATIVES |
| |
| The user function can also compute function derivatives, which are |
| used in the minimization process. This can be useful to save time, or |
| when the derivative is tricky to evaluate numerically. |
| |
| Users should pass the "pars" parameter (see below), and for the 'side' |
| structure field, the value of 3 should be passed. This indicates to |
| mpfit() that it should request the user function will compute the |
| derivative. NOTE: this setting is specified on a parameter by |
| parameter basis. This means that users are free to choose which |
| parameters' derivatives are computed explicitly and which |
| numerically. A combination of both is allowed. However, since the |
| default value of 'side' is 0, derivative estimates will be numerical |
| unless explicitly requested using the 'side' structure field. |
| |
| The user function should only compute derivatives when the 'derivs' |
| parameter is non-null. Note that even when user-computed derivatives |
| are enabled, mpfit() may call the user function with derivs set to |
| null, so the user function must check before accessing this pointer. |
| |
| The 'derivs' parameter is an array of pointers, one pointer for each |
| parameter. If derivs[i] is non-null, then derivatives are requested |
| for the ith parameter. Note that if some parameters are fixed, or |
| pars[i].side is not 3 for some parameters, then derivs[i] will be null |
| and the derivatives do not need to be computed for those parameters. |
| |
| derivs[j][i] is the derivative of the ith deviate with respect to the |
| jth parameter (for 0<=i<m, 0<=j<n). Storage has already been |
| allocated for this array, and the user is not responsible for freeing |
| it. The user function may assume that derivs[j][i] are initialized to |
| zero. |
| |
| The function prototype for user-computed derivatives is: |
| |
| int myfunct_with_derivs(int m, int n, double *p, double *deviates, |
| double **derivs, void *private) |
| { |
| int i; |
| /* Compute function deviates */ |
| for (i=0; i<m; i++) { |
| deviates[i] = {function of x[i], p and private data}; |
| } |
| |
| /* If derivs is non-zero then user-computed derivatives are |
| requested */ |
| if (derivs) { |
| int j; |
| for (j=0; j<n; j++) if (derivs[j]) { |
| /* It is only required to compute derivatives when |
| derivs[ipar] is non-null */ |
| for (i=0; i<m; i++) { |
| derivs[j][i] = {derivative of the ith deviate with respect to |
| the jth parameter = d(deviates[i])/d(par[j])} |
| } |
| } |
| } |
| |
| return 0; |
| } |
| |
| TESTING EXPLICIT DERIVATIVES |
| |
| In principle, the process of computing explicit derivatives should be |
| straightforward. In practice, the computation can be error prone, |
| often being wrong by a sign or a scale factor. |
| |
| In order to be sure that the explicit derivatives are correct, the |
| user can set pars[i].deriv_debug = 1 for parameter i (see below for a |
| description of the "pars" structure). This will cause mpfit() to |
| print *both* explicit derivatives and numerical derivatives to the |
| console so that the user can compare the results. This would |
| typically be used during development and debugging to be sure the |
| calculated derivatives are correct, than then deriv_debug would be set |
| to zero for production use. |
| |
| If you want debugging derivatives, it is important to set pars[i].side |
| to the kind of numerical derivative you want to compare with. |
| pars[i].side should be set to 0, 1, -1, or 2, and *not* set to 3. |
| When pars[i].deriv_debug is set, then mpfit() automatically |
| understands to request user-computed derivatives. |
| |
| The console output will be sent to the standard output, and will |
| appear as a block of ASCII text like this: |
| FJAC DEBUG BEGIN |
| # IPNT FUNC DERIV_U DERIV_N DIFF_ABS DIFF_REL |
| FJAC PARM 1 |
| .... derivative data for parameter 1 .... |
| FJAC PARM 2 |
| .... derivative data for parameter 2 .... |
| .... and so on .... |
| FJAC DEBUG END |
| |
| which is to say, debugging data will be bracketed by pairs of "FJAC |
| DEBUG" BEGIN/END phrases. Derivative data for individual parameter i |
| will be labeled by "FJAC PARM i". The columns are, in order, |
| |
| IPNT - data point number j |
| FUNC - user function evaluated at X[j] |
| DERIV_U - user-calculated derivative d(FUNC(X[j]))/d(P[i]) |
| DERIV_N - numerically calculated derivative according to pars[i].side value |
| DIFF_ABS - difference between DERIV_U and DERIV_N = fabs(DERIV_U-DERIV_N) |
| DIFF_REL - relative difference = fabs(DERIV_U-DERIV_N)/DERIV_U |
| |
| Since individual numerical derivative values may contain significant |
| round-off errors, it is up to the user to critically compare DERIV_U |
| and DERIV_N, using DIFF_ABS and DIFF_REL as a guide. |
| |
| |
| CONSTRAINING PARAMETER VALUES WITH THE PARS PARAMETER |
| |
| The behavior of MPFIT can be modified with respect to each |
| parameter to be fitted. A parameter value can be fixed; simple |
| boundary constraints can be imposed; and limitations on the |
| parameter changes can be imposed. |
| |
| If fitting constraints are to be supplied, then the user should pass |
| an array of mp_par structures to mpfit() in the pars parameter. If |
| pars is set to 0, then the fitting parameters are asssumed to be |
| unconstrained. |
| |
| pars should be an array of structures, one for each parameter. |
| Each parameter is associated with one element of the array, in |
| numerical order. The structure is declared to have the following |
| fields: |
| |
| .fixed - a boolean value, whether the parameter is to be held |
| fixed or not. Fixed parameters are not varied by |
| MPFIT, but are passed on to MYFUNCT for evaluation. |
| |
| .limited - a two-element boolean array. If the first/second |
| element is set, then the parameter is bounded on the |
| lower/upper side. A parameter can be bounded on both |
| sides. |
| |
| |
| .limits - a two-element float or double array. Gives the |
| parameter limits on the lower and upper sides, |
| respectively. Zero, one or two of these values can be |
| set, depending on the values of LIMITED. |
| |
| .parname - a string, giving the name of the parameter. The |
| fitting code of MPFIT does not use this tag in any |
| way. However, it may be used for output purposes. |
| |
| .step - the step size to be used in calculating the numerical |
| derivatives. If set to zero, then the step size is |
| computed automatically. |
| This value is superceded by the RELSTEP value. |
| |
| .relstep - the *relative* step size to be used in calculating |
| the numerical derivatives. This number is the |
| fractional size of the step, compared to the |
| parameter value. This value supercedes the STEP |
| setting. If the parameter is zero, then a default |
| step size is chosen. |
| |
| .side - the sidedness of the finite difference when computing |
| numerical derivatives. This field can take four |
| values: |
| |
| 0 - one-sided derivative computed automatically |
| 1 - one-sided derivative (f(x+h) - f(x) )/h |
| -1 - one-sided derivative (f(x) - f(x-h))/h |
| 2 - two-sided derivative (f(x+h) - f(x-h))/(2*h) |
| 3 - user-computed explicit derivatives |
| |
| Where H is the STEP parameter described above. The |
| "automatic" one-sided derivative method will chose a |
| direction for the finite difference which does not |
| violate any constraints. The other methods do not |
| perform this check. The two-sided method is in |
| principle more precise, but requires twice as many |
| function evaluations. Default: 0. |
| |
| .deriv_debug - flag to enable/disable console debug logging of |
| user-computed derivatives, as described above. 1=enable |
| debugging; 0=disable debugging. Note that when |
| pars[i].deriv_debug is set, then pars[i].side should be |
| set to 0, 1, -1 or 2, depending on which numerical |
| derivative you wish to compare to. |
| Default: 0. |
| |
| RETURN VALUE of MPFIT() |
| |
| mpfit() returns an integer status code. In principle, any positive return |
| value from mpfit() indicates success or partial success. |
| |
| The possible error codes are integer constants: |
| MP_ERR_INPUT /* General input parameter error */ |
| MP_ERR_NAN /* User function produced non-finite values */ |
| MP_ERR_FUNC /* No user function was supplied */ |
| MP_ERR_NPOINTS /* No user data points were supplied */ |
| MP_ERR_NFREE /* No free parameters */ |
| MP_ERR_MEMORY /* Memory allocation error */ |
| MP_ERR_INITBOUNDS /* Initial values inconsistent w constraints*/ |
| MP_ERR_BOUNDS /* Initial constraints inconsistent */ |
| MP_ERR_PARAM /* General input parameter error */ |
| MP_ERR_DOF /* Not enough degrees of freedom */ |
| |
| The possible success status codes are: |
| MP_OK_CHI /* Convergence in chi-square value */ |
| MP_OK_PAR /* Convergence in parameter value */ |
| MP_OK_BOTH /* Both MP_OK_PAR and MP_OK_CHI hold */ |
| MP_OK_DIR /* Convergence in orthogonality */ |
| MP_MAXITER /* Maximum number of iterations reached */ |
| MP_FTOL /* ftol is too small; no further improvement*/ |
| MP_XTOL /* xtol is too small; no further improvement*/ |
| MP_GTOL /* gtol is too small; no further improvement*/ |
| |
| |
| CONFIGURING MPFIT() |
| |
| mpfit() is primarily configured using the config parameter, which in |
| turn is defined as a pointer to the following structure: |
| |
| struct mp_config_struct { |
| /* NOTE: the user may set the value explicitly; OR, if the passed |
| value is zero, then the "Default" value will be substituted by |
| mpfit(). */ |
| double ftol; /* Relative chi-square convergence criterium Default: 1e-10 */ |
| double xtol; /* Relative parameter convergence criterium Default: 1e-10 */ |
| double gtol; /* Orthogonality convergence criterium Default: 1e-10 */ |
| double epsfcn; /* Finite derivative step size Default: MP_MACHEP0 */ |
| double stepfactor; /* Initial step bound Default: 100.0 */ |
| double covtol; /* Range tolerance for covariance calculation Default: 1e-14 */ |
| int maxiter; /* Maximum number of iterations. If maxiter == MP_NO_ITER, |
| then basic error checking is done, and parameter |
| errors/covariances are estimated based on input |
| parameter values, but no fitting iterations are done. |
| Default: 200 |
| */ |
| int maxfev; /* Maximum number of function evaluations, or 0 for no limit |
| Default: 0 (no limit) */ |
| int nprint; /* Default: 1 */ |
| int douserscale;/* Scale variables by user values? |
| 1 = yes, user scale values in diag; |
| 0 = no, variables scaled internally (Default) */ |
| int nofinitecheck; /* Disable check for infinite quantities from user? |
| 0 = do not perform check (Default) |
| 1 = perform check |
| */ |
| mp_iterproc iterproc; /* Placeholder pointer - must set to 0 */ |
| }; |
| |
| Generally speaking, a value of 0 for a field in the structure above |
| indicates that a default value should be used, indicated in |
| parentheses. Therefore, a user should zero this structure before |
| passing it. |
| |
| |
| EXTRACTING RESULTS FROM MPFIT() |
| |
| The basic result of mpfit() is the set of updated parameters, xall. |
| However, there are other auxiliary quantities that the user can |
| extract by using the results parameter. This is a structure defined |
| like this: |
| |
| /* Definition of results structure, for when fit completes */ |
| struct mp_result_struct { |
| double bestnorm; /* Final chi^2 */ |
| double orignorm; /* Starting value of chi^2 */ |
| int niter; /* Number of iterations */ |
| int nfev; /* Number of function evaluations */ |
| int status; /* Fitting status code */ |
| |
| int npar; /* Total number of parameters */ |
| int nfree; /* Number of free parameters */ |
| int npegged; /* Number of pegged parameters */ |
| int nfunc; /* Number of residuals (= num. of data points) */ |
| |
| double *resid; /* Final residuals |
| nfunc-vector, or 0 if not desired */ |
| double *xerror; /* Final parameter uncertainties (1-sigma) |
| npar-vector, or 0 if not desired */ |
| double *covar; /* Final parameter covariance matrix |
| npar x npar array, or 0 if not desired */ |
| char version[20]; /* MPFIT version string */ |
| }; |
| |
| All of the scalar numeric quantities are filled when mpfit() returns, |
| and any incoming value will be overwritten. |
| |
| If the user would like the final residuals, parameter errors, or the |
| covariance matrix, then they should allocate the required storage, and |
| pass a pointer to it in the corresponding structure field. A pointer |
| value of zero indicates that the array should not be returned. Thus, |
| by default, the user should zero this structure before passing it. |
| The user is responsible for allocating and freeing resid, xerror and |
| covar storage. |
| |
| The 'version' string can be used to interpret the behavior of MPFIT, |
| in case the behavior changes over time. Version numbers will be of |
| the form "i.j" or "i.j.k" where i, j and k are integers. |
| |
| |
| EXAMPLE 2 - Basic call |
| |
| This example does the basics, which is to perform the optimization |
| with no constraints, and returns the scalar results in the result |
| structure. |
| |
| mp_result result; |
| |
| memset(&result, 0, sizeof(result)); |
| status = mpfit(myfunct, m, n, p, 0, 0, 0, &result); |
| |
| |
| EXAMPLE 3 - Requesting Parameter Errors |
| |
| The only modification needed to retrieve parameter uncertainties is to |
| allocate storage for it, and place a pointer to it in result. |
| |
| mp_result result; |
| double perror[n]; |
| |
| memset(&result, 0, sizeof(result)); |
| result.xerror = perror; |
| status = mpfit(myfunct, m, n, p, 0, 0, 0, &result); |
| |
| |
| EXAMPLE 4 - Fixing a parameter |
| |
| This example fixes the third parameter (i.e. p[2]) at its starting |
| value |
| |
| mp_result result; |
| mp_par pars[n]; |
| |
| memset(&pars[0], 0, sizeof(pars)); |
| pars[2].fixed = 1; |
| |
| memset(&result, 0, sizeof(result)); |
| status = mpfit(myfunct, m, n, p, pars, 0, 0, &result); |
| |
| EXAMPLE 5 - Applying parameter constraints |
| |
| This example ensures that the third parameter (i.e. p[2]) is always |
| greater or equal to ten. |
| |
| mp_result result; |
| mp_par pars[n]; |
| |
| memset(&pars[0], 0, sizeof(pars)); |
| pars[2].limited[0] = 1; /* limited[0] indicates lower bound */ |
| pars[2].limits[0] = 10.0; /* Actual constraint p[2] >= 10 */ |
| |
| memset(&result, 0, sizeof(result)); |
| status = mpfit(myfunct, m, n, p, pars, 0, 0, &result); |
| |
| |
| EXAMPLE 6 - Increasing maximum number of iterations |
| |
| This example changes the maximum number of iterations from its default |
| to 1000. |
| |
| mp_config config; |
| mp_result result; |
| |
| memset(&config, 0, sizeof(config)); |
| config.maxiter = 1000; |
| memset(&result, 0, sizeof(result)); |
| status = mpfit(myfunct, m, n, p, 0, &config, 0, &result); |
| |
| |
| EXAMPLE 7 - Passing private data to user function |
| |
| This example passes "private" data to its user function using the |
| private parameter. It assumes that three variables, x, y, and ey, |
| already exist, and that the user function will know what to do with |
| them. |
| |
| mp_result result; |
| struct mystruct { |
| double *x; |
| double *y; |
| double *ey; |
| }; |
| |
| struct mystruct mydata; |
| |
| mydata.x = x; |
| mydata.y = y; |
| mydata.ey = ey; |
| |
| memset(&result, 0, sizeof(result)); |
| status = mpfit(myfunct, m, n, p, 0, 0, (void *)&mydata, &result); |
| |
| |
| EXAMPLE 8 - Complete example |
| |
| |
| This example shows how to fit sample data to a linear function |
| y = f(x) = a - b*x, where a and b are fitted parameters. There is sample |
| data included within the program. The parameters are initialized to |
| a = 1.0, b = 1.0, and then the fitting occurs. The function |
| residuals; within linfunc(), the value of f(x) is computed. The main |
| routine, <b>main()</b> initializes the variables and calls mpfit(). |
| |
| |
| #include "mpfit.h" |
| #include <stdio.h> |
| #include <stdlib.h> |
| #include <math.h> |
| |
| /* This is the private data structure which contains the data points |
| and their uncertainties */ |
| struct vars_struct { |
| double *x; |
| double *y; |
| double *ey; |
| }; |
| |
| /* |
| * linear fit function |
| * |
| * m - number of data points |
| * n - number of parameters (2) |
| * p - array of fit parameters |
| * dy - array of residuals to be returned |
| * vars - private data (struct vars_struct *) |
| * |
| * RETURNS: error code (0 = success) |
| */ |
| int linfunc(int m, int n, double *p, double *dy, double **dvec, void *vars) |
| { |
| int i; |
| struct vars_struct *v = (struct vars_struct *) vars; |
| double *x, *y, *ey, f; |
| |
| x = v->x; |
| y = v->y; |
| ey = v->ey; |
| |
| for (i=0; i<m; i++) { |
| f = p[0] - p[1]*x[i]; /* Linear fit function; note f = a - b*x */ |
| dy[i] = (y[i] - f)/ey[i]; |
| } |
| |
| return 0; |
| } |
| |
| /* Test harness routine, which contains test data, invokes mpfit() */ |
| int main(int argc, char *argv[]) |
| { |
| /* X - independent variable */ |
| double x[] = {-1.7237128E+00,1.8712276E+00,-9.6608055E-01, |
| -2.8394297E-01,1.3416969E+00,1.3757038E+00, |
| -1.3703436E+00,4.2581975E-02,-1.4970151E-01, |
| 8.2065094E-01}; |
| /* Y - measured value of dependent quantity */ |
| double y[] = {1.9000429E-01,6.5807428E+00,1.4582725E+00, |
| 2.7270851E+00,5.5969253E+00,5.6249280E+00, |
| 0.787615,3.2599759E+00,2.9771762E+00, |
| 4.5936475E+00}; |
| double ey[10]; /* Measurement uncertainty - initialized below */ |
| |
| double p[2] = {1.0, 1.0}; /* Initial conditions */ |
| double pactual[2] = {3.20, 1.78}; /* Actual values used to make data */ |
| double perror[2]; /* Returned parameter errors */ |
| int i; |
| struct vars_struct v; /* Private data structure */ |
| int status; |
| mp_result result; |
| |
| memset(&result,0,sizeof(result)); /* Zero results structure */ |
| result.xerror = perror; |
| for (i=0; i<10; i++) ey[i] = 0.07; /* Data errors */ |
| |
| /* Fill private data structure */ |
| v.x = x; |
| v.y = y; |
| v.ey = ey; |
| |
| /* Call fitting function for 10 data points and 2 parameters */ |
| status = mpfit(linfunc, 10, 2, p, 0, 0, (void *) &v, &result); |
| |
| printf("*** testlinfit status = %d\n", status); |
| /* ... print or use the results of the fitted parametres p[] here! ... */ |
| |
| return 0; |
| } |
| |
| |