1 /* The MIT License
2
3 Copyright (c) 2013-2015 Genome Research Ltd.
4
5 Author: Petr Danecek <pd3@sanger.ac.uk>
6
7 Permission is hereby granted, free of charge, to any person obtaining a copy
8 of this software and associated documentation files (the "Software"), to deal
9 in the Software without restriction, including without limitation the rights
10 to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
11 copies of the Software, and to permit persons to whom the Software is
12 furnished to do so, subject to the following conditions:
13
14 The above copyright notice and this permission notice shall be included in
15 all copies or substantial portions of the Software.
16
17 THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
18 IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
19 FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
20 AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
21 LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
22 OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
23 THE SOFTWARE.
24
25 */
26
27 #include "peakfit.h"
28 #include <stdio.h>
29 #include <gsl/gsl_version.h>
30 #include <gsl/gsl_vector.h>
31 #include <gsl/gsl_multifit_nlin.h>
32 #include <htslib/hts.h>
33 #include <htslib/kstring.h>
34 #include <assert.h>
35 #include <math.h>
36
37 #define NPARAMS 5
38
39 // gauss params: sqrt(scale), center, sigma
40 typedef struct _peak_t
41 {
42 int fit_mask;
43 double params[NPARAMS], ori_params[NPARAMS]; // current and input parameters
44 struct { int scan; double min, max, best; } mc[NPARAMS]; // monte-carlo settings and best parameter
45 void (*calc_f) (int nvals, double *xvals, double *yvals, void *args);
46 void (*calc_df) (int nvals, double *xvals, double *yvals, double *dfvals, int idf, void *args);
47 void (*print_func) (struct _peak_t *pk, kstring_t *str);
48 void (*convert_get) (struct _peak_t *pk, double *params);
49 double (*convert_set) (struct _peak_t *pk, int iparam, double value);
50 }
51 peak_t;
52
53 struct _peakfit_t
54 {
55 int npeaks, mpeaks, nparams, mparams;
56 peak_t *peaks;
57 double *params;
58 int nvals, mvals;
59 double *xvals, *yvals, *vals;
60 kstring_t str;
61 int verbose, nmc_iter;
62 };
63
64
65 /*
66 Gaussian peak with the center bound in the interval <d,e>:
67 yi = scale^2 * exp(-(xi-z)^2/sigma^2)
68
69 dy/dscale = 2*scale * EXP
70 dy/dcenter = -scale^2 * sin(center) * (e-d) * (xi - z) * EXP / sigma^2
71 dy/dsigma = 2*scale^2 * (xi - z)^2 * EXP / sigma^3
72
73 where
74 z = 0.5*(cos(center)+1)*(e-d) + d
75 EXP = exp(-(xi-z)^2/sigma^2)
76 */
bounded_gaussian_calc_f(int nvals,double * xvals,double * yvals,void * args)77 void bounded_gaussian_calc_f(int nvals, double *xvals, double *yvals, void *args)
78 {
79 peak_t *pk = (peak_t*) args;
80
81 double scale2 = pk->params[0] * pk->params[0];
82 double center = pk->params[1];
83 double sigma = pk->params[2];
84 double d = pk->params[3];
85 double e = pk->params[4];
86 double z = 0.5*(cos(center)+1)*(e-d) + d;
87
88 int i;
89 for (i=0; i<nvals; i++)
90 {
91 double tmp = (xvals[i] - z)/sigma;
92 yvals[i] += scale2 * exp(-tmp*tmp);
93 }
94 }
bounded_gaussian_calc_df(int nvals,double * xvals,double * yvals,double * dfvals,int idf,void * args)95 void bounded_gaussian_calc_df(int nvals, double *xvals, double *yvals, double *dfvals, int idf, void *args)
96 {
97 peak_t *pk = (peak_t*) args;
98
99 double scale = pk->params[0];
100 double center = pk->params[1];
101 double sigma = pk->params[2];
102 double d = pk->params[3];
103 double e = pk->params[4];
104 double z = 0.5*(cos(center)+1)*(e-d) + d;
105
106 int i;
107 for (i=0; i<nvals; i++)
108 {
109 double EXP = exp(-(xvals[i]-z)*(xvals[i]-z)/sigma/sigma);
110 double zi = xvals[i] - z;
111 if ( idf==0 ) // dscale
112 dfvals[i] += 2*scale*EXP;
113 else if ( idf==1 ) // dcenter
114 dfvals[i] -= scale*scale*sin(center)*(e-d)*zi*EXP/sigma/sigma;
115 else if ( idf==2 ) // dsigma
116 dfvals[i] += 2*scale*scale*zi*zi*EXP/sigma/sigma/sigma;
117 }
118 }
bounded_gaussian_sprint_func(peak_t * pk,kstring_t * str)119 void bounded_gaussian_sprint_func(peak_t *pk, kstring_t *str)
120 {
121 double center = pk->params[1];
122 double d = pk->params[3];
123 double e = pk->params[4];
124 double z = 0.5*(cos(center)+1)*(e-d) + d;
125 ksprintf(str,"%f**2 * exp(-(x-%f)**2/%f**2)",fabs(pk->params[0]),z,fabs(pk->params[2]));
126 }
bounded_gaussian_convert_set(peak_t * pk,int iparam,double value)127 double bounded_gaussian_convert_set(peak_t *pk, int iparam, double value)
128 {
129 if ( iparam!=1 ) return value;
130 double d = pk->ori_params[3];
131 double e = pk->ori_params[4];
132 if ( value<d ) value = d;
133 else if ( value>e ) value = e;
134 return acos(2*(value-d)/(e-d) - 1);
135 }
bounded_gaussian_convert_get(peak_t * pk,double * params)136 void bounded_gaussian_convert_get(peak_t *pk, double *params)
137 {
138 params[0] = fabs(pk->params[0]);
139 params[2] = fabs(pk->params[2]);
140 double center = pk->params[1];
141 double d = pk->params[3];
142 double e = pk->params[4];
143 params[1] = 0.5*(cos(center)+1)*(e-d) + d;
144 }
145
peakfit_add_bounded_gaussian(peakfit_t * pkf,double a,double b,double c,double d,double e,int fit_mask)146 void peakfit_add_bounded_gaussian(peakfit_t *pkf, double a, double b, double c, double d, double e, int fit_mask)
147 {
148 pkf->npeaks++;
149 hts_expand0(peak_t,pkf->npeaks,pkf->mpeaks,pkf->peaks);
150
151 int i, nfit = 0;
152 for (i=0; i<NPARAMS; i++)
153 if ( fit_mask & (1<<i) ) nfit++;
154
155 assert( d<e );
156
157 pkf->nparams += nfit;
158 hts_expand0(double,pkf->nparams,pkf->mparams,pkf->params);
159
160 peak_t *pk = &pkf->peaks[pkf->npeaks-1];
161 memset(pk, 0, sizeof(peak_t));
162
163 pk->calc_f = bounded_gaussian_calc_f;
164 pk->calc_df = bounded_gaussian_calc_df;
165 pk->print_func = bounded_gaussian_sprint_func;
166 pk->convert_set = bounded_gaussian_convert_set;
167 pk->convert_get = bounded_gaussian_convert_get;
168 pk->fit_mask = fit_mask;
169 pk->ori_params[0] = a;
170 pk->ori_params[2] = c;
171 pk->ori_params[3] = d;
172 pk->ori_params[4] = e;
173 pk->ori_params[1] = pk->convert_set(pk, 1, b);
174 }
175
176
177 /*
178 Gaussian peak:
179 yi = scale^2 * exp(-(x-center)^2/sigma^2)
180
181 dy/dscale = 2 * scale * EXP
182 dy/dcenter = 2 * scale^2 * (x-center) * EXP / sigma^2
183 dy/dsigma = 2 * scale^2 * (x-center)^2 * EXP / sigma^3
184
185 where
186 EXP = exp(-(x-center)^2/sigma^2)
187 */
gaussian_calc_f(int nvals,double * xvals,double * yvals,void * args)188 void gaussian_calc_f(int nvals, double *xvals, double *yvals, void *args)
189 {
190 peak_t *pk = (peak_t*) args;
191
192 double scale2 = pk->params[0] * pk->params[0];
193 double center = pk->params[1];
194 double sigma = pk->params[2];
195
196 int i;
197 for (i=0; i<nvals; i++)
198 {
199 double tmp = (xvals[i] - center)/sigma;
200 yvals[i] += scale2 * exp(-tmp*tmp);
201 }
202 }
gaussian_calc_df(int nvals,double * xvals,double * yvals,double * dfvals,int idf,void * args)203 void gaussian_calc_df(int nvals, double *xvals, double *yvals, double *dfvals, int idf, void *args)
204 {
205 peak_t *pk = (peak_t*) args;
206
207 double scale = pk->params[0];
208 double center = pk->params[1];
209 double sigma = pk->params[2];
210
211 int i;
212 for (i=0; i<nvals; i++)
213 {
214 double zi = xvals[i] - center;
215 double EXP = exp(-zi*zi/(sigma*sigma));
216 if ( idf==0 ) // dscale
217 dfvals[i] += 2*scale*EXP;
218 else if ( idf==1 ) // dcenter
219 dfvals[i] += 2*scale*scale*zi*EXP/(sigma*sigma);
220 else if ( idf==2 ) // dsigma
221 dfvals[i] += 2*scale*scale*zi*zi*EXP/(sigma*sigma*sigma);
222 }
223 }
gaussian_sprint_func(struct _peak_t * pk,kstring_t * str)224 void gaussian_sprint_func(struct _peak_t *pk, kstring_t *str)
225 {
226 ksprintf(str,"%f**2 * exp(-(x-%f)**2/%f**2)",fabs(pk->params[0]),pk->params[1],fabs(pk->params[2]));
227 }
gaussian_convert_get(peak_t * pk,double * params)228 void gaussian_convert_get(peak_t *pk, double *params)
229 {
230 params[0] = fabs(pk->params[0]);
231 params[1] = fabs(pk->params[1]);
232 params[2] = fabs(pk->params[2]);
233 }
234
235
peakfit_add_gaussian(peakfit_t * pkf,double a,double b,double c,int fit_mask)236 void peakfit_add_gaussian(peakfit_t *pkf, double a, double b, double c, int fit_mask)
237 {
238 pkf->npeaks++;
239 hts_expand0(peak_t,pkf->npeaks,pkf->mpeaks,pkf->peaks);
240
241 int i, nfit = 0;
242 for (i=0; i<NPARAMS; i++)
243 if ( fit_mask & (1<<i) ) nfit++;
244
245 pkf->nparams += nfit;
246 hts_expand0(double,pkf->nparams,pkf->mparams,pkf->params);
247
248 peak_t *pk = &pkf->peaks[pkf->npeaks-1];
249 memset(pk, 0, sizeof(peak_t));
250
251 pk->calc_f = gaussian_calc_f;
252 pk->calc_df = gaussian_calc_df;
253 pk->print_func = gaussian_sprint_func;
254 pk->convert_get = gaussian_convert_get;
255 pk->fit_mask = fit_mask;
256 pk->ori_params[0] = a;
257 pk->ori_params[1] = b;
258 pk->ori_params[2] = c;
259 }
260
261
262 /*
263 exp peak:
264 yi = scale^2 * exp((x-center)/sigma^2)
265
266 dy/dscale = 2 * scale * EXP
267 dy/dcenter = -scale^2 * EXP / sigma^2
268 dy/dsigma = -2 * scale^2 * (x-center) * EXP / sigma^3
269
270 where
271 EXP = exp((x-center)/sigma^2)
272 */
exp_calc_f(int nvals,double * xvals,double * yvals,void * args)273 void exp_calc_f(int nvals, double *xvals, double *yvals, void *args)
274 {
275 peak_t *pk = (peak_t*) args;
276
277 double scale2 = pk->params[0] * pk->params[0];
278 double center = pk->params[1];
279 double sigma = pk->params[2];
280
281 int i;
282 for (i=0; i<nvals; i++)
283 {
284 yvals[i] += scale2 * exp((xvals[i]-center)/sigma/sigma);
285 }
286 }
exp_calc_df(int nvals,double * xvals,double * yvals,double * dfvals,int idf,void * args)287 void exp_calc_df(int nvals, double *xvals, double *yvals, double *dfvals, int idf, void *args)
288 {
289 peak_t *pk = (peak_t*) args;
290
291 double scale = pk->params[0];
292 double center = pk->params[1];
293 double sigma = pk->params[2];
294
295 int i;
296 for (i=0; i<nvals; i++)
297 {
298 double EXP = exp((xvals[i]-center)/sigma/sigma);
299 if ( idf==0 ) // dscale
300 dfvals[i] += 2*scale*EXP;
301 else if ( idf==2 ) // dsigma
302 dfvals[i] -= 2*scale*scale*(xvals[i]-center)*EXP/sigma/sigma/sigma;
303 }
304 }
exp_sprint_func(struct _peak_t * pk,kstring_t * str)305 void exp_sprint_func(struct _peak_t *pk, kstring_t *str)
306 {
307 ksprintf(str,"%f**2 * exp((x-%f)/%f**2)",fabs(pk->params[0]),pk->params[1],fabs(pk->params[2]));
308 }
exp_convert_get(peak_t * pk,double * params)309 void exp_convert_get(peak_t *pk, double *params)
310 {
311 params[0] = fabs(pk->params[0]);
312 params[1] = fabs(pk->params[1]);
313 params[2] = fabs(pk->params[2]);
314 }
315
peakfit_add_exp(peakfit_t * pkf,double a,double b,double c,int fit_mask)316 void peakfit_add_exp(peakfit_t *pkf, double a, double b, double c, int fit_mask)
317 {
318 pkf->npeaks++;
319 hts_expand0(peak_t,pkf->npeaks,pkf->mpeaks,pkf->peaks);
320
321 int i, nfit = 0;
322 for (i=0; i<NPARAMS; i++)
323 if ( fit_mask & (1<<i) ) nfit++;
324
325 assert( !(fit_mask&2) );
326
327 pkf->nparams += nfit;
328 hts_expand0(double,pkf->nparams,pkf->mparams,pkf->params);
329
330 peak_t *pk = &pkf->peaks[pkf->npeaks-1];
331 memset(pk, 0, sizeof(peak_t));
332
333 pk->calc_f = exp_calc_f;
334 pk->calc_df = exp_calc_df;
335 pk->print_func = exp_sprint_func;
336 pk->convert_get = exp_convert_get;
337 pk->fit_mask = fit_mask;
338 pk->ori_params[0] = a;
339 pk->ori_params[1] = b;
340 pk->ori_params[2] = c;
341 }
342
343
peakfit_set_params(peakfit_t * pkf,int ipk,double * params,int nparams)344 void peakfit_set_params(peakfit_t *pkf, int ipk, double *params, int nparams)
345 {
346 peak_t *pk = &pkf->peaks[ipk];
347 int i;
348 if ( pk->convert_set )
349 for (i=0; i<nparams; i++) pk->params[i] = pk->convert_set(pk, i, params[i]);
350 else
351 for (i=0; i<nparams; i++) pk->params[i] = params[i];
352 }
353
peakfit_get_params(peakfit_t * pkf,int ipk,double * params,int nparams)354 void peakfit_get_params(peakfit_t *pkf, int ipk, double *params, int nparams)
355 {
356 peak_t *pk = &pkf->peaks[ipk];
357 if ( pk->convert_get ) pk->convert_get(pk, params);
358 else
359 {
360 int i;
361 for (i=0; i<nparams; i++) params[i] = pk->params[i];
362 }
363 }
364
peakfit_init(void)365 peakfit_t *peakfit_init(void)
366 {
367 return (peakfit_t*)calloc(1,sizeof(peakfit_t));
368 }
369
peakfit_reset(peakfit_t * pkf)370 void peakfit_reset(peakfit_t *pkf)
371 {
372 pkf->npeaks = pkf->nparams = 0;
373 memset(pkf->peaks,0,sizeof(peak_t)*pkf->mpeaks);
374 }
375
peakfit_destroy(peakfit_t * pkf)376 void peakfit_destroy(peakfit_t *pkf)
377 {
378 free(pkf->str.s);
379 free(pkf->vals);
380 free(pkf->params);
381 free(pkf->peaks);
382 free(pkf);
383 }
384
peakfit_calc_f(const gsl_vector * params,void * data,gsl_vector * yvals)385 int peakfit_calc_f(const gsl_vector *params, void *data, gsl_vector *yvals)
386 {
387 peakfit_t *pkf = (peakfit_t *) data;
388
389 int i,j;
390 for (i=0; i<pkf->nvals; i++)
391 pkf->vals[i] = 0;
392
393 int iparam = 0;
394 for (i=0; i<pkf->npeaks; i++)
395 {
396 peak_t *pk = &pkf->peaks[i];
397 for (j=0; j<NPARAMS; j++)
398 {
399 if ( !(pk->fit_mask & (1<<j)) ) continue;
400 pk->params[j] = gsl_vector_get(params,iparam);
401 iparam++;
402 }
403 pk->calc_f(pkf->nvals, pkf->xvals, pkf->vals, pk);
404 }
405
406 for (i=0; i<pkf->nvals; i++)
407 gsl_vector_set(yvals, i, (pkf->vals[i] - pkf->yvals[i])/0.01);
408
409 return GSL_SUCCESS;
410 }
peakfit_calc_df(const gsl_vector * params,void * data,gsl_matrix * jacobian)411 int peakfit_calc_df(const gsl_vector *params, void *data, gsl_matrix *jacobian)
412 {
413 peakfit_t *pkf = (peakfit_t *) data;
414
415 int i,j,k,iparam = 0;
416 for (i=0; i<pkf->npeaks; i++)
417 {
418 peak_t *pk = &pkf->peaks[i];
419 int iparam_prev = iparam;
420 for (j=0; j<NPARAMS; j++)
421 {
422 if ( !(pk->fit_mask & (1<<j)) ) continue;
423 pk->params[j] = gsl_vector_get(params,iparam);
424 iparam++;
425 }
426 iparam = iparam_prev;
427 for (j=0; j<NPARAMS; j++)
428 {
429 if ( !(pk->fit_mask & (1<<j)) ) continue;
430 for (k=0; k<pkf->nvals; k++) pkf->vals[k] = 0;
431 pk->calc_df(pkf->nvals, pkf->xvals, pkf->yvals, pkf->vals, j, pk);
432 for (k=0; k<pkf->nvals; k++) gsl_matrix_set(jacobian, k, iparam, pkf->vals[k]);
433 iparam++;
434 }
435 }
436 return GSL_SUCCESS;
437 }
peakfit_calc_fdf(const gsl_vector * params,void * data,gsl_vector * yvals,gsl_matrix * jacobian)438 int peakfit_calc_fdf(const gsl_vector *params, void *data, gsl_vector *yvals, gsl_matrix *jacobian)
439 {
440 peakfit_calc_f(params, data, yvals);
441 peakfit_calc_df(params, data, jacobian);
442 return GSL_SUCCESS;
443 }
444
peakfit_evaluate(peakfit_t * pkf)445 double peakfit_evaluate(peakfit_t *pkf)
446 {
447 int i;
448 for (i=0; i<pkf->nvals; i++)
449 pkf->vals[i] = 0;
450
451 for (i=0; i<pkf->npeaks; i++)
452 pkf->peaks[i].calc_f(pkf->nvals, pkf->xvals, pkf->vals, &pkf->peaks[i]);
453
454 double sum = 0;
455 for (i=0; i<pkf->nvals; i++)
456 sum += fabs(pkf->vals[i] - pkf->yvals[i]);
457
458 return sum;
459 }
460
peakfit_sprint_func(peakfit_t * pkf)461 const char *peakfit_sprint_func(peakfit_t *pkf)
462 {
463 pkf->str.l = 0;
464 int i;
465 for (i=0; i<pkf->npeaks; i++)
466 {
467 if ( i>0 ) kputs(" + ", &pkf->str);
468 pkf->peaks[i].print_func(&pkf->peaks[i], &pkf->str);
469 }
470 return (const char*)pkf->str.s;
471 }
472
peakfit_verbose(peakfit_t * pkf,int level)473 void peakfit_verbose(peakfit_t *pkf, int level)
474 {
475 pkf->verbose = level;
476 }
477
peakfit_set_mc(peakfit_t * pkf,double xmin,double xmax,int iparam,int niter)478 void peakfit_set_mc(peakfit_t *pkf, double xmin, double xmax, int iparam, int niter)
479 {
480 peak_t *pk = &pkf->peaks[ pkf->npeaks-1 ];
481 pk->mc[iparam].scan = 1;
482 pk->mc[iparam].min = xmin;
483 pk->mc[iparam].max = xmax;
484 pkf->nmc_iter = niter;
485 }
486
peakfit_run(peakfit_t * pkf,int nvals,double * xvals,double * yvals)487 double peakfit_run(peakfit_t *pkf, int nvals, double *xvals, double *yvals)
488 {
489 srand(0); // for reproducibility
490
491 pkf->nvals = nvals;
492 pkf->xvals = xvals;
493 pkf->yvals = yvals;
494 hts_expand0(double,pkf->nvals,pkf->mvals,pkf->vals);
495 if ( !pkf->nparams ) return peakfit_evaluate(pkf);
496
497 gsl_multifit_function_fdf mfunc;
498 mfunc.f = &peakfit_calc_f;
499 mfunc.df = &peakfit_calc_df;
500 mfunc.fdf = &peakfit_calc_fdf;
501 mfunc.n = nvals;
502 mfunc.p = pkf->nparams;
503 mfunc.params = pkf;
504
505 const gsl_multifit_fdfsolver_type *solver_type;
506 gsl_multifit_fdfsolver *solver;
507 solver_type = gsl_multifit_fdfsolver_lmsder;
508 solver = gsl_multifit_fdfsolver_alloc(solver_type, nvals, mfunc.p);
509 gsl_vector *grad = gsl_vector_alloc(pkf->nparams);
510
511 int imc_iter, i,j, iparam;
512 double best_fit = HUGE_VAL;
513 for (imc_iter=0; imc_iter<=pkf->nmc_iter; imc_iter++) // possibly multiple monte-carlo iterations
514 {
515 // set GSL parameters
516 iparam = 0;
517 for (i=0; i<pkf->npeaks; i++)
518 {
519 peak_t *pk = &pkf->peaks[i];
520 for (j=0; j<NPARAMS; j++)
521 {
522 pk->params[j] = pk->ori_params[j];
523 if ( pk->mc[j].scan )
524 {
525 pk->params[j] = rand()*(pk->mc[j].max - pk->mc[j].min)/RAND_MAX + pk->mc[j].min;
526 if ( pk->convert_set ) pk->params[j] = pk->convert_set(pk, j, pk->params[j]);
527 }
528 if ( !(pk->fit_mask & (1<<j)) ) continue;
529 pkf->params[iparam] = pk->params[j];
530 iparam++;
531 }
532 }
533
534 gsl_vector_view vview = gsl_vector_view_array(pkf->params, mfunc.p);
535 gsl_multifit_fdfsolver_set(solver, &mfunc, &vview.vector);
536
537 // iterate until convergence (or lack of it)
538 int ret, test1 = 0, test2 = 0, niter = 0, niter_max = 500;
539 do
540 {
541 ret = gsl_multifit_fdfsolver_iterate(solver);
542 if ( pkf->verbose >1 )
543 {
544 fprintf(stderr, "%d: ", niter);
545 for (i=0; i<pkf->npeaks; i++)
546 {
547 peak_t *pk = &pkf->peaks[i];
548 fprintf(stderr,"\t%f %f %f", pk->params[0],pk->params[1],pk->params[2]);
549 }
550 fprintf(stderr, "\t.. %s\n", gsl_strerror(ret));
551 }
552 if ( ret ) break;
553
554 #if GSL_MAJOR_VERSION >= 2
555 int info;
556 test1 = gsl_multifit_fdfsolver_test(solver, 1e-8,1e-8, 0.0, &info);
557 #else
558 gsl_multifit_gradient(solver->J, solver->f, grad);
559 test1 = gsl_multifit_test_gradient(grad, 1e-8);
560 test2 = gsl_multifit_test_delta(solver->dx, solver->x, 1e-8, 1e-8);
561 #endif
562 }
563 while ((test1==GSL_CONTINUE || test2==GSL_CONTINUE) && ++niter<niter_max);
564 if ( pkf->verbose >1 )
565 {
566 fprintf(stderr,"test1=%s\n", gsl_strerror(test1));
567 fprintf(stderr,"test2=%s\n", gsl_strerror(test2));
568 }
569
570 // recover parameters
571 iparam = 0;
572 for (i=0; i<pkf->npeaks; i++)
573 {
574 peak_t *pk = &pkf->peaks[i];
575 for (j=0; j<NPARAMS; j++)
576 {
577 if ( !(pk->fit_mask & (1<<j)) ) continue;
578 pk->params[j] = gsl_vector_get(solver->x, iparam++);
579 }
580 }
581
582 // evaluate fit, update best parameters
583 double fit = peakfit_evaluate(pkf);
584 if ( fit<best_fit )
585 {
586 for (i=0; i<pkf->npeaks; i++)
587 {
588 peak_t *pk = &pkf->peaks[i];
589 for (j=0; j<NPARAMS; j++) pk->mc[j].best = pk->params[j];
590 }
591 }
592 if ( fit<best_fit ) best_fit = fit;
593 }
594 gsl_multifit_fdfsolver_free(solver);
595 gsl_vector_free(grad);
596
597 for (i=0; i<pkf->npeaks; i++)
598 {
599 peak_t *pk = &pkf->peaks[i];
600 for (j=0; j<NPARAMS; j++) pk->params[j] = pk->mc[j].best;
601 }
602 return best_fit;
603 }
604
605
606