1 /* multifit/fdfsolver.c
2 *
3 * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 Brian Gough
4 *
5 * This program is free software; you can redistribute it and/or modify
6 * it under the terms of the GNU General Public License as published by
7 * the Free Software Foundation; either version 3 of the License, or (at
8 * your option) any later version.
9 *
10 * This program is distributed in the hope that it will be useful, but
11 * WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
13 * General Public License for more details.
14 *
15 * You should have received a copy of the GNU General Public License
16 * along with this program; if not, write to the Free Software
17 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
18 */
19
20 #include <config.h>
21 #include <stdlib.h>
22 #include <string.h>
23 #include <gsl/gsl_math.h>
24 #include <gsl/gsl_errno.h>
25 #include <gsl/gsl_multifit_nlin.h>
26
27 gsl_multifit_fdfsolver *
gsl_multifit_fdfsolver_alloc(const gsl_multifit_fdfsolver_type * T,size_t n,size_t p)28 gsl_multifit_fdfsolver_alloc (const gsl_multifit_fdfsolver_type * T,
29 size_t n, size_t p)
30 {
31 int status;
32
33 gsl_multifit_fdfsolver * s;
34
35 if (n < p)
36 {
37 GSL_ERROR_VAL ("insufficient data points, n < p", GSL_EINVAL, 0);
38 }
39
40 s = (gsl_multifit_fdfsolver *) calloc (1, sizeof (gsl_multifit_fdfsolver));
41 if (s == 0)
42 {
43 GSL_ERROR_VAL ("failed to allocate space for multifit solver struct",
44 GSL_ENOMEM, 0);
45 }
46
47 s->x = gsl_vector_calloc (p);
48
49 if (s->x == 0)
50 {
51 gsl_multifit_fdfsolver_free (s);
52 GSL_ERROR_VAL ("failed to allocate space for x", GSL_ENOMEM, 0);
53 }
54
55 s->f = gsl_vector_calloc (n);
56
57 if (s->f == 0)
58 {
59 gsl_multifit_fdfsolver_free (s);
60 GSL_ERROR_VAL ("failed to allocate space for f", GSL_ENOMEM, 0);
61 }
62
63 s->dx = gsl_vector_calloc (p);
64
65 if (s->dx == 0)
66 {
67 gsl_multifit_fdfsolver_free (s);
68 GSL_ERROR_VAL ("failed to allocate space for dx", GSL_ENOMEM, 0);
69 }
70
71 s->g = gsl_vector_alloc (p);
72
73 if (s->g == 0)
74 {
75 gsl_multifit_fdfsolver_free (s);
76 GSL_ERROR_VAL ("failed to allocate space for g", GSL_ENOMEM, 0);
77 }
78
79 s->sqrt_wts = gsl_vector_calloc (n);
80
81 if (s->sqrt_wts == 0)
82 {
83 gsl_multifit_fdfsolver_free (s);
84 GSL_ERROR_VAL ("failed to allocate space for sqrt_wts", GSL_ENOMEM, 0);
85 }
86
87 s->state = calloc (1, T->size);
88
89 if (s->state == 0)
90 {
91 gsl_multifit_fdfsolver_free (s);
92 GSL_ERROR_VAL ("failed to allocate space for multifit solver state",
93 GSL_ENOMEM, 0);
94 }
95
96 s->type = T ;
97
98 status = (s->type->alloc)(s->state, n, p);
99
100 if (status != GSL_SUCCESS)
101 {
102 gsl_multifit_fdfsolver_free (s);
103 GSL_ERROR_VAL ("failed to set solver", status, 0);
104 }
105
106 s->fdf = NULL;
107
108 s->niter = 0;
109
110 return s;
111 }
112
113 int
gsl_multifit_fdfsolver_set(gsl_multifit_fdfsolver * s,gsl_multifit_function_fdf * f,const gsl_vector * x)114 gsl_multifit_fdfsolver_set (gsl_multifit_fdfsolver * s,
115 gsl_multifit_function_fdf * f,
116 const gsl_vector * x)
117 {
118 return gsl_multifit_fdfsolver_wset(s, f, x, NULL);
119 }
120
121 int
gsl_multifit_fdfsolver_wset(gsl_multifit_fdfsolver * s,gsl_multifit_function_fdf * f,const gsl_vector * x,const gsl_vector * wts)122 gsl_multifit_fdfsolver_wset (gsl_multifit_fdfsolver * s,
123 gsl_multifit_function_fdf * f,
124 const gsl_vector * x,
125 const gsl_vector * wts)
126 {
127 const size_t n = s->f->size;
128
129 if (n != f->n)
130 {
131 GSL_ERROR ("function size does not match solver", GSL_EBADLEN);
132 }
133 else if (s->x->size != x->size)
134 {
135 GSL_ERROR ("vector length does not match solver", GSL_EBADLEN);
136 }
137 else if (wts != NULL && n != wts->size)
138 {
139 GSL_ERROR ("weight vector length does not match solver", GSL_EBADLEN);
140 }
141 else
142 {
143 size_t i;
144
145 s->fdf = f;
146 gsl_vector_memcpy(s->x, x);
147 s->niter = 0;
148
149 if (wts)
150 {
151 for (i = 0; i < n; ++i)
152 {
153 double wi = gsl_vector_get(wts, i);
154 gsl_vector_set(s->sqrt_wts, i, sqrt(wi));
155 }
156 }
157 else
158 gsl_vector_set_all(s->sqrt_wts, 1.0);
159
160 return (s->type->set) (s->state, s->sqrt_wts, s->fdf, s->x, s->f, s->dx);
161 }
162 }
163
164 int
gsl_multifit_fdfsolver_iterate(gsl_multifit_fdfsolver * s)165 gsl_multifit_fdfsolver_iterate (gsl_multifit_fdfsolver * s)
166 {
167 int status =
168 (s->type->iterate) (s->state, s->sqrt_wts, s->fdf, s->x, s->f, s->dx);
169
170 s->niter++;
171
172 return status;
173 }
174
175 /*
176 gsl_multifit_fdfsolver_driver()
177 Iterate the nonlinear least squares solver until completion
178
179 Inputs: s - fdfsolver
180 maxiter - maximum iterations to allow
181 xtol - tolerance in step x
182 gtol - tolerance in gradient
183 ftol - tolerance in ||f||
184 info - (output) info flag on why iteration terminated
185 1 = stopped due to small step size ||dx|
186 2 = stopped due to small gradient
187 3 = stopped due to small change in f
188 GSL_ETOLX = ||dx|| has converged to within machine
189 precision (and xtol is too small)
190 GSL_ETOLG = ||g||_inf is smaller than machine
191 precision (gtol is too small)
192 GSL_ETOLF = change in ||f|| is smaller than machine
193 precision (ftol is too small)
194
195 Return: GSL_SUCCESS if converged, GSL_MAXITER if maxiter exceeded without
196 converging
197 */
198 int
gsl_multifit_fdfsolver_driver(gsl_multifit_fdfsolver * s,const size_t maxiter,const double xtol,const double gtol,const double ftol,int * info)199 gsl_multifit_fdfsolver_driver (gsl_multifit_fdfsolver * s,
200 const size_t maxiter,
201 const double xtol,
202 const double gtol,
203 const double ftol,
204 int *info)
205 {
206 int status;
207 size_t iter = 0;
208
209 do
210 {
211 status = gsl_multifit_fdfsolver_iterate (s);
212
213 /*
214 * if status is GSL_ENOPROG or GSL_SUCCESS, continue iterating,
215 * otherwise the method has converged with a GSL_ETOLx flag
216 */
217 if (status != GSL_SUCCESS && status != GSL_ENOPROG)
218 break;
219
220 /* test for convergence */
221 status = gsl_multifit_fdfsolver_test(s, xtol, gtol, ftol, info);
222 }
223 while (status == GSL_CONTINUE && ++iter < maxiter);
224
225 /*
226 * the following error codes mean that the solution has converged
227 * to within machine precision, so record the error code in info
228 * and return success
229 */
230 if (status == GSL_ETOLF || status == GSL_ETOLX || status == GSL_ETOLG)
231 {
232 *info = status;
233 status = GSL_SUCCESS;
234 }
235
236 /* check if max iterations reached */
237 if (iter >= maxiter && status != GSL_SUCCESS)
238 status = GSL_EMAXITER;
239
240 return status;
241 } /* gsl_multifit_fdfsolver_driver() */
242
243 int
gsl_multifit_fdfsolver_jac(gsl_multifit_fdfsolver * s,gsl_matrix * J)244 gsl_multifit_fdfsolver_jac (gsl_multifit_fdfsolver * s, gsl_matrix * J)
245 {
246 const size_t n = s->f->size;
247 const size_t p = s->x->size;
248
249 if (n != J->size1 || p != J->size2)
250 {
251 GSL_ERROR ("Jacobian dimensions do not match solver", GSL_EBADLEN);
252 }
253 else
254 {
255 return (s->type->jac) (s->state, J);
256 }
257 } /* gsl_multifit_fdfsolver_jac() */
258
259 void
gsl_multifit_fdfsolver_free(gsl_multifit_fdfsolver * s)260 gsl_multifit_fdfsolver_free (gsl_multifit_fdfsolver * s)
261 {
262 RETURN_IF_NULL (s);
263
264 if (s->state)
265 {
266 (s->type->free) (s->state);
267 free (s->state);
268 }
269
270 if (s->dx)
271 gsl_vector_free (s->dx);
272
273 if (s->x)
274 gsl_vector_free (s->x);
275
276 if (s->f)
277 gsl_vector_free (s->f);
278
279 if (s->sqrt_wts)
280 gsl_vector_free (s->sqrt_wts);
281
282 if (s->g)
283 gsl_vector_free (s->g);
284
285 free (s);
286 }
287
288 const char *
gsl_multifit_fdfsolver_name(const gsl_multifit_fdfsolver * s)289 gsl_multifit_fdfsolver_name (const gsl_multifit_fdfsolver * s)
290 {
291 return s->type->name;
292 }
293
294 gsl_vector *
gsl_multifit_fdfsolver_position(const gsl_multifit_fdfsolver * s)295 gsl_multifit_fdfsolver_position (const gsl_multifit_fdfsolver * s)
296 {
297 return s->x;
298 }
299
300 gsl_vector *
gsl_multifit_fdfsolver_residual(const gsl_multifit_fdfsolver * s)301 gsl_multifit_fdfsolver_residual (const gsl_multifit_fdfsolver * s)
302 {
303 return s->f;
304 }
305
306 size_t
gsl_multifit_fdfsolver_niter(const gsl_multifit_fdfsolver * s)307 gsl_multifit_fdfsolver_niter (const gsl_multifit_fdfsolver * s)
308 {
309 return s->niter;
310 }
311
312 /*
313 gsl_multifit_eval_wf()
314 Compute residual vector y with user callback function, and apply
315 weighting transform if given:
316
317 y~ = sqrt(W) y
318
319 Inputs: fdf - callback function
320 x - model parameters
321 swts - weight matrix sqrt(W) = sqrt(diag(w1,w2,...,wn))
322 set to NULL for unweighted fit
323 y - (output) (weighted) residual vector
324 y_i = sqrt(w_i) f_i where f_i is unweighted residual
325 */
326
327 int
gsl_multifit_eval_wf(gsl_multifit_function_fdf * fdf,const gsl_vector * x,const gsl_vector * swts,gsl_vector * y)328 gsl_multifit_eval_wf(gsl_multifit_function_fdf *fdf, const gsl_vector *x,
329 const gsl_vector *swts, gsl_vector *y)
330 {
331 int s = ((*((fdf)->f)) (x, fdf->params, y));
332 ++(fdf->nevalf);
333
334 /* y <- sqrt(W) y */
335 if (swts)
336 gsl_vector_mul(y, swts);
337
338 return s;
339 }
340
341 /*
342 gsl_multifit_eval_wdf()
343 Compute Jacobian matrix J with user callback function, and apply
344 weighting transform if given:
345
346 J~ = sqrt(W) J
347
348 Inputs: fdf - callback function
349 x - model parameters
350 swts - weight matrix W = diag(w1,w2,...,wn)
351 set to NULL for unweighted fit
352 dy - (output) (weighted) Jacobian matrix
353 dy = sqrt(W) dy where dy is unweighted Jacobian
354 */
355
356 int
gsl_multifit_eval_wdf(gsl_multifit_function_fdf * fdf,const gsl_vector * x,const gsl_vector * swts,gsl_matrix * dy)357 gsl_multifit_eval_wdf(gsl_multifit_function_fdf *fdf, const gsl_vector *x,
358 const gsl_vector *swts, gsl_matrix *dy)
359 {
360 int status = ((*((fdf)->df)) (x, fdf->params, dy));
361
362 ++(fdf->nevaldf);
363
364 /* J <- sqrt(W) J */
365 if (swts)
366 {
367 const size_t n = swts->size;
368 size_t i;
369
370 for (i = 0; i < n; ++i)
371 {
372 double swi = gsl_vector_get(swts, i);
373 gsl_vector_view v = gsl_matrix_row(dy, i);
374
375 gsl_vector_scale(&v.vector, swi);
376 }
377 }
378
379 return status;
380 }
381