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

Revision 7781, 8.9 KB checked in by paulbarber, 8 years ago (diff)

Added GPL license to all source code.

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        double sig[],
33        double *z,
34        double *a,
35        double *tau,
36        double fitted[],
37        double *chi_square,
38        double chi_square_target
39        ) {
40
41    int noise = 0;
42    int n_data = fit_end + 1;
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 + 1, // want exclusive, not inclusive 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        double sig[],
122        double param[],
123        int param_free[],
124        int n_param,
125        double fitted[],
126        double *chi_square,
127        double chi_square_target
128        ) {
129
130    int noise = 0;
131    int restrain = 0;
132    float chi_square_delta = 0;
133    float chi_square_percent = 500;
134    int n_data = fit_end + 1;
135    float *y_float = (float *)malloc(n_data * sizeof(float));
136    float *sig_float = (float *)malloc(n_data * sizeof(float));
137    float *instr_float = 0;
138    float *fitted_float = 0;
139    float *param_float;
140    float *residuals_float = (float *)malloc(n_data * sizeof(float));
141    float **covar_float    = GCI_ecf_matrix(n_data, n_data);
142    float **alpha_float    = GCI_ecf_matrix(n_data, n_data);
143    float **err_axes_float = GCI_ecf_matrix(n_data, n_data);
144    float x_inc_float      = (float) x_inc;
145    float chi_square_float = (float) *chi_square;
146    void (*fitfunc)(float, float [], float *, float[], int) = NULL;
147    int return_value;
148    int i;
149
150    for (i = 0; i < n_data; ++i) {
151       y_float[i] = (float) y[i];
152       sig_float[i] = ( sig ? (float) sig[i] : 1.0f );
153    }
154
155    if (instr) {
156        int i;
157        instr_float = (float *)malloc(n_instr * sizeof(float));
158        for (i = 0; i < n_instr; ++i) {
159            instr_float[i] = (float) instr[i];
160        }
161    }
162
163    if (fitted) {
164        fitted_float = (float*)malloc(n_data * sizeof(float));
165    }
166    param_float = (float *)malloc(n_param * sizeof(float));
167    switch (n_param) {
168        case 3:
169            // single exponential fit
170            param_float[0] = (float) param[1]; // z
171            param_float[1] = (float) param[2]; // a
172            param_float[2] = (float) param[3]; // tau
173            break;
174        case 4:
175            // stretched exponential fit
176            param_float[0] = (float) param[1]; // z
177            param_float[1] = (float) param[2]; // a
178            param_float[2] = (float) param[3]; // tau
179            param_float[3] = (float) param[4]; // h
180            break;
181        case 5:
182            // double exponential fit
183            param_float[0] = (float) param[1]; // z
184            param_float[1] = (float) param[2]; // a1
185            param_float[2] = (float) param[3]; // tau1
186            param_float[3] = (float) param[4]; // a2
187            param_float[4] = (float) param[5]; // tau2
188            break;
189        case 7:
190            // triple exponential fit
191            param_float[0] = (float) param[1]; // z
192            param_float[1] = (float) param[2]; // a1
193            param_float[2] = (float) param[3]; // tau1
194            param_float[3] = (float) param[4]; // a2
195            param_float[4] = (float) param[5]; // tau2
196            param_float[5] = (float) param[6]; // a3
197            param_float[6] = (float) param[7]; // tau3
198            break;
199        default:
200            break;
201    }
202
203    // choose appropriate fitting function
204    if (4 == n_param) {
205        fitfunc = GCI_stretchedexp;
206    }
207    else {
208        fitfunc = GCI_multiexp_tau;
209    }
210
211    return_value = GCI_marquardt_fitting_engine(
212            x_inc_float, //TODO not necessary, just cast it, I think
213            y_float,
214            n_data,
215            fit_start,
216            fit_end + 1, // want exclusive, not inclusive end
217            instr_float,
218            n_instr,
219            noise,
220            sig_float,
221            param_float,
222            param_free,
223            n_param,
224            restrain,
225            fitfunc,
226            fitted_float,
227            residuals_float,
228            &chi_square_float,
229            covar_float,
230            alpha_float,
231            err_axes_float,
232            chi_square_target,
233            chi_square_delta,
234            chi_square_percent);
235
236    *chi_square = (double) chi_square_float;
237
238    free(y_float);
239    free(sig_float);
240
241    if (instr_float) {
242        free(instr_float);
243    }
244    if (fitted_float) {
245        int i;
246        for (i = 0; i < n_data; ++i) {
247            fitted[i] = (double) fitted_float[i];
248        }
249        free(fitted_float);
250    }
251    switch (n_param) {
252        case 3:
253            // single exponential fit
254            param[0] = (double) chi_square_float;
255            param[1] = (double) param_float[0]; // z
256            param[2] = (double) param_float[1]; // a
257            param[3] = (double) param_float[2]; // tau
258            break;
259        case 4:
260            // stretched exponential fit
261            param[0] = (double) chi_square_float;
262            param[1] = (double) param_float[0]; // z
263            param[2] = (double) param_float[1]; // a
264            param[3] = (double) param_float[2]; // tau
265            param[4] = (double) param_float[3]; // h
266            break;
267        case 5:
268            // double exponential fit
269            param[0] = (double) chi_square_float;
270            param[1] = (double) param_float[0]; // z
271            param[2] = (double) param_float[1]; // a1
272            param[3] = (double) param_float[2]; // tau1
273            param[4] = (double) param_float[3]; // a2
274            param[5] = (double) param_float[4]; // tau2
275            break;
276        case 7:
277            // triple exponential fit
278            param[0] = (double) chi_square_float;
279            param[1] = (double) param_float[0]; // z
280            param[2] = (double) param_float[1]; // a1
281            param[3] = (double) param_float[2]; // tau1
282            param[4] = (double) param_float[3]; // a2
283            param[5] = (double) param_float[4]; // tau2
284            param[6] = (double) param_float[5]; // a3
285            param[7] = (double) param_float[6]; // tau3
286            break;
287    }
288    free(param_float);
289   
290    free(residuals_float);
291    GCI_ecf_free_matrix(covar_float);
292    GCI_ecf_free_matrix(alpha_float);
293    GCI_ecf_free_matrix(err_axes_float);
294
295    //if (returnValue) {
296    //    printf("returnValue is %d\n", returnValue);
297    //}
298
299    return return_value;
300}
301
Note: See TracBrowser for help on using the repository browser.