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

Revision 7723, 8.1 KB checked in by curtis, 9 years ago (diff)

Fix capitalization of Ecf.h when including.

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