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

Revision 7387, 8.0 KB checked in by aivar, 9 years ago (diff)

A version of Paul Barber's ECF library with functions used in SLIM Plugin. These have been rewritten to remove NR dependencies. Functions not used and other code that still uses NR all removed.

Line 
1#include "EcfWrapper.h"
2#include "ecf.h"
3#include <stdio.h>
4#include <stdlib.h>
5
6int RLD_fit(
7        double x_inc,
8        double y[],
9        int fit_start,
10        int fit_end,
11        double instr[],
12        int n_instr,
13        double sig[],
14        double *z,
15        double *a,
16        double *tau,
17        double fitted[],
18        double *chi_square,
19        double chi_square_target
20        ) {
21
22    int noise = 0;
23    int n_data = fit_end + 1;
24
25    float *y_float = (float *)malloc(n_data * sizeof(float));
26    float *sig_float = (float *)malloc(n_data * sizeof(float));
27    int i;
28    for (i = 0; i < n_data; ++i) {
29       y_float[i] = (float) y[i];
30       sig_float[i] = ( sig ? (float) sig[i] : 1.0f );
31    }
32
33    float *instr_float = 0;
34    if (instr) {
35        instr_float = (float *)malloc(n_instr * sizeof(float));
36        int i;
37        for (i = 0; i < n_instr; ++i) {
38            instr_float[i] = (float) instr[i];
39        }
40    }
41    float *fitted_float = 0;
42    if (fitted) {
43        fitted_float = (float*)malloc(n_data * sizeof(float));
44    }
45    float *residuals_float = (float *)malloc(n_data * sizeof(float));
46
47    float x_inc_float      = (float) x_inc;
48    float z_float          = (float) *z;
49    float a_float          = (float) *a;
50    float tau_float        = (float) *tau;
51    float chi_square_float = (float) *chi_square;
52
53    int returnValue =  GCI_triple_integral_fitting_engine(
54            x_inc_float,
55            y_float,
56            fit_start,
57            fit_end + 1, // want exclusive, not inclusive end
58            instr_float,
59            n_instr,
60            noise,
61            sig_float,
62            &z_float,
63            &a_float,
64            &tau_float,
65            fitted_float,
66            residuals_float,
67            &chi_square_float,
68            chi_square_target
69            );
70
71    *z          = (double) z_float;
72    *a          = (double) a_float;
73    *tau        = (double) tau_float;
74    *chi_square = (double) chi_square_float;
75   
76    free(y_float);
77    free(sig_float);
78
79    if (instr_float) {
80        free(instr_float);
81    }
82    if (fitted_float) {
83        int i;
84        for (i = 0; i < n_data; ++i) {
85            fitted[i] = (double) fitted_float[i];
86        }
87        free(fitted_float);
88    }
89    free(residuals_float);
90
91    return returnValue;
92}
93
94int LMA_fit(
95        double x_inc,
96        double y[],
97        int fit_start,
98        int fit_end,
99        double instr[],
100        int n_instr,
101        double sig[],
102        double param[],
103        int param_free[],
104        int n_param,
105        double fitted[],
106        double *chi_square,
107        double chi_square_target
108        ) {
109
110    int noise = 0;
111    int restrain = 0;
112    float chi_square_delta = 0;
113    float chi_square_percent = 500;
114    int n_data = fit_end + 1;
115
116    float *y_float = (float *)malloc(n_data * sizeof(float));
117    float *sig_float = (float *)malloc(n_data * sizeof(float));
118    int i;
119    for (i = 0; i < n_data; ++i) {
120       y_float[i] = (float) y[i];
121       sig_float[i] = ( sig ? (float) sig[i] : 1.0f );
122    }
123
124    float *instr_float = 0;
125    if (instr) {
126        instr_float = (float *)malloc(n_instr * sizeof(float));
127        int i;
128        for (i = 0; i < n_instr; ++i) {
129            instr_float[i] = (float) instr[i];
130        }
131    }
132    float *fitted_float = 0;
133    if (fitted) {
134        fitted_float = (float*)malloc(n_data * sizeof(float));
135    }
136    float *param_float = (float *)malloc(n_param * sizeof(float));
137    switch (n_param) {
138        case 3:
139            // single exponential fit
140            param_float[0] = (float) param[1]; // z
141            param_float[1] = (float) param[2]; // a
142            param_float[2] = (float) param[3]; // tau
143            break;
144        case 4:
145            // stretched exponential fit
146            param_float[0] = (float) param[1]; // z
147            param_float[1] = (float) param[2]; // a
148            param_float[2] = (float) param[3]; // tau
149            param_float[3] = (float) param[4]; // h
150            break;
151        case 5:
152            // double exponential fit
153            param_float[0] = (float) param[1]; // z
154            param_float[1] = (float) param[2]; // a1
155            param_float[2] = (float) param[3]; // tau1
156            param_float[3] = (float) param[4]; // a2
157            param_float[4] = (float) param[5]; // tau2
158            break;
159        case 7:
160            // triple exponential fit
161            param_float[0] = (float) param[1]; // z
162            param_float[1] = (float) param[2]; // a1
163            param_float[2] = (float) param[3]; // tau1
164            param_float[3] = (float) param[4]; // a2
165            param_float[4] = (float) param[5]; // tau2
166            param_float[5] = (float) param[6]; // a3
167            param_float[6] = (float) param[7]; // tau3
168            break;
169        default:
170            break;
171    }
172
173    float *residuals_float = (float *)malloc(n_data * sizeof(float));
174    float **covar_float    = GCI_ecf_matrix(n_data, n_data);
175    float **alpha_float    = GCI_ecf_matrix(n_data, n_data);
176    float **err_axes_float = GCI_ecf_matrix(n_data, n_data);
177
178    float x_inc_float      = (float) x_inc;
179    float chi_square_float = (float) *chi_square;
180
181    // choose appropriate fitting function
182    void (*fitfunc)(float, float [], float *, float[], int) = NULL;
183    if (4 == n_param) {
184        fitfunc = GCI_stretchedexp;
185    }
186    else {
187        fitfunc = GCI_multiexp_tau;
188    }
189
190    int returnValue = GCI_marquardt_fitting_engine(
191            x_inc_float, //TODO not necessary, just cast it, I think
192            y_float,
193            n_data,
194            fit_start,
195            fit_end + 1, // want exclusive, not inclusive end
196            instr_float,
197            n_instr,
198            noise,
199            sig_float,
200            param_float,
201            param_free,
202            n_param,
203            restrain,
204            fitfunc,
205            fitted_float,
206            residuals_float,
207            &chi_square_float,
208            covar_float,
209            alpha_float,
210            err_axes_float,
211            chi_square_target,
212            chi_square_delta,
213            chi_square_percent);
214
215    *chi_square = (double) chi_square_float;
216
217    free(y_float);
218    free(sig_float);
219
220    if (instr_float) {
221        free(instr_float);
222    }
223    if (fitted_float) {
224        int i;
225        for (i = 0; i < n_data; ++i) {
226            fitted[i] = (double) fitted_float[i];
227        }
228        free(fitted_float);
229    }
230    switch (n_param) {
231        case 3:
232            // single exponential fit
233            param[0] = (double) chi_square_float;
234            param[1] = (double) param_float[0]; // z
235            param[2] = (double) param_float[1]; // a
236            param[3] = (double) param_float[2]; // tau
237            break;
238        case 4:
239            // stretched exponential fit
240            param[0] = (double) chi_square_float;
241            param[1] = (double) param_float[0]; // z
242            param[2] = (double) param_float[1]; // a
243            param[3] = (double) param_float[2]; // tau
244            param[4] = (double) param_float[3]; // h
245            break;
246        case 5:
247            // double exponential fit
248            param[0] = (double) chi_square_float;
249            param[1] = (double) param_float[0]; // z
250            param[2] = (double) param_float[1]; // a1
251            param[3] = (double) param_float[2]; // tau1
252            param[4] = (double) param_float[3]; // a2
253            param[5] = (double) param_float[4]; // tau2
254            break;
255        case 7:
256            // triple exponential fit
257            param[0] = (double) chi_square_float;
258            param[1] = (double) param_float[0]; // z
259            param[2] = (double) param_float[1]; // a1
260            param[3] = (double) param_float[2]; // tau1
261            param[4] = (double) param_float[3]; // a2
262            param[5] = (double) param_float[4]; // tau2
263            param[6] = (double) param_float[5]; // a3
264            param[7] = (double) param_float[6]; // tau3
265            break;
266    }
267    free(param_float);
268   
269    free(residuals_float);
270    GCI_ecf_free_matrix(covar_float);
271    GCI_ecf_free_matrix(alpha_float);
272    GCI_ecf_free_matrix(err_axes_float);
273
274    //if (returnValue) {
275    //    printf("returnValue is %d\n", returnValue);
276    //}
277
278    return returnValue;
279}
280
Note: See TracBrowser for help on using the repository browser.