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 */ 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 */ 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 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