1 | ////////////////////////////////////////////////////////////////////////////////////
|
---|
2 | // Example program that shows how to use levmar in order to fit the three-
|
---|
3 | // parameter exponential model x_i = p[0]*exp(-p[1]*i) + p[2] to a set of
|
---|
4 | // data measurements; example is based on a similar one from GSL.
|
---|
5 | //
|
---|
6 | // Copyright (C) 2008 Manolis Lourakis (lourakis at ics forth gr)
|
---|
7 | // Institute of Computer Science, Foundation for Research & Technology - Hellas
|
---|
8 | // Heraklion, Crete, Greece.
|
---|
9 | //
|
---|
10 | // This program is free software; you can redistribute it and/or modify
|
---|
11 | // it under the terms of the GNU General Public License as published by
|
---|
12 | // the Free Software Foundation; either version 2 of the License, or
|
---|
13 | // (at your option) any later version.
|
---|
14 | //
|
---|
15 | // This program is distributed in the hope that it will be useful,
|
---|
16 | // but WITHOUT ANY WARRANTY; without even the implied warranty of
|
---|
17 | // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
---|
18 | // GNU General Public License for more details.
|
---|
19 | //
|
---|
20 | ////////////////////////////////////////////////////////////////////////////////////
|
---|
21 |
|
---|
22 | #include <stdio.h>
|
---|
23 | #include <stdlib.h>
|
---|
24 | #include <math.h>
|
---|
25 |
|
---|
26 | #include <levmar.h>
|
---|
27 |
|
---|
28 | #ifndef LM_DBL_PREC
|
---|
29 | #error Example program assumes that levmar has been compiled with double precision, see LM_DBL_PREC!
|
---|
30 | #endif
|
---|
31 |
|
---|
32 |
|
---|
33 | /* the following macros concern the initialization of a random number generator for adding noise */
|
---|
34 | #undef REPEATABLE_RANDOM
|
---|
35 | #define DBL_RAND_MAX (double)(RAND_MAX)
|
---|
36 |
|
---|
37 | #ifdef _MSC_VER // MSVC
|
---|
38 | #include <process.h>
|
---|
39 | #define GETPID _getpid
|
---|
40 | #elif defined(__GNUC__) // GCC
|
---|
41 | #include <sys/types.h>
|
---|
42 | #include <unistd.h>
|
---|
43 | #define GETPID getpid
|
---|
44 | #else
|
---|
45 | #warning Do not know the name of the function returning the process id for your OS/compiler combination
|
---|
46 | #define GETPID 0
|
---|
47 | #endif /* _MSC_VER */
|
---|
48 |
|
---|
49 | #ifdef REPEATABLE_RANDOM
|
---|
50 | #define INIT_RANDOM(seed) srandom(seed)
|
---|
51 | #else
|
---|
52 | #define INIT_RANDOM(seed) srandom((int)GETPID()) // seed unused
|
---|
53 | #endif
|
---|
54 |
|
---|
55 | /* Gaussian noise with mean m and variance s, uses the Box-Muller transformation */
|
---|
56 | double gNoise(double m, double s)
|
---|
57 | {
|
---|
58 | double r1, r2, val;
|
---|
59 |
|
---|
60 | r1=((double)random())/DBL_RAND_MAX;
|
---|
61 | r2=((double)random())/DBL_RAND_MAX;
|
---|
62 |
|
---|
63 | val=sqrt(-2.0*log(r1))*cos(2.0*M_PI*r2);
|
---|
64 |
|
---|
65 | val=s*val+m;
|
---|
66 |
|
---|
67 | return val;
|
---|
68 | }
|
---|
69 |
|
---|
70 | /* model to be fitted to measurements: x_i = p[0]*exp(-p[1]*i) + p[2], i=0...n-1 */
|
---|
71 | void expfunc(double *p, double *x, int m, int n, void *data)
|
---|
72 | {
|
---|
73 | register int i;
|
---|
74 |
|
---|
75 | for(i=0; i<n; ++i){
|
---|
76 | x[i]=p[0]*exp(-p[1]*i) + p[2];
|
---|
77 | }
|
---|
78 | }
|
---|
79 |
|
---|
80 | /* Jacobian of expfunc() */
|
---|
81 | void jacexpfunc(double *p, double *jac, int m, int n, void *data)
|
---|
82 | {
|
---|
83 | register int i, j;
|
---|
84 |
|
---|
85 | /* fill Jacobian row by row */
|
---|
86 | for(i=j=0; i<n; ++i){
|
---|
87 | jac[j++]=exp(-p[1]*i);
|
---|
88 | jac[j++]=-p[0]*i*exp(-p[1]*i);
|
---|
89 | jac[j++]=1.0;
|
---|
90 | }
|
---|
91 | }
|
---|
92 |
|
---|
93 | int main()
|
---|
94 | {
|
---|
95 | const int n=40, m=3; // 40 measurements, 3 parameters
|
---|
96 | double p[m], x[n], opts[LM_OPTS_SZ], info[LM_INFO_SZ];
|
---|
97 | register int i;
|
---|
98 | int ret;
|
---|
99 |
|
---|
100 | /* generate some measurement using the exponential model with
|
---|
101 | * parameters (5.0, 0.1, 1.0), corrupted with zero-mean
|
---|
102 | * Gaussian noise of s=0.1
|
---|
103 | */
|
---|
104 | INIT_RANDOM(0);
|
---|
105 | for(i=0; i<n; ++i)
|
---|
106 | x[i]=(5.0*exp(-0.1*i) + 1.0) + gNoise(0.0, 0.1);
|
---|
107 |
|
---|
108 | /* initial parameters estimate: (1.0, 0.0, 0.0) */
|
---|
109 | p[0]=1.0; p[1]=0.0; p[2]=0.0;
|
---|
110 |
|
---|
111 | /* optimization control parameters; passing to levmar NULL instead of opts reverts to defaults */
|
---|
112 | opts[0]=LM_INIT_MU; opts[1]=1E-15; opts[2]=1E-15; opts[3]=1E-20;
|
---|
113 | opts[4]=LM_DIFF_DELTA; // relevant only if the finite difference Jacobian version is used
|
---|
114 |
|
---|
115 | /* invoke the optimization function */
|
---|
116 | ret=dlevmar_der(expfunc, jacexpfunc, p, x, m, n, 1000, opts, info, NULL, NULL, NULL); // with analytic Jacobian
|
---|
117 | //ret=dlevmar_dif(expfunc, p, x, m, n, 1000, opts, info, NULL, NULL, NULL); // without Jacobian
|
---|
118 | printf("Levenberg-Marquardt returned in %g iter, reason %g, sumsq %g [%g]\n", info[5], info[6], info[1], info[0]);
|
---|
119 | printf("Best fit parameters: %.7g %.7g %.7g\n", p[0], p[1], p[2]);
|
---|
120 |
|
---|
121 | exit(0);
|
---|
122 | }
|
---|