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