source: trunk/projects/slim-curve/src/main/c/EcfWrapper.c @ 7901

Revision 7901, 8.8 KB checked in by aivar, 8 years ago (diff)

SLIM Curve: Work in progress; tweak to agree with TRI2.

Line 
1/*
2This file is part of the SLIM-curve package for exponential curve fitting of spectral lifetime data.
3
4Copyright (c) 2010, 2011, Gray Institute University of Oxford & UW-Madison LOCI.
5
6    This program is free software: you can redistribute it and/or modify
7    it under the terms of the GNU General Public License as published by
8    the Free Software Foundation, either version 3 of the License, or
9    (at your option) any later version.
10
11    This program is distributed in the hope that it will be useful,
12    but WITHOUT ANY WARRANTY; without even the implied warranty of
13    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
14    GNU General Public License for more details.
15
16    You should have received a copy of the GNU General Public License
17    along with this program.  If not, see <http://www.gnu.org/licenses/>.
18 */
19
20#include "EcfWrapper.h"
21#include "Ecf.h"
22#include <stdio.h>
23#include <stdlib.h>
24
25int RLD_fit(
26        double x_inc,
27        double y[],
28        int fit_start,
29        int fit_end,
30        double instr[],
31        int n_instr,
32        int noise,
33        double sig[],
34        double *z,
35        double *a,
36        double *tau,
37        double fitted[],
38        double *chi_square,
39        double chi_square_target
40        ) {
41
42    int n_data = fit_end;
43    float *y_float = (float *)malloc(n_data * sizeof(float));
44    float *sig_float = (float *)malloc(n_data * sizeof(float));
45    float *instr_float = 0;
46    float *fitted_float = 0;
47    float *residuals_float = (float *)malloc(n_data * sizeof(float));
48    float x_inc_float      = (float) x_inc;
49    float z_float          = (float) *z;
50    float a_float          = (float) *a;
51    float tau_float        = (float) *tau;
52    float chi_square_float = (float) *chi_square;
53    int return_value;
54    int i;
55
56    for (i = 0; i < n_data; ++i) {
57       y_float[i] = (float) y[i];
58       sig_float[i] = ( sig ? (float) sig[i] : 1.0f );
59    }
60
61    if (instr) {
62        int i;
63        instr_float = (float *)malloc(n_instr * sizeof(float));
64        for (i = 0; i < n_instr; ++i) {
65            instr_float[i] = (float) instr[i];
66        }
67    }
68
69    if (fitted) {
70        fitted_float = (float*)malloc(n_data * sizeof(float));
71    }
72
73    return_value =  GCI_triple_integral_fitting_engine(
74            x_inc_float,
75            y_float,
76            fit_start,
77            fit_end,
78            instr_float,
79            n_instr,
80            noise,
81            sig_float,
82            &z_float,
83            &a_float,
84            &tau_float,
85            fitted_float,
86            residuals_float,
87            &chi_square_float,
88            chi_square_target
89            );
90
91    *z          = (double) z_float;
92    *a          = (double) a_float;
93    *tau        = (double) tau_float;
94    *chi_square = (double) chi_square_float;
95   
96    free(y_float);
97    free(sig_float);
98
99    if (instr_float) {
100        free(instr_float);
101    }
102    if (fitted_float) {
103        int i;
104        for (i = 0; i < n_data; ++i) {
105            fitted[i] = (double) fitted_float[i];
106        }
107        free(fitted_float);
108    }
109    free(residuals_float);
110
111    return return_value;
112}
113
114int LMA_fit(
115        double x_inc,
116        double y[],
117        int fit_start,
118        int fit_end,
119        double instr[],
120        int n_instr,
121        int noise,
122        double sig[],
123        double param[],
124        int param_free[],
125        int n_param,
126        double fitted[],
127        double *chi_square,
128        double chi_square_target,
129        double chi_square_delta
130        ) {
131
132    int restrain = 0;
133    float chi_square_percent = 95;
134
135    int n_data = fit_end;
136    float *y_float = (float *)malloc(n_data * sizeof(float));
137    float *sig_float = (float *)malloc(n_data * sizeof(float));
138    float *instr_float = 0;
139    float *fitted_float = 0;
140    float *param_float;
141    float *residuals_float = (float *)malloc(n_data * sizeof(float));
142    float **covar_float    = GCI_ecf_matrix(n_data, n_data);
143    float **alpha_float    = GCI_ecf_matrix(n_data, n_data);
144    float **err_axes_float = GCI_ecf_matrix(n_data, n_data);
145    float x_inc_float      = (float) x_inc;
146    float chi_square_float = (float) *chi_square;
147    void (*fitfunc)(float, float [], float *, float[], int) = NULL;
148    int return_value;
149    int i;
150
151    for (i = 0; i < n_data; ++i) {
152       y_float[i] = (float) y[i];
153       sig_float[i] = ( sig ? (float) sig[i] : 1.0f );
154    }
155
156    if (instr) {
157        int i;
158        instr_float = (float *)malloc(n_instr * sizeof(float));
159        for (i = 0; i < n_instr; ++i) {
160            instr_float[i] = (float) instr[i];
161        }
162    }
163
164    if (fitted) {
165        fitted_float = (float*)malloc(n_data * sizeof(float));
166    }
167    param_float = (float *)malloc(n_param * sizeof(float));
168    switch (n_param) {
169        case 3:
170            // single exponential fit
171            param_float[0] = (float) param[1]; // z
172            param_float[1] = (float) param[2]; // a
173            param_float[2] = (float) param[3]; // tau
174            break;
175        case 4:
176            // stretched exponential fit
177            param_float[0] = (float) param[1]; // z
178            param_float[1] = (float) param[2]; // a
179            param_float[2] = (float) param[3]; // tau
180            param_float[3] = (float) param[4]; // h
181            break;
182        case 5:
183            // double exponential fit
184            param_float[0] = (float) param[1]; // z
185            param_float[1] = (float) param[2]; // a1
186            param_float[2] = (float) param[3]; // tau1
187            param_float[3] = (float) param[4]; // a2
188            param_float[4] = (float) param[5]; // tau2
189            break;
190        case 7:
191            // triple exponential fit
192            param_float[0] = (float) param[1]; // z
193            param_float[1] = (float) param[2]; // a1
194            param_float[2] = (float) param[3]; // tau1
195            param_float[3] = (float) param[4]; // a2
196            param_float[4] = (float) param[5]; // tau2
197            param_float[5] = (float) param[6]; // a3
198            param_float[6] = (float) param[7]; // tau3
199            break;
200        default:
201            break;
202    }
203
204    // choose appropriate fitting function
205    if (4 == n_param) {
206        fitfunc = GCI_stretchedexp;
207    }
208    else {
209        fitfunc = GCI_multiexp_tau;
210    }
211
212    return_value = GCI_marquardt_fitting_engine(
213            x_inc_float, //TODO not necessary, just cast it, I think
214            y_float,
215            n_data,
216            fit_start,
217            fit_end,
218            instr_float,
219            n_instr,
220            noise,
221            sig_float,
222            param_float,
223            param_free,
224            n_param,
225            restrain,
226            fitfunc,
227            fitted_float,
228            residuals_float,
229            &chi_square_float,
230            covar_float,
231            alpha_float,
232            err_axes_float,
233            chi_square_target,
234            chi_square_delta,
235            chi_square_percent);
236
237    *chi_square = (double) chi_square_float;
238
239    free(y_float);
240    free(sig_float);
241
242    if (instr_float) {
243        free(instr_float);
244    }
245    if (fitted_float) {
246        int i;
247        for (i = 0; i < n_data; ++i) {
248            fitted[i] = (double) fitted_float[i];
249        }
250        free(fitted_float);
251    }
252    switch (n_param) {
253        case 3:
254            // single exponential fit
255            param[0] = (double) chi_square_float;
256            param[1] = (double) param_float[0]; // z
257            param[2] = (double) param_float[1]; // a
258            param[3] = (double) param_float[2]; // tau
259            break;
260        case 4:
261            // stretched exponential fit
262            param[0] = (double) chi_square_float;
263            param[1] = (double) param_float[0]; // z
264            param[2] = (double) param_float[1]; // a
265            param[3] = (double) param_float[2]; // tau
266            param[4] = (double) param_float[3]; // h
267            break;
268        case 5:
269            // double exponential fit
270            param[0] = (double) chi_square_float;
271            param[1] = (double) param_float[0]; // z
272            param[2] = (double) param_float[1]; // a1
273            param[3] = (double) param_float[2]; // tau1
274            param[4] = (double) param_float[3]; // a2
275            param[5] = (double) param_float[4]; // tau2
276            break;
277        case 7:
278            // triple exponential fit
279            param[0] = (double) chi_square_float;
280            param[1] = (double) param_float[0]; // z
281            param[2] = (double) param_float[1]; // a1
282            param[3] = (double) param_float[2]; // tau1
283            param[4] = (double) param_float[3]; // a2
284            param[5] = (double) param_float[4]; // tau2
285            param[6] = (double) param_float[5]; // a3
286            param[7] = (double) param_float[6]; // tau3
287            break;
288    }
289    free(param_float);
290   
291    free(residuals_float);
292    GCI_ecf_free_matrix(covar_float);
293    GCI_ecf_free_matrix(alpha_float);
294    GCI_ecf_free_matrix(err_axes_float);
295
296    return return_value;
297}
298
Note: See TracBrowser for help on using the repository browser.