#
source:
trunk/projects/slim-curve/src/main/c/EcfGlobal.c
@
7611

Revision 7611, 99.7 KB checked in by aivar, 9 years ago (diff) |
---|

Rev | Line | |
---|---|---|

[7606] | 1 | /* The 2003 version of the ECF library. This has basically been |

2 | completely rewritten using the Numerical Recipes code (modified as | |

3 | appropriate). Also, this takes account of the fact that we may be | |

4 | handling Poisson noise. | |

5 | ||

6 | This file contains functions for global analysis, referring to a | |

7 | little code from the single transient analysis functions. Utility | |

8 | code is found in EcfUtil.c and single transient analysis code in | |

9 | EcfSingle.c. | |

10 | */ | |

11 | ||

12 | #include <stdio.h> | |

13 | #include <stdlib.h> | |

14 | #include <math.h> | |

15 | #include "EcfInternal.h" | |

16 | ||

17 | typedef struct { | |

18 | float **P, **Q, ***S; | |

19 | } global_matrix; | |

20 | ||

21 | typedef struct { | |

22 | float *global, *local; | |

23 | } global_vector; | |

24 | ||

25 | /* Predeclarations */ | |

26 | ||

27 | int GCI_alloc_global_matrix(global_matrix *m, | |

28 | int global, int local, int ntrans); | |

29 | void GCI_free_global_matrix(global_matrix *m); | |

30 | void GCI_copy_global_matrix(global_matrix dest, global_matrix src, | |

31 | int global, int local, int ntrans); | |

32 | int GCI_alloc_global_vector(global_vector *v, | |

33 | int global, int local, int ntrans); | |

34 | void GCI_free_global_vector(global_vector *v); | |

35 | void GCI_copy_global_vector(global_vector dest, global_vector src, | |

36 | int global, int local, int ntrans); | |

37 | int GCI_marquardt_global_exps_est_globals_instr( | |

38 | float xincr, float **trans, int ndata, int ntrans, | |

39 | int fit_start, int fit_end, float instr[], int ninstr, | |

40 | noise_type noise, float sig[], int ftype, | |

41 | float **param, int paramfree[], int nparam, float gparam[], | |

42 | restrain_type restrain, float chisq_delta, | |

43 | float fitted[], float residuals[], | |

44 | float **covar, float **alpha, float *chisq_global); | |

45 | int GCI_marquardt_global_exps_est_params_instr( | |

46 | float xincr, float **trans, int ndata, int ntrans, | |

47 | int fit_start, int fit_end, float instr[], int ninstr, | |

48 | noise_type noise, float sig[], int ftype, | |

49 | float **param, int paramfree[], int nparam, restrain_type restrain, float chisq_delta, | |

50 | float exp_pure[], float *exp_conv[], | |

51 | float **fitted, float **residuals, | |

52 | float **covar, float **alpha, float chisq_trans[], | |

53 | int drop_bad_transients); | |

54 | int GCI_marquardt_global_exps_calculate_exps_instr( | |

55 | float xincr, int ndata, float instr[], int ninstr, | |

56 | int ftype, float param[], int nparam, | |

57 | float exp_pure[], float *exp_conv[]); | |

58 | int GCI_marquardt_global_exps_do_fit_single( | |

59 | float xincr, float y[], int ndata, int fit_start, int fit_end, | |

60 | noise_type noise, float sig[], int ftype, | |

61 | float param[], int paramfree[], int nparam, restrain_type restrain, float chisq_delta, | |

62 | float *exp_conv[], float *fitted, float *residuals, | |

63 | float **covar, float **alpha, float *chisq_trans); | |

64 | int GCI_marquardt_global_exps_single_step( | |

65 | float xincr, float y[], | |

66 | int ndata, int fit_start, int fit_end, | |

67 | noise_type noise, float sig[], int ftype, | |

68 | float param[], int paramfree[], int nparam, | |

69 | restrain_type restrain, | |

70 | float *exp_conv[], | |

71 | float yfit[], float dy[], | |

72 | float **covar, float **alpha, float *chisq, | |

73 | float *alambda); | |

74 | int GCI_marquardt_global_compute_exps_fn( | |

75 | float xincr, float y[], | |

76 | int ndata, int fit_start, int fit_end, | |

77 | noise_type noise, float sig[], int ftype, | |

78 | float param[], int paramfree[], int nparam, | |

79 | float *exp_conv[], | |

80 | float yfit[], float dy[], | |

81 | float **alpha, float *beta, float *chisq); | |

82 | int GCI_marquardt_global_compute_exps_fn_final( | |

83 | float xincr, float y[], | |

84 | int ndata, int fit_start, int fit_end, | |

85 | noise_type noise, float sig[], int ftype, | |

86 | float param[], int paramfree[], int nparam, | |

87 | float *exp_conv[], | |

88 | float yfit[], float dy[], float *chisq); | |

89 | int GCI_marquardt_global_exps_do_fit_instr( | |

90 | float xincr, float **trans, int ndata, int ntrans, | |

91 | int fit_start, int fit_end, float instr[], int ninstr, | |

92 | noise_type noise, float sig[], int ftype, | |

93 | float **param, int paramfree[], int nparam, restrain_type restrain, float chisq_delta, | |

94 | float exp_pure[], float *exp_conv[], | |

95 | float **fitted, float **residuals, | |

96 | float **covar_scratch, float **alpha_scratch, | |

97 | float *chisq_trans, float *chisq_global, | |

98 | int drop_bad_transients); | |

99 | int GCI_marquardt_global_exps_global_step( | |

100 | float xincr, float **trans, | |

101 | int ndata, int ntrans, int fit_start, int fit_end, | |

102 | float instr[], int ninstr, | |

103 | noise_type noise, float sig[], int ftype, | |

104 | float **param, int paramfree[], int nparam, | |

105 | restrain_type restrain, | |

106 | float exp_pure[], float *exp_conv[], | |

107 | float **yfit, float **dy, | |

108 | float *chisq_trans, float *chisq_global, | |

109 | float **alpha_scratch, float *alambda, | |

110 | int drop_bad_transients); | |

111 | int GCI_marquardt_global_compute_global_exps_fn( | |

112 | float xincr, float **trans, int ndata, int ntrans, | |

113 | int fit_start, int fit_end, float instr[], int ninstr, | |

114 | noise_type noise, float sig[], int ftype, | |

115 | float **param, int paramfree[], int nparam, | |

116 | int mfit_global, int mfit_local, int gindex[], int lindex[], | |

117 | float exp_pure[], float *exp_conv[], | |

118 | float **yfit, float **dy, global_matrix alpha, global_vector beta, | |

119 | float **alpha_scratch, float *chisq_trans, float *chisq_global, | |

120 | int drop_bad_transients); | |

121 | int GCI_marquardt_global_compute_global_exps_fn_final( | |

122 | float xincr, float **trans, int ndata, int ntrans, | |

123 | int fit_start, int fit_end, float instr[], int ninstr, | |

124 | noise_type noise, float sig[], int ftype, | |

125 | float **param, int paramfree[], int nparam, | |

126 | int mfit_global, int mfit_local, int gindex[], int lindex[], | |

127 | float exp_pure[], float *exp_conv[], | |

128 | float **yfit, float **dy, | |

129 | float *chisq_trans, float *chisq_global, int drop_bad_transients); | |

130 | ||

131 | int GCI_marquardt_global_solve_eqn(global_matrix A, global_vector b, | |

132 | int mfit_global, int mfit_local, int ntrans); | |

133 | ||

134 | int GCI_marquardt_global_generic_do_fit_instr( | |

135 | float xincr, float **trans, int ndata, int ntrans, | |

136 | int fit_start, int fit_end, float instr[], int ninstr, | |

137 | noise_type noise, float sig[], | |

138 | float **param, int paramfree[], int nparam, int gparam[], | |

139 | restrain_type restrain, float chisq_delta, | |

140 | void (*fitfunc)(float, float [], float *, float [], int), | |

141 | float **fitted, float **residuals, | |

142 | float **covar_scratch, float **alpha_scratch, | |

143 | float *chisq_trans, float *chisq_global); | |

144 | int GCI_marquardt_global_generic_global_step( | |

145 | float xincr, float **trans, | |

146 | int ndata, int ntrans, int fit_start, int fit_end, | |

147 | float instr[], int ninstr, | |

148 | noise_type noise, float sig[], | |

149 | float **param, int paramfree[], int nparam, int gparam[], | |

150 | restrain_type restrain, float chisq_delta, | |

151 | void (*fitfunc)(float, float [], float *, float [], int), | |

152 | float **yfit, float **dy, | |

153 | float *chisq_trans, float *chisq_global, | |

154 | float **alpha_scratch, float *alambda); | |

155 | int GCI_marquardt_global_compute_global_generic_fn( | |

156 | float xincr, float **trans, int ndata, int ntrans, | |

157 | int fit_start, int fit_end, float instr[], int ninstr, | |

158 | noise_type noise, float sig[], | |

159 | float **param, int paramfree[], int nparam, int gparam[], | |

160 | int mfit_global, int mfit_local, int gindex[], int lindex[], | |

161 | void (*fitfunc)(float, float [], float *, float [], int), | |

162 | float **yfit, float **dy, global_matrix alpha, global_vector beta, | |

163 | float **alpha_scratch, float *chisq_trans, float *chisq_global, | |

164 | float alambda, | |

165 | float **pfnvals, float ***pdy_dparam_pure, float ***pdy_dparam_conv, | |

166 | int *pfnvals_len, int *pdy_dparam_nparam_size); | |

167 | int GCI_marquardt_global_compute_global_generic_fn_final( | |

168 | float xincr, float **trans, int ndata, int ntrans, | |

169 | int fit_start, int fit_end, float instr[], int ninstr, | |

170 | noise_type noise, float sig[], | |

171 | float **param, int paramfree[], int nparam, int gparam[], | |

172 | int mfit_global, int mfit_local, int gindex[], int lindex[], | |

173 | void (*fitfunc)(float, float [], float *, float [], int), | |

174 | float **yfit, float **dy, | |

175 | float *chisq_trans, float *chisq_global, | |

176 | float **pfnvals, float ***pdy_dparam_pure, float ***pdy_dparam_conv, | |

177 | int *pfnvals_len, int *pdy_dparam_nparam_size); | |

178 | ||

179 | ||

180 | /******************************************************************** | |

181 | ||

182 | GLOBAL ANALYSIS | |

183 | ||

184 | ********************************************************************/ | |

185 | ||

186 | /* We only work with the case of multiple transients, all with the | |

187 | same instrument/prompt response, the same xincr, the same number of | |

188 | points, and so on. The recommended fitting algorithm is a | |

189 | three-step process: | |

190 | ||

191 | (1) Sum the transients and use this to get initial estimates for | |

192 | the global parameters we are estimating. | |

193 | ||

194 | (2) Fixing these parameters, perform a Marquardt fit on each of the | |

195 | transients. In our cases, this will be fairly efficient, as we | |

196 | will not need to repeatedly calculate the exponential decay. | |

197 | ||

198 | (3) Finally, perform the global fit. There's lots of interesting | |

199 | maths here which we'll discuss when we get to it. Note that | |

200 | again we will not need to repeatedly calculate the exponential | |

201 | decays for each transient, which will hopefully make the | |

202 | process significantly faster. | |

203 | ||

204 | We provide special code to perform these steps in the case of | |

205 | multiexponential and stretched exponential fits where we are aiming | |

206 | to globally fit all of the taus and the h parameter (in the | |

207 | stretched exponential case); these can be performed far more | |

208 | efficiently than the general case. We also provide code to perform | |

209 | step (3) in the general case; steps (1) and (2) will have to be | |

210 | handled on a case-by-case basis by the calling code. We also | |

211 | provide a version of step (3) to handle the case of arbitrary | |

212 | x data with no instrument response. | |

213 | */ | |

214 | ||

215 | /* ***** UTILITY CODE ***** */ | |

216 | ||

217 | /* First, we define functions to handle the data structures we will be | |

218 | using later on to store the alpha and covar matrices and the beta | |

219 | etc. vectors for global analysis */ | |

220 | ||

221 | int GCI_alloc_global_matrix(global_matrix *m, | |

222 | int global, int local, int ntrans) | |

223 | { | |

224 | if (global <= 0 || local < 0 || ntrans <= 0) | |

225 | return -2; | |

226 | ||

227 | if ((m->P = GCI_ecf_matrix(global, global)) == NULL) | |

228 | return -1; | |

229 | if (local > 0) { | |

230 | if ((m->Q = GCI_ecf_matrix(global, ntrans*local)) == NULL) { | |

231 | GCI_ecf_free_matrix(m->P); | |

232 | return -1; | |

233 | } | |

234 | if ((m->S = GCI_ecf_matrix_array(ntrans, local, local)) == NULL) { | |

235 | GCI_ecf_free_matrix(m->P); | |

236 | GCI_ecf_free_matrix(m->Q); | |

237 | return -1; | |

238 | } | |

239 | } else { | |

240 | m->Q = NULL; | |

241 | m->S = NULL; | |

242 | } | |

243 | ||

244 | return 0; | |

245 | } | |

246 | ||

247 | void GCI_free_global_matrix(global_matrix *m) | |

248 | { | |

249 | GCI_ecf_free_matrix(m->P); | |

250 | if (m->Q != NULL) { | |

251 | GCI_ecf_free_matrix(m->Q); | |

252 | GCI_ecf_free_matrix_array(m->S); | |

253 | } | |

254 | } | |

255 | ||

256 | void GCI_copy_global_matrix(global_matrix dest, global_matrix src, | |

257 | int global, int local, int ntrans) | |

258 | { | |

259 | int i, j, k; | |

260 | ||

261 | for (i=0; i<global; i++) | |

262 | for (j=0; j<global; j++) | |

263 | dest.P[i][j] = src.P[i][j]; | |

264 | ||

265 | if (local > 0) { | |

266 | for (i=0; i<global; i++) | |

267 | for (j=0; j<ntrans*local; j++) | |

268 | dest.Q[i][j] = src.Q[i][j]; | |

269 | ||

270 | for (i=0; i<ntrans; i++) | |

271 | for (j=0; j<local; j++) | |

272 | for (k=0; k<local; k++) | |

273 | dest.S[i][j][k] = src.S[i][j][k]; | |

274 | } | |

275 | } | |

276 | ||

277 | ||

278 | int GCI_alloc_global_vector(global_vector *v, | |

279 | int global, int local, int ntrans) | |

280 | { | |

281 | if (global <= 0 || local < 0 || ntrans <= 0) | |

282 | return -2; | |

283 | ||

284 | if ((v->global = (float *) malloc(global * sizeof(float))) == NULL) | |

285 | return -1; | |

286 | if (local > 0) { | |

287 | if ((v->local = | |

288 | (float *) malloc(ntrans * local * sizeof(float))) == NULL) { | |

289 | free(v->global); | |

290 | return -1; | |

291 | } | |

292 | } else { | |

293 | v->local = NULL; | |

294 | } | |

295 | ||

296 | return 0; | |

297 | } | |

298 | ||

299 | void GCI_free_global_vector(global_vector *v) | |

300 | { | |

301 | free(v->global); | |

302 | if (v->local != NULL) | |

303 | free(v->local); | |

304 | } | |

305 | ||

306 | void GCI_copy_global_vector(global_vector dest, global_vector src, | |

307 | int global, int local, int ntrans) | |

308 | { | |

309 | int i; | |

310 | ||

311 | for (i=0; i<global; i++) | |

312 | dest.global[i] = src.global[i]; | |

313 | ||

314 | if (local > 0) { | |

315 | for (i=0; i<ntrans*local; i++) | |

316 | dest.local[i] = src.local[i]; | |

317 | } | |

318 | } | |

319 | ||

320 | ||

321 | /* ***** EXPONENTIALS CODE ***** */ | |

322 | ||

323 | /* Now the code for performing a global fit for multiexponential taus | |

324 | and stretched exponentials. This is the function which is called | |

325 | from external programs. */ | |

326 | ||

327 | /* I don't even want to _contemplate_ error axes for this! It would | |

328 | be computationally very messy, as there would be very large numbers | |

329 | of large vectors involved. */ | |

330 | ||

331 | int GCI_marquardt_global_exps_instr(float xincr, float **trans, | |

332 | int ndata, int ntrans, int fit_start, int fit_end, | |

333 | float instr[], int ninstr, | |

334 | noise_type noise, float sig[], int ftype, | |

335 | float **param, int paramfree[], int nparam, | |

336 | restrain_type restrain, float chisq_delta, | |

337 | float **fitted, float **residuals, | |

338 | float chisq_trans[], float *chisq_global, int *df, | |

339 | int drop_bad_transients) | |

340 | { | |

341 | float **covar, **alpha, *scaled_instr, instrsum; | |

342 | int i, j, ret; | |

343 | int mlocal, mglobal; | |

344 | float gparam[MAXFIT]; | |

345 | float *exp_pure, *exp_conv[MAXFIT]; | |

346 | // double time; | |

347 | ||

348 | // time=Timer(); | |

349 | ||

350 | /* Some basic parameter checks */ | |

351 | if (xincr <= 0) return -1; | |

352 | if (ntrans < 1) return -1; | |

353 | if (ndata < 1) return -1; | |

354 | if (fit_start < 0 || fit_end > ndata) return -1; | |

355 | // if (ninstr < 1) return -1; | |

356 | if (nparam < 1) return -1; | |

357 | ||

358 | switch (ftype) { | |

359 | case FIT_GLOBAL_MULTIEXP: | |

360 | if (nparam % 2 != 1) { | |

361 | dbgprintf(1, "global fitting: " | |

362 | "multiexp needs odd number of parameters\n"); | |

363 | return -1; | |

364 | } | |

365 | break; | |

366 | ||

367 | case FIT_GLOBAL_STRETCHEDEXP: | |

368 | if (nparam != 4) { | |

369 | dbgprintf(1, "global fitting: " | |

370 | "stretched exp needs precisely 4 parameters\n"); | |

371 | return -1; | |

372 | } | |

373 | break; | |

374 | ||

375 | default: | |

376 | dbgprintf(1, "global fitting: unknown fitting type\n"); | |

377 | return -1; | |

378 | } | |

379 | ||

380 | if ((covar = GCI_ecf_matrix(nparam, nparam)) == NULL) | |

381 | return -2; | |

382 | ||

383 | if ((alpha = GCI_ecf_matrix(nparam, nparam)) == NULL) { | |

384 | GCI_ecf_free_matrix(covar); | |

385 | return -3; | |

386 | } | |

387 | ||

388 | if ((scaled_instr = (float *) malloc(ninstr * sizeof(float))) == NULL && ninstr>1) { | |

389 | GCI_ecf_free_matrix(covar); | |

390 | GCI_ecf_free_matrix(alpha); | |

391 | return -4; | |

392 | } | |

393 | ||

394 | /* Also allocate space for the exp_pure and exp_conv arrays */ | |

395 | if ((exp_conv[0] = (float *) malloc(nparam * ndata * sizeof(float))) | |

396 | == NULL) { | |

397 | GCI_ecf_free_matrix(covar); | |

398 | GCI_ecf_free_matrix(alpha); | |

399 | free(scaled_instr); | |

400 | return -5; | |

401 | } | |

402 | ||

403 | exp_pure = exp_conv[0]; | |

404 | for (i=1; i<nparam; i++) | |

405 | exp_conv[i] = exp_conv[0] + i * ndata; | |

406 | ||

407 | /* Scale the instrument response */ | |

408 | for (i=0, instrsum=0; i<ninstr; i++) | |

409 | instrsum += instr[i]; | |

410 | if (instrsum == 0) instrsum=1.0; //return -6; | |

411 | for (i=0; i<ninstr; i++) | |

412 | scaled_instr[i] = instr[i] / instrsum; | |

413 | ||

414 | // printf("that took: %f secs.\n", Timer()-time); time=Timer(); | |

415 | dbgprintf(2, "About to enter step (1)\n"); | |

416 | ||

417 | /* Step (1): estimate the global taus */ | |

418 | ret = GCI_marquardt_global_exps_est_globals_instr( | |

419 | xincr, trans, ndata, ntrans, fit_start, fit_end, | |

420 | scaled_instr, ninstr, noise, sig, | |

421 | ftype, param, paramfree, nparam, gparam, restrain, chisq_delta, | |

422 | fitted[0], residuals[0], covar, alpha, chisq_global); | |

423 | ||

424 | if (ret != 0) { | |

425 | dbgprintf(1, "Step (1) failed, ret = %d\n", ret); | |

426 | GCI_ecf_free_matrix(covar); | |

427 | GCI_ecf_free_matrix(alpha); | |

428 | free(scaled_instr); | |

429 | free(exp_conv[0]); | |

430 | return -10 + ret; | |

431 | } | |

432 | ||

433 | // printf("that took: %f secs.\n", Timer()-time); time=Timer(); | |

434 | /* Copy the estimated global taus to the parameters array */ | |

435 | ||

436 | switch (ftype) | |

437 | { | |

438 | case FIT_GLOBAL_MULTIEXP: | |

439 | for (i=2; i<nparam; i+=2) | |

440 | for (j=0; j<ntrans; j++) | |

441 | param[j][i] = gparam[i]; | |

442 | break; | |

443 | ||

444 | case FIT_GLOBAL_STRETCHEDEXP: | |

445 | for (i=2; i<nparam; i++) /* param 2 = tau, param 3 = h */ | |

446 | for (j=0; j<ntrans; j++) | |

447 | param[j][i] = gparam[i]; | |

448 | break; | |

449 | ||

450 | default: | |

451 | dbgprintf(1, "global_exps_instr: please update me!\n"); | |

452 | return -1; | |

453 | } | |

454 | ||

455 | // printf("that took: %f secs.\n", Timer()-time); time=Timer(); | |

456 | dbgprintf(2, "About to enter step (2)\n"); | |

457 | ||

458 | /* Step (2): use these taus to estimate initial values for all of | |

459 | the individual transient parameters */ | |

460 | ret = GCI_marquardt_global_exps_est_params_instr( | |

461 | xincr, trans, ndata, ntrans, fit_start, fit_end, | |

462 | scaled_instr, ninstr, noise, sig, ftype, | |

463 | param, paramfree, nparam, restrain, chisq_delta, | |

464 | exp_pure, exp_conv, fitted, residuals, covar, alpha, | |

465 | chisq_trans, drop_bad_transients); | |

466 | ||

467 | if (ret != 0) { | |

468 | dbgprintf(1, "Step (2) failed, ret = %d\n", ret); | |

469 | GCI_ecf_free_matrix(covar); | |

470 | GCI_ecf_free_matrix(alpha); | |

471 | free(scaled_instr); | |

472 | free(exp_conv[0]); | |

473 | return -20 + ret; | |

474 | } | |

475 | ||

476 | // printf("that took: %f secs.\n", Timer()-time); time=Timer(); | |

477 | dbgprintf(2, "About to enter step (3)\n"); | |

478 | ||

479 | /* Step (3): now that we have estimates for initial values for all | |

480 | parameters, we can do the global Marquardt fitting. Note that | |

481 | covar and alpha are only provided as scratch space. */ | |

482 | ret = GCI_marquardt_global_exps_do_fit_instr( | |

483 | xincr, trans, ndata, ntrans, fit_start, fit_end, | |

484 | scaled_instr, ninstr, noise, sig, ftype, | |

485 | param, paramfree, nparam, restrain, chisq_delta, | |

486 | exp_pure, exp_conv, fitted, residuals, covar, alpha, | |

487 | chisq_trans, chisq_global, drop_bad_transients); | |

488 | ||

489 | GCI_ecf_free_matrix(covar); | |

490 | GCI_ecf_free_matrix(alpha); | |

491 | free(scaled_instr); | |

492 | free(exp_conv[0]); | |

493 | ||

494 | if (ret < 0) { | |

495 | dbgprintf(1, "Step (3) failed, ret = %d\n", ret); | |

496 | return -30 + ret; | |

497 | } | |

498 | ||

499 | // printf("that took: %f secs.\n", Timer()-time); time=Timer(); | |

500 | dbgprintf(2, "Step (3) succeeded, ret = %d\n", ret); | |

501 | ||

502 | /* Before we return, calculate the number of degrees of freedom */ | |

503 | /* The number of degrees of freedom is given by: | |

504 | d.f. = ntrans * ((fit_end - fit_start) - # free local parameters) | |

505 | - # free global parameters | |

506 | */ | |

507 | ||

508 | if (drop_bad_transients) { | |

509 | *df = 0; | |

510 | for (i=0; i<ntrans; i++) { | |

511 | if (chisq_trans[i] > 0) | |

512 | (*df)++; | |

513 | } | |

514 | } else | |

515 | *df = ntrans; | |

516 | ||

517 | mglobal = mlocal = 0; | |

518 | ||

519 | switch (ftype) | |

520 | { | |

521 | case FIT_GLOBAL_MULTIEXP: | |

522 | for (i=2; i<nparam; i+=2) | |

523 | if (paramfree[i]) mglobal++; | |

524 | for (i=1; i<nparam; i+=2) | |

525 | if (paramfree[i]) mlocal++; | |

526 | if (paramfree[0]) mlocal++; | |

527 | break; | |

528 | ||

529 | case FIT_GLOBAL_STRETCHEDEXP: | |

530 | if (paramfree[0]) mlocal++; | |

531 | if (paramfree[1]) mlocal++; | |

532 | if (paramfree[2]) mglobal++; | |

533 | if (paramfree[3]) mglobal++; | |

534 | break; | |

535 | ||

536 | default: | |

537 | dbgprintf(1, "global_exps_instr: please update me!\n"); | |

538 | return -1; | |

539 | } | |

540 | ||

541 | *df *= ((fit_end - fit_start) - mlocal); | |

542 | *df -= mglobal; | |

543 | ||

544 | // printf("mlocal %d\nmglobal %d\ndf %d\n", mlocal, mglobal, *df); | |

545 | ||

546 | // printf("that took: %f secs.\n", Timer()-time); time=Timer(); | |

547 | return ret; | |

548 | } | |

549 | ||

550 | ||

551 | int GCI_marquardt_global_exps_est_globals_instr( | |

552 | float xincr, float **trans, int ndata, int ntrans, | |

553 | int fit_start, int fit_end, float instr[], int ninstr, | |

554 | noise_type noise, float sig[], int ftype, | |

555 | float **param, int paramfree[], int nparam, float gparam[], | |

556 | restrain_type restrain, float chisq_delta, | |

557 | float fitted[], float residuals[], | |

558 | float **covar, float **alpha, float *chisq_global) | |

559 | { | |

560 | int i, j, ret, nparamfree; | |

561 | float *summed, *tptr; | |

562 | int data_start; | |

563 | float Z, A, tau; | |

564 | void (*fitfunc)(float, float [], float *, float [], int); | |

565 | ||

566 | if ((summed = (float *) calloc(ndata, sizeof(float))) == NULL) | |

567 | return -1; | |

568 | ||

569 | for (i=0; i<ntrans; i++) { | |

570 | tptr = trans[i]; | |

571 | for (j=0; j<ndata; j++) | |

572 | summed[j] += tptr[j]; | |

573 | } | |

574 | ||

575 | /* This code is now lifted from fitting.c, appropriately modified */ | |

576 | data_start = fit_start + ECF_Find_Float_Max(&summed[fit_start], | |

577 | fit_end - fit_start, &A); | |

578 | // ret = GCI_triple_integral_instr(xincr, summed, data_start, fit_end, | |

579 | // instr, ninstr, noise, sig, | |

580 | // &Z, &A, &tau, NULL, NULL, NULL); | |

581 | ||

582 | ret = GCI_triple_integral_fitting_engine(xincr, summed, data_start, fit_end, | |

583 | instr, ninstr, noise, sig, | |

584 | &Z, &A, &tau, NULL, NULL, NULL, 1.5*(fit_end-fit_start-3)); | |

585 | ||

586 | dbgprintf(3, "In est_globals_instr, triple integral ret = %d\n", ret); | |

587 | ||

588 | if (ret < 0) { | |

589 | Z = 0; | |

590 | ECF_Find_Float_Max(&summed[fit_start], fit_end - fit_start, &A); | |

591 | tau = 2.0; | |

592 | } | |

593 | ||

594 | ||

595 | /* We set gparam[] to be an array which holds initial estimates of | |

596 | the parameters for the _sum_ of all the transients, for those | |

597 | parameters which are fixed and therefore "known"; the rest will | |

598 | be estimated later. It doesn't matter if we set a few other | |

599 | values, as these will be overwritten later. We could also | |

600 | merge this switch() with the next one, but then the code would | |

601 | possibly be a little less easy to follow, so we won't. */ | |

602 | ||

603 | switch (ftype) { | |

604 | case FIT_GLOBAL_MULTIEXP: | |

605 | /* Z */ | |

606 | if (! paramfree[0]) | |

607 | { | |

608 | gparam[0] = 0; | |

609 | for (j=0; j<ntrans; j++) gparam[0] += param[j][0]; | |

610 | } | |

611 | ||

612 | /* A's */ | |

613 | for (i=1; i<nparam; i++) | |

614 | if (! paramfree[i]) | |

615 | { | |

616 | gparam[i] = 0; | |

617 | for (j=0; j<ntrans; j++) gparam[i] += param[j][i]; | |

618 | } | |

619 | ||

620 | /* taus last (was first) */ | |

621 | for (i=2; i<nparam; i+=2) | |

622 | gparam[i] = param[0][i]; | |

623 | ||

624 | break; | |

625 | ||

626 | case FIT_GLOBAL_STRETCHEDEXP: | |

627 | /* Z */ | |

628 | if (! paramfree[0]) { | |

629 | gparam[0] = 0; | |

630 | for (j=0; j<ntrans; j++) gparam[0] += param[j][0]; | |

631 | } | |

632 | ||

633 | /* A */ | |

634 | if (! paramfree[1]) { | |

635 | gparam[1] = 0; | |

636 | for (j=0; j<ntrans; j++) gparam[1] += param[j][1]; | |

637 | } | |

638 | ||

639 | /* tau and h last (were first) */ | |

640 | for (i=2; i<nparam; i++) | |

641 | gparam[i] = param[0][i]; | |

642 | ||

643 | break; | |

644 | ||

645 | default: | |

646 | dbgprintf(1, "global_exps_est_globals_instr: please update me!\n"); | |

647 | free(summed); | |

648 | return -1; | |

649 | } | |

650 | ||

651 | ||

652 | /* Now we can set any non-fixed parameters to meaningful initial | |

653 | estimates */ | |

654 | switch (ftype) { | |

655 | case FIT_GLOBAL_MULTIEXP: | |

656 | fitfunc = GCI_multiexp_tau; | |

657 | ||

658 | switch (nparam) { | |

659 | case 3: | |

660 | if (paramfree[0]) gparam[0] = Z; | |

661 | if (paramfree[1]) gparam[1] = A; | |

662 | if (paramfree[2]) gparam[2] = tau; | |

663 | break; | |

664 | ||

665 | case 5: | |

666 | if (paramfree[0]) gparam[0] = Z; | |

667 | if (paramfree[1]) gparam[1] = A*3/4; | |

668 | if (paramfree[2]) gparam[2] = tau; | |

669 | if (paramfree[3]) gparam[3] = A*1/4; | |

670 | if (paramfree[4]) gparam[4] = tau*2/3; | |

671 | break; | |

672 | ||

673 | default: | |

674 | if (nparam<7) { | |

675 | free(summed); | |

676 | return -2; | |

677 | } | |

678 | if (paramfree[0]) gparam[0] = Z; | |

679 | if (paramfree[1]) gparam[1] = A*3/4; | |

680 | if (paramfree[2]) gparam[2] = tau; | |

681 | if (paramfree[3]) gparam[3] = A*1/6; | |

682 | if (paramfree[4]) gparam[4] = tau*2/3; | |

683 | if (paramfree[5]) gparam[5] = A*1/6; | |

684 | if (paramfree[6]) gparam[6] = tau/3; | |

685 | for (i=7; i<nparam; i+=2) { | |

686 | if (paramfree[i]) gparam[i] = A/i; | |

687 | if (paramfree[i+1]) gparam[i+1] = tau/i; | |

688 | } | |

689 | break; | |

690 | } | |

691 | break; | |

692 | ||

693 | case FIT_GLOBAL_STRETCHEDEXP: | |

694 | fitfunc = GCI_stretchedexp; | |

695 | ||

696 | if (paramfree[0]) gparam[0] = Z; | |

697 | if (paramfree[1]) gparam[1] = A; | |

698 | if (paramfree[2]) gparam[2] = tau; | |

699 | if (paramfree[3]) gparam[3] = 1.5; /* h */ | |

700 | break; | |

701 | ||

702 | default: | |

703 | dbgprintf(1, "est_globals_instr: please update me!\n"); | |

704 | free(summed); | |

705 | return -1; | |

706 | } | |

707 | ||

708 | for (i=0, nparamfree=0; i<nparam; i++) if (paramfree[i]) nparamfree++; | |

709 | ||

710 | /* Note that the only values in the gparam array which are of | |

711 | interest are the taus and h for stretched exp, so we don't need | |

712 | to rescale Z and the A's back again */ | |

713 | // ret = GCI_marquardt_instr(xincr, summed, ndata, fit_start, fit_end, | |

714 | // instr, ninstr, noise, sig, | |

715 | // gparam, paramfree, nparam, restrain, fitfunc, | |

716 | // fitted, residuals, covar, alpha, chisq_global, | |

717 | // 0, NULL); | |

718 | ||

719 | ret = GCI_marquardt_fitting_engine(xincr, summed, ndata, fit_start, fit_end, | |

720 | instr, ninstr, noise, sig, | |

721 | gparam, paramfree, nparam, restrain, fitfunc, | |

722 | fitted, residuals, chisq_global, covar, alpha, | |

723 | NULL, 1.5*(fit_end-fit_start-nparamfree), chisq_delta, 0); | |

724 | ||

725 | dbgprintf(3, "In est_globals_instr, marquardt ret = %d\n", ret); | |

726 | ||

727 | free(summed); | |

728 | ||

729 | if (ret < 0) | |

730 | return -3; | |

731 | else | |

732 | return 0; | |

733 | } | |

734 | ||

735 | ||

736 | int GCI_marquardt_global_exps_est_params_instr( | |

737 | float xincr, float **trans, int ndata, int ntrans, | |

738 | int fit_start, int fit_end, float instr[], int ninstr, | |

739 | noise_type noise, float sig[], int ftype, | |

740 | float **param, int paramfree[], int nparam, restrain_type restrain, float chisq_delta, | |

741 | float exp_pure[], float *exp_conv[], | |

742 | float **fitted, float **residuals, | |

743 | float **covar, float **alpha, float chisq_trans[], | |

744 | int drop_bad_transients) | |

745 | { | |

746 | int i, j, sortkey[MAXFIT], paramfree_local[MAXFIT], tempi, ret; | |

747 | int data_start; | |

748 | float Z, A, tau; | |

749 | ||

750 | /* We begin by estimating the non-tau parameters for each | |

751 | transient. We do this by performing a triple-integral fit and | |

752 | using the Z and A resulting as a basis for our initial | |

753 | parameters. In the case of multiexponential fitting, we also | |

754 | sort the taus into decreasing order, and assume that the | |

755 | largest tau is the most significant component. */ | |

756 | ||

757 | // PRB 03/07 Although **fitted and **residuals are provided only one "transient" is required and used, fitted[0] and residuals[0] | |

758 | ||

759 | if (ftype == FIT_GLOBAL_MULTIEXP) { | |

760 | /* Initialise */ | |

761 | for (i=2; i<nparam; i+=2) | |

762 | sortkey[i] = i; | |

763 | ||

764 | /* Bubblesort :-) */ | |

765 | for (i=2; i<nparam; i+=2) | |

766 | for (j=2*(nparam/2); j>i; j-=2) /* nparam is odd */ | |

767 | if (param[0][sortkey[j]] > param[0][sortkey[j-2]]) { | |

768 | tempi = sortkey[j]; | |

769 | sortkey[j] = sortkey[j-2]; | |

770 | sortkey[j-2] = tempi; | |

771 | } | |

772 | } | |

773 | ||

774 | dbgprintf(3, "In est_params_instr, parameters initialised to:\n"); | |

775 | ||

776 | for (i=0; i<ntrans; i++) { | |

777 | /* This code is now lifted from fitting.c, appropriately modified */ | |

778 | data_start = fit_start + ECF_Find_Float_Max(&trans[i][fit_start], | |

779 | fit_end - fit_start, &A); | |

780 | // ret = GCI_triple_integral_instr(xincr, trans[i], data_start, fit_end, | |

781 | // instr, ninstr, noise, sig, | |

782 | // &Z, &A, &tau, NULL, NULL, NULL); | |

783 | ||

784 | ret = GCI_triple_integral_fitting_engine(xincr, trans[i], data_start, fit_end, | |

785 | instr, ninstr, noise, sig, | |

786 | &Z, &A, &tau, NULL, NULL, NULL, 1.5*(fit_end-fit_start-3)); | |

787 | if (ret < 0) { | |

788 | Z = 0; | |

789 | ECF_Find_Float_Max(&trans[i][fit_start], fit_end - fit_start, &A); | |

790 | } | |

791 | ||

792 | switch (ftype) { | |

793 | case FIT_GLOBAL_MULTIEXP: | |

794 | switch (nparam) { | |

795 | case 3: | |

796 | if (paramfree[0]) param[i][0] = Z; | |

797 | if (paramfree[1]) param[i][1] = A; | |

798 | break; | |

799 | ||

800 | case 5: | |

801 | if (paramfree[0]) param[i][0] = Z; | |

802 | if (paramfree[sortkey[2]-1]) param[i][sortkey[2]-1] = A*3/4; | |

803 | if (paramfree[sortkey[4]-1]) param[i][sortkey[4]-1] = A*1/4; | |

804 | break; | |

805 | ||

806 | default: | |

807 | if (nparam<7) { /* only actually need to do this once */ | |

808 | return -1; | |

809 | } | |

810 | if (paramfree[0]) param[i][0] = Z; | |

811 | if (paramfree[sortkey[2]-1]) param[i][sortkey[2]-1] = A*3/4; | |

812 | if (paramfree[sortkey[4]-1]) param[i][sortkey[4]-1] = A*1/6; | |

813 | if (paramfree[sortkey[6]-1]) param[i][sortkey[6]-1] = A*1/6; | |

814 | /* this is all pretty meaningless from here on, anyway */ | |

815 | for (j=8; j<nparam; j+=2) { | |

816 | if (paramfree[sortkey[j]-1]) param[i][sortkey[j]-1] = A/j; | |

817 | } | |

818 | break; | |

819 | } | |

820 | break; | |

821 | ||

822 | case FIT_GLOBAL_STRETCHEDEXP: | |

823 | if (paramfree[0]) param[i][0] = Z; | |

824 | if (paramfree[1]) param[i][1] = A; | |

825 | break; | |

826 | ||

827 | default: | |

828 | dbgprintf(1, "est_params_instr: please update me!\n"); | |

829 | return -1; | |

830 | } | |

831 | ||

832 | if (ECF_debug >= 3) { | |

833 | for (j=0; j<nparam; j++) | |

834 | dbgprintf(3, "param[%d][%d] = %.4g\n", i, j, param[i][j]); | |

835 | } | |

836 | } | |

837 | ||

838 | ||

839 | /* OK, the initial parameters are set up, we now do a Marquardt | |

840 | fit on all of the transients simultaneously to get more decent | |

841 | initial A and Z values. But we do this manually without | |

842 | recalculating the exponentials repeatedly. Furthermore, since | |

843 | the instrument response is convolved linearly with the | |

844 | exponentials, we can get by doing the convolution once only as | |

845 | well, making further major time savings. */ | |

846 | ||

847 | if (GCI_marquardt_global_exps_calculate_exps_instr( | |

848 | xincr, ndata, instr, ninstr, ftype, param[0], nparam, | |

849 | exp_pure, exp_conv) != 0) | |

850 | return -2; | |

851 | ||

852 | /* Now we can do a Marquardt fit on each of the transients */ | |

853 | ||

854 | /* Create a paramfree[] array which fixes all of the taus */ | |

855 | switch (ftype) { | |

856 | case FIT_GLOBAL_MULTIEXP: | |

857 | paramfree_local[0] = paramfree[0]; /* Z */ | |

858 | for (i=1; i<nparam; i+=2) | |

859 | paramfree_local[i] = paramfree[i]; /* the A's */ | |

860 | for (i=2; i<nparam; i+=2) | |

861 | paramfree_local[i] = 0; /* the taus */ | |

862 | break; | |

863 | ||

864 | case FIT_GLOBAL_STRETCHEDEXP: | |

865 | paramfree_local[0] = paramfree[0]; /* Z */ | |

866 | paramfree_local[1] = paramfree[1]; /* A */ | |

867 | paramfree_local[2] = 0; /* tau */ | |

868 | paramfree_local[3] = 0; /* h */ | |

869 | break; | |

870 | ||

871 | default: | |

872 | dbgprintf(1, "est_params_instr: please update me!\n"); | |

873 | return -1; | |

874 | } | |

875 | ||

876 | dbgprintf(3, "In est_params_instr, after do_fit_single, " | |

877 | "parameters initialised to:\n"); | |

878 | for (i=0; i<ntrans; i++) { | |

879 | ret = GCI_marquardt_global_exps_do_fit_single( | |

880 | xincr, trans[i], ndata, fit_start, fit_end, | |

881 | noise, sig, ftype, | |

882 | param[i], paramfree_local, nparam, restrain, chisq_delta, exp_conv, | |

883 | fitted[0], residuals[0], covar, alpha, &chisq_trans[i]); | |

884 | ||

885 | if (ret < 0) { | |

886 | if (drop_bad_transients) { | |

887 | dbgprintf(2, "In est_params_instr, transient %d gave " | |

888 | "do_fit_single return val %d; dropping it\n", | |

889 | i, ret); | |

890 | chisq_trans[i] = -1; | |

891 | continue; | |

892 | } else { | |

893 | dbgprintf(1, "In est_params_instr, transient %d gave " | |

894 | "do_fit_single return val %d\n", i, ret); | |

895 | return -10 + ret; | |

896 | } | |

897 | } | |

898 | ||

899 | /* We try a second time with these parameters if we got | |

900 | nonsense */ | |

901 | if (chisq_trans[i] > 20 * (fit_end - fit_start)) { | |

902 | ret = GCI_marquardt_global_exps_do_fit_single( | |

903 | xincr, trans[i], ndata, fit_start, fit_end, | |

904 | noise, sig, ftype, | |

905 | param[i], paramfree_local, nparam, restrain, chisq_delta, exp_conv, | |

906 | fitted[0], residuals[0], covar, alpha, &chisq_trans[i]); | |

907 | ||

908 | /* Improved? */ | |

909 | if (ret < 0 || chisq_trans[i] > 20 * (fit_end - fit_start)) { | |

910 | if (drop_bad_transients) { | |

911 | dbgprintf(2, "In est_params_instr, transient %d gave " | |

912 | "do_fit_single return val %d, chisq value %.3f; " | |

913 | "dropping it\n", i, ret, chisq_trans[i]); | |

914 | chisq_trans[i] = -1; | |

915 | continue; | |

916 | } else { | |

917 | dbgprintf(1, "In est_params_instr, transient %d gave " | |

918 | "do_fit_single return val %d, " | |

919 | "chisq value %.3f\n", i, ret, chisq_trans[i]); | |

920 | return -10 + ret; | |

921 | } | |

922 | } | |

923 | ||

924 | if (ECF_debug >= 3) { | |

925 | for (j=0; j<nparam; j++) | |

926 | dbgprintf(3, "param[%d][%d] = %.4g\n", i, j, param[i][j]); | |

927 | } | |

928 | } | |

929 | } | |

930 | ||

931 | return 0; | |

932 | } | |

933 | ||

934 | ||

935 | /* This finds values of exp(-x/tau) and x*exp(-x/tau)/tau^2, which are | |

936 | needed later for the multiexponential case, and finds the | |

937 | equivalents in the stretched exponential case. */ | |

938 | // should also now handle no instrument response, i.e. instr=NULL | |

939 | ||

940 | int GCI_marquardt_global_exps_calculate_exps_instr( | |

941 | float xincr, int ndata, float instr[], int ninstr, | |

942 | int ftype, float param[], int nparam, | |

943 | float exp_pure[], float *exp_conv[]) | |

944 | { | |

945 | int i, j, k; | |

946 | int convpts; | |

947 | double excur; /* exp(-x/tau) */ | |

948 | double exincr; /* exp(-xincr/tau) */ | |

949 | float *expi; | |

950 | float ex, lxa, xah, a2inv, a3inv; /* for stetched exp */ | |

951 | double xa, xaincr; | |

952 | ||

953 | switch (ftype) { | |

954 | case FIT_GLOBAL_MULTIEXP: | |

955 | /* Not quite the most efficient way to do this, but not bad */ | |

956 | /* First we calculate the exp(-x/tau) array */ | |

957 | for (i=2; i<nparam; i+=2) { | |

958 | expi = exp_conv[i]; | |

959 | ||

960 | excur = 1.0; | |

961 | exincr = exp(-xincr/(double)param[i]); | |

962 | for (j=0; j<ndata; j++) { | |

963 | exp_pure[j] = excur; | |

964 | excur *= exincr; | |

965 | ||

966 | /* And convolve the exponentials with the instrument response if possible */ | |

967 | expi[j] = 0; | |

968 | convpts = (ninstr <= j) ? ninstr-1 : j; | |

969 | if (convpts<=0 || instr==NULL) expi[j] = exp_pure[j]; | |

970 | else for (k=0; k<=convpts; k++) expi[j] += exp_pure[j-k]*instr[k]; | |

971 | } | |

972 | } | |

973 | ||

974 | /* Now we repeat the exercise for x*exp(-x/tau) / tau^2 */ | |

975 | for (i=2; i<nparam; i+=2) { | |

976 | expi = exp_conv[i-1]; | |

977 | ||

978 | excur = 1.0 / (param[i]*param[i]); /* 1/tau^2 */ | |

979 | exincr = exp(-xincr/(double)param[i]); | |

980 | for (j=0; j<ndata; j++) { | |

981 | exp_pure[j] = (xincr*i) * excur; /* x*exp(-x/tau) / tau^2 */ | |

982 | excur *= exincr; | |

983 | ||

984 | /* And convolve the exponentials with the instrument response if possible */ | |

985 | expi[j] = 0; | |

986 | convpts = (ninstr <= j) ? ninstr-1 : j; | |

987 | if (convpts<=0 || instr==NULL) expi[j] = exp_pure[j]; | |

988 | else for (k=0; k<=convpts; k++) expi[j] += exp_pure[j-k]*instr[k]; | |

989 | } | |

990 | } | |

991 | break; | |

992 | ||

993 | case FIT_GLOBAL_STRETCHEDEXP: | |

994 | /* We start by essentially repeating the stretchedexp_array | |

995 | function */ | |

996 | xa=0; | |

997 | xaincr = xincr / param[2]; | |

998 | a2inv = 1/param[2]; | |

999 | a3inv = 1/param[3]; | |

1000 | ||

1001 | /* When x=0 */ | |

1002 | exp_conv[1][0] = 1.0; | |

1003 | exp_conv[2][0] = exp_conv[3][0] = 0.0; | |

1004 | ||

1005 | for (i=1; i<ndata; i++) { | |

1006 | xa += xaincr; /* xa = (xincr*i)/param[2] */ | |

1007 | lxa = log(xa); /* lxa = log(x/param[2]) */ | |

1008 | xah = exp(lxa * a3inv); /* xah = exp(log(x/param[2])/param[3]) | |

1009 | = (x/param[2])^(1/param[3]) */ | |

1010 | exp_conv[1][i] = ex = exp(-xah); | |

1011 | /* ex = exp(-(x/param[2])^(1/param[3])) */ | |

1012 | ex *= xah * a3inv; /* ex = exp(...) * (x/param[2])^(1/param[3]) * | |

1013 | 1/param[3] */ | |

1014 | exp_conv[2][i] = ex * a2inv; | |

1015 | exp_conv[3][i] = ex * lxa * a3inv; | |

1016 | } | |

1017 | ||

1018 | if (ninstr>0 && instr!=NULL) // else exp_conv already contains the pure data | |

1019 | { | |

1020 | /* Now convolve these with the instrument response */ | |

1021 | for (i=1; i<4; i++) | |

1022 | { | |

1023 | expi = exp_conv[i]; | |

1024 | ||

1025 | for (j=0; j<ndata; j++) | |

1026 | exp_pure[j] = expi[j]; /* save the array in temp storage */ | |

1027 | ||

1028 | for (j=0; j<ndata; j++) | |

1029 | { | |

1030 | expi[j] = 0; | |

1031 | convpts = (ninstr <= j) ? ninstr-1 : j; | |

1032 | for (k=0; k<=convpts; k++) | |

1033 | expi[j] += exp_pure[j-k]*instr[k]; | |

1034 | } | |

1035 | } | |

1036 | } | |

1037 | break; | |

1038 | ||

1039 | default: | |

1040 | dbgprintf(1, "calculate_exps_instr: please update me!\n"); | |

1041 | return -1; | |

1042 | } | |

1043 | ||

1044 | return 0; | |

1045 | } | |

1046 | ||

1047 | ||

1048 | /* This is just like the normal GCI_marquardt_instr, except that it | |

1049 | is designed for the multiexp case where we provide the exponential | |

1050 | decays in advance, and where we don't care about error axes */ | |

1051 | int GCI_marquardt_global_exps_do_fit_single( | |

1052 | float xincr, float y[], int ndata, int fit_start, int fit_end, | |

1053 | noise_type noise, float sig[], int ftype, | |

1054 | float param[], int paramfree[], int nparam, restrain_type restrain, | |

1055 | float chisq_delta, float *exp_conv[], float *fitted, float *residuals, | |

1056 | float **covar, float **alpha, float *chisq) | |

1057 | { | |

1058 | float alambda, ochisq; | |

1059 | int mfit; | |

1060 | int i, k, itst, itst_max; | |

1061 | ||

1062 | itst_max = (restrain == ECF_RESTRAIN_DEFAULT) ? 4 : 6; | |

1063 | ||

1064 | mfit = 0; | |

1065 | for (i=0; i<nparam; i++) { | |

1066 | if (paramfree[i]) | |

1067 | mfit++; | |

1068 | } | |

1069 | ||

1070 | alambda = -1; | |

1071 | if (GCI_marquardt_global_exps_single_step( | |

1072 | xincr, y, ndata, fit_start, fit_end, | |

1073 | noise, sig, ftype, param, paramfree, nparam, restrain, | |

1074 | exp_conv, fitted, residuals, covar, alpha, | |

1075 | chisq, &alambda) != 0) { | |

1076 | return -1; | |

1077 | } | |

1078 | ||

1079 | k = 1; /* Iteration counter */ | |

1080 | itst = 0; | |

1081 | for (;;) { | |

1082 | k++; | |

1083 | if (k > MAXITERS) { | |

1084 | return -2; | |

1085 | } | |

1086 | ||

1087 | ochisq = *chisq; | |

1088 | if (GCI_marquardt_global_exps_single_step( | |

1089 | xincr, y, ndata, fit_start, fit_end, | |

1090 | noise, sig, ftype, param, paramfree, nparam, restrain, | |

1091 | exp_conv, fitted, residuals, covar, alpha, | |

1092 | chisq, &alambda) != 0) { | |

1093 | return -3; | |

1094 | } | |

1095 | ||

1096 | if (*chisq > ochisq) | |

1097 | itst = 0; | |

1098 | else if (ochisq - *chisq < chisq_delta) | |

1099 | itst++; | |

1100 | ||

1101 | if (itst < itst_max) continue; | |

1102 | ||

1103 | /* Endgame; this also handles correcting the chi-squared values */ | |

1104 | alambda=0.0; | |

1105 | if (GCI_marquardt_global_exps_single_step( | |

1106 | xincr, y, ndata, fit_start, fit_end, | |

1107 | noise, sig, ftype, param, paramfree, nparam, restrain, | |

1108 | exp_conv, fitted, residuals, covar, alpha, | |

1109 | chisq, &alambda) != 0) { | |

1110 | return -4; | |

1111 | } | |

1112 | ||

1113 | return k; /* We're done now */ | |

1114 | } | |

1115 | } | |

1116 | ||

1117 | ||

1118 | /* And this one is basically a specialised GCI_marquardt_instr_step */ | |

1119 | int GCI_marquardt_global_exps_single_step( | |

1120 | float xincr, float y[], | |

1121 | int ndata, int fit_start, int fit_end, | |

1122 | noise_type noise, float sig[], int ftype, | |

1123 | float param[], int paramfree[], int nparam, | |

1124 | restrain_type restrain, | |

1125 | float *exp_conv[], | |

1126 | float yfit[], float dy[], | |

1127 | float **covar, float **alpha, float *chisq, | |

1128 | float *alambda) | |

1129 | { | |

1130 | int j, k, l, ret; | |

1131 | static int mfit; | |

1132 | static float ochisq, paramtry[MAXFIT], beta[MAXFIT], dparam[MAXFIT]; | |

1133 | static void (*fitfunc)(float, float [], float *, float [], int); | |

1134 | ||

1135 | if (nparam > MAXFIT) | |

1136 | return -10; | |

1137 | if (xincr <= 0) | |

1138 | return -11; | |

1139 | if (fit_start < 0 || fit_start > fit_end || fit_end > ndata) | |

1140 | return -12; | |

1141 | ||

1142 | /* Initialisation */ | |

1143 | /* We assume we're given sensible starting values for param[] */ | |

1144 | if (*alambda < 0.0) { | |

1145 | switch (ftype) { | |

1146 | case FIT_GLOBAL_MULTIEXP: | |

1147 | fitfunc = GCI_multiexp_tau; | |

1148 | break; | |

1149 | ||

1150 | case FIT_GLOBAL_STRETCHEDEXP: | |

1151 | fitfunc = GCI_stretchedexp; | |

1152 | break; | |

1153 | ||

1154 | default: | |

1155 | dbgprintf(1, "exps_single_step: please update me!\n"); | |

1156 | return -1; | |

1157 | } | |

1158 | ||

1159 | for (mfit=0, j=0; j<nparam; j++) | |

1160 | if (paramfree[j]) | |

1161 | mfit++; | |

1162 | ||

1163 | if (GCI_marquardt_global_compute_exps_fn( | |

1164 | xincr, y, ndata, fit_start, fit_end, noise, sig, | |

1165 | ftype, param, paramfree, nparam, exp_conv, | |

1166 | yfit, dy, alpha, beta, chisq) != 0) | |

1167 | return -2; | |

1168 | ||

1169 | *alambda = 0.001; | |

1170 | ochisq = *chisq; | |

1171 | for (j=0; j<nparam; j++) | |

1172 | paramtry[j] = param[j]; | |

1173 | } | |

1174 | ||

1175 | /* Alter linearised fitting matrix by augmenting diagonal elements */ | |

1176 | for (j=0; j<mfit; j++) { | |

1177 | for (k=0; k<mfit; k++) | |

1178 | covar[j][k] = alpha[j][k]; | |

1179 | covar[j][j] = alpha[j][j] * (1.0 + (*alambda)); | |

1180 | dparam[j] = beta[j]; | |

1181 | } | |

1182 | ||

1183 | /* Matrix solution; GCI_gauss_jordan solves Ax=b rather than AX=B */ | |

1184 | if (GCI_gauss_jordan(covar, mfit, dparam) != 0) | |

1185 | return -1; | |

1186 | ||

1187 | /* Once converged, calculate corrected chi-squared values */ | |

1188 | if (*alambda == 0) { | |

1189 | if (GCI_marquardt_global_compute_exps_fn_final( | |

1190 | xincr, y, ndata, fit_start, fit_end, noise, sig, | |

1191 | ftype, param, paramfree, nparam, exp_conv, | |

1192 | yfit, dy, chisq) != 0) | |

1193 | return -4; | |

1194 | /* Irrelevant */ | |

1195 | // if (mfit < nparam) { /* no need to do this otherwise */ | |

1196 | // GCI_covar_sort(covar, nparam, paramfree, mfit); | |

1197 | // GCI_covar_sort(alpha, nparam, paramfree, mfit); | |

1198 | // } | |

1199 | return 0; | |

1200 | } | |

1201 | ||

1202 | /* Did the trial succeed? */ | |

1203 | for (j=0, l=0; l<nparam; l++) | |

1204 | if (paramfree[l]) | |

1205 | paramtry[l] = param[l] + dparam[j++]; | |

1206 | ||

1207 | if (restrain == ECF_RESTRAIN_DEFAULT) | |

1208 | ret = check_ecf_params (paramtry, nparam, fitfunc); | |

1209 | else | |

1210 | ret = check_ecf_user_params (paramtry, nparam, fitfunc); | |

1211 | ||

1212 | if (ret != 0) { | |

1213 | /* Bad parameters, increase alambda and return */ | |

1214 | *alambda *= 10.0; | |

1215 | return 0; | |

1216 | } | |

1217 | ||

1218 | if (GCI_marquardt_global_compute_exps_fn( | |

1219 | xincr, y, ndata, fit_start, fit_end, noise, sig, | |

1220 | ftype, paramtry, paramfree, nparam, exp_conv, | |

1221 | yfit, dy, covar, dparam, chisq) != 0) | |

1222 | return -2; | |

1223 | ||

1224 | /* Success, accept the new solution */ | |

1225 | if (*chisq < ochisq) { | |

1226 | *alambda *= 0.1; | |

1227 | ochisq = *chisq; | |

1228 | for (j=0; j<mfit; j++) { | |

1229 | for (k=0; k<mfit; k++) | |

1230 | alpha[j][k] = covar[j][k]; | |

1231 | beta[j] = dparam[j]; | |

1232 | } | |

1233 | for (l=0; l<nparam; l++) | |

1234 | param[l] = paramtry[l]; | |

1235 | } else { /* Failure, increase alambda and return */ | |

1236 | *alambda *= 10.0; | |

1237 | *chisq = ochisq; | |

1238 | } | |

1239 | ||

1240 | return 0; | |

1241 | } | |

1242 | ||

1243 | ||

1244 | /* This is a streamlined GCI_marquardt_compute_fn_instr */ | |

1245 | int GCI_marquardt_global_compute_exps_fn( | |

1246 | float xincr, float y[], | |

1247 | int ndata, int fit_start, int fit_end, | |

1248 | noise_type noise, float sig[], int ftype, | |

1249 | float param[], int paramfree[], int nparam, | |

1250 | float *exp_conv[], | |

1251 | float yfit[], float dy[], | |

1252 | float **alpha, float *beta, float *chisq) | |

1253 | { | |

1254 | int i, j, k, l, m, mfit; | |

1255 | float wt, sig2i, y_ymod, dy_dparam[MAXFIT]; | |

1256 | ||

1257 | for (j=0, mfit=0; j<nparam; j++) | |

1258 | if (paramfree[j]) | |

1259 | mfit++; | |

1260 | ||

1261 | /* Initialise (symmetric) alpha, beta */ | |

1262 | for (j=0; j<mfit; j++) { | |

1263 | for (k=0; k<=j; k++) | |

1264 | alpha[j][k] = 0.0; | |

1265 | beta[j] = 0.0; | |

1266 | } | |

1267 | ||

1268 | /* Calculation of the fitting data will depend upon the type of | |

1269 | noise. Since there's no convolution involved here, this is | |

1270 | very easy. */ | |

1271 | ||

1272 | switch (ftype) { | |

1273 | case FIT_GLOBAL_MULTIEXP: | |

1274 | switch (noise) { | |

1275 | case NOISE_CONST: | |

1276 | *chisq = 0.0; | |

1277 | /* Summation loop over all data */ | |

1278 | for (i=fit_start; i<fit_end; i++) { | |

1279 | yfit[i] = param[0]; /* Z */ | |

1280 | dy_dparam[0] = 1.0; | |

1281 | ||

1282 | for (j=1; j<nparam; j+=2) { | |

1283 | yfit[i] += param[j] * exp_conv[j+1][i]; | |

1284 | /* A_j . exp(-x/tau_j) */ | |

1285 | dy_dparam[j] = exp_conv[j+1][i]; | |

1286 | /* exp(-x/tau_j) */ | |

1287 | dy_dparam[j+1] = param[j] * exp_conv[j][i]; | |

1288 | /* A_j * x * exp(-x/tau_j) / tau_j^2 */ | |

1289 | } | |

1290 | dy[i] = y[i] - yfit[i]; | |

1291 | ||

1292 | for (j=0, l=0; l<nparam; l++) { | |

1293 | if (paramfree[l]) { | |

1294 | wt = dy_dparam[l]; | |

1295 | for (k=0, m=0; m<=l; m++) | |

1296 | if (paramfree[m]) | |

1297 | alpha[j][k++] += wt * dy_dparam[m]; | |

1298 | beta[j] += wt * dy[i]; | |

1299 | j++; | |

1300 | } | |

1301 | } | |

1302 | /* And find chi^2 */ | |

1303 | *chisq += dy[i] * dy[i]; | |

1304 | } | |

1305 | ||

1306 | /* Now divide everything by sigma^2 */ | |

1307 | sig2i = 1.0 / (sig[0] * sig[0]); | |

1308 | *chisq *= sig2i; | |

1309 | for (j=0; j<mfit; j++) { | |

1310 | for (k=0; k<=j; k++) | |

1311 | alpha[j][k] *= sig2i; | |

1312 | beta[j] *= sig2i; | |

1313 | } | |

1314 | break; | |

1315 | ||

1316 | case NOISE_GIVEN: /* This is essentially the NR version */ | |

1317 | *chisq = 0.0; | |

1318 | /* Summation loop over all data */ | |

1319 | for (i=fit_start; i<fit_end; i++) { | |

1320 | yfit[i] = param[0]; /* Z */ | |

1321 | dy_dparam[0] = 1.0; | |

1322 | ||

1323 | for (j=1; j<nparam; j+=2) { | |

1324 | yfit[i] += param[j] * exp_conv[j+1][i]; | |

1325 | /* A_j . exp(-x/tau_j) */ | |

1326 | dy_dparam[j] = exp_conv[j+1][i]; | |

1327 | /* exp(-x/tau_j) */ | |

1328 | dy_dparam[j+1] = param[j] * exp_conv[j][i]; | |

1329 | /* A_j * x * exp(-x/tau_j) / tau_j^2 */ | |

1330 | } | |

1331 | sig2i = 1.0 / (sig[i] * sig[i]); | |

1332 | dy[i] = y[i] - yfit[i]; | |

1333 | ||

1334 | for (j=0, l=0; l<nparam; l++) { | |

1335 | if (paramfree[l]) { | |

1336 | wt = dy_dparam[l] * sig2i; | |

1337 | for (k=0, m=0; m<=l; m++) | |

1338 | if (paramfree[m]) | |

1339 | alpha[j][k++] += wt * dy_dparam[m]; | |

1340 | beta[j] += wt * dy[i]; | |

1341 | j++; | |

1342 | } | |

1343 | } | |

1344 | /* And find chi^2 */ | |

1345 | *chisq += dy[i] * dy[i] * sig2i; | |

1346 | } | |

1347 | break; | |

1348 | ||

1349 | case NOISE_POISSON_DATA: | |

1350 | *chisq = 0.0; | |

1351 | /* Summation loop over all data */ | |

1352 | for (i=fit_start; i<fit_end; i++) { | |

1353 | yfit[i] = param[0]; /* Z */ | |

1354 | dy_dparam[0] = 1.0; | |

1355 | ||

1356 | for (j=1; j<nparam; j+=2) { | |

1357 | yfit[i] += param[j] * exp_conv[j+1][i]; | |

1358 | /* A_j . exp(-x/tau_j) */ | |

1359 | dy_dparam[j] = exp_conv[j+1][i]; | |

1360 | /* exp(-x/tau_j) */ | |

1361 | dy_dparam[j+1] = param[j] * exp_conv[j][i]; | |

1362 | /* A_j * x * exp(-x/tau_j) / tau_j^2 */ | |

1363 | } | |

1364 | sig2i = (y[i] > 15 ? 1.0/y[i] : 1.0/15); | |

1365 | dy[i] = y[i] - yfit[i]; | |

1366 | ||

1367 | for (j=0, l=0; l<nparam; l++) { | |

1368 | if (paramfree[l]) { | |

1369 | wt = dy_dparam[l] * sig2i; | |

1370 | for (k=0, m=0; m<=l; m++) | |

1371 | if (paramfree[m]) | |

1372 | alpha[j][k++] += wt * dy_dparam[m]; | |

1373 | beta[j] += wt * dy[i]; | |

1374 | j++; | |

1375 | } | |

1376 | } | |

1377 | /* And find chi^2 */ | |

1378 | *chisq += dy[i] * dy[i] * sig2i; | |

1379 | } | |

1380 | break; | |

1381 | ||

1382 | case NOISE_POISSON_FIT: | |

1383 | *chisq = 0.0; | |

1384 | /* Summation loop over all data */ | |

1385 | for (i=fit_start; i<fit_end; i++) { | |

1386 | yfit[i] = param[0]; /* Z */ | |

1387 | dy_dparam[0] = 1.0; | |

1388 | ||

1389 | for (j=1; j<nparam; j+=2) { | |

1390 | yfit[i] += param[j] * exp_conv[j+1][i]; | |

1391 | /* A_j . exp(-x/tau_j) */ | |

1392 | dy_dparam[j] = exp_conv[j+1][i]; | |

1393 | /* exp(-x/tau_j) */ | |

1394 | dy_dparam[j+1] = param[j] * exp_conv[j][i]; | |

1395 | /* A_j * x * exp(-x/tau_j) / tau_j^2 */ | |

1396 | } | |

1397 | sig2i = (yfit[i] > 15 ? 1.0/yfit[i] : 1.0/15); | |

1398 | dy[i] = y[i] - yfit[i]; | |

1399 | ||

1400 | for (j=0, l=0; l<nparam; l++) { | |

1401 | if (paramfree[l]) { | |

1402 | wt = dy_dparam[l] * sig2i; | |

1403 | for (k=0, m=0; m<=l; m++) | |

1404 | if (paramfree[m]) | |

1405 | alpha[j][k++] += wt * dy_dparam[m]; | |

1406 | beta[j] += wt * dy[i]; | |

1407 | j++; | |

1408 | } | |

1409 | } | |

1410 | /* And find chi^2 */ | |

1411 | *chisq += dy[i] * dy[i] * sig2i; | |

1412 | } | |

1413 | break; | |

1414 | ||

1415 | case NOISE_GAUSSIAN_FIT: | |

1416 | *chisq = 0.0; | |

1417 | /* Summation loop over all data */ | |

1418 | for (i=fit_start; i<fit_end; i++) { | |

1419 | yfit[i] = param[0]; /* Z */ | |

1420 | dy_dparam[0] = 1.0; | |

1421 | ||

1422 | for (j=1; j<nparam; j+=2) { | |

1423 | yfit[i] += param[j] * exp_conv[j+1][i]; | |

1424 | /* A_j . exp(-x/tau_j) */ | |

1425 | dy_dparam[j] = exp_conv[j+1][i]; | |

1426 | /* exp(-x/tau_j) */ | |

1427 | dy_dparam[j+1] = param[j] * exp_conv[j][i]; | |

1428 | /* A_j * x * exp(-x/tau_j) / tau_j^2 */ | |

1429 | } | |

1430 | sig2i = (yfit[i] > 1 ? 1.0/yfit[i] : 1.0); | |

1431 | dy[i] = y[i] - yfit[i]; | |

1432 | ||

1433 | for (j=0, l=0; l<nparam; l++) { | |

1434 | if (paramfree[l]) { | |

1435 | wt = dy_dparam[l] * sig2i; | |

1436 | for (k=0, m=0; m<=l; m++) | |

1437 | if (paramfree[m]) | |

1438 | alpha[j][k++] += wt * dy_dparam[m]; | |

1439 | beta[j] += wt * dy[i]; | |

1440 | j++; | |

1441 | } | |

1442 | } | |

1443 | /* And find chi^2 */ | |

1444 | *chisq += dy[i] * dy[i] * sig2i; | |

1445 | } | |

1446 | break; | |

1447 | ||

1448 | case NOISE_MLE: | |

1449 | *chisq = 0.0; | |

1450 | /* Summation loop over all data */ | |

1451 | for (i=fit_start; i<fit_end; i++) { | |

1452 | yfit[i] = param[0]; /* Z */ | |

1453 | dy_dparam[0] = 1.0; | |

1454 | ||

1455 | for (j=1; j<nparam; j+=2) { | |

1456 | yfit[i] += param[j] * exp_conv[j+1][i]; | |

1457 | /* A_j . exp(-x/tau_j) */ | |

1458 | dy_dparam[j] = exp_conv[j+1][i]; | |

1459 | /* exp(-x/tau_j) */ | |

1460 | dy_dparam[j+1] = param[j] * exp_conv[j][i]; | |

1461 | /* A_j * x * exp(-x/tau_j) / tau_j^2 */ | |

1462 | } | |

1463 | sig2i = (yfit[i] > 1 ? 1.0/yfit[i] : 1.0); | |

1464 | dy[i] = y[i] - yfit[i]; | |

1465 | y_ymod=y[i]/yfit[i]; | |

1466 | ||

1467 | for (j=0, l=0; l<nparam; l++) { | |

1468 | if (paramfree[l]) { | |

1469 | wt = dy_dparam[l] * sig2i; | |

1470 | for (k=0, m=0; m<=l; m++) | |

1471 | if (paramfree[m]) | |

1472 | alpha[j][k++] += wt*dy_dparam[m]*y_ymod; // wt * dy_dparam[m]; | |

1473 | beta[j] += wt * dy[i]; | |

1474 | j++; | |

1475 | } | |

1476 | } | |

1477 | // And find chi^2 | |

1478 | if (yfit[i]<=0.0) | |

1479 | ; // do nothing | |

1480 | else if (y[i]==0.0) | |

1481 | *chisq += 2.0*yfit[i]; // to avoid NaN from log | |

1482 | else | |

1483 | *chisq += 2.0*(yfit[i]-y[i]) - 2.0*y[i]*log(yfit[i]/y[i]); // was dy[i] * dy[i] * sig2i; | |

1484 | } | |

1485 | if (*chisq <= 0.0) *chisq = 1.0e308; // don't let chisq=0 through yfit being all -ve | |

1486 | break; | |

1487 | ||

1488 | default: | |

1489 | return -3; | |

1490 | /* break; */ // (unreachable code) | |

1491 | } | |

1492 | break; | |

1493 | ||

1494 | case FIT_GLOBAL_STRETCHEDEXP: | |

1495 | switch (noise) { | |

1496 | case NOISE_CONST: | |

1497 | *chisq = 0.0; | |

1498 | /* Summation loop over all data */ | |

1499 | for (i=fit_start; i<fit_end; i++) { | |

1500 | yfit[i] = param[0] + param[1] * exp_conv[1][i]; | |

1501 | dy[i] = y[i] - yfit[i]; | |

1502 | ||

1503 | dy_dparam[0] = 1.0; | |

1504 | dy_dparam[1] = exp_conv[1][i]; | |

1505 | dy_dparam[2] = param[1] * exp_conv[2][i]; | |

1506 | dy_dparam[3] = param[1] * exp_conv[3][i]; | |

1507 | ||

1508 | for (j=0, l=0; l<nparam; l++) { | |

1509 | if (paramfree[l]) { | |

1510 | wt = dy_dparam[l]; | |

1511 | for (k=0, m=0; m<=l; m++) | |

1512 | if (paramfree[m]) | |

1513 | alpha[j][k++] += wt * dy_dparam[m]; | |

1514 | beta[j] += wt * dy[i]; | |

1515 | j++; | |

1516 | } | |

1517 | } | |

1518 | /* And find chi^2 */ | |

1519 | *chisq += dy[i] * dy[i]; | |

1520 | } | |

1521 | ||

1522 | /* Now divide everything by sigma^2 */ | |

1523 | sig2i = 1.0 / (sig[0] * sig[0]); | |

1524 | *chisq *= sig2i; | |

1525 | for (j=0; j<mfit; j++) { | |

1526 | for (k=0; k<=j; k++) | |

1527 | alpha[j][k] *= sig2i; | |

1528 | beta[j] *= sig2i; | |

1529 | } | |

1530 | break; | |

1531 | ||

1532 | case NOISE_GIVEN: /* This is essentially the NR version */ | |

1533 | *chisq = 0.0; | |

1534 | /* Summation loop over all data */ | |

1535 | for (i=fit_start; i<fit_end; i++) { | |

1536 | yfit[i] = param[0] + param[1] * exp_conv[1][i]; | |

1537 | dy[i] = y[i] - yfit[i]; | |

1538 | ||

1539 | dy_dparam[0] = 1.0; | |

1540 | dy_dparam[1] = exp_conv[1][i]; | |

1541 | dy_dparam[2] = param[1] * exp_conv[2][i]; | |

1542 | dy_dparam[3] = param[1] * exp_conv[3][i]; | |

1543 | ||

1544 | sig2i = 1.0 / (sig[i] * sig[i]); | |

1545 | ||

1546 | for (j=0, l=0; l<nparam; l++) { | |

1547 | if (paramfree[l]) { | |

1548 | wt = dy_dparam[l] * sig2i; | |

1549 | for (k=0, m=0; m<=l; m++) | |

1550 | if (paramfree[m]) | |

1551 | alpha[j][k++] += wt * dy_dparam[m]; | |

1552 | beta[j] += wt * dy[i]; | |

1553 | j++; | |

1554 | } | |

1555 | } | |

1556 | /* And find chi^2 */ | |

1557 | *chisq += dy[i] * dy[i] * sig2i; | |

1558 | } | |

1559 | break; | |

1560 | ||

1561 | case NOISE_POISSON_DATA: | |

1562 | *chisq = 0.0; | |

1563 | /* Summation loop over all data */ | |

1564 | for (i=fit_start; i<fit_end; i++) { | |

1565 | yfit[i] = param[0] + param[1] * exp_conv[1][i]; | |

1566 | dy[i] = y[i] - yfit[i]; | |

1567 | ||

1568 | dy_dparam[0] = 1.0; | |

1569 | dy_dparam[1] = exp_conv[1][i]; | |

1570 | dy_dparam[2] = param[1] * exp_conv[2][i]; | |

1571 | dy_dparam[3] = param[1] * exp_conv[3][i]; | |

1572 | ||

1573 | sig2i = (y[i] > 15 ? 1.0/y[i] : 1.0/15); | |

1574 | ||

1575 | for (j=0, l=0; l<nparam; l++) { | |

1576 | if (paramfree[l]) { | |

1577 | wt = dy_dparam[l] * sig2i; | |

1578 | for (k=0, m=0; m<=l; m++) | |

1579 | if (paramfree[m]) | |

1580 | alpha[j][k++] += wt * dy_dparam[m]; | |

1581 | beta[j] += wt * dy[i]; | |

1582 | j++; | |

1583 | } | |

1584 | } | |

1585 | /* And find chi^2 */ | |

1586 | *chisq += dy[i] * dy[i] * sig2i; | |

1587 | } | |

1588 | break; | |

1589 | ||

1590 | case NOISE_POISSON_FIT: | |

1591 | *chisq = 0.0; | |

1592 | /* Summation loop over all data */ | |

1593 | for (i=fit_start; i<fit_end; i++) { | |

1594 | yfit[i] = param[0] + param[1] * exp_conv[1][i]; | |

1595 | dy[i] = y[i] - yfit[i]; | |

1596 | ||

1597 | dy_dparam[0] = 1.0; | |

1598 | dy_dparam[1] = exp_conv[1][i]; | |

1599 | dy_dparam[2] = param[1] * exp_conv[2][i]; | |

1600 | dy_dparam[3] = param[1] * exp_conv[3][i]; | |

1601 | ||

1602 | sig2i = (yfit[i] > 15 ? 1.0/yfit[i] : 1.0/15); | |

1603 | ||

1604 | for (j=0, l=0; l<nparam; l++) { | |

1605 | if (paramfree[l]) { | |

1606 | wt = dy_dparam[l] * sig2i; | |

1607 | for (k=0, m=0; m<=l; m++) | |

1608 | if (paramfree[m]) | |

1609 | alpha[j][k++] += wt * dy_dparam[m]; | |

1610 | beta[j] += wt * dy[i]; | |

1611 | j++; | |

1612 | } | |

1613 | } | |

1614 | /* And find chi^2 */ | |

1615 | *chisq += dy[i] * dy[i] * sig2i; | |

1616 | } | |

1617 | break; | |

1618 | ||

1619 | case NOISE_GAUSSIAN_FIT: | |

1620 | *chisq = 0.0; | |

1621 | /* Summation loop over all data */ | |

1622 | for (i=fit_start; i<fit_end; i++) { | |

1623 | yfit[i] = param[0] + param[1] * exp_conv[1][i]; | |

1624 | dy[i] = y[i] - yfit[i]; | |

1625 | ||

1626 | dy_dparam[0] = 1.0; | |

1627 | dy_dparam[1] = exp_conv[1][i]; | |

1628 | dy_dparam[2] = param[1] * exp_conv[2][i]; | |

1629 | dy_dparam[3] = param[1] * exp_conv[3][i]; | |

1630 | ||

1631 | sig2i = (yfit[i] > 1 ? 1.0/yfit[i] : 1.0); | |

1632 | ||

1633 | for (j=0, l=0; l<nparam; l++) { | |

1634 | if (paramfree[l]) { | |

1635 | wt = dy_dparam[l] * sig2i; | |

1636 | for (k=0, m=0; m<=l; m++) | |

1637 | if (paramfree[m]) | |

1638 | alpha[j][k++] += wt * dy_dparam[m]; | |

1639 | beta[j] += wt * dy[i]; | |

1640 | j++; | |

1641 | } | |

1642 | } | |

1643 | /* And find chi^2 */ | |

1644 | *chisq += dy[i] * dy[i] * sig2i; | |

1645 | } | |

1646 | break; | |

1647 | ||

1648 | case NOISE_MLE: | |

1649 | *chisq = 0.0; | |

1650 | /* Summation loop over all data */ | |

1651 | for (i=fit_start; i<fit_end; i++) { | |

1652 | yfit[i] = param[0] + param[1] * exp_conv[1][i]; | |

1653 | dy[i] = y[i] - yfit[i]; | |

1654 | ||

1655 | dy_dparam[0] = 1.0; | |

1656 | dy_dparam[1] = exp_conv[1][i]; | |

1657 | dy_dparam[2] = param[1] * exp_conv[2][i]; | |

1658 | dy_dparam[3] = param[1] * exp_conv[3][i]; | |

1659 | ||

1660 | sig2i = (yfit[i] > 1 ? 1.0/yfit[i] : 1.0); | |

1661 | y_ymod=y[i]/yfit[i]; | |

1662 | ||

1663 | for (j=0, l=0; l<nparam; l++) { | |

1664 | if (paramfree[l]) { | |

1665 | wt = dy_dparam[l] * sig2i; | |

1666 | for (k=0, m=0; m<=l; m++) | |

1667 | if (paramfree[m]) | |

1668 | alpha[j][k++] += wt*dy_dparam[m]*y_ymod; //wt * dy_dparam[m]; | |

1669 | beta[j] += wt * dy[i]; | |

1670 | j++; | |

1671 | } | |

1672 | } | |

1673 | // And find chi^2 | |

1674 | if (yfit[i]<=0.0) | |

1675 | ; // do nothing | |

1676 | else if (y[i]==0.0) | |

1677 | *chisq += 2.0*yfit[i]; // to avoid NaN from log | |

1678 | else | |

1679 | *chisq += 2.0*(yfit[i]-y[i]) - 2.0*y[i]*log(yfit[i]/y[i]); // was dy[i] * dy[i] * sig2i; | |

1680 | } | |

1681 | if (*chisq <= 0.0) *chisq = 1.0e308; // don't let chisq=0 through yfit being all -ve | |

1682 | break; | |

1683 | ||

1684 | default: | |

1685 | return -3; | |

1686 | /* break; */ // (unreachable code) | |

1687 | } | |

1688 | break; | |

1689 | ||

1690 | default: | |

1691 | dbgprintf(1, "compute_exps_fn: please update me!\n"); | |

1692 | return -1; | |

1693 | } | |

1694 | ||

1695 | /* Fill in the symmetric side */ | |

1696 | for (j=1; j<mfit; j++) | |

1697 | for (k=0; k<j; k++) | |

1698 | alpha[k][j] = alpha[j][k]; | |

1699 | ||

1700 | return 0; | |

1701 | } | |

1702 | ||

1703 | ||

1704 | /* And this is a final variant which computes the true chi-squared | |

1705 | values and the full fit, as in EcfSingle.c */ | |

1706 | int GCI_marquardt_global_compute_exps_fn_final( | |

1707 | float xincr, float y[], | |

1708 | int ndata, int fit_start, int fit_end, | |

1709 | noise_type noise, float sig[], int ftype, | |

1710 | float param[], int paramfree[], int nparam, | |

1711 | float *exp_conv[], | |

1712 | float yfit[], float dy[], float *chisq) | |

1713 | { | |

1714 | int i, j, mfit; | |

1715 | float sig2i; | |

1716 | ||

1717 | for (j=0, mfit=0; j<nparam; j++) | |

1718 | if (paramfree[j]) | |

1719 | mfit++; | |

1720 | ||

1721 | /* Calculation of the fitting data will depend upon the type of | |

1722 | noise. Since there's no convolution involved here, this is | |

1723 | very easy. */ | |

1724 | ||

1725 | switch (ftype) { | |

1726 | case FIT_GLOBAL_MULTIEXP: | |

1727 | switch (noise) { | |

1728 | case NOISE_CONST: | |

1729 | *chisq = 0.0; | |

1730 | /* Summation loop over all data */ | |

1731 | for (i=0; i<ndata; i++) { | |

1732 | yfit[i] = param[0]; /* Z */ | |

1733 | for (j=1; j<nparam; j+=2) { | |

1734 | yfit[i] += param[j] * exp_conv[j+1][i]; | |

1735 | /* A_j . exp(-x/tau_j) */ | |

1736 | } | |

1737 | dy[i] = y[i] - yfit[i]; | |

1738 | ||

1739 | /* And find chi^2 */ | |

1740 | if (i >= fit_start && i < fit_end) | |

1741 | *chisq += dy[i] * dy[i]; | |

1742 | } | |

1743 | ||

1744 | /* Now divide by sigma^2 */ | |

1745 | sig2i = 1.0 / (sig[0] * sig[0]); | |

1746 | *chisq *= sig2i; | |

1747 | break; | |

1748 | ||

1749 | case NOISE_GIVEN: /* This is essentially the NR version */ | |

1750 | *chisq = 0.0; | |

1751 | /* Summation loop over all data */ | |

1752 | for (i=0; i<ndata; i++) { | |

1753 | yfit[i] = param[0]; /* Z */ | |

1754 | for (j=1; j<nparam; j+=2) { | |

1755 | yfit[i] += param[j] * exp_conv[j+1][i]; | |

1756 | /* A_j . exp(-x/tau_j) */ | |

1757 | } | |

1758 | dy[i] = y[i] - yfit[i]; | |

1759 | ||

1760 | /* And find chi^2 */ | |

1761 | if (i >= fit_start && i < fit_end) { | |

1762 | sig2i = 1.0 / (sig[i] * sig[i]); | |

1763 | *chisq += dy[i] * dy[i] * sig2i; | |

1764 | } | |

1765 | } | |

1766 | break; | |

1767 | ||

1768 | case NOISE_POISSON_DATA: | |

1769 | *chisq = 0.0; | |

1770 | /* Summation loop over all data */ | |

1771 | for (i=0; i<ndata; i++) { | |

1772 | yfit[i] = param[0]; /* Z */ | |

1773 | for (j=1; j<nparam; j+=2) { | |

1774 | yfit[i] += param[j] * exp_conv[j+1][i]; | |

1775 | /* A_j . exp(-x/tau_j) */ | |

1776 | } | |

1777 | dy[i] = y[i] - yfit[i]; | |

1778 | ||

1779 | /* And find chi^2 */ | |

1780 | if (i >= fit_start && i < fit_end) { | |

1781 | /* we still don't let the variance drop below 1 */ | |

1782 | sig2i = (y[i] > 1 ? 1.0/y[i] : 1.0); | |

1783 | *chisq += dy[i] * dy[i] * sig2i; | |

1784 | } | |

1785 | } | |

1786 | break; | |

1787 | ||

1788 | case NOISE_POISSON_FIT: | |

1789 | *chisq = 0.0; | |

1790 | /* Summation loop over all data */ | |

1791 | for (i=0; i<ndata; i++) { | |

1792 | yfit[i] = param[0]; /* Z */ | |

1793 | for (j=1; j<nparam; j+=2) { | |

1794 | yfit[i] += param[j] * exp_conv[j+1][i]; | |

1795 | /* A_j . exp(-x/tau_j) */ | |

1796 | } | |

1797 | dy[i] = y[i] - yfit[i]; | |

1798 | ||

1799 | /* And find chi^2 */ | |

1800 | if (i >= fit_start && i < fit_end) { | |

1801 | /* we still don't let the variance drop below 1 */ | |

1802 | sig2i = (yfit[i] > 1 ? 1.0/yfit[i] : 1.0); | |

1803 | *chisq += dy[i] * dy[i] * sig2i; | |

1804 | } | |

1805 | } | |

1806 | break; | |

1807 | ||

1808 | case NOISE_MLE: | |

1809 | *chisq = 0.0; | |

1810 | /* Summation loop over all data */ | |

1811 | for (i=0; i<ndata; i++) { | |

1812 | yfit[i] = param[0]; /* Z */ | |

1813 | for (j=1; j<nparam; j+=2) { | |

1814 | yfit[i] += param[j] * exp_conv[j+1][i]; | |

1815 | /* A_j . exp(-x/tau_j) */ | |

1816 | } | |

1817 | dy[i] = y[i] - yfit[i]; | |

1818 | ||

1819 | /* And find chi^2 */ | |

1820 | if (i >= fit_start && i < fit_end) { | |

1821 | // sig2i = (yfit[i] > 1 ? 1.0/yfit[i] : 1.0); | |

1822 | // *chisq += dy[i] * dy[i] * sig2i; | |

1823 | if (yfit[i]<=0.0) | |

1824 | ; // do nothing | |

1825 | else if (y[i]==0.0) | |

1826 | *chisq += 2.0*yfit[i]; // to avoid NaN from log | |

1827 | else | |

1828 | *chisq += 2.0*(yfit[i]-y[i]) - 2.0*y[i]*log(yfit[i]/y[i]); // was dy[i] * dy[i] * sig2i; | |

1829 | } | |

1830 | } | |

1831 | if (*chisq <= 0.0) *chisq = 1.0e308; // don't let chisq=0 through yfit being all -ve | |

1832 | break; | |

1833 | ||

1834 | ||

1835 | case NOISE_GAUSSIAN_FIT: | |

1836 | *chisq = 0.0; | |

1837 | /* Summation loop over all data */ | |

1838 | for (i=0; i<ndata; i++) { | |

1839 | yfit[i] = param[0]; /* Z */ | |

1840 | for (j=1; j<nparam; j+=2) { | |

1841 | yfit[i] += param[j] * exp_conv[j+1][i]; | |

1842 | /* A_j . exp(-x/tau_j) */ | |

1843 | } | |

1844 | dy[i] = y[i] - yfit[i]; | |

1845 | ||

1846 | /* And find chi^2 */ | |

1847 | if (i >= fit_start && i < fit_end) { | |

1848 | sig2i = (yfit[i] > 1 ? 1.0/yfit[i] : 1.0); | |

1849 | *chisq += dy[i] * dy[i] * sig2i; | |

1850 | } | |

1851 | } | |

1852 | break; | |

1853 | ||

1854 | default: | |

1855 | return -3; | |

1856 | /* break; */ // (unreachable code) | |

1857 | } | |

1858 | break; | |

1859 | ||

1860 | case FIT_GLOBAL_STRETCHEDEXP: | |

1861 | switch (noise) { | |

1862 | case NOISE_CONST: | |

1863 | *chisq = 0.0; | |

1864 | /* Summation loop over all data */ | |

1865 | for (i=0; i<ndata; i++) { | |

1866 | yfit[i] = param[0] + param[1] * exp_conv[1][i]; | |

1867 | dy[i] = y[i] - yfit[i]; | |

1868 | ||

1869 | /* And find chi^2 */ | |

1870 | if (i >= fit_start && i < fit_end) | |

1871 | *chisq += dy[i] * dy[i]; | |

1872 | } | |

1873 | ||

1874 | /* Now divide by sigma^2 */ | |

1875 | sig2i = 1.0 / (sig[0] * sig[0]); | |

1876 | *chisq *= sig2i; | |

1877 | break; | |

1878 | ||

1879 | case NOISE_GIVEN: /* This is essentially the NR version */ | |

1880 | *chisq = 0.0; | |

1881 | /* Summation loop over all data */ | |

1882 | for (i=0; i<ndata; i++) { | |

1883 | yfit[i] = param[0] + param[1] * exp_conv[1][i]; | |

1884 | dy[i] = y[i] - yfit[i]; | |

1885 | ||

1886 | /* And find chi^2 */ | |

1887 | if (i >= fit_start && i < fit_end) { | |

1888 | sig2i = 1.0 / (sig[i] * sig[i]); | |

1889 | *chisq += dy[i] * dy[i] * sig2i; | |

1890 | } | |

1891 | } | |

1892 | break; | |

1893 | ||

1894 | case NOISE_POISSON_DATA: | |

1895 | *chisq = 0.0; | |

1896 | /* Summation loop over all data */ | |

1897 | for (i=0; i<ndata; i++) { | |

1898 | yfit[i] = param[0] + param[1] * exp_conv[1][i]; | |

1899 | dy[i] = y[i] - yfit[i]; | |

1900 | ||

1901 | /* And find chi^2 */ | |

1902 | if (i >= fit_start && i < fit_end) { | |

1903 | /* we still don't let the variance drop below 1 */ | |

1904 | sig2i = (y[i] > 1 ? 1.0/y[i] : 1.0); | |

1905 | *chisq += dy[i] * dy[i] * sig2i; | |

1906 | } | |

1907 | } | |

1908 | break; | |

1909 | ||

1910 | case NOISE_POISSON_FIT: | |

1911 | *chisq = 0.0; | |

1912 | /* Summation loop over all data */ | |

1913 | for (i=0; i<ndata; i++) { | |

1914 | yfit[i] = param[0] + param[1] * exp_conv[1][i]; | |

1915 | dy[i] = y[i] - yfit[i]; | |

1916 | ||

1917 | /* And find chi^2 */ | |

1918 | if (i >= fit_start && i < fit_end) { | |

1919 | /* we still don't let the variance drop below 1 */ | |

1920 | sig2i = (yfit[i] > 1 ? 1.0/yfit[i] : 1.0); | |

1921 | *chisq += dy[i] * dy[i] * sig2i; | |

1922 | } | |

1923 | } | |

1924 | break; | |

1925 | ||

1926 | case NOISE_MLE: | |

1927 | *chisq = 0.0; | |

1928 | /* Summation loop over all data */ | |

1929 | for (i=0; i<ndata; i++) { | |

1930 | yfit[i] = param[0] + param[1] * exp_conv[1][i]; | |

1931 | dy[i] = y[i] - yfit[i]; | |

1932 | ||

1933 | /* And find chi^2 */ | |

1934 | if (i >= fit_start && i < fit_end) { | |

1935 | // sig2i = (yfit[i] > 1 ? 1.0/yfit[i] : 1.0); | |

1936 | // *chisq += dy[i] * dy[i] * sig2i; | |

1937 | if (yfit[i]<=0.0) | |

1938 | ; // do nothing | |

1939 | else if (y[i]==0.0) | |

1940 | *chisq += 2.0*yfit[i]; // to avoid NaN from log | |

1941 | else | |

1942 | *chisq += 2.0*(yfit[i]-y[i]) - 2.0*y[i]*log(yfit[i]/y[i]); // was dy[i] * dy[i] * sig2i; | |

1943 | } | |

1944 | } | |

1945 | if (*chisq <= 0.0) *chisq = 1.0e308; // don't let chisq=0 through yfit being all -ve | |

1946 | break; | |

1947 | ||

1948 | case NOISE_GAUSSIAN_FIT: | |

1949 | *chisq = 0.0; | |

1950 | /* Summation loop over all data */ | |

1951 | for (i=0; i<ndata; i++) { | |

1952 | yfit[i] = param[0] + param[1] * exp_conv[1][i]; | |

1953 | dy[i] = y[i] - yfit[i]; | |

1954 | ||

1955 | /* And find chi^2 */ | |

1956 | if (i >= fit_start && i < fit_end) { | |

1957 | /* we still don't let the variance drop below 1 */ | |

1958 | sig2i = (yfit[i] > 1 ? 1.0/yfit[i] : 1.0); | |

1959 | *chisq += dy[i] * dy[i] * sig2i; | |

1960 | } | |

1961 | } | |

1962 | break; | |

1963 | ||

1964 | default: | |

1965 | return -3; | |

1966 | /* break; */ // (unreachable code) | |

1967 | } | |

1968 | break; | |

1969 | ||

1970 | default: | |

1971 | dbgprintf(1, "compute_exps_fn: please update me!\n"); | |

1972 | return -1; | |

1973 | } | |

1974 | ||

1975 | return 0; | |

1976 | } | |

1977 | ||

1978 | ||

1979 | /* This one does the actual global fitting for multiexponential taus. | |

1980 | It is basically similar to the above do_fit_single function, except | |

1981 | that now we handle the extra intricacies involved in global | |

1982 | fitting, in particular, the much larger alpha matrix is handled in | |

1983 | a special way. */ | |

1984 | ||

1985 | int GCI_marquardt_global_exps_do_fit_instr( | |

1986 | float xincr, float **trans, int ndata, int ntrans, | |

1987 | int fit_start, int fit_end, float instr[], int ninstr, | |

1988 | noise_type noise, float sig[], int ftype, | |

1989 | float **param, int paramfree[], int nparam, restrain_type restrain, | |

1990 | float chisq_delta, float exp_pure[], float *exp_conv[], | |

1991 | float **fitted, float **residuals, | |

1992 | float **covar_scratch, float **alpha_scratch, | |

1993 | float *chisq_trans, float *chisq_global, | |

1994 | int drop_bad_transients) | |

1995 | { | |

1996 | float alambda, ochisq_global, *ochisq_trans; | |

1997 | int i, k, itst, itst_max; | |

1998 | int ret; | |

1999 | ||

2000 | itst_max = (restrain == ECF_RESTRAIN_DEFAULT) ? 4 : 6; | |

2001 | ||

2002 | /* If there are no global parameters being fitted, we simply fit | |

2003 | each local set. */ | |

2004 | switch (ftype) { | |

2005 | case FIT_GLOBAL_MULTIEXP: | |

2006 | for (i=2; i<nparam; i+=2) { | |

2007 | if (paramfree[i]) { | |

2008 | i = 0; /* sentinel value */ | |

2009 | break; | |

2010 | } | |

2011 | } | |

2012 | break; | |

2013 | ||

2014 | case FIT_GLOBAL_STRETCHEDEXP: | |

2015 | for (i=2; i<nparam; i++) { | |

2016 | if (paramfree[i]) { | |

2017 | i = 0; /* sentinel value */ | |

2018 | break; | |

2019 | } | |

2020 | } | |

2021 | break; | |

2022 | ||

2023 | default: | |

2024 | dbgprintf(1, "exps_do_fit_instr: please update me!\n"); | |

2025 | return -1; | |

2026 | } | |

2027 | ||

2028 | if (i > 0) { /* no globals to fit */ | |

2029 | if (GCI_marquardt_global_exps_calculate_exps_instr( | |

2030 | xincr, ndata, instr, ninstr, ftype, param[0], nparam, | |

2031 | exp_pure, exp_conv) != 0) | |

2032 | return -2; | |

2033 | ||

2034 | *chisq_global = 0; | |

2035 | ||

2036 | for (i=0; i<ntrans; i++) { | |

2037 | if (drop_bad_transients && chisq_trans[i] < 0) | |

2038 | continue; | |

2039 | ||

2040 | ret = GCI_marquardt_global_exps_do_fit_single( | |

2041 | xincr, trans[i], ndata, fit_start, fit_end, | |

2042 | noise, sig, ftype, param[i], paramfree, nparam, restrain, | |

2043 | //// exp_conv, fitted[i], residuals[i], | |

2044 | chisq_delta, exp_conv, fitted[0], residuals[0], | |

2045 | covar_scratch, alpha_scratch, &chisq_trans[i]); | |

2046 | if (ret < 0) { | |

2047 | if (drop_bad_transients) { | |

2048 | dbgprintf(1, "In do_fit_instr, do_fit_single returned %d " | |

2049 | "for transient %d, dropping it\n", ret, i); | |

2050 | chisq_trans[i] = -1; | |

2051 | } else { | |

2052 | dbgprintf(1, "In do_fit_instr, do_fit_single returned %d " | |

2053 | "for transient %d\n", ret, i); | |

2054 | return -10 + ret; | |

2055 | } | |

2056 | } else { | |

2057 | *chisq_global += chisq_trans[i]; | |

2058 | } | |

2059 | } | |

2060 | return 0; | |

2061 | } | |

2062 | ||

2063 | /* If there are no free local variables to fit, we still do the | |

2064 | global fitting, but we have to be a little careful in some of | |

2065 | the later routines */ | |

2066 | ||

2067 | /* Now allocate all of the arrays we will need. */ | |

2068 | ||

2069 | if ((ochisq_trans = (float *) malloc(ntrans * sizeof(float))) == NULL) | |

2070 | return -1; | |

2071 | ||

2072 | /* We now begin our standard Marquardt loop, with several | |

2073 | modifications */ | |

2074 | alambda = -1; | |

2075 | ret = GCI_marquardt_global_exps_global_step( | |

2076 | xincr, trans, ndata, ntrans, fit_start, fit_end, | |

2077 | instr, ninstr, noise, sig, ftype, | |

2078 | param, paramfree, nparam, restrain, exp_pure, exp_conv, | |

2079 | fitted, residuals, chisq_trans, chisq_global, | |

2080 | alpha_scratch, &alambda, drop_bad_transients); | |

2081 | if (ret != 0) { | |

2082 | dbgprintf(1, "In do_fit_instr, first global_step returned %d\n", ret); | |

2083 | if (ret != -1) { | |

2084 | /* Wasn't a memory error, so unallocate arrays */ | |

2085 | alambda = 0.0; | |

2086 | GCI_marquardt_global_exps_global_step( | |

2087 | xincr, trans, ndata, ntrans, fit_start, fit_end, | |

2088 | instr, ninstr, noise, sig, ftype, | |

2089 | param, paramfree, nparam, restrain, exp_pure, exp_conv, | |

2090 | fitted, residuals, chisq_trans, chisq_global, | |

2091 | alpha_scratch, &alambda, drop_bad_transients); | |

2092 | } | |

2093 | free(ochisq_trans); | |

2094 | return ret; | |

2095 | } | |

2096 | ||

2097 | k = 1; /* Iteration counter */ | |

2098 | itst = 0; | |

2099 | for (;;) { | |

2100 | dbgprintf(3, "In do_fit_instr, beginning iteration %d:\n", k); | |

2101 | dbgprintf(3, " itst = %d, chisq_global = %.4f\n", itst, *chisq_global); | |

2102 | ||

2103 | k++; | |

2104 | if (k > MAXITERS) { | |

2105 | return -2; | |

2106 | } | |

2107 | ||

2108 | ochisq_global = *chisq_global; | |

2109 | for (i=0; i<ntrans; i++) | |

2110 | ochisq_trans[i] = chisq_trans[i]; | |

2111 | ||

2112 | ret = GCI_marquardt_global_exps_global_step( | |

2113 | xincr, trans, ndata, ntrans, fit_start, fit_end, | |

2114 | instr, ninstr, noise, sig, ftype, | |

2115 | param, paramfree, nparam, restrain, exp_pure, exp_conv, | |

2116 | fitted, residuals, chisq_trans, chisq_global, | |

2117 | alpha_scratch, &alambda, drop_bad_transients); | |

2118 | if (ret != 0) { | |

2119 | dbgprintf(1, "In do_fit_instr, second global_step returned %d\n", | |

2120 | ret); | |

2121 | /* Unallocate arrays */ | |

2122 | alambda = 0.0; | |

2123 | GCI_marquardt_global_exps_global_step( | |

2124 | xincr, trans, ndata, ntrans, fit_start, fit_end, | |

2125 | instr, ninstr, noise, sig, ftype, | |

2126 | param, paramfree, nparam, restrain, exp_pure, exp_conv, | |

2127 | fitted, residuals, chisq_trans, chisq_global, | |

2128 | alpha_scratch, &alambda, drop_bad_transients); | |

2129 | free(ochisq_trans); | |

2130 | return ret; | |

2131 | } | |

2132 | ||

2133 | if (*chisq_global > ochisq_global) | |

2134 | itst = 0; | |

2135 | else { | |

2136 | /* Let's try this approach; I really don't know what will | |

2137 | be best */ | |

2138 | float maxdiff; | |

2139 | ||

2140 | maxdiff = 0.0; | |

2141 | for (i=0; i<ntrans; i++) | |

2142 | if (ochisq_trans[i] - chisq_trans[i] > maxdiff) | |

2143 | maxdiff = ochisq_trans[i] - chisq_trans[i]; | |

2144 | ||

2145 | if (maxdiff < chisq_delta) | |

2146 | itst++; | |

2147 | dbgprintf(3, "In do_fit_instr, maxdiff = %.3f:\n", maxdiff); | |

2148 | } | |

2149 | ||

2150 | if (itst < itst_max) continue; | |

2151 | ||

2152 | /* Endgame */ | |

2153 | alambda = 0.0; | |

2154 | ret = GCI_marquardt_global_exps_global_step( | |

2155 | xincr, trans, ndata, ntrans, fit_start, fit_end, | |

2156 | instr, ninstr, noise, sig, ftype, | |

2157 | param, paramfree, nparam, restrain, exp_pure, exp_conv, | |

2158 | fitted, residuals, chisq_trans, chisq_global, | |

2159 | alpha_scratch, &alambda, drop_bad_transients); | |

2160 | if (ret != 0) { | |

2161 | dbgprintf(1, "In do_fit_instr, final global_step returned %d\n", | |

2162 | ret); | |

2163 | free(ochisq_trans); | |

2164 | return ret; | |

2165 | } | |

2166 | ||

2167 | free(ochisq_trans); | |

2168 | return k; /* We're done now */ | |

2169 | } | |

2170 | } | |

2171 | ||

2172 | ||

2173 | /* And this one is basically a specialised GCI_marquardt_instr_step | |

2174 | for the global fitting setup. */ | |

2175 | int GCI_marquardt_global_exps_global_step( | |

2176 | float xincr, float **trans, | |

2177 | int ndata, int ntrans, int fit_start, int fit_end, | |

2178 | float instr[], int ninstr, | |

2179 | noise_type noise, float sig[], int ftype, | |

2180 | float **param, int paramfree[], int nparam, | |

2181 | restrain_type restrain, | |

2182 | float exp_pure[], float *exp_conv[], | |

2183 | float **yfit, float **dy, | |

2184 | float *chisq_trans, float *chisq_global, | |

2185 | float **alpha_scratch, float *alambda, | |

2186 | int drop_bad_transients) | |

2187 | { | |

2188 | int i, j, ret; | |

2189 | static global_matrix alpha, covar; | |

2190 | static global_vector beta, dparam; | |

2191 | static float **paramtry; | |

2192 | static int mfit_local, mfit_global; | |

2193 | static int gindex[MAXFIT], lindex[MAXFIT]; | |

2194 | static float ochisq_global, *ochisq_trans; | |

2195 | static void (*fitfunc)(float, float [], float *, float [], int); | |

2196 | static int initialised=0; | |

2197 | ||

2198 | if (nparam > MAXFIT) | |

2199 | return -10; | |

2200 | if (xincr <= 0) | |

2201 | return -11; | |

2202 | if (fit_start < 0 || fit_start > fit_end || fit_end > ndata) | |

2203 | return -12; | |

2204 | ||

2205 | /* Initialisation */ | |

2206 | /* We assume we're given sensible starting values for param[] */ | |

2207 | if (*alambda < 0.0) { | |

2208 | /* Start by allocating lots of variables we will need */ | |

2209 | mfit_local = mfit_global = 0; | |

2210 | ||

2211 | switch (ftype) { | |

2212 | case FIT_GLOBAL_MULTIEXP: | |

2213 | fitfunc = GCI_multiexp_tau; | |

2214 | ||

2215 | /* We know that all of param[2i], the taus, are the global | |

2216 | variables, and that the param[2i+1], the A's, are the | |

2217 | local variables, along with param[0] = Z. We stored | |

2218 | the indices of the local and global free variables in | |

2219 | lindex and gindex respectively. */ | |

2220 | if (paramfree[0]) { | |

2221 | lindex[mfit_local++] = 0; | |

2222 | } | |

2223 | for (i=1; i<nparam; i+=2) { | |

2224 | if (paramfree[i]) | |

2225 | lindex[mfit_local++] = i; | |

2226 | } | |

2227 | for (i=2; i<nparam; i+=2) { | |

2228 | if (paramfree[i]) | |

2229 | gindex[mfit_global++] = i; | |

2230 | } | |

2231 | break; | |

2232 | ||

2233 | case FIT_GLOBAL_STRETCHEDEXP: | |

2234 | fitfunc = GCI_stretchedexp; | |

2235 | ||

2236 | if (paramfree[0]) | |

2237 | lindex[mfit_local++] = 0; /* Z */ | |

2238 | if (paramfree[1]) | |

2239 | lindex[mfit_local++] = 1; /* A */ | |

2240 | if (paramfree[2]) | |

2241 | gindex[mfit_global++] = 2; /* tau */ | |

2242 | if (paramfree[3]) | |

2243 | gindex[mfit_global++] = 3; /* h */ | |

2244 | break; | |

2245 | ||

2246 | default: | |

2247 | dbgprintf(1, "exps_global_step: please update me!\n"); | |

2248 | return -1; | |

2249 | } | |

2250 | ||

2251 | if (initialised) { | |

2252 | GCI_free_global_matrix(&alpha); GCI_free_global_matrix(&covar); | |

2253 | GCI_ecf_free_matrix(paramtry); GCI_free_global_vector(&beta); | |

2254 | GCI_free_global_vector(&dparam); free(ochisq_trans); | |

2255 | initialised = 0; | |

2256 | } | |

2257 | ||

2258 | if (GCI_alloc_global_matrix(&alpha, mfit_global, mfit_local, ntrans) | |

2259 | != 0) | |

2260 | return -1; | |

2261 | ||

2262 | if (GCI_alloc_global_matrix(&covar, mfit_global, mfit_local, ntrans) | |

2263 | != 0) { | |

2264 | GCI_free_global_matrix(&alpha); | |

2265 | return -1; | |

2266 | } | |

2267 | ||

2268 | if ((paramtry = GCI_ecf_matrix(ntrans, nparam)) == NULL) { | |

2269 | GCI_free_global_matrix(&alpha); GCI_free_global_matrix(&covar); | |

2270 | return -1; | |

2271 | } | |

2272 | ||

2273 | if (GCI_alloc_global_vector(&beta, mfit_global, mfit_local, ntrans) | |

2274 | != 0) { | |

2275 | GCI_free_global_matrix(&alpha); GCI_free_global_matrix(&covar); | |

2276 | GCI_ecf_free_matrix(paramtry); | |

2277 | return -1; | |

2278 | } | |

2279 | ||

2280 | if (GCI_alloc_global_vector(&dparam, mfit_global, mfit_local, ntrans) | |

2281 | != 0) { | |

2282 | GCI_free_global_matrix(&alpha); GCI_free_global_matrix(&covar); | |

2283 | GCI_ecf_free_matrix(paramtry); GCI_free_global_vector(&beta); | |

2284 | return -1; | |

2285 | } | |

2286 | ||

2287 | if ((ochisq_trans = (float *) malloc(ntrans * sizeof(float))) | |

2288 | == NULL) { | |

2289 | GCI_free_global_matrix(&alpha); GCI_free_global_matrix(&covar); | |

2290 | GCI_ecf_free_matrix(paramtry); GCI_free_global_vector(&beta); | |

2291 | GCI_free_global_vector(&dparam); | |

2292 | return -1; | |

2293 | } | |

2294 | ||

2295 | initialised = 1; | |

2296 | ||

2297 | if (GCI_marquardt_global_compute_global_exps_fn( | |

2298 | xincr, trans, ndata, ntrans, | |

2299 | fit_start, fit_end, instr, ninstr, noise, sig, ftype, | |

2300 | param, paramfree, nparam, mfit_global, mfit_local, | |

2301 | gindex, lindex, exp_pure, exp_conv, | |

2302 | yfit, dy, alpha, beta, alpha_scratch, | |

2303 | chisq_trans, chisq_global, drop_bad_transients) != 0) | |

2304 | return -2; | |

2305 | ||

2306 | *alambda = 0.001; | |

2307 | ochisq_global = *chisq_global; | |

2308 | for (i=0; i<ntrans; i++) | |

2309 | ochisq_trans[i] = chisq_trans[i]; | |

2310 | ||

2311 | /* Initialise paramtry to param */ | |

2312 | for (i=0; i<ntrans; i++) { | |

2313 | for (j=0; j<nparam; j++) | |

2314 | paramtry[i][j] = param[i][j]; | |

2315 | } | |

2316 | } | |

2317 | ||

2318 | /* Once converged, evaluate covariance matrix */ | |

2319 | if (*alambda == 0) { | |

2320 | if (GCI_marquardt_global_compute_global_exps_fn_final( | |

2321 | xincr, trans, ndata, ntrans, | |

2322 | fit_start, fit_end, instr, ninstr, noise, sig, ftype, | |

2323 | param, paramfree, nparam, mfit_global, mfit_local, | |

2324 | gindex, lindex, exp_pure, exp_conv, yfit, dy, | |

2325 | chisq_trans, chisq_global, drop_bad_transients) != 0) | |

2326 | return -3; | |

2327 | /* Don't need to do this here; if we wished to, we'd have to | |

2328 | move this code (the "if (*alambda == 0)" block) to after | |

2329 | the Gauss-Jordan call. We'd also need to rewrite it for | |

2330 | our situation.... */ | |

2331 | // if (mfit < nparam) { /* no need to do this otherwise */ | |

2332 | // GCI_covar_sort(covar, nparam, paramfree, mfit); | |

2333 | // GCI_covar_sort(alpha, nparam, paramfree, mfit); | |

2334 | // } | |

2335 | GCI_free_global_matrix(&alpha); GCI_free_global_matrix(&covar); | |

2336 | GCI_ecf_free_matrix(paramtry); GCI_free_global_vector(&beta); | |

2337 | GCI_free_global_vector(&dparam); free(ochisq_trans); | |

2338 | initialised = 0; | |

2339 | return 0; | |

2340 | } | |

2341 | ||

2342 | /* Alter linearised fitting matrix by augmenting diagonal | |

2343 | elements. */ | |

2344 | GCI_copy_global_matrix(covar, alpha, mfit_global, mfit_local, ntrans); | |

2345 | GCI_copy_global_vector(dparam, beta, mfit_global, mfit_local, ntrans); | |

2346 | for (j=0; j<mfit_global; j++) | |

2347 | covar.P[j][j] *= 1.0 + (*alambda); | |

2348 | for (i=0; i<ntrans; i++) | |

2349 | for (j=0; j<mfit_local; j++) | |

2350 | covar.S[i][j][j] *= 1.0 + (*alambda); | |

2351 | ||

2352 | /* Matrix solution; GCI_gauss_jordan solves Ax=b rather than AX=B */ | |

2353 | if (GCI_marquardt_global_solve_eqn(covar, dparam, | |

2354 | mfit_global, mfit_local, ntrans) != 0) | |

2355 | return -3; | |

2356 | ||

2357 | /* Did the trial succeed? Modify param by dparam... */ | |

2358 | for (i=0; i<ntrans; i++) { | |

2359 | for (j=0; j<mfit_global; j++) | |

2360 | paramtry[i][gindex[j]] = param[i][gindex[j]] + dparam.global[j]; | |

2361 | for (j=0; j<mfit_local; j++) | |

2362 | paramtry[i][lindex[j]] = | |

2363 | param[i][lindex[j]] + dparam.local[i*mfit_local + j]; | |

2364 | } | |

2365 | ||

2366 | for (i=0; i<ntrans; i++) { | |

2367 | if (drop_bad_transients && chisq_trans[i] < 0) | |

2368 | continue; | |

2369 | ||

2370 | if (restrain == ECF_RESTRAIN_DEFAULT) | |

2371 | ret = check_ecf_params (paramtry[i], nparam, fitfunc); | |

2372 | else | |

2373 | ret = check_ecf_user_params (paramtry[i], nparam, fitfunc); | |

2374 | ||

2375 | if (ret != 0) { | |

2376 | /* Bad parameters, increase alambda and return */ | |

2377 | *alambda *= 10.0; | |

2378 | return 0; | |

2379 | } | |

2380 | } | |

2381 | ||

2382 | if (GCI_marquardt_global_compute_global_exps_fn( | |

2383 | xincr, trans, ndata, ntrans, | |

2384 | fit_start, fit_end, instr, ninstr, noise, sig, ftype, | |

2385 | paramtry, paramfree, nparam, mfit_global, mfit_local, | |

2386 | gindex, lindex, exp_pure, exp_conv, | |

2387 | yfit, dy, covar, dparam, alpha_scratch, | |

2388 | chisq_trans, chisq_global, drop_bad_transients) != 0) | |

2389 | return -2; | |

2390 | ||

2391 | /* Success, accept the new solution */ | |

2392 | if (*chisq_global < ochisq_global) { | |

2393 | *alambda *= 0.1; | |

2394 | ochisq_global = *chisq_global; | |

2395 | for (i=0; i<ntrans; i++) | |

2396 | ochisq_trans[i] = chisq_trans[i]; | |

2397 | GCI_copy_global_matrix(alpha, covar, mfit_global, mfit_local, ntrans); | |

2398 | GCI_copy_global_vector(beta, dparam, mfit_global, mfit_local, ntrans); | |

2399 | for (i=0; i<ntrans; i++) { | |

2400 | for (j=0; j<nparam; j++) | |

2401 | param[i][j] = paramtry[i][j]; | |

2402 | } | |

2403 | } else { /* Failure, increase alambda and return */ | |

2404 | *alambda *= 10.0; | |

2405 | *chisq_global = ochisq_global; | |

2406 | for (i=0; i<ntrans; i++) | |

2407 | chisq_trans[i] = ochisq_trans[i]; | |

2408 | } | |

2409 | ||

2410 | return 0; | |

2411 | } | |

2412 | ||

2413 | ||

2414 | /* Here we use alpha only for scratch space */ | |

2415 | int GCI_marquardt_global_compute_global_exps_fn( | |

2416 | float xincr, float **trans, int ndata, int ntrans, | |

2417 | int fit_start, int fit_end, float instr[], int ninstr, | |

2418 | noise_type noise, float sig[], int ftype, | |

2419 | float **param, int paramfree[], int nparam, | |

2420 | int mfit_global, int mfit_local, int gindex[], int lindex[], | |

2421 | float exp_pure[], float *exp_conv[], | |

2422 | float **yfit, float **dy, global_matrix alpha, global_vector beta, | |

2423 | float **alpha_scratch, float *chisq_trans, float *chisq_global, | |

2424 | int drop_bad_transients) | |

2425 | { | |

2426 | int i, j, k, ret; | |

2427 | float beta_scratch[MAXFIT]; /* scratch space */ | |

2428 | ||

2429 | /* Calculate the exponential array once only */ | |

2430 | if (GCI_marquardt_global_exps_calculate_exps_instr( | |

2431 | xincr, ndata, instr, ninstr, ftype, param[0], nparam, | |

2432 | exp_pure, exp_conv) != 0) | |

2433 | return -1; | |

2434 | ||

2435 | /* We initialise P and beta_global to zero; the others don't | |

2436 | matter, as they will be totally overwritten */ | |

2437 | for (i=0; i<mfit_global; i++) { | |

2438 | for (j=0; j<mfit_global; j++) | |

2439 | alpha.P[i][j] = 0; | |

2440 | beta.global[i] = 0; | |

2441 | } | |

2442 | *chisq_global = 0.0; | |

2443 | ||

2444 | for (i=0; i<ntrans; i++) { | |

2445 | if (drop_bad_transients && chisq_trans[i] < 0) { | |

2446 | for (j=0; j<mfit_global; j++) | |

2447 | for (k=0; k<mfit_local; k++) | |

2448 | alpha.Q[j][i*mfit_local + k] = 0.0; | |

2449 | for (j=0; j<mfit_local; j++) { | |

2450 | /* Make this component of S an identity matrix and of | |

2451 | beta zero */ | |

2452 | for (k=0; k<mfit_local; k++) | |

2453 | alpha.S[i][j][k] = (j == k) ? 1.0 : 0.0; | |

2454 | beta.local[i*mfit_local + j] = 0; | |

2455 | } | |

2456 | continue; | |

2457 | } | |

2458 | ||

2459 | /* This transient is fine! */ | |

2460 | ret = GCI_marquardt_global_compute_exps_fn( | |

2461 | xincr, trans[i], ndata, fit_start, fit_end, noise, sig, | |

2462 | ftype, param[i], paramfree, nparam, | |

2463 | //// exp_conv, yfit[i], dy[i], | |

2464 | exp_conv, yfit[0], dy[0], | |

2465 | alpha_scratch, beta_scratch, &chisq_trans[i]); | |

2466 | ||

2467 | if (ret != 0) { | |

2468 | if (drop_bad_transients) { | |

2469 | dbgprintf(3, "In compute_global_exps_fn, " | |

2470 | "compute_exps_fn returned %d for transient %d\n", | |

2471 | ret, i); | |

2472 | chisq_trans[i] = -1; | |

2473 | continue; | |

2474 | } else { | |

2475 | dbgprintf(1, "In compute_global_exps_fn, " | |

2476 | "compute_exps_fn returned %d for transient %d\n", | |

2477 | ret, i); | |

2478 | return -2; | |

2479 | } | |

2480 | } | |

2481 | ||

2482 | /* So now have to populate alpha and beta with the contents of | |

2483 | alpha_scratch and beta_scratch. */ | |

2484 | ||

2485 | for (j=0; j<mfit_global; j++) { | |

2486 | for (k=0; k<mfit_global; k++) | |

2487 | alpha.P[j][k] += alpha_scratch[gindex[j]][gindex[k]]; | |

2488 | for (k=0; k<mfit_local; k++) | |

2489 | alpha.Q[j][i*mfit_local + k] = | |

2490 | alpha_scratch[gindex[j]][lindex[k]]; | |

2491 | beta.global[j] += beta_scratch[gindex[j]]; | |

2492 | } | |

2493 | for (j=0; j<mfit_local; j++) { | |

2494 | for (k=0; k<mfit_local; k++) | |

2495 | alpha.S[i][j][k] = alpha_scratch[lindex[j]][lindex[k]]; | |

2496 | beta.local[i*mfit_local + j] = beta_scratch[lindex[j]]; | |

2497 | } | |

2498 | ||

2499 | *chisq_global += chisq_trans[i]; | |

2500 | } | |

2501 | ||

2502 | return 0; | |

2503 | } | |

2504 | ||

2505 | ||

2506 | /* The final variant */ | |

2507 | int GCI_marquardt_global_compute_global_exps_fn_final( | |

2508 | float xincr, float **trans, int ndata, int ntrans, | |

2509 | int fit_start, int fit_end, float instr[], int ninstr, | |

2510 | noise_type noise, float sig[], int ftype, | |

2511 | float **param, int paramfree[], int nparam, | |

2512 | int mfit_global, int mfit_local, int gindex[], int lindex[], | |

2513 | float exp_pure[], float *exp_conv[], | |

2514 | float **yfit, float **dy, | |

2515 | float *chisq_trans, float *chisq_global, int drop_bad_transients) | |

2516 | { | |

2517 | int i, ret; | |

2518 | ||

2519 | /* Calculate the exponential array once only */ | |

2520 | if (GCI_marquardt_global_exps_calculate_exps_instr( | |

2521 | xincr, ndata, instr, ninstr, ftype, param[0], nparam, | |

2522 | exp_pure, exp_conv) != 0) | |

2523 | return -1; | |

2524 | ||

2525 | *chisq_global = 0.0; | |

2526 | ||

2527 | for (i=0; i<ntrans; i++) { | |

2528 | if (drop_bad_transients && chisq_trans[i] < 0) | |

2529 | continue; | |

2530 | ||

2531 | /* This transient is fine! */ | |

2532 | ret = GCI_marquardt_global_compute_exps_fn_final( | |

2533 | xincr, trans[i], ndata, fit_start, fit_end, noise, sig, | |

2534 | ftype, param[i], paramfree, nparam, | |

2535 | //// exp_conv, yfit[i], dy[i], &chisq_trans[i]); | |

2536 | exp_conv, yfit[0], dy[0], &chisq_trans[i]); | |

2537 | ||

2538 | if (ret != 0) { | |

2539 | if (drop_bad_transients) { | |

2540 | dbgprintf(3, "In compute_global_exps_fn_final, " | |

2541 | "compute_exps_fn_final returned %d " | |

2542 | "for transient %d\n", | |

2543 | ret, i); | |

2544 | chisq_trans[i] = -1; | |

2545 | continue; | |

2546 | } else { | |

2547 | dbgprintf(1, "In compute_global_exps_fn_final, " | |

2548 | "compute_exps_fn_final returned %d " | |

2549 | "for transient %d\n", | |

2550 | ret, i); | |

2551 | return -2; | |

2552 | } | |

2553 | } | |

2554 | ||

2555 | *chisq_global += chisq_trans[i]; | |

2556 | } | |

2557 | ||

2558 | return 0; | |

2559 | } | |

2560 | ||

2561 | ||

2562 | /* This function solves the equation Ax=b, where A is the alpha | |

2563 | matrix, which has the form: | |

2564 | ||

2565 | A = (P Q) | |

2566 | (R S) | |

2567 | ||

2568 | Here P is an mfit_global x mfit_global square matrix, S is a block | |

2569 | diagonal matrix with ntrans blocks, each of size mfit_local x | |

2570 | mfit_local, and Q and R are the right sizes to make the whole | |

2571 | matrix square. We solve it by inverting the matrix A using the | |

2572 | formulae given in Numerical Recipes section 2.7, then multiplying | |

2573 | the inverse by b to get x. We are not too concerned about | |

2574 | accuracy, as this does not affect the solution found, only the | |

2575 | route taken to find it. | |

2576 | ||

2577 | Numerical Recipes, section 2.7, notes that A^{-1} is given by: | |

2578 | ||

2579 | (P' Q') | |

2580 | (R' S') | |

2581 | ||

2582 | with: | |

2583 | ||

2584 | P' = (P - Q.S^{-1}.R)^{-1} | |

2585 | Q' = - (P - Q.S^{-1}.R)^{-1} . (Q.S^{-1}) | |

2586 | R' = - (S^{-1}.R) . (P - Q.S^{-1}.R)^{-1} | |

2587 | S' = S^{-1} + (S^{-1}.R) . (P - Q.S^{-1}.R)^{-1} . (Q.S^{-1}) | |

2588 | ||

2589 | We also make use of the fact that A is symmetric, so in particular, | |

2590 | (S^{-1}.R) = (Q.S^{-1})^T and R' = Q'^T. | |

2591 | ||

2592 | We are given A as a global_matrix and b as a global_vector. This | |

2593 | function destroys the original matrix and returns the solution in | |

2594 | place of b. | |

2595 | */ | |

2596 | ||

2597 | int GCI_marquardt_global_solve_eqn(global_matrix A, global_vector b, | |

2598 | int mfit_global, int mfit_local, int ntrans) | |

2599 | { | |

2600 | int row, col, block, i, j; | |

2601 | float x_temp[MAXFIT], x_temp2[MAXFIT]; | |

2602 | static float **QS; | |

2603 | static global_vector x; | |

2604 | static int saved_global=0, saved_local=0, saved_ntrans=0; | |

2605 | ||

2606 | /* If no local parameters, just do a straight Gauss-Jordan method */ | |

2607 | if (mfit_local == 0) { | |

2608 | if (GCI_gauss_jordan(A.P, mfit_global, b.global) != 0) | |

2609 | return -2; | |

2610 | return 0; | |

2611 | } | |

2612 | ||

2613 | /* Allocate arrays if necessary */ | |

2614 | if ((saved_global != mfit_global) || (saved_local != mfit_local) || | |

2615 | (saved_ntrans != ntrans)) { | |

2616 | if (saved_global > 0) { | |

2617 | GCI_ecf_free_matrix(QS); | |

2618 | GCI_free_global_vector(&x); | |

2619 | saved_global = 0; | |

2620 | } | |

2621 | if ((QS = GCI_ecf_matrix(mfit_global, mfit_local*ntrans)) == NULL) | |

2622 | return -1; | |

2623 | if (GCI_alloc_global_vector(&x, mfit_global, mfit_local, ntrans) | |

2624 | != 0) { | |

2625 | GCI_ecf_free_matrix(QS); | |

2626 | return -1; | |

2627 | } | |

2628 | saved_global = mfit_global; | |

2629 | saved_local = mfit_local; | |

2630 | saved_ntrans = ntrans; | |

2631 | } | |

2632 | ||

2633 | /* Start by inverting S */ | |

2634 | for (block=0; block<ntrans; block++) | |

2635 | if (GCI_gauss_jordan(A.S[block], mfit_local, NULL) != 0) | |

2636 | return -2; | |

2637 | ||

2638 | /* Calculate Q.S^{-1} */ | |

2639 | for (row=0; row<mfit_global; row++) | |

2640 | for (block=0; block<ntrans; block++) | |

2641 | for (i=0; i<mfit_local; i++) { | |

2642 | QS[row][block*mfit_local + i] = 0; | |

2643 | for (j=0; j<mfit_local; j++) | |

2644 | QS[row][block*mfit_local + i] += | |

2645 | A.Q[row][block*mfit_local + j] * A.S[block][j][i]; | |

2646 | } | |

2647 | ||

2648 | /* Now find P - Q.S^{-1}.R */ | |

2649 | for (row=0; row<mfit_global; row++) | |

2650 | for (col=0; col<mfit_global; col++) | |

2651 | for (i=0; i<ntrans*mfit_local; i++) | |

2652 | A.P[row][col] -= QS[row][i] * A.Q[col][i]; /* Q = R^T */ | |

2653 | ||

2654 | /* And invert it to get P' */ | |

2655 | if (GCI_gauss_jordan(A.P, mfit_global, NULL) != 0) | |

2656 | return -3; | |

2657 | ||

2658 | /* Now overwrite Q with Q' */ | |

2659 | for (row=0; row<mfit_global; row++) | |

2660 | for (col=0; col<ntrans*mfit_local; col++) { | |

2661 | A.Q[row][col] = 0; | |

2662 | for (i=0; i<mfit_global; i++) | |

2663 | A.Q[row][col] -= A.P[row][i] * QS[i][col]; /* P contains P' */ | |

2664 | } | |

2665 | ||

2666 | /* Finally, we can solve to find x */ | |

2667 | /* We have x.global = P'.(b.global) + Q'.(b.local) | |

2668 | and x.local = R'.(b.global) + S'.(b.local) | |

2669 | We do x.global first. */ | |

2670 | for (row=0; row<mfit_global; row++) { | |

2671 | x.global[row] = 0; | |

2672 | for (i=0; i<mfit_global; i++) | |

2673 | x.global[row] += A.P[row][i] * b.global[i]; | |

2674 | /* Recall that Q now contains Q' */ | |

2675 | for (i=0; i < ntrans * mfit_local; i++) | |

2676 | x.global[row] += A.Q[row][i] * b.local[i]; | |

2677 | } | |

2678 | ||

2679 | /* Now x_local; the R'.b_global component first, recalling that R' | |

2680 | is Q' transposed and that Q' is stored in Q. */ | |

2681 | for (row=0; row < ntrans * mfit_local; row++) { | |

2682 | x.local[row] = 0; | |

2683 | for (i=0; i<mfit_global; i++) | |

2684 | x.local[row] += A.Q[i][row] * b.global[i]; | |

2685 | } | |

2686 | ||

2687 | /* Now S' = S^{-1} + (S^{-1}.R).(P-Q.S^{-1}.R)^{-1}.(Q.S^{-1}). | |

2688 | We first handle the S^{-1} term, then the remaining term. */ | |

2689 | for (block=0; block<ntrans; block++) | |

2690 | for (row=0; row<mfit_local; row++) { | |

2691 | for (j=0; j<mfit_local; j++) | |

2692 | x.local[block*mfit_local + row] += | |

2693 | A.S[block][row][j] * b.local[block*mfit_local + j]; | |

2694 | } | |

2695 | ||

2696 | /* For the remaining term, we have an x.local[row] contribution of | |

2697 | ||

2698 | sum_j sum_k sum_l | |

2699 | (S^{-1}.R)_{row,j} . (P-Q.S^{-1}.R)^{-1}_{j,k} . | |

2700 | (Q.S^{-1})_{k,l} . b.local_{l} | |

2701 | ||

2702 | In order to save computations, we calculate the matrices once | |

2703 | only. We start with (Q.S^{-1}) . b_local, which is an | |

2704 | mfit_global x 1 array, and store this in x_temp; premultiplying | |

2705 | this by the square matrix P' again gives an mfit_global x 1 | |

2706 | array, which goes in x_temp2, then premultiplying this by | |

2707 | (S^{-1}.R) gives an (ntrans * mfit_local) x 1 array which is | |

2708 | added directly onto x_local. Recall also that S^{-1}.R is the | |

2709 | transpose of Q.S^{-1}, which is currently stored in QS, and | |

2710 | that the middle term is currently stored in P. */ | |

2711 | ||

2712 | for (row=0; row<mfit_global; row++) { | |

2713 | x_temp[row] = 0; | |

2714 | for (i=0; i < ntrans*mfit_local; i++) | |

2715 | x_temp[row] += QS[row][i] * b.local[i]; | |

2716 | } | |

2717 | for (row=0; row<mfit_global; row++) { | |

2718 | x_temp2[row] = 0; | |

2719 | for (i=0; i<mfit_global; i++) | |

2720 | x_temp2[row] += A.P[row][i] * x_temp[i]; | |

2721 | } | |

2722 | /* Again, S^{-1}.R is the transpose of Q.S^{-1} */ | |

2723 | for (row=0; row < ntrans * mfit_local; row++) { | |

2724 | for (i=0; i<mfit_global; i++) | |

2725 | x.local[row] += QS[i][row] * x_temp2[i]; | |

2726 | } | |

2727 | ||

2728 | /* And we're done, once we've copied x into b */ | |

2729 | GCI_copy_global_vector(b, x, mfit_global, mfit_local, ntrans); | |

2730 | ||

2731 | return 0; | |

2732 | } | |

2733 | ||

2734 | ||

2735 | /* ***** GENERIC FUNCTION CODE ***** */ | |

2736 | ||

2737 | /* These functions are essentially the same as the above functions | |

2738 | GCI_marquardt_global_exps_instr and | |

2739 | GCI_marquardt_global_exps_do_fit_instr and the latter's dependents, | |

2740 | except that this version takes an arbitrary function and a list of | |

2741 | global parameters. This function is designed to be called from | |

2742 | external code. It is nowhere near as efficient as the streamlined | |

2743 | code above for globally fitting taus for multi-exponential models. | |

2744 | Also, this function must be provided with meaningful starting | |

2745 | estimates for all parameters. | |

2746 | */ | |

2747 | ||

2748 | int GCI_marquardt_global_generic_instr(float xincr, float **trans, | |

2749 | int ndata, int ntrans, int fit_start, int fit_end, | |

2750 | float instr[], int ninstr, | |

2751 | noise_type noise, float sig[], | |

2752 | float **param, int paramfree[], int nparam, int gparam[], | |

2753 | restrain_type restrain, float chisq_delta, | |

2754 | void (*fitfunc)(float, float [], float *, float [], int), | |

2755 | float **fitted, float **residuals, | |

2756 | float chisq_trans[], float *chisq_global, int *df) | |

2757 | { | |

2758 | float **covar, **alpha, *scaled_instr, instrsum; | |

2759 | int i, ret; | |

2760 | int mlocal, mglobal; | |

2761 | ||

2762 | /* Some basic parameter checks */ | |

2763 | if (xincr <= 0) return -1; | |

2764 | if (ntrans < 1) return -1; | |

2765 | if (ndata < 1) return -1; | |

2766 | if (fit_start < 0 || fit_end > ndata) return -1; | |

2767 | if (ninstr < 1) return -1; | |

2768 | if (nparam < 1) return -1; | |

2769 | ||

2770 | if ((covar = GCI_ecf_matrix(nparam, nparam)) == NULL) | |

2771 | return -2; | |

2772 | ||

2773 | if ((alpha = GCI_ecf_matrix(nparam, nparam)) == NULL) { | |

2774 | GCI_ecf_free_matrix(covar); | |

2775 | return -3; | |

2776 | } | |

2777 | ||

2778 | if ((scaled_instr = (float *) malloc(ninstr * sizeof(float))) == NULL) { | |

2779 | GCI_ecf_free_matrix(covar); | |

2780 | GCI_ecf_free_matrix(alpha); | |

2781 | return -4; | |

2782 | } | |

2783 | ||

2784 | /* Scale the instrument response */ | |

2785 | for (i=0, instrsum=0; i<ninstr; i++) | |

2786 | instrsum += instr[i]; | |

2787 | if (instrsum == 0) return -6; | |

2788 | for (i=0; i<ninstr; i++) | |

2789 | scaled_instr[i] = instr[i] / instrsum; | |

2790 | ||

2791 | /* Now call the global fitting function. */ | |

2792 | ret = GCI_marquardt_global_generic_do_fit_instr( | |

2793 | xincr, trans, ndata, ntrans, fit_start, fit_end, | |

2794 | scaled_instr, ninstr, noise, sig, | |

2795 | param, paramfree, nparam, gparam, restrain, chisq_delta, | |

2796 | fitfunc, fitted, residuals, covar, alpha, | |

2797 | chisq_trans, chisq_global); | |

2798 | ||

2799 | GCI_ecf_free_matrix(covar); | |

2800 | GCI_ecf_free_matrix(alpha); | |

2801 | free(scaled_instr); | |

2802 | GCI_marquardt_cleanup(); | |

2803 | ||

2804 | if (ret < 0) { | |

2805 | dbgprintf(1, "Fit failed, ret = %d\n", ret); | |

2806 | return -10 + ret; | |

2807 | } | |

2808 | ||

2809 | dbgprintf(2, "Fit succeeded, ret = %d\n", ret); | |

2810 | ||

2811 | /* Before we return, calculate the number of degrees of freedom */ | |

2812 | /* The number of degrees of freedom is given by: | |

2813 | d.f. = ntrans * ((fit_end - fit_start) - # free local parameters) | |

2814 | - # free global parameters | |

2815 | */ | |

2816 | ||

2817 | mglobal = mlocal = 0; | |

2818 | for (i=0; i<nparam; i++) | |

2819 | if (paramfree[i]) { | |

2820 | if (gparam[i]) mglobal++; | |

2821 | else mlocal++; | |

2822 | } | |

2823 | ||

2824 | *df = ntrans * ((fit_end - fit_start) - mlocal) - mglobal; | |

2825 | ||

2826 | return ret; | |

2827 | } | |

2828 | ||

2829 | ||

2830 | int GCI_marquardt_global_generic_do_fit_instr( | |

2831 | float xincr, float **trans, int ndata, int ntrans, | |

2832 | int fit_start, int fit_end, float instr[], int ninstr, | |

2833 | noise_type noise, float sig[], | |

2834 | float **param, int paramfree[], int nparam, int gparam[], | |

2835 | restrain_type restrain, float chisq_delta, | |

2836 | void (*fitfunc)(float, float [], float *, float [], int), | |

2837 | float **fitted, float **residuals, | |

2838 | float **covar_scratch, float **alpha_scratch, | |

2839 | float *chisq_trans, float *chisq_global) | |

2840 | { | |

2841 | // PRB 03/07 Although **fitted and **residuals are provided only one "transient" is required and used, fitted[0] and residuals[0] | |

2842 | ||

2843 | float alambda, ochisq_global, *ochisq_trans; | |

2844 | int i, k, itst, itst_max; | |

2845 | int ret; | |

2846 | ||

2847 | itst_max = (restrain == ECF_RESTRAIN_DEFAULT) ? 4 : 6; | |

2848 | ||

2849 | /* If there are no global parameters being fitted, we simply fit | |

2850 | each local set. */ | |

2851 | for (i=0; i<nparam; i++) { | |

2852 | if (gparam[i] && paramfree[i]) { | |

2853 | i = -1; /* sentinel value */ | |

2854 | break; | |

2855 | } | |

2856 | } | |

2857 | ||

2858 | if (i >= 0) { /* no globals to fit */ | |

2859 | *chisq_global = 0; | |

2860 | ||

2861 | for (i=0; i<ntrans; i++) { | |

2862 | ret = GCI_marquardt_instr(xincr, trans[i], | |

2863 | ndata, fit_start, fit_end, instr, ninstr, noise, sig, | |

2864 | param[i], paramfree, nparam, restrain, fitfunc, | |

2865 | fitted[0], residuals[0], covar_scratch, alpha_scratch, | |

2866 | &chisq_trans[i], chisq_delta, 0, NULL); | |

2867 | if (ret < 0) { | |

2868 | dbgprintf(1, "In do_fit_instr, marquardt_instr returned %d " | |

2869 | "for transient %d\n", ret, i); | |

2870 | return -10 + ret; | |

2871 | } else { | |

2872 | *chisq_global += chisq_trans[i]; | |

2873 | } | |

2874 | } | |

2875 | return 0; | |

2876 | } | |

2877 | ||

2878 | /* If there are no free local variables to fit, we still do the | |

2879 | global fitting, but we have to be a little careful in some of | |

2880 | the later routines */ | |

2881 | ||

2882 | /* Now allocate all of the arrays we will need. */ | |

2883 | ||

2884 | if ((ochisq_trans = (float *) malloc(ntrans * sizeof(float))) == NULL) | |

2885 | return -1; | |

2886 | ||

2887 | /* We now begin our standard Marquardt loop, with several | |

2888 | modifications */ | |

2889 | alambda = -1; | |

2890 | ret = GCI_marquardt_global_generic_global_step( | |

2891 | xincr, trans, ndata, ntrans, fit_start, fit_end, | |

2892 | instr, ninstr, noise, sig, | |

2893 | param, paramfree, nparam, gparam, restrain, chisq_delta, fitfunc, | |

2894 | fitted, residuals, chisq_trans, chisq_global, | |

2895 | alpha_scratch, &alambda); | |

2896 | if (ret != 0) { | |

2897 | dbgprintf(1, "In do_fit_instr, first global_step returned %d\n", ret); | |

2898 | if (ret != -1) { | |

2899 | /* Wasn't a memory error, so unallocate arrays */ | |

2900 | alambda = 0.0; | |

2901 | GCI_marquardt_global_generic_global_step( | |

2902 | xincr, trans, ndata, ntrans, fit_start, fit_end, | |

2903 | instr, ninstr, noise, sig, | |

2904 | param, paramfree, nparam, gparam, restrain, chisq_delta, fitfunc, | |

2905 | fitted, residuals, chisq_trans, chisq_global, | |

2906 | alpha_scratch, &alambda); | |

2907 | } | |

2908 | free(ochisq_trans); | |

2909 | return ret; | |

2910 | } | |

2911 | ||

2912 | k = 1; /* Iteration counter */ | |

2913 | itst = 0; | |

2914 | for (;;) { | |

2915 | dbgprintf(3, "In do_fit_instr, beginning iteration %d:\n", k); | |

2916 | dbgprintf(3, " itst = %d, chisq_global = %.4f\n", itst, *chisq_global); | |

2917 | ||

2918 | k++; | |

2919 | if (k > MAXITERS) { | |

2920 | return -2; | |

2921 | } | |

2922 | ||

2923 | ochisq_global = *chisq_global; | |

2924 | for (i=0; i<ntrans; i++) | |

2925 | ochisq_trans[i] = chisq_trans[i]; | |

2926 | ||

2927 | ret = GCI_marquardt_global_generic_global_step( | |

2928 | xincr, trans, ndata, ntrans, fit_start, fit_end, | |

2929 | instr, ninstr, noise, sig, | |

2930 | param, paramfree, nparam, gparam, restrain, chisq_delta, fitfunc, | |

2931 | fitted, residuals, chisq_trans, chisq_global, | |

2932 | alpha_scratch, &alambda); | |

2933 | if (ret != 0) { | |

2934 | dbgprintf(1, "In do_fit_instr, second global_step returned %d\n", | |

2935 | ret); | |

2936 | /* Unallocate arrays */ | |

2937 | alambda = 0.0; | |

2938 | GCI_marquardt_global_generic_global_step( | |

2939 | xincr, trans, ndata, ntrans, fit_start, fit_end, | |

2940 | instr, ninstr, noise, sig, | |

2941 | param, paramfree, nparam, gparam, restrain, chisq_delta, fitfunc, | |

2942 | fitted, residuals, chisq_trans, chisq_global, | |

2943 | alpha_scratch, &alambda); | |

2944 | free(ochisq_trans); | |

2945 | return ret; | |

2946 | } | |

2947 | ||

2948 | if (*chisq_global > ochisq_global) | |

2949 | itst = 0; | |

2950 | else { | |

2951 | /* Let's try this approach; I really don't know what will | |

2952 | be best */ | |

2953 | float maxdiff; | |

2954 | ||

2955 | maxdiff = 0.0; | |

2956 | for (i=0; i<ntrans; i++) | |

2957 | if (ochisq_trans[i] - chisq_trans[i] > maxdiff) | |

2958 | maxdiff = ochisq_trans[i] - chisq_trans[i]; | |

2959 | ||

2960 | if (maxdiff < chisq_delta) | |

2961 | itst++; | |

2962 | dbgprintf(3, "In do_fit_instr, maxdiff = %.3f:\n", maxdiff); | |

2963 | } | |

2964 | ||

2965 | if (itst < itst_max) continue; | |

2966 | ||

2967 | /* Endgame */ | |

2968 | alambda = 0.0; | |

2969 | ret = GCI_marquardt_global_generic_global_step( | |

2970 | xincr, trans, ndata, ntrans, fit_start, fit_end, | |

2971 | instr, ninstr, noise, sig, | |

2972 | param, paramfree, nparam, gparam, restrain, chisq_delta, fitfunc, | |

2973 | fitted, residuals, chisq_trans, chisq_global, | |

2974 | alpha_scratch, &alambda); | |

2975 | if (ret != 0) { | |

2976 | dbgprintf(1, "In do_fit_instr, final global_step returned %d\n", | |

2977 | ret); | |

2978 | free(ochisq_trans); | |

2979 | return ret; | |

2980 | } | |

2981 | ||

2982 | free(ochisq_trans); | |

2983 | return k; /* We're done now */ | |

2984 | } | |

2985 | } | |

2986 | ||

2987 | ||

2988 | /* And this one is basically a specialised GCI_marquardt_instr_step | |

2989 | for the global fitting setup. */ | |

2990 | ||

2991 | #define do_frees \ | |

2992 | if (fnvals) free(fnvals);\ | |

2993 | if (dy_dparam_pure) GCI_ecf_free_matrix(dy_dparam_pure);\ | |

2994 | if (dy_dparam_conv) GCI_ecf_free_matrix(dy_dparam_conv); | |

2995 | ||

2996 | int GCI_marquardt_global_generic_global_step( | |

2997 | float xincr, float **trans, | |

2998 | int ndata, int ntrans, int fit_start, int fit_end, | |

2999 | float instr[], int ninstr, | |

3000 | noise_type noise, float sig[], | |

3001 | float **param, int paramfree[], int nparam, int gparam[], | |

[7611] | 3002 | restrain_type restrain, float chisq_delta, |

[7606] | 3003 | void (*fitfunc)(float, float [], float *, float [], int), |

3004 | float **yfit, float **dy, | |

3005 | float *chisq_trans, float *chisq_global, | |

3006 | float **alpha_scratch, float *alambda) | |

3007 | { | |

3008 | int i, j, ret; | |

3009 | static global_matrix alpha, covar; | |

3010 | static global_vector beta, dparam; | |

3011 | static float **paramtry; | |

3012 | static int mfit_local, mfit_global; | |

3013 | static int gindex[MAXFIT], lindex[MAXFIT]; | |

3014 | static float ochisq_global, *ochisq_trans; | |

3015 | static int initialised=0; | |

3016 | ||

3017 | // The following are declared here to retain some optimisation by not repeatedly mallocing | |

3018 | // (only once per transient), but to remain thread safe. | |

3019 | // They are malloced by lower fns but at the end, freed by this fn. | |

3020 | // These vars were global or static before thread safety was introduced. | |

3021 | float *fnvals=NULL, **dy_dparam_pure=NULL, **dy_dparam_conv=NULL; | |

3022 | int fnvals_len=0, dy_dparam_nparam_size=0; | |

3023 | ||

3024 | if (nparam > MAXFIT) | |

3025 | return -10; | |

3026 | if (xincr <= 0) | |

3027 | return -11; | |

3028 | if (fit_start < 0 || fit_start > fit_end || fit_end > ndata) | |

3029 | return -12; | |

3030 | ||

3031 | /* Initialisation */ | |

3032 | /* We assume we're given sensible starting values for param[] */ | |

3033 | if (*alambda < 0.0) { | |

3034 | /* Start by allocating lots of variables we will need */ | |

3035 | mfit_local = mfit_global = 0; | |

3036 | ||

3037 | for (i=0; i<nparam; i++) { | |

3038 | if (paramfree[i]) { | |

3039 | if (gparam[i]) | |

3040 | gindex[mfit_global++] = i; | |

3041 | else | |

3042 | lindex[mfit_local++] = i; | |

3043 | } | |

3044 | } | |

3045 | ||

3046 | if (initialised) { | |

3047 | GCI_free_global_matrix(&alpha); GCI_free_global_matrix(&covar); | |

3048 | GCI_ecf_free_matrix(paramtry); GCI_free_global_vector(&beta); | |

3049 | GCI_free_global_vector(&dparam); free(ochisq_trans); | |

3050 | initialised = 0; | |

3051 | } | |

3052 | ||

3053 | if (GCI_alloc_global_matrix(&alpha, mfit_global, mfit_local, ntrans) | |

3054 | != 0) | |

3055 | return -1; | |

3056 | ||

3057 | if (GCI_alloc_global_matrix(&covar, mfit_global, mfit_local, ntrans) | |

3058 | != 0) { | |

3059 | GCI_free_global_matrix(&alpha); | |

3060 | return -1; | |

3061 | } | |

3062 | ||

3063 | if ((paramtry = GCI_ecf_matrix(ntrans, nparam)) == NULL) { | |

3064 | GCI_free_global_matrix(&alpha); GCI_free_global_matrix(&covar); | |

3065 | return -1; | |

3066 | } | |

3067 | ||

3068 | if (GCI_alloc_global_vector(&beta, mfit_global, mfit_local, ntrans) | |

3069 | != 0) { | |

3070 | GCI_free_global_matrix(&alpha); GCI_free_global_matrix(&covar); | |

3071 | GCI_ecf_free_matrix(paramtry); | |

3072 | return -1; | |

3073 | } | |

3074 | ||

3075 | if (GCI_alloc_global_vector(&dparam, mfit_global, mfit_local, ntrans) | |

3076 | != 0) { | |

3077 | GCI_free_global_matrix(&alpha); GCI_free_global_matrix(&covar); | |

3078 | GCI_ecf_free_matrix(paramtry); GCI_free_global_vector(&beta); | |

3079 | return -1; | |

3080 | } | |

3081 | ||

3082 | if ((ochisq_trans = (float *) malloc(ntrans * sizeof(float))) | |

3083 | == NULL) { | |

3084 | GCI_free_global_matrix(&alpha); GCI_free_global_matrix(&covar); | |

3085 | GCI_ecf_free_matrix(paramtry); GCI_free_global_vector(&beta); | |

3086 | GCI_free_global_vector(&dparam); | |

3087 | return -1; | |

3088 | } | |

3089 | ||

3090 | initialised = 1; | |

3091 | ||

3092 | if (GCI_marquardt_global_compute_global_generic_fn( | |

3093 | xincr, trans, ndata, ntrans, | |

3094 | fit_start, fit_end, instr, ninstr, noise, sig, | |

3095 | param, paramfree, nparam, gparam, | |

3096 | mfit_global, mfit_local, gindex, lindex, fitfunc, | |

3097 | yfit, dy, alpha, beta, alpha_scratch, | |

3098 | chisq_trans, chisq_global, *alambda, | |

3099 | &fnvals, &dy_dparam_pure, &dy_dparam_conv, | |

3100 | &fnvals_len, &dy_dparam_nparam_size) != 0) | |

3101 | return -2; | |

3102 | ||

3103 | *alambda = 0.001; | |

3104 | ochisq_global = *chisq_global; | |

3105 | for (i=0; i<ntrans; i++) | |

3106 | ochisq_trans[i] = chisq_trans[i]; | |

3107 | ||

3108 | /* Initialise paramtry to param */ | |

3109 | for (i=0; i<ntrans; i++) { | |

3110 | for (j=0; j<nparam; j++) | |

3111 | paramtry[i][j] = param[i][j]; | |

3112 | } | |

3113 | } | |

3114 | ||

3115 | /* Once converged, evaluate covariance matrix */ | |

3116 | if (*alambda == 0) { | |

3117 | if (GCI_marquardt_global_compute_global_generic_fn_final( | |

3118 | xincr, trans, ndata, ntrans, | |

3119 | fit_start, fit_end, instr, ninstr, noise, sig, | |

3120 | param, paramfree, nparam, gparam, | |

3121 | mfit_global, mfit_local, gindex, lindex, fitfunc, | |

3122 | yfit, dy, chisq_trans, chisq_global, | |

3123 | &fnvals, &dy_dparam_pure, &dy_dparam_conv, | |

3124 | &fnvals_len, &dy_dparam_nparam_size) != 0) | |

3125 | return -3; | |

3126 | /* Don't need to do this here; if we wished to, we'd have to | |

3127 | move this code (the "if (*alambda == 0)" block) to after | |

3128 | the Gauss-Jordan call. We'd also need to rewrite it for | |

3129 | our situation.... */ | |

3130 | // if (mfit < nparam) { /* no need to do this otherwise */ | |

3131 | // GCI_covar_sort(covar, nparam, paramfree, mfit); | |

3132 | // GCI_covar_sort(alpha, nparam, paramfree, mfit); | |

3133 | // } | |

3134 | GCI_free_global_matrix(&alpha); GCI_free_global_matrix(&covar); | |

3135 | GCI_ecf_free_matrix(paramtry); GCI_free_global_vector(&beta); | |

3136 | GCI_free_global_vector(&dparam); free(ochisq_trans); | |

3137 | initialised = 0; | |

3138 | return 0; | |

3139 | } | |

3140 | ||

3141 | /* Alter linearised fitting matrix by augmenting diagonal | |

3142 | elements. */ | |

3143 | GCI_copy_global_matrix(covar, alpha, mfit_global, mfit_local, ntrans); | |

3144 | GCI_copy_global_vector(dparam, beta, mfit_global, mfit_local, ntrans); | |

3145 | for (j=0; j<mfit_global; j++) | |

3146 | covar.P[j][j] *= 1.0 + (*alambda); | |

3147 | for (i=0; i<ntrans; i++) | |

3148 | for (j=0; j<mfit_local; j++) | |

3149 | covar.S[i][j][j] *= 1.0 + (*alambda); | |

3150 | ||

3151 | /* Matrix solution; GCI_gauss_jordan solves Ax=b rather than AX=B */ | |

3152 | if (GCI_marquardt_global_solve_eqn(covar, dparam, | |

3153 | mfit_global, mfit_local, ntrans) != 0) | |

3154 | return -3; | |

3155 | ||

3156 | /* Did the trial succeed? Modify param by dparam... */ | |

3157 | for (i=0; i<ntrans; i++) { | |

3158 | for (j=0; j<mfit_global; j++) | |

3159 | paramtry[i][gindex[j]] = param[i][gindex[j]] + dparam.global[j]; | |

3160 | for (j=0; j<mfit_local; j++) | |

3161 | paramtry[i][lindex[j]] = | |

3162 | param[i][lindex[j]] + dparam.local[i*mfit_local + j]; | |

3163 | } | |

3164 | ||

3165 | for (i=0; i<ntrans; i++) { | |

3166 | if (restrain == ECF_RESTRAIN_DEFAULT) | |

3167 | ret = check_ecf_params (paramtry[i], nparam, fitfunc); | |

3168 | else | |

3169 | ret = check_ecf_user_params (paramtry[i], nparam, fitfunc); | |

3170 | ||

3171 | if (ret != 0) { | |

3172 | /* Bad parameters, increase alambda and return */ | |

3173 | *alambda *= 10.0; | |

3174 | return 0; | |

3175 | } | |

3176 | } | |

3177 | ||

3178 | if (GCI_marquardt_global_compute_global_generic_fn( | |

3179 | xincr, trans, ndata, ntrans, | |

3180 | fit_start, fit_end, instr, ninstr, noise, sig, | |

3181 | paramtry, paramfree, nparam, gparam, | |

3182 | mfit_global, mfit_local, gindex, lindex, fitfunc, | |

3183 | yfit, dy, covar, dparam, alpha_scratch, | |

3184 | chisq_trans, chisq_global, *alambda, | |

3185 | &fnvals, &dy_dparam_pure, &dy_dparam_conv, | |

3186 | &fnvals_len, &dy_dparam_nparam_size) != 0) | |

3187 | return -2; | |

3188 | ||

3189 | /* Success, accept the new solution */ | |

3190 | if (*chisq_global < ochisq_global) { | |

3191 | *alambda *= 0.1; | |

3192 | ochisq_global = *chisq_global; | |

3193 | for (i=0; i<ntrans; i++) | |

3194 | ochisq_trans[i] = chisq_trans[i]; | |

3195 | GCI_copy_global_matrix(alpha, covar, mfit_global, mfit_local, ntrans); | |

3196 | GCI_copy_global_vector(beta, dparam, mfit_global, mfit_local, ntrans); | |

3197 | for (i=0; i<ntrans; i++) { | |

3198 | for (j=0; j<nparam; j++) | |

3199 | param[i][j] = paramtry[i][j]; | |

3200 | } | |

3201 | } else { /* Failure, increase alambda and return */ | |

3202 | *alambda *= 10.0; | |

3203 | *chisq_global = ochisq_global; | |

3204 | for (i=0; i<ntrans; i++) | |

3205 | chisq_trans[i] = ochisq_trans[i]; | |

3206 | } | |

3207 | ||

3208 | return 0; | |

3209 | } | |

3210 | ||

3211 | ||

3212 | /* Here we use alpha only for scratch space */ | |

3213 | int GCI_marquardt_global_compute_global_generic_fn( | |

3214 | float xincr, float **trans, int ndata, int ntrans, | |

3215 | int fit_start, int fit_end, float instr[], int ninstr, | |

3216 | noise_type noise, float sig[], | |

3217 | float **param, int paramfree[], int nparam, int gparam[], | |

3218 | int mfit_global, int mfit_local, int gindex[], int lindex[], | |

3219 | void (*fitfunc)(float, float [], float *, float [], int), | |

3220 | float **yfit, float **dy, global_matrix alpha, global_vector beta, | |

3221 | float **alpha_scratch, float *chisq_trans, float *chisq_global, | |

3222 | float alambda, | |

3223 | float **pfnvals, float ***pdy_dparam_pure, float ***pdy_dparam_conv, | |

3224 | int *pfnvals_len, int *pdy_dparam_nparam_size) | |

3225 | { | |

3226 | int i, j, k, ret; | |

3227 | float beta_scratch[MAXFIT]; /* scratch space */ | |

3228 | ||

3229 | /* We initialise P and beta_global to zero; the others don't | |

3230 | matter, as they will be totally overwritten */ | |

3231 | for (i=0; i<mfit_global; i++) { | |

3232 | for (j=0; j<mfit_global; j++) | |

3233 | alpha.P[i][j] = 0; | |

3234 | beta.global[i] = 0; | |

3235 | } | |

3236 | *chisq_global = 0.0; | |

3237 | ||

3238 | for (i=0; i<ntrans; i++) { | |

3239 | /* Only pass the true alambda, used for initialisation, for | |

3240 | the first transient */ | |

3241 | ret = GCI_marquardt_compute_fn_instr( | |

3242 | xincr, trans[i], ndata, fit_start, fit_end, | |

3243 | instr, ninstr, noise, sig, | |

3244 | param[i], paramfree, nparam, fitfunc, | |

3245 | //// yfit[i], dy[i], alpha_scratch, beta_scratch, | |

3246 | yfit[0], dy[0], alpha_scratch, beta_scratch, | |

3247 | &chisq_trans[i], (i == 0) ? alambda : 0.0, | |

3248 | pfnvals, pdy_dparam_pure, pdy_dparam_conv, | |

3249 | pfnvals_len, pdy_dparam_nparam_size); | |

3250 | ||

3251 | if (ret != 0) { | |

3252 | dbgprintf(1, "In compute_global_generic_fn, " | |

3253 | "compute_fn_instr returned %d for transient %d\n", | |

3254 | ret, i); | |

3255 | return -2; | |

3256 | } | |

3257 | ||

3258 | /* So now have to populate alpha and beta with the contents of | |

3259 | alpha_scratch and beta_scratch. */ | |

3260 | ||

3261 | for (j=0; j<mfit_global; j++) { | |

3262 | for (k=0; k<mfit_global; k++) | |

3263 | alpha.P[j][k] += alpha_scratch[gindex[j]][gindex[k]]; | |

3264 | for (k=0; k<mfit_local; k++) | |

3265 | alpha.Q[j][i*mfit_local + k] = | |

3266 | alpha_scratch[gindex[j]][lindex[k]]; | |

3267 | beta.global[j] += beta_scratch[gindex[j]]; | |

3268 | } | |

3269 | for (j=0; j<mfit_local; j++) { | |

3270 | for (k=0; k<mfit_local; k++) | |

3271 | alpha.S[i][j][k] = alpha_scratch[lindex[j]][lindex[k]]; | |

3272 | beta.local[i*mfit_local + j] = beta_scratch[lindex[j]]; | |

3273 | } | |

3274 | ||

3275 | *chisq_global += chisq_trans[i]; | |

3276 | } | |

3277 | ||

3278 | return 0; | |

3279 | } | |

3280 | ||

3281 | ||

3282 | /* And the final variant */ | |

3283 | int GCI_marquardt_global_compute_global_generic_fn_final( | |

3284 | float xincr, float **trans, int ndata, int ntrans, | |

3285 | int fit_start, int fit_end, float instr[], int ninstr, | |

3286 | noise_type noise, float sig[], | |

3287 | float **param, int paramfree[], int nparam, int gparam[], | |

3288 | int mfit_global, int mfit_local, int gindex[], int lindex[], | |

3289 | void (*fitfunc)(float, float [], float *, float [], int), | |

3290 | float **yfit, float **dy, | |

3291 | float *chisq_trans, float *chisq_global, | |

3292 | float **pfnvals, float ***pdy_dparam_pure, float ***pdy_dparam_conv, | |

3293 | int *pfnvals_len, int *pdy_dparam_nparam_size) | |

3294 | { | |

3295 | int i, ret; | |

3296 | ||

3297 | *chisq_global = 0.0; | |

3298 | ||

3299 | for (i=0; i<ntrans; i++) { | |

3300 | /* Only pass the true alambda, used for initialisation, for | |

3301 | the first transient */ | |

3302 | ret = GCI_marquardt_compute_fn_final_instr( | |

3303 | xincr, trans[i], ndata, fit_start, fit_end, | |

3304 | instr, ninstr, noise, sig, | |

3305 | param[i], paramfree, nparam, fitfunc, | |

3306 | // yfit[i], dy[i], &chisq_trans[i]); | |

3307 | yfit[0], dy[0], &chisq_trans[i], | |

3308 | pfnvals, pdy_dparam_pure, pdy_dparam_conv, | |

3309 | pfnvals_len, pdy_dparam_nparam_size); | |

3310 | ||

3311 | if (ret != 0) { | |

3312 | dbgprintf(1, "In compute_global_generic_fn_final, " | |

3313 | "compute_fn_final_instr returned %d for transient %d\n", | |

3314 | ret, i); | |

3315 | return -2; | |

3316 | } | |

3317 | ||

3318 | *chisq_global += chisq_trans[i]; | |

3319 | } | |

3320 | ||

3321 | return 0; | |

3322 | } | |

3323 | ||

3324 | ||

3325 | // Emacs settings: | |

3326 | // Local variables: | |

3327 | // mode: c | |

3328 | // c-basic-offset: 4 | |

3329 | // tab-width: 4 | |

3330 | // End: |

**Note:**See TracBrowser for help on using the repository browser.