1 #include <stdlib.h>
2 #include <stdio.h>
3 #include <gsl/gsl_vector.h>
4 #include <gsl/gsl_matrix.h>
5 #include <gsl/gsl_blas.h>
6 #include <gsl/gsl_multifit_nlinear.h>
7 #include <gsl/gsl_rng.h>
8 #include <gsl/gsl_randist.h>
9
10 struct data
11 {
12 double *t;
13 double *y;
14 size_t n;
15 };
16
17 /* model function: a * exp( -1/2 * [ (t - b) / c ]^2 ) */
18 double
gaussian(const double a,const double b,const double c,const double t)19 gaussian(const double a, const double b, const double c, const double t)
20 {
21 const double z = (t - b) / c;
22 return (a * exp(-0.5 * z * z));
23 }
24
25 int
func_f(const gsl_vector * x,void * params,gsl_vector * f)26 func_f (const gsl_vector * x, void *params, gsl_vector * f)
27 {
28 struct data *d = (struct data *) params;
29 double a = gsl_vector_get(x, 0);
30 double b = gsl_vector_get(x, 1);
31 double c = gsl_vector_get(x, 2);
32 size_t i;
33
34 for (i = 0; i < d->n; ++i)
35 {
36 double ti = d->t[i];
37 double yi = d->y[i];
38 double y = gaussian(a, b, c, ti);
39
40 gsl_vector_set(f, i, yi - y);
41 }
42
43 return GSL_SUCCESS;
44 }
45
46 int
func_df(const gsl_vector * x,void * params,gsl_matrix * J)47 func_df (const gsl_vector * x, void *params, gsl_matrix * J)
48 {
49 struct data *d = (struct data *) params;
50 double a = gsl_vector_get(x, 0);
51 double b = gsl_vector_get(x, 1);
52 double c = gsl_vector_get(x, 2);
53 size_t i;
54
55 for (i = 0; i < d->n; ++i)
56 {
57 double ti = d->t[i];
58 double zi = (ti - b) / c;
59 double ei = exp(-0.5 * zi * zi);
60
61 gsl_matrix_set(J, i, 0, -ei);
62 gsl_matrix_set(J, i, 1, -(a / c) * ei * zi);
63 gsl_matrix_set(J, i, 2, -(a / c) * ei * zi * zi);
64 }
65
66 return GSL_SUCCESS;
67 }
68
69 int
func_fvv(const gsl_vector * x,const gsl_vector * v,void * params,gsl_vector * fvv)70 func_fvv (const gsl_vector * x, const gsl_vector * v,
71 void *params, gsl_vector * fvv)
72 {
73 struct data *d = (struct data *) params;
74 double a = gsl_vector_get(x, 0);
75 double b = gsl_vector_get(x, 1);
76 double c = gsl_vector_get(x, 2);
77 double va = gsl_vector_get(v, 0);
78 double vb = gsl_vector_get(v, 1);
79 double vc = gsl_vector_get(v, 2);
80 size_t i;
81
82 for (i = 0; i < d->n; ++i)
83 {
84 double ti = d->t[i];
85 double zi = (ti - b) / c;
86 double ei = exp(-0.5 * zi * zi);
87 double Dab = -zi * ei / c;
88 double Dac = -zi * zi * ei / c;
89 double Dbb = a * ei / (c * c) * (1.0 - zi*zi);
90 double Dbc = a * zi * ei / (c * c) * (2.0 - zi*zi);
91 double Dcc = a * zi * zi * ei / (c * c) * (3.0 - zi*zi);
92 double sum;
93
94 sum = 2.0 * va * vb * Dab +
95 2.0 * va * vc * Dac +
96 vb * vb * Dbb +
97 2.0 * vb * vc * Dbc +
98 vc * vc * Dcc;
99
100 gsl_vector_set(fvv, i, sum);
101 }
102
103 return GSL_SUCCESS;
104 }
105
106 void
callback(const size_t iter,void * params,const gsl_multifit_nlinear_workspace * w)107 callback(const size_t iter, void *params,
108 const gsl_multifit_nlinear_workspace *w)
109 {
110 gsl_vector *f = gsl_multifit_nlinear_residual(w);
111 gsl_vector *x = gsl_multifit_nlinear_position(w);
112 double avratio = gsl_multifit_nlinear_avratio(w);
113 double rcond;
114
115 (void) params; /* not used */
116
117 /* compute reciprocal condition number of J(x) */
118 gsl_multifit_nlinear_rcond(&rcond, w);
119
120 fprintf(stderr, "iter %2zu: a = %.4f, b = %.4f, c = %.4f, |a|/|v| = %.4f cond(J) = %8.4f, |f(x)| = %.4f\n",
121 iter,
122 gsl_vector_get(x, 0),
123 gsl_vector_get(x, 1),
124 gsl_vector_get(x, 2),
125 avratio,
126 1.0 / rcond,
127 gsl_blas_dnrm2(f));
128 }
129
130 void
solve_system(gsl_vector * x,gsl_multifit_nlinear_fdf * fdf,gsl_multifit_nlinear_parameters * params)131 solve_system(gsl_vector *x, gsl_multifit_nlinear_fdf *fdf,
132 gsl_multifit_nlinear_parameters *params)
133 {
134 const gsl_multifit_nlinear_type *T = gsl_multifit_nlinear_trust;
135 const size_t max_iter = 200;
136 const double xtol = 1.0e-8;
137 const double gtol = 1.0e-8;
138 const double ftol = 1.0e-8;
139 const size_t n = fdf->n;
140 const size_t p = fdf->p;
141 gsl_multifit_nlinear_workspace *work =
142 gsl_multifit_nlinear_alloc(T, params, n, p);
143 gsl_vector * f = gsl_multifit_nlinear_residual(work);
144 gsl_vector * y = gsl_multifit_nlinear_position(work);
145 int info;
146 double chisq0, chisq, rcond;
147
148 /* initialize solver */
149 gsl_multifit_nlinear_init(x, fdf, work);
150
151 /* store initial cost */
152 gsl_blas_ddot(f, f, &chisq0);
153
154 /* iterate until convergence */
155 gsl_multifit_nlinear_driver(max_iter, xtol, gtol, ftol,
156 callback, NULL, &info, work);
157
158 /* store final cost */
159 gsl_blas_ddot(f, f, &chisq);
160
161 /* store cond(J(x)) */
162 gsl_multifit_nlinear_rcond(&rcond, work);
163
164 gsl_vector_memcpy(x, y);
165
166 /* print summary */
167
168 fprintf(stderr, "NITER = %zu\n", gsl_multifit_nlinear_niter(work));
169 fprintf(stderr, "NFEV = %zu\n", fdf->nevalf);
170 fprintf(stderr, "NJEV = %zu\n", fdf->nevaldf);
171 fprintf(stderr, "NAEV = %zu\n", fdf->nevalfvv);
172 fprintf(stderr, "initial cost = %.12e\n", chisq0);
173 fprintf(stderr, "final cost = %.12e\n", chisq);
174 fprintf(stderr, "final x = (%.12e, %.12e, %12e)\n",
175 gsl_vector_get(x, 0), gsl_vector_get(x, 1), gsl_vector_get(x, 2));
176 fprintf(stderr, "final cond(J) = %.12e\n", 1.0 / rcond);
177
178 gsl_multifit_nlinear_free(work);
179 }
180
181 int
main(void)182 main (void)
183 {
184 const size_t n = 300; /* number of data points to fit */
185 const size_t p = 3; /* number of model parameters */
186 const double a = 5.0; /* amplitude */
187 const double b = 0.4; /* center */
188 const double c = 0.15; /* width */
189 const gsl_rng_type * T = gsl_rng_default;
190 gsl_vector *f = gsl_vector_alloc(n);
191 gsl_vector *x = gsl_vector_alloc(p);
192 gsl_multifit_nlinear_fdf fdf;
193 gsl_multifit_nlinear_parameters fdf_params =
194 gsl_multifit_nlinear_default_parameters();
195 struct data fit_data;
196 gsl_rng * r;
197 size_t i;
198
199 gsl_rng_env_setup ();
200 r = gsl_rng_alloc (T);
201
202 fit_data.t = malloc(n * sizeof(double));
203 fit_data.y = malloc(n * sizeof(double));
204 fit_data.n = n;
205
206 /* generate synthetic data with noise */
207 for (i = 0; i < n; ++i)
208 {
209 double t = (double)i / (double) n;
210 double y0 = gaussian(a, b, c, t);
211 double dy = gsl_ran_gaussian (r, 0.1 * y0);
212
213 fit_data.t[i] = t;
214 fit_data.y[i] = y0 + dy;
215 }
216
217 /* define function to be minimized */
218 fdf.f = func_f;
219 fdf.df = func_df;
220 fdf.fvv = func_fvv;
221 fdf.n = n;
222 fdf.p = p;
223 fdf.params = &fit_data;
224
225 /* starting point */
226 gsl_vector_set(x, 0, 1.0);
227 gsl_vector_set(x, 1, 0.0);
228 gsl_vector_set(x, 2, 1.0);
229
230 fdf_params.trs = gsl_multifit_nlinear_trs_lmaccel;
231 solve_system(x, &fdf, &fdf_params);
232
233 /* print data and model */
234 {
235 double A = gsl_vector_get(x, 0);
236 double B = gsl_vector_get(x, 1);
237 double C = gsl_vector_get(x, 2);
238
239 for (i = 0; i < n; ++i)
240 {
241 double ti = fit_data.t[i];
242 double yi = fit_data.y[i];
243 double fi = gaussian(A, B, C, ti);
244
245 printf("%f %f %f\n", ti, yi, fi);
246 }
247 }
248
249 gsl_vector_free(f);
250 gsl_vector_free(x);
251 gsl_rng_free(r);
252
253 return 0;
254 }
255