1 /* multifit/fdfridge.c
2 *
3 * Copyright (C) 2014 Patrick Alken
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 <gsl/gsl_math.h>
23 #include <gsl/gsl_errno.h>
24 #include <gsl/gsl_vector.h>
25 #include <gsl/gsl_matrix.h>
26 #include <gsl/gsl_multifit_nlin.h>
27 #include <gsl/gsl_blas.h>
28
29 static int fdfridge_f(const gsl_vector * x, void * params, gsl_vector * f);
30 static int fdfridge_df(const gsl_vector * x, void * params, gsl_matrix * J);
31
32 gsl_multifit_fdfridge *
gsl_multifit_fdfridge_alloc(const gsl_multifit_fdfsolver_type * T,const size_t n,const size_t p)33 gsl_multifit_fdfridge_alloc (const gsl_multifit_fdfsolver_type * T,
34 const size_t n, const size_t p)
35 {
36 gsl_multifit_fdfridge * work;
37
38 work = calloc(1, sizeof(gsl_multifit_fdfridge));
39 if (work == NULL)
40 {
41 GSL_ERROR_VAL("failed to allocate workspace",
42 GSL_ENOMEM, 0);
43 }
44
45 work->s = gsl_multifit_fdfsolver_alloc (T, n + p, p);
46 if (work->s == NULL)
47 {
48 gsl_multifit_fdfridge_free(work);
49 GSL_ERROR_VAL("failed to allocate space for fdfsolver",
50 GSL_ENOMEM, 0);
51 }
52
53 work->wts = gsl_vector_alloc(n + p);
54 if (work->wts == NULL)
55 {
56 gsl_multifit_fdfridge_free(work);
57 GSL_ERROR_VAL("failed to allocate space for weight vector",
58 GSL_ENOMEM, 0);
59 }
60
61 work->f = gsl_vector_alloc(n);
62 if (work->f == NULL)
63 {
64 gsl_multifit_fdfridge_free(work);
65 GSL_ERROR_VAL("failed to allocate space for f vector",
66 GSL_ENOMEM, 0);
67 }
68
69 work->n = n;
70 work->p = p;
71 work->lambda = 0.0;
72
73 /* initialize weights to 1 (for augmented portion of vector) */
74 gsl_vector_set_all(work->wts, 1.0);
75
76 return work;
77 } /* gsl_multifit_fdfridge_alloc() */
78
79 void
gsl_multifit_fdfridge_free(gsl_multifit_fdfridge * work)80 gsl_multifit_fdfridge_free(gsl_multifit_fdfridge *work)
81 {
82 if (work->s)
83 gsl_multifit_fdfsolver_free(work->s);
84
85 if (work->wts)
86 gsl_vector_free(work->wts);
87
88 if (work->f)
89 gsl_vector_free(work->f);
90
91 free(work);
92 }
93
94 const char *
gsl_multifit_fdfridge_name(const gsl_multifit_fdfridge * w)95 gsl_multifit_fdfridge_name(const gsl_multifit_fdfridge * w)
96 {
97 return gsl_multifit_fdfsolver_name(w->s);
98 }
99
100 gsl_vector *
gsl_multifit_fdfridge_position(const gsl_multifit_fdfridge * w)101 gsl_multifit_fdfridge_position (const gsl_multifit_fdfridge * w)
102 {
103 return gsl_multifit_fdfsolver_position(w->s);
104 }
105
106 gsl_vector *
gsl_multifit_fdfridge_residual(const gsl_multifit_fdfridge * w)107 gsl_multifit_fdfridge_residual (const gsl_multifit_fdfridge * w)
108 {
109 return gsl_multifit_fdfsolver_residual(w->s);
110 }
111
112 size_t
gsl_multifit_fdfridge_niter(const gsl_multifit_fdfridge * w)113 gsl_multifit_fdfridge_niter (const gsl_multifit_fdfridge * w)
114 {
115 return w->s->niter;
116 }
117
118 int
gsl_multifit_fdfridge_set(gsl_multifit_fdfridge * w,gsl_multifit_function_fdf * f,const gsl_vector * x,const double lambda)119 gsl_multifit_fdfridge_set (gsl_multifit_fdfridge * w,
120 gsl_multifit_function_fdf * f,
121 const gsl_vector * x,
122 const double lambda)
123 {
124 return gsl_multifit_fdfridge_wset(w, f, x, lambda, NULL);
125 } /* gsl_multifit_fdfridge_set() */
126
127 int
gsl_multifit_fdfridge_wset(gsl_multifit_fdfridge * w,gsl_multifit_function_fdf * f,const gsl_vector * x,const double lambda,const gsl_vector * wts)128 gsl_multifit_fdfridge_wset (gsl_multifit_fdfridge * w,
129 gsl_multifit_function_fdf * f,
130 const gsl_vector * x,
131 const double lambda,
132 const gsl_vector * wts)
133 {
134 const size_t n = w->n;
135 const size_t p = w->p;
136
137 if (n != f->n || p != f->p)
138 {
139 GSL_ERROR ("function size does not match solver", GSL_EBADLEN);
140 }
141 else if (p != x->size)
142 {
143 GSL_ERROR ("vector length does not match solver", GSL_EBADLEN);
144 }
145 else if (wts != NULL && n != wts->size)
146 {
147 GSL_ERROR ("weight vector length does not match solver", GSL_EBADLEN);
148 }
149 else
150 {
151 int status;
152 gsl_vector_view wv = gsl_vector_subvector(w->wts, 0, n);
153
154 /* save user defined fdf */
155 w->fdf = f;
156
157 /* build modified fdf for Tikhonov terms */
158 w->fdftik.f = &fdfridge_f;
159 w->fdftik.df = &fdfridge_df;
160 w->fdftik.n = n + p; /* add p for Tikhonov terms */
161 w->fdftik.p = p;
162 w->fdftik.params = (void *) w;
163
164 /* store damping parameter */
165 w->lambda = lambda;
166 w->L = NULL;
167
168 if (wts)
169 {
170 /* copy weight vector into user portion of w->wts */
171 gsl_vector_memcpy(&wv.vector, wts);
172 status = gsl_multifit_fdfsolver_wset(w->s, &(w->fdftik), x, w->wts);
173 }
174 else
175 {
176 status = gsl_multifit_fdfsolver_wset(w->s, &(w->fdftik), x, NULL);
177 }
178
179 /* update function/Jacobian evaluations */
180 f->nevalf = w->fdftik.nevalf;
181 f->nevaldf = w->fdftik.nevaldf;
182
183 return status;
184 }
185 } /* gsl_multifit_fdfridge_wset() */
186
187 int
gsl_multifit_fdfridge_set2(gsl_multifit_fdfridge * w,gsl_multifit_function_fdf * f,const gsl_vector * x,const gsl_vector * lambda)188 gsl_multifit_fdfridge_set2 (gsl_multifit_fdfridge * w,
189 gsl_multifit_function_fdf * f,
190 const gsl_vector * x,
191 const gsl_vector * lambda)
192 {
193 return gsl_multifit_fdfridge_wset2(w, f, x, lambda, NULL);
194 } /* gsl_multifit_fdfridge_set2() */
195
196 int
gsl_multifit_fdfridge_wset2(gsl_multifit_fdfridge * w,gsl_multifit_function_fdf * f,const gsl_vector * x,const gsl_vector * lambda,const gsl_vector * wts)197 gsl_multifit_fdfridge_wset2 (gsl_multifit_fdfridge * w,
198 gsl_multifit_function_fdf * f,
199 const gsl_vector * x,
200 const gsl_vector * lambda,
201 const gsl_vector * wts)
202 {
203 const size_t n = w->n;
204 const size_t p = w->p;
205
206 if (n != f->n || p != f->p)
207 {
208 GSL_ERROR ("function size does not match solver", GSL_EBADLEN);
209 }
210 else if (p != x->size)
211 {
212 GSL_ERROR ("vector length does not match solver", GSL_EBADLEN);
213 }
214 else if (lambda->size != p)
215 {
216 GSL_ERROR ("lambda vector size does not match solver", GSL_EBADLEN);
217 }
218 else if (wts != NULL && n != wts->size)
219 {
220 GSL_ERROR ("weight vector length does not match solver", GSL_EBADLEN);
221 }
222 else
223 {
224 int status;
225 gsl_vector_view wv = gsl_vector_subvector(w->wts, 0, n);
226
227 /* save user defined fdf */
228 w->fdf = f;
229 w->fdf->nevalf = 0;
230 w->fdf->nevaldf = 0;
231
232 /* build modified fdf for Tikhonov terms */
233 w->fdftik.f = &fdfridge_f;
234 w->fdftik.df = &fdfridge_df;
235 w->fdftik.n = n + p; /* add p for Tikhonov terms */
236 w->fdftik.p = p;
237 w->fdftik.params = (void *) w;
238
239 /* store damping matrix */
240 w->lambda = 0.0;
241 w->L_diag = lambda;
242 w->L = NULL;
243
244 if (wts)
245 {
246 /* copy weight vector into user portion */
247 gsl_vector_memcpy(&wv.vector, wts);
248 status = gsl_multifit_fdfsolver_wset(w->s, &(w->fdftik), x, w->wts);
249 }
250 else
251 {
252 status = gsl_multifit_fdfsolver_wset(w->s, &(w->fdftik), x, NULL);
253 }
254
255 /* update function/Jacobian evaluations */
256 f->nevalf = w->fdftik.nevalf;
257 f->nevaldf = w->fdftik.nevaldf;
258
259 return status;
260 }
261 } /* gsl_multifit_fdfridge_wset2() */
262
263 int
gsl_multifit_fdfridge_set3(gsl_multifit_fdfridge * w,gsl_multifit_function_fdf * f,const gsl_vector * x,const gsl_matrix * L)264 gsl_multifit_fdfridge_set3 (gsl_multifit_fdfridge * w,
265 gsl_multifit_function_fdf * f,
266 const gsl_vector * x,
267 const gsl_matrix * L)
268 {
269 return gsl_multifit_fdfridge_wset3(w, f, x, L, NULL);
270 } /* gsl_multifit_fdfridge_set3() */
271
272 int
gsl_multifit_fdfridge_wset3(gsl_multifit_fdfridge * w,gsl_multifit_function_fdf * f,const gsl_vector * x,const gsl_matrix * L,const gsl_vector * wts)273 gsl_multifit_fdfridge_wset3 (gsl_multifit_fdfridge * w,
274 gsl_multifit_function_fdf * f,
275 const gsl_vector * x,
276 const gsl_matrix * L,
277 const gsl_vector * wts)
278 {
279 const size_t n = w->n;
280 const size_t p = w->p;
281
282 if (n != f->n || p != f->p)
283 {
284 GSL_ERROR ("function size does not match solver", GSL_EBADLEN);
285 }
286 else if (p != x->size)
287 {
288 GSL_ERROR ("vector length does not match solver", GSL_EBADLEN);
289 }
290 else if (L->size2 != p)
291 {
292 GSL_ERROR ("L matrix size2 does not match solver", GSL_EBADLEN);
293 }
294 else if (wts != NULL && n != wts->size)
295 {
296 GSL_ERROR ("weight vector length does not match solver", GSL_EBADLEN);
297 }
298 else
299 {
300 int status;
301 gsl_vector_view wv = gsl_vector_subvector(w->wts, 0, n);
302
303 /* save user defined fdf */
304 w->fdf = f;
305 w->fdf->nevalf = 0;
306 w->fdf->nevaldf = 0;
307
308 /* build modified fdf for Tikhonov terms */
309 w->fdftik.f = &fdfridge_f;
310 w->fdftik.df = &fdfridge_df;
311 w->fdftik.n = n + p; /* add p for Tikhonov terms */
312 w->fdftik.p = p;
313 w->fdftik.params = (void *) w;
314
315 /* store damping matrix */
316 w->lambda = 0.0;
317 w->L_diag = NULL;
318 w->L = L;
319
320 if (wts)
321 {
322 /* copy weight vector into user portion */
323 gsl_vector_memcpy(&wv.vector, wts);
324 status = gsl_multifit_fdfsolver_wset(w->s, &(w->fdftik), x, w->wts);
325 }
326 else
327 {
328 status = gsl_multifit_fdfsolver_wset(w->s, &(w->fdftik), x, NULL);
329 }
330
331 /* update function/Jacobian evaluations */
332 f->nevalf = w->fdftik.nevalf;
333 f->nevaldf = w->fdftik.nevaldf;
334
335 return status;
336 }
337 } /* gsl_multifit_fdfridge_wset3() */
338
339 int
gsl_multifit_fdfridge_iterate(gsl_multifit_fdfridge * w)340 gsl_multifit_fdfridge_iterate (gsl_multifit_fdfridge * w)
341 {
342 int status = gsl_multifit_fdfsolver_iterate(w->s);
343
344 /* update function/Jacobian evaluations */
345 w->fdf->nevalf = w->fdftik.nevalf;
346 w->fdf->nevaldf = w->fdftik.nevaldf;
347
348 return status;
349 }
350
351 int
gsl_multifit_fdfridge_driver(gsl_multifit_fdfridge * w,const size_t maxiter,const double xtol,const double gtol,const double ftol,int * info)352 gsl_multifit_fdfridge_driver (gsl_multifit_fdfridge * w,
353 const size_t maxiter,
354 const double xtol,
355 const double gtol,
356 const double ftol,
357 int *info)
358 {
359 int status = gsl_multifit_fdfsolver_driver(w->s, maxiter, xtol,
360 gtol, ftol, info);
361 return status;
362 } /* gsl_multifit_fdfridge_driver() */
363
364 /*
365 fdfridge_f()
366 Callback function to provide residuals, including extra p
367 Tikhonov terms. The residual vector will have the form:
368
369 f~ = [ f ]
370 [ \lambda x ]
371
372 where f is the user supplied residuals, x are the model
373 parameters, and \lambda is the Tikhonov damping parameter
374
375 Inputs: x - model parameters (size p)
376 params - pointer to fdfridge workspace
377 f - (output) (n+p) vector to store f~
378 */
379
380 static int
fdfridge_f(const gsl_vector * x,void * params,gsl_vector * f)381 fdfridge_f(const gsl_vector * x, void * params, gsl_vector * f)
382 {
383 int status;
384 gsl_multifit_fdfridge *w = (gsl_multifit_fdfridge *) params;
385 const size_t n = w->n;
386 const size_t p = w->p;
387 gsl_vector_view f_user = gsl_vector_subvector(f, 0, n);
388 gsl_vector_view f_tik = gsl_vector_subvector(f, n, p);
389
390 /* call user callback function to get residual vector f */
391 status = gsl_multifit_eval_wf(w->fdf, x, NULL, &f_user.vector);
392 if (status)
393 return status;
394
395 if (w->L_diag)
396 {
397 /* store diag(L_diag) x in Tikhonov portion of f~ */
398 gsl_vector_memcpy(&f_tik.vector, x);
399 gsl_vector_mul(&f_tik.vector, w->L_diag);
400 }
401 else if (w->L)
402 {
403 /* store Lx in Tikhonov portion of f~ */
404 gsl_blas_dgemv(CblasNoTrans, 1.0, w->L, x, 0.0, &f_tik.vector);
405 }
406 else
407 {
408 /* store \lambda x in Tikhonov portion of f~ */
409 gsl_vector_memcpy(&f_tik.vector, x);
410 gsl_vector_scale(&f_tik.vector, w->lambda);
411 }
412
413 return GSL_SUCCESS;
414 } /* fdfridge_f() */
415
416 static int
fdfridge_df(const gsl_vector * x,void * params,gsl_matrix * J)417 fdfridge_df(const gsl_vector * x, void * params, gsl_matrix * J)
418 {
419 int status;
420 gsl_multifit_fdfridge *w = (gsl_multifit_fdfridge *) params;
421 const size_t n = w->n;
422 const size_t p = w->p;
423 gsl_matrix_view J_user = gsl_matrix_submatrix(J, 0, 0, n, p);
424 gsl_matrix_view J_tik = gsl_matrix_submatrix(J, n, 0, p, p);
425 gsl_vector_view diag = gsl_matrix_diagonal(&J_tik.matrix);
426
427 /* compute Jacobian */
428 if (w->fdf->df)
429 status = gsl_multifit_eval_wdf(w->fdf, x, NULL, &J_user.matrix);
430 else
431 {
432 /* compute f(x) and then finite difference Jacobian */
433 status = gsl_multifit_eval_wf(w->fdf, x, NULL, w->f);
434 status += gsl_multifit_fdfsolver_dif_df(x, NULL, w->fdf, w->f,
435 &J_user.matrix);
436 }
437
438 if (status)
439 return status;
440
441 if (w->L_diag)
442 {
443 /* store diag(L_diag) in Tikhonov portion of J */
444 gsl_matrix_set_zero(&J_tik.matrix);
445 gsl_vector_memcpy(&diag.vector, w->L_diag);
446 }
447 else if (w->L)
448 {
449 /* store L in Tikhonov portion of J */
450 gsl_matrix_memcpy(&J_tik.matrix, w->L);
451 }
452 else
453 {
454 /* store \lambda I_p in Tikhonov portion of J */
455 gsl_matrix_set_zero(&J_tik.matrix);
456 gsl_vector_set_all(&diag.vector, w->lambda);
457 }
458
459 return GSL_SUCCESS;
460 } /* fdfridge_df() */
461