1*5ddc57e5SXin Li /*
2*5ddc57e5SXin Li * Library: lmfit (Levenberg-Marquardt least squares fitting)
3*5ddc57e5SXin Li *
4*5ddc57e5SXin Li * File: surface1.c
5*5ddc57e5SXin Li *
6*5ddc57e5SXin Li * Contents: Example for generic minimization with lmmin():
7*5ddc57e5SXin Li fit data y(t) by a function f(t;p), where t is a 2-vector.
8*5ddc57e5SXin Li *
9*5ddc57e5SXin Li * Note: Any modification of this example should be copied to
10*5ddc57e5SXin Li * the manual page source lmmin.pod and to the wiki.
11*5ddc57e5SXin Li *
12*5ddc57e5SXin Li * Author: Joachim Wuttke 2010, following a suggestion by Mario Rudolphi
13*5ddc57e5SXin Li *
14*5ddc57e5SXin Li * Licence: see ../COPYING (FreeBSD)
15*5ddc57e5SXin Li *
16*5ddc57e5SXin Li * Homepage: apps.jcns.fz-juelich.de/lmfit
17*5ddc57e5SXin Li */
18*5ddc57e5SXin Li
19*5ddc57e5SXin Li #include "lmmin.h"
20*5ddc57e5SXin Li #include <stdio.h>
21*5ddc57e5SXin Li
22*5ddc57e5SXin Li /* fit model: a plane p0 + p1*tx + p2*tz */
f(double tx,double tz,const double * p)23*5ddc57e5SXin Li double f( double tx, double tz, const double *p )
24*5ddc57e5SXin Li {
25*5ddc57e5SXin Li return p[0] + p[1]*tx + p[2]*tz;
26*5ddc57e5SXin Li }
27*5ddc57e5SXin Li
28*5ddc57e5SXin Li /* data structure to transmit arrays and fit model */
29*5ddc57e5SXin Li typedef struct {
30*5ddc57e5SXin Li double *tx, *tz;
31*5ddc57e5SXin Li double *y;
32*5ddc57e5SXin Li double (*f)( double tx, double tz, const double *p );
33*5ddc57e5SXin Li } data_struct;
34*5ddc57e5SXin Li
35*5ddc57e5SXin Li /* function evaluation, determination of residues */
evaluate_surface(const double * par,int m_dat,const void * data,double * fvec,int * info)36*5ddc57e5SXin Li void evaluate_surface( const double *par, int m_dat, const void *data,
37*5ddc57e5SXin Li double *fvec, int *info )
38*5ddc57e5SXin Li {
39*5ddc57e5SXin Li /* for readability, explicit type conversion */
40*5ddc57e5SXin Li data_struct *D;
41*5ddc57e5SXin Li D = (data_struct*)data;
42*5ddc57e5SXin Li
43*5ddc57e5SXin Li int i;
44*5ddc57e5SXin Li for ( i = 0; i < m_dat; i++ )
45*5ddc57e5SXin Li fvec[i] = D->y[i] - D->f( D->tx[i], D->tz[i], par );
46*5ddc57e5SXin Li }
47*5ddc57e5SXin Li
main()48*5ddc57e5SXin Li int main()
49*5ddc57e5SXin Li {
50*5ddc57e5SXin Li /* parameter vector */
51*5ddc57e5SXin Li int n_par = 3; /* number of parameters in model function f */
52*5ddc57e5SXin Li double par[3] = { -1, 0, 1 }; /* arbitrary starting value */
53*5ddc57e5SXin Li
54*5ddc57e5SXin Li /* data points */
55*5ddc57e5SXin Li int m_dat = 4;
56*5ddc57e5SXin Li double tx[4] = { -1, -1, 1, 1 };
57*5ddc57e5SXin Li double tz[4] = { -1, 1, -1, 1 };
58*5ddc57e5SXin Li double y[4] = { 0, 1, 1, 2 };
59*5ddc57e5SXin Li
60*5ddc57e5SXin Li data_struct data = { tx, tz, y, f };
61*5ddc57e5SXin Li
62*5ddc57e5SXin Li /* auxiliary parameters */
63*5ddc57e5SXin Li lm_status_struct status;
64*5ddc57e5SXin Li lm_control_struct control = lm_control_double;
65*5ddc57e5SXin Li control.verbosity = 9;
66*5ddc57e5SXin Li
67*5ddc57e5SXin Li /* perform the fit */
68*5ddc57e5SXin Li printf( "Fitting:\n" );
69*5ddc57e5SXin Li lmmin( n_par, par, m_dat, (const void*) &data,
70*5ddc57e5SXin Li evaluate_surface, &control, &status );
71*5ddc57e5SXin Li
72*5ddc57e5SXin Li /* print results */
73*5ddc57e5SXin Li printf( "\nResults:\n" );
74*5ddc57e5SXin Li printf( "status after %d function evaluations:\n %s\n",
75*5ddc57e5SXin Li status.nfev, lm_infmsg[status.outcome] );
76*5ddc57e5SXin Li
77*5ddc57e5SXin Li printf("obtained parameters:\n");
78*5ddc57e5SXin Li int i;
79*5ddc57e5SXin Li for ( i=0; i<n_par; ++i )
80*5ddc57e5SXin Li printf(" par[%i] = %12g\n", i, par[i]);
81*5ddc57e5SXin Li printf("obtained norm:\n %12g\n", status.fnorm );
82*5ddc57e5SXin Li
83*5ddc57e5SXin Li printf("fitting data as follows:\n");
84*5ddc57e5SXin Li double ff;
85*5ddc57e5SXin Li for ( i=0; i<m_dat; ++i ){
86*5ddc57e5SXin Li ff = f(tx[i], tz[i], par);
87*5ddc57e5SXin Li printf( " t[%2d]=%12g,%12g y=%12g fit=%12g residue=%12g\n",
88*5ddc57e5SXin Li i, tx[i], tz[i], y[i], ff, y[i] - ff );
89*5ddc57e5SXin Li }
90*5ddc57e5SXin Li
91*5ddc57e5SXin Li return 0;
92*5ddc57e5SXin Li }
93