diff options
author | Alberto Ruiz <aruiz@um.es> | 2014-05-21 10:30:55 +0200 |
---|---|---|
committer | Alberto Ruiz <aruiz@um.es> | 2014-05-21 10:30:55 +0200 |
commit | 197e88c3b56d28840217010a2871c6ea3a4dd1a4 (patch) | |
tree | 825be9d6c9d87d23f7e5497c0133d11d52c63535 /packages/gsl/src/Numeric/GSL/gsl-aux.c | |
parent | e07c3dee7235496b71a89233106d93f6cc94ada1 (diff) |
update dependencies, move examples etc
Diffstat (limited to 'packages/gsl/src/Numeric/GSL/gsl-aux.c')
-rw-r--r-- | packages/gsl/src/Numeric/GSL/gsl-aux.c | 945 |
1 files changed, 945 insertions, 0 deletions
diff --git a/packages/gsl/src/Numeric/GSL/gsl-aux.c b/packages/gsl/src/Numeric/GSL/gsl-aux.c new file mode 100644 index 0000000..ffc5c20 --- /dev/null +++ b/packages/gsl/src/Numeric/GSL/gsl-aux.c | |||
@@ -0,0 +1,945 @@ | |||
1 | #include <gsl/gsl_complex.h> | ||
2 | |||
3 | #define RVEC(A) int A##n, double*A##p | ||
4 | #define RMAT(A) int A##r, int A##c, double* A##p | ||
5 | #define KRVEC(A) int A##n, const double*A##p | ||
6 | #define KRMAT(A) int A##r, int A##c, const double* A##p | ||
7 | |||
8 | #define CVEC(A) int A##n, gsl_complex*A##p | ||
9 | #define CMAT(A) int A##r, int A##c, gsl_complex* A##p | ||
10 | #define KCVEC(A) int A##n, const gsl_complex*A##p | ||
11 | #define KCMAT(A) int A##r, int A##c, const gsl_complex* A##p | ||
12 | |||
13 | #define FVEC(A) int A##n, float*A##p | ||
14 | #define FMAT(A) int A##r, int A##c, float* A##p | ||
15 | #define KFVEC(A) int A##n, const float*A##p | ||
16 | #define KFMAT(A) int A##r, int A##c, const float* A##p | ||
17 | |||
18 | #define QVEC(A) int A##n, gsl_complex_float*A##p | ||
19 | #define QMAT(A) int A##r, int A##c, gsl_complex_float* A##p | ||
20 | #define KQVEC(A) int A##n, const gsl_complex_float*A##p | ||
21 | #define KQMAT(A) int A##r, int A##c, const gsl_complex_float* A##p | ||
22 | |||
23 | #include <gsl/gsl_blas.h> | ||
24 | #include <gsl/gsl_math.h> | ||
25 | #include <gsl/gsl_errno.h> | ||
26 | #include <gsl/gsl_fft_complex.h> | ||
27 | #include <gsl/gsl_integration.h> | ||
28 | #include <gsl/gsl_deriv.h> | ||
29 | #include <gsl/gsl_poly.h> | ||
30 | #include <gsl/gsl_multimin.h> | ||
31 | #include <gsl/gsl_multiroots.h> | ||
32 | #include <gsl/gsl_min.h> | ||
33 | #include <gsl/gsl_complex_math.h> | ||
34 | #include <gsl/gsl_rng.h> | ||
35 | #include <gsl/gsl_randist.h> | ||
36 | #include <gsl/gsl_roots.h> | ||
37 | #include <gsl/gsl_multifit_nlin.h> | ||
38 | #include <string.h> | ||
39 | #include <stdio.h> | ||
40 | |||
41 | #define MACRO(B) do {B} while (0) | ||
42 | #define ERROR(CODE) MACRO(return CODE;) | ||
43 | #define REQUIRES(COND, CODE) MACRO(if(!(COND)) {ERROR(CODE);}) | ||
44 | #define OK return 0; | ||
45 | |||
46 | #define MIN(A,B) ((A)<(B)?(A):(B)) | ||
47 | #define MAX(A,B) ((A)>(B)?(A):(B)) | ||
48 | |||
49 | #ifdef DBG | ||
50 | #define DEBUGMSG(M) printf("*** calling aux C function: %s\n",M); | ||
51 | #else | ||
52 | #define DEBUGMSG(M) | ||
53 | #endif | ||
54 | |||
55 | #define CHECK(RES,CODE) MACRO(if(RES) return CODE;) | ||
56 | |||
57 | #ifdef DBG | ||
58 | #define DEBUGMAT(MSG,X) printf(MSG" = \n"); gsl_matrix_fprintf(stdout,X,"%f"); printf("\n"); | ||
59 | #else | ||
60 | #define DEBUGMAT(MSG,X) | ||
61 | #endif | ||
62 | |||
63 | #ifdef DBG | ||
64 | #define DEBUGVEC(MSG,X) printf(MSG" = \n"); gsl_vector_fprintf(stdout,X,"%f"); printf("\n"); | ||
65 | #else | ||
66 | #define DEBUGVEC(MSG,X) | ||
67 | #endif | ||
68 | |||
69 | #define DVVIEW(A) gsl_vector_view A = gsl_vector_view_array(A##p,A##n) | ||
70 | #define DMVIEW(A) gsl_matrix_view A = gsl_matrix_view_array(A##p,A##r,A##c) | ||
71 | #define CVVIEW(A) gsl_vector_complex_view A = gsl_vector_complex_view_array((double*)A##p,A##n) | ||
72 | #define CMVIEW(A) gsl_matrix_complex_view A = gsl_matrix_complex_view_array((double*)A##p,A##r,A##c) | ||
73 | #define KDVVIEW(A) gsl_vector_const_view A = gsl_vector_const_view_array(A##p,A##n) | ||
74 | #define KDMVIEW(A) gsl_matrix_const_view A = gsl_matrix_const_view_array(A##p,A##r,A##c) | ||
75 | #define KCVVIEW(A) gsl_vector_complex_const_view A = gsl_vector_complex_const_view_array((double*)A##p,A##n) | ||
76 | #define KCMVIEW(A) gsl_matrix_complex_const_view A = gsl_matrix_complex_const_view_array((double*)A##p,A##r,A##c) | ||
77 | |||
78 | #define FVVIEW(A) gsl_vector_float_view A = gsl_vector_float_view_array(A##p,A##n) | ||
79 | #define FMVIEW(A) gsl_matrix_float_view A = gsl_matrix_float_view_array(A##p,A##r,A##c) | ||
80 | #define QVVIEW(A) gsl_vector_complex_float_view A = gsl_vector_float_complex_view_array((float*)A##p,A##n) | ||
81 | #define QMVIEW(A) gsl_matrix_complex_float_view A = gsl_matrix_float_complex_view_array((float*)A##p,A##r,A##c) | ||
82 | #define KFVVIEW(A) gsl_vector_float_const_view A = gsl_vector_float_const_view_array(A##p,A##n) | ||
83 | #define KFMVIEW(A) gsl_matrix_float_const_view A = gsl_matrix_float_const_view_array(A##p,A##r,A##c) | ||
84 | #define KQVVIEW(A) gsl_vector_complex_float_const_view A = gsl_vector_complex_float_const_view_array((float*)A##p,A##n) | ||
85 | #define KQMVIEW(A) gsl_matrix_complex_float_const_view A = gsl_matrix_complex_float_const_view_array((float*)A##p,A##r,A##c) | ||
86 | |||
87 | #define V(a) (&a.vector) | ||
88 | #define M(a) (&a.matrix) | ||
89 | |||
90 | #define GCVEC(A) int A##n, gsl_complex*A##p | ||
91 | #define KGCVEC(A) int A##n, const gsl_complex*A##p | ||
92 | |||
93 | #define GQVEC(A) int A##n, gsl_complex_float*A##p | ||
94 | #define KGQVEC(A) int A##n, const gsl_complex_float*A##p | ||
95 | |||
96 | #define BAD_SIZE 2000 | ||
97 | #define BAD_CODE 2001 | ||
98 | #define MEM 2002 | ||
99 | #define BAD_FILE 2003 | ||
100 | |||
101 | |||
102 | void no_abort_on_error() { | ||
103 | gsl_set_error_handler_off(); | ||
104 | } | ||
105 | |||
106 | |||
107 | |||
108 | int fft(int code, KCVEC(X), CVEC(R)) { | ||
109 | REQUIRES(Xn == Rn,BAD_SIZE); | ||
110 | DEBUGMSG("fft"); | ||
111 | int s = Xn; | ||
112 | gsl_fft_complex_wavetable * wavetable = gsl_fft_complex_wavetable_alloc (s); | ||
113 | gsl_fft_complex_workspace * workspace = gsl_fft_complex_workspace_alloc (s); | ||
114 | gsl_vector_const_view X = gsl_vector_const_view_array((double*)Xp, 2*Xn); | ||
115 | gsl_vector_view R = gsl_vector_view_array((double*)Rp, 2*Rn); | ||
116 | gsl_blas_dcopy(&X.vector,&R.vector); | ||
117 | if(code==0) { | ||
118 | gsl_fft_complex_forward ((double*)Rp, 1, s, wavetable, workspace); | ||
119 | } else { | ||
120 | gsl_fft_complex_inverse ((double*)Rp, 1, s, wavetable, workspace); | ||
121 | } | ||
122 | gsl_fft_complex_wavetable_free (wavetable); | ||
123 | gsl_fft_complex_workspace_free (workspace); | ||
124 | OK | ||
125 | } | ||
126 | |||
127 | |||
128 | int deriv(int code, double f(double, void*), double x, double h, double * result, double * abserr) | ||
129 | { | ||
130 | gsl_function F; | ||
131 | F.function = f; | ||
132 | F.params = 0; | ||
133 | |||
134 | if(code==0) return gsl_deriv_central (&F, x, h, result, abserr); | ||
135 | |||
136 | if(code==1) return gsl_deriv_forward (&F, x, h, result, abserr); | ||
137 | |||
138 | if(code==2) return gsl_deriv_backward (&F, x, h, result, abserr); | ||
139 | |||
140 | return 0; | ||
141 | } | ||
142 | |||
143 | |||
144 | int integrate_qng(double f(double, void*), double a, double b, double aprec, double prec, | ||
145 | double *result, double*error) { | ||
146 | DEBUGMSG("integrate_qng"); | ||
147 | gsl_function F; | ||
148 | F.function = f; | ||
149 | F.params = NULL; | ||
150 | size_t neval; | ||
151 | int res = gsl_integration_qng (&F, a,b, aprec, prec, result, error, &neval); | ||
152 | CHECK(res,res); | ||
153 | OK | ||
154 | } | ||
155 | |||
156 | int integrate_qags(double f(double,void*), double a, double b, double aprec, double prec, int w, | ||
157 | double *result, double* error) { | ||
158 | DEBUGMSG("integrate_qags"); | ||
159 | gsl_integration_workspace * wk = gsl_integration_workspace_alloc (w); | ||
160 | gsl_function F; | ||
161 | F.function = f; | ||
162 | F.params = NULL; | ||
163 | int res = gsl_integration_qags (&F, a,b, aprec, prec, w,wk, result, error); | ||
164 | CHECK(res,res); | ||
165 | gsl_integration_workspace_free (wk); | ||
166 | OK | ||
167 | } | ||
168 | |||
169 | int integrate_qagi(double f(double,void*), double aprec, double prec, int w, | ||
170 | double *result, double* error) { | ||
171 | DEBUGMSG("integrate_qagi"); | ||
172 | gsl_integration_workspace * wk = gsl_integration_workspace_alloc (w); | ||
173 | gsl_function F; | ||
174 | F.function = f; | ||
175 | F.params = NULL; | ||
176 | int res = gsl_integration_qagi (&F, aprec, prec, w,wk, result, error); | ||
177 | CHECK(res,res); | ||
178 | gsl_integration_workspace_free (wk); | ||
179 | OK | ||
180 | } | ||
181 | |||
182 | |||
183 | int integrate_qagiu(double f(double,void*), double a, double aprec, double prec, int w, | ||
184 | double *result, double* error) { | ||
185 | DEBUGMSG("integrate_qagiu"); | ||
186 | gsl_integration_workspace * wk = gsl_integration_workspace_alloc (w); | ||
187 | gsl_function F; | ||
188 | F.function = f; | ||
189 | F.params = NULL; | ||
190 | int res = gsl_integration_qagiu (&F, a, aprec, prec, w,wk, result, error); | ||
191 | CHECK(res,res); | ||
192 | gsl_integration_workspace_free (wk); | ||
193 | OK | ||
194 | } | ||
195 | |||
196 | |||
197 | int integrate_qagil(double f(double,void*), double b, double aprec, double prec, int w, | ||
198 | double *result, double* error) { | ||
199 | DEBUGMSG("integrate_qagil"); | ||
200 | gsl_integration_workspace * wk = gsl_integration_workspace_alloc (w); | ||
201 | gsl_function F; | ||
202 | F.function = f; | ||
203 | F.params = NULL; | ||
204 | int res = gsl_integration_qagil (&F, b, aprec, prec, w,wk, result, error); | ||
205 | CHECK(res,res); | ||
206 | gsl_integration_workspace_free (wk); | ||
207 | OK | ||
208 | } | ||
209 | |||
210 | int integrate_cquad(double f(double,void*), double a, double b, double aprec, double prec, | ||
211 | int w, double *result, double* error, int *neval) { | ||
212 | DEBUGMSG("integrate_cquad"); | ||
213 | gsl_integration_cquad_workspace * wk = gsl_integration_cquad_workspace_alloc (w); | ||
214 | gsl_function F; | ||
215 | F.function = f; | ||
216 | F.params = NULL; | ||
217 | size_t * sneval = NULL; | ||
218 | int res = gsl_integration_cquad (&F, a, b, aprec, prec, wk, result, error, sneval); | ||
219 | *neval = *sneval; | ||
220 | CHECK(res,res); | ||
221 | gsl_integration_cquad_workspace_free (wk); | ||
222 | OK | ||
223 | } | ||
224 | |||
225 | |||
226 | int polySolve(KRVEC(a), CVEC(z)) { | ||
227 | DEBUGMSG("polySolve"); | ||
228 | REQUIRES(an>1,BAD_SIZE); | ||
229 | gsl_poly_complex_workspace * w = gsl_poly_complex_workspace_alloc (an); | ||
230 | int res = gsl_poly_complex_solve ((double*)ap, an, w, (double*)zp); | ||
231 | CHECK(res,res); | ||
232 | gsl_poly_complex_workspace_free (w); | ||
233 | OK; | ||
234 | } | ||
235 | |||
236 | int vector_fscanf(char*filename, RVEC(a)) { | ||
237 | DEBUGMSG("gsl_vector_fscanf"); | ||
238 | DVVIEW(a); | ||
239 | FILE * f = fopen(filename,"r"); | ||
240 | CHECK(!f,BAD_FILE); | ||
241 | int res = gsl_vector_fscanf(f,V(a)); | ||
242 | CHECK(res,res); | ||
243 | fclose (f); | ||
244 | OK | ||
245 | } | ||
246 | |||
247 | int vector_fprintf(char*filename, char*fmt, RVEC(a)) { | ||
248 | DEBUGMSG("gsl_vector_fprintf"); | ||
249 | DVVIEW(a); | ||
250 | FILE * f = fopen(filename,"w"); | ||
251 | CHECK(!f,BAD_FILE); | ||
252 | int res = gsl_vector_fprintf(f,V(a),fmt); | ||
253 | CHECK(res,res); | ||
254 | fclose (f); | ||
255 | OK | ||
256 | } | ||
257 | |||
258 | int vector_fread(char*filename, RVEC(a)) { | ||
259 | DEBUGMSG("gsl_vector_fread"); | ||
260 | DVVIEW(a); | ||
261 | FILE * f = fopen(filename,"r"); | ||
262 | CHECK(!f,BAD_FILE); | ||
263 | int res = gsl_vector_fread(f,V(a)); | ||
264 | CHECK(res,res); | ||
265 | fclose (f); | ||
266 | OK | ||
267 | } | ||
268 | |||
269 | int vector_fwrite(char*filename, RVEC(a)) { | ||
270 | DEBUGMSG("gsl_vector_fwrite"); | ||
271 | DVVIEW(a); | ||
272 | FILE * f = fopen(filename,"w"); | ||
273 | CHECK(!f,BAD_FILE); | ||
274 | int res = gsl_vector_fwrite(f,V(a)); | ||
275 | CHECK(res,res); | ||
276 | fclose (f); | ||
277 | OK | ||
278 | } | ||
279 | |||
280 | int matrix_fprintf(char*filename, char*fmt, int ro, RMAT(m)) { | ||
281 | DEBUGMSG("matrix_fprintf"); | ||
282 | FILE * f = fopen(filename,"w"); | ||
283 | CHECK(!f,BAD_FILE); | ||
284 | int i,j,sr,sc; | ||
285 | if (ro==1) { sr = mc; sc = 1;} else { sr = 1; sc = mr;} | ||
286 | #define AT(M,r,c) (M##p[(r)*sr+(c)*sc]) | ||
287 | for (i=0; i<mr; i++) { | ||
288 | for (j=0; j<mc-1; j++) { | ||
289 | fprintf(f,fmt,AT(m,i,j)); | ||
290 | fprintf(f," "); | ||
291 | } | ||
292 | fprintf(f,fmt,AT(m,i,j)); | ||
293 | fprintf(f,"\n"); | ||
294 | } | ||
295 | fclose (f); | ||
296 | OK | ||
297 | } | ||
298 | |||
299 | //--------------------------------------------------------------- | ||
300 | |||
301 | typedef double Trawfun(int, double*); | ||
302 | |||
303 | double only_f_aux_min(const gsl_vector*x, void *pars) { | ||
304 | Trawfun * f = (Trawfun*) pars; | ||
305 | double* p = (double*)calloc(x->size,sizeof(double)); | ||
306 | int k; | ||
307 | for(k=0;k<x->size;k++) { | ||
308 | p[k] = gsl_vector_get(x,k); | ||
309 | } | ||
310 | double res = f(x->size,p); | ||
311 | free(p); | ||
312 | return res; | ||
313 | } | ||
314 | |||
315 | double only_f_aux_root(double x, void *pars); | ||
316 | int uniMinimize(int method, double f(double), | ||
317 | double epsrel, int maxit, double min, | ||
318 | double xl, double xu, RMAT(sol)) { | ||
319 | REQUIRES(solr == maxit && solc == 4,BAD_SIZE); | ||
320 | DEBUGMSG("minimize_only_f"); | ||
321 | gsl_function my_func; | ||
322 | my_func.function = only_f_aux_root; | ||
323 | my_func.params = f; | ||
324 | size_t iter = 0; | ||
325 | int status; | ||
326 | const gsl_min_fminimizer_type *T; | ||
327 | gsl_min_fminimizer *s; | ||
328 | // Starting point | ||
329 | switch(method) { | ||
330 | case 0 : {T = gsl_min_fminimizer_goldensection; break; } | ||
331 | case 1 : {T = gsl_min_fminimizer_brent; break; } | ||
332 | case 2 : {T = gsl_min_fminimizer_quad_golden; break; } | ||
333 | default: ERROR(BAD_CODE); | ||
334 | } | ||
335 | s = gsl_min_fminimizer_alloc (T); | ||
336 | gsl_min_fminimizer_set (s, &my_func, min, xl, xu); | ||
337 | do { | ||
338 | double current_min, current_lo, current_hi; | ||
339 | status = gsl_min_fminimizer_iterate (s); | ||
340 | current_min = gsl_min_fminimizer_x_minimum (s); | ||
341 | current_lo = gsl_min_fminimizer_x_lower (s); | ||
342 | current_hi = gsl_min_fminimizer_x_upper (s); | ||
343 | solp[iter*solc] = iter + 1; | ||
344 | solp[iter*solc+1] = current_min; | ||
345 | solp[iter*solc+2] = current_lo; | ||
346 | solp[iter*solc+3] = current_hi; | ||
347 | iter++; | ||
348 | if (status) /* check if solver is stuck */ | ||
349 | break; | ||
350 | |||
351 | status = | ||
352 | gsl_min_test_interval (current_lo, current_hi, 0, epsrel); | ||
353 | } | ||
354 | while (status == GSL_CONTINUE && iter < maxit); | ||
355 | int i; | ||
356 | for (i=iter; i<solr; i++) { | ||
357 | solp[i*solc+0] = iter; | ||
358 | solp[i*solc+1]=0.; | ||
359 | solp[i*solc+2]=0.; | ||
360 | solp[i*solc+3]=0.; | ||
361 | } | ||
362 | gsl_min_fminimizer_free(s); | ||
363 | OK | ||
364 | } | ||
365 | |||
366 | |||
367 | |||
368 | // this version returns info about intermediate steps | ||
369 | int minimize(int method, double f(int, double*), double tolsize, int maxit, | ||
370 | KRVEC(xi), KRVEC(sz), RMAT(sol)) { | ||
371 | REQUIRES(xin==szn && solr == maxit && solc == 3+xin,BAD_SIZE); | ||
372 | DEBUGMSG("minimizeList (nmsimplex)"); | ||
373 | gsl_multimin_function my_func; | ||
374 | // extract function from pars | ||
375 | my_func.f = only_f_aux_min; | ||
376 | my_func.n = xin; | ||
377 | my_func.params = f; | ||
378 | size_t iter = 0; | ||
379 | int status; | ||
380 | double size; | ||
381 | const gsl_multimin_fminimizer_type *T; | ||
382 | gsl_multimin_fminimizer *s = NULL; | ||
383 | // Initial vertex size vector | ||
384 | KDVVIEW(sz); | ||
385 | // Starting point | ||
386 | KDVVIEW(xi); | ||
387 | // Minimizer nmsimplex, without derivatives | ||
388 | switch(method) { | ||
389 | case 0 : {T = gsl_multimin_fminimizer_nmsimplex; break; } | ||
390 | #ifdef GSL110 | ||
391 | case 1 : {T = gsl_multimin_fminimizer_nmsimplex; break; } | ||
392 | #else | ||
393 | case 1 : {T = gsl_multimin_fminimizer_nmsimplex2; break; } | ||
394 | #endif | ||
395 | default: ERROR(BAD_CODE); | ||
396 | } | ||
397 | s = gsl_multimin_fminimizer_alloc (T, my_func.n); | ||
398 | gsl_multimin_fminimizer_set (s, &my_func, V(xi), V(sz)); | ||
399 | do { | ||
400 | status = gsl_multimin_fminimizer_iterate (s); | ||
401 | size = gsl_multimin_fminimizer_size (s); | ||
402 | |||
403 | solp[iter*solc+0] = iter+1; | ||
404 | solp[iter*solc+1] = s->fval; | ||
405 | solp[iter*solc+2] = size; | ||
406 | |||
407 | int k; | ||
408 | for(k=0;k<xin;k++) { | ||
409 | solp[iter*solc+k+3] = gsl_vector_get(s->x,k); | ||
410 | } | ||
411 | iter++; | ||
412 | if (status) break; | ||
413 | status = gsl_multimin_test_size (size, tolsize); | ||
414 | } while (status == GSL_CONTINUE && iter < maxit); | ||
415 | int i,j; | ||
416 | for (i=iter; i<solr; i++) { | ||
417 | solp[i*solc+0] = iter; | ||
418 | for(j=1;j<solc;j++) { | ||
419 | solp[i*solc+j]=0.; | ||
420 | } | ||
421 | } | ||
422 | gsl_multimin_fminimizer_free(s); | ||
423 | OK | ||
424 | } | ||
425 | |||
426 | // working with the gradient | ||
427 | |||
428 | typedef struct {double (*f)(int, double*); int (*df)(int, double*, int, double*);} Tfdf; | ||
429 | |||
430 | double f_aux_min(const gsl_vector*x, void *pars) { | ||
431 | Tfdf * fdf = ((Tfdf*) pars); | ||
432 | double* p = (double*)calloc(x->size,sizeof(double)); | ||
433 | int k; | ||
434 | for(k=0;k<x->size;k++) { | ||
435 | p[k] = gsl_vector_get(x,k); | ||
436 | } | ||
437 | double res = fdf->f(x->size,p); | ||
438 | free(p); | ||
439 | return res; | ||
440 | } | ||
441 | |||
442 | |||
443 | void df_aux_min(const gsl_vector * x, void * pars, gsl_vector * g) { | ||
444 | Tfdf * fdf = ((Tfdf*) pars); | ||
445 | double* p = (double*)calloc(x->size,sizeof(double)); | ||
446 | double* q = (double*)calloc(g->size,sizeof(double)); | ||
447 | int k; | ||
448 | for(k=0;k<x->size;k++) { | ||
449 | p[k] = gsl_vector_get(x,k); | ||
450 | } | ||
451 | |||
452 | fdf->df(x->size,p,g->size,q); | ||
453 | |||
454 | for(k=0;k<x->size;k++) { | ||
455 | gsl_vector_set(g,k,q[k]); | ||
456 | } | ||
457 | free(p); | ||
458 | free(q); | ||
459 | } | ||
460 | |||
461 | void fdf_aux_min(const gsl_vector * x, void * pars, double * f, gsl_vector * g) { | ||
462 | *f = f_aux_min(x,pars); | ||
463 | df_aux_min(x,pars,g); | ||
464 | } | ||
465 | |||
466 | |||
467 | int minimizeD(int method, double f(int, double*), int df(int, double*, int, double*), | ||
468 | double initstep, double minimpar, double tolgrad, int maxit, | ||
469 | KRVEC(xi), RMAT(sol)) { | ||
470 | REQUIRES(solr == maxit && solc == 2+xin,BAD_SIZE); | ||
471 | DEBUGMSG("minimizeWithDeriv (conjugate_fr)"); | ||
472 | gsl_multimin_function_fdf my_func; | ||
473 | // extract function from pars | ||
474 | my_func.f = f_aux_min; | ||
475 | my_func.df = df_aux_min; | ||
476 | my_func.fdf = fdf_aux_min; | ||
477 | my_func.n = xin; | ||
478 | Tfdf stfdf; | ||
479 | stfdf.f = f; | ||
480 | stfdf.df = df; | ||
481 | my_func.params = &stfdf; | ||
482 | size_t iter = 0; | ||
483 | int status; | ||
484 | const gsl_multimin_fdfminimizer_type *T; | ||
485 | gsl_multimin_fdfminimizer *s = NULL; | ||
486 | // Starting point | ||
487 | KDVVIEW(xi); | ||
488 | // conjugate gradient fr | ||
489 | switch(method) { | ||
490 | case 0 : {T = gsl_multimin_fdfminimizer_conjugate_fr; break; } | ||
491 | case 1 : {T = gsl_multimin_fdfminimizer_conjugate_pr; break; } | ||
492 | case 2 : {T = gsl_multimin_fdfminimizer_vector_bfgs; break; } | ||
493 | case 3 : {T = gsl_multimin_fdfminimizer_vector_bfgs2; break; } | ||
494 | case 4 : {T = gsl_multimin_fdfminimizer_steepest_descent; break; } | ||
495 | default: ERROR(BAD_CODE); | ||
496 | } | ||
497 | s = gsl_multimin_fdfminimizer_alloc (T, my_func.n); | ||
498 | gsl_multimin_fdfminimizer_set (s, &my_func, V(xi), initstep, minimpar); | ||
499 | do { | ||
500 | status = gsl_multimin_fdfminimizer_iterate (s); | ||
501 | solp[iter*solc+0] = iter+1; | ||
502 | solp[iter*solc+1] = s->f; | ||
503 | int k; | ||
504 | for(k=0;k<xin;k++) { | ||
505 | solp[iter*solc+k+2] = gsl_vector_get(s->x,k); | ||
506 | } | ||
507 | iter++; | ||
508 | if (status) break; | ||
509 | status = gsl_multimin_test_gradient (s->gradient, tolgrad); | ||
510 | } while (status == GSL_CONTINUE && iter < maxit); | ||
511 | int i,j; | ||
512 | for (i=iter; i<solr; i++) { | ||
513 | solp[i*solc+0] = iter; | ||
514 | for(j=1;j<solc;j++) { | ||
515 | solp[i*solc+j]=0.; | ||
516 | } | ||
517 | } | ||
518 | gsl_multimin_fdfminimizer_free(s); | ||
519 | OK | ||
520 | } | ||
521 | |||
522 | //--------------------------------------------------------------- | ||
523 | |||
524 | double only_f_aux_root(double x, void *pars) { | ||
525 | double (*f)(double) = (double (*)(double)) pars; | ||
526 | return f(x); | ||
527 | } | ||
528 | |||
529 | int root(int method, double f(double), | ||
530 | double epsrel, int maxit, | ||
531 | double xl, double xu, RMAT(sol)) { | ||
532 | REQUIRES(solr == maxit && solc == 4,BAD_SIZE); | ||
533 | DEBUGMSG("root_only_f"); | ||
534 | gsl_function my_func; | ||
535 | // extract function from pars | ||
536 | my_func.function = only_f_aux_root; | ||
537 | my_func.params = f; | ||
538 | size_t iter = 0; | ||
539 | int status; | ||
540 | const gsl_root_fsolver_type *T; | ||
541 | gsl_root_fsolver *s; | ||
542 | // Starting point | ||
543 | switch(method) { | ||
544 | case 0 : {T = gsl_root_fsolver_bisection; printf("7\n"); break; } | ||
545 | case 1 : {T = gsl_root_fsolver_falsepos; break; } | ||
546 | case 2 : {T = gsl_root_fsolver_brent; break; } | ||
547 | default: ERROR(BAD_CODE); | ||
548 | } | ||
549 | s = gsl_root_fsolver_alloc (T); | ||
550 | gsl_root_fsolver_set (s, &my_func, xl, xu); | ||
551 | do { | ||
552 | double best, current_lo, current_hi; | ||
553 | status = gsl_root_fsolver_iterate (s); | ||
554 | best = gsl_root_fsolver_root (s); | ||
555 | current_lo = gsl_root_fsolver_x_lower (s); | ||
556 | current_hi = gsl_root_fsolver_x_upper (s); | ||
557 | solp[iter*solc] = iter + 1; | ||
558 | solp[iter*solc+1] = best; | ||
559 | solp[iter*solc+2] = current_lo; | ||
560 | solp[iter*solc+3] = current_hi; | ||
561 | iter++; | ||
562 | if (status) /* check if solver is stuck */ | ||
563 | break; | ||
564 | |||
565 | status = | ||
566 | gsl_root_test_interval (current_lo, current_hi, 0, epsrel); | ||
567 | } | ||
568 | while (status == GSL_CONTINUE && iter < maxit); | ||
569 | int i; | ||
570 | for (i=iter; i<solr; i++) { | ||
571 | solp[i*solc+0] = iter; | ||
572 | solp[i*solc+1]=0.; | ||
573 | solp[i*solc+2]=0.; | ||
574 | solp[i*solc+3]=0.; | ||
575 | } | ||
576 | gsl_root_fsolver_free(s); | ||
577 | OK | ||
578 | } | ||
579 | |||
580 | typedef struct { | ||
581 | double (*f)(double); | ||
582 | double (*jf)(double); | ||
583 | } uniTfjf; | ||
584 | |||
585 | double f_aux_uni(double x, void *pars) { | ||
586 | uniTfjf * fjf = ((uniTfjf*) pars); | ||
587 | return (fjf->f)(x); | ||
588 | } | ||
589 | |||
590 | double jf_aux_uni(double x, void * pars) { | ||
591 | uniTfjf * fjf = ((uniTfjf*) pars); | ||
592 | return (fjf->jf)(x); | ||
593 | } | ||
594 | |||
595 | void fjf_aux_uni(double x, void * pars, double * f, double * g) { | ||
596 | *f = f_aux_uni(x,pars); | ||
597 | *g = jf_aux_uni(x,pars); | ||
598 | } | ||
599 | |||
600 | int rootj(int method, double f(double), | ||
601 | double df(double), | ||
602 | double epsrel, int maxit, | ||
603 | double x, RMAT(sol)) { | ||
604 | REQUIRES(solr == maxit && solc == 2,BAD_SIZE); | ||
605 | DEBUGMSG("root_fjf"); | ||
606 | gsl_function_fdf my_func; | ||
607 | // extract function from pars | ||
608 | my_func.f = f_aux_uni; | ||
609 | my_func.df = jf_aux_uni; | ||
610 | my_func.fdf = fjf_aux_uni; | ||
611 | uniTfjf stfjf; | ||
612 | stfjf.f = f; | ||
613 | stfjf.jf = df; | ||
614 | my_func.params = &stfjf; | ||
615 | size_t iter = 0; | ||
616 | int status; | ||
617 | const gsl_root_fdfsolver_type *T; | ||
618 | gsl_root_fdfsolver *s; | ||
619 | // Starting point | ||
620 | switch(method) { | ||
621 | case 0 : {T = gsl_root_fdfsolver_newton;; break; } | ||
622 | case 1 : {T = gsl_root_fdfsolver_secant; break; } | ||
623 | case 2 : {T = gsl_root_fdfsolver_steffenson; break; } | ||
624 | default: ERROR(BAD_CODE); | ||
625 | } | ||
626 | s = gsl_root_fdfsolver_alloc (T); | ||
627 | |||
628 | gsl_root_fdfsolver_set (s, &my_func, x); | ||
629 | |||
630 | do { | ||
631 | double x0; | ||
632 | status = gsl_root_fdfsolver_iterate (s); | ||
633 | x0 = x; | ||
634 | x = gsl_root_fdfsolver_root(s); | ||
635 | solp[iter*solc+0] = iter+1; | ||
636 | solp[iter*solc+1] = x; | ||
637 | |||
638 | iter++; | ||
639 | if (status) /* check if solver is stuck */ | ||
640 | break; | ||
641 | |||
642 | status = | ||
643 | gsl_root_test_delta (x, x0, 0, epsrel); | ||
644 | } | ||
645 | while (status == GSL_CONTINUE && iter < maxit); | ||
646 | |||
647 | int i; | ||
648 | for (i=iter; i<solr; i++) { | ||
649 | solp[i*solc+0] = iter; | ||
650 | solp[i*solc+1]=0.; | ||
651 | } | ||
652 | gsl_root_fdfsolver_free(s); | ||
653 | OK | ||
654 | } | ||
655 | |||
656 | |||
657 | //--------------------------------------------------------------- | ||
658 | |||
659 | typedef void TrawfunV(int, double*, int, double*); | ||
660 | |||
661 | int only_f_aux_multiroot(const gsl_vector*x, void *pars, gsl_vector*y) { | ||
662 | TrawfunV * f = (TrawfunV*) pars; | ||
663 | double* p = (double*)calloc(x->size,sizeof(double)); | ||
664 | double* q = (double*)calloc(y->size,sizeof(double)); | ||
665 | int k; | ||
666 | for(k=0;k<x->size;k++) { | ||
667 | p[k] = gsl_vector_get(x,k); | ||
668 | } | ||
669 | f(x->size,p,y->size,q); | ||
670 | for(k=0;k<y->size;k++) { | ||
671 | gsl_vector_set(y,k,q[k]); | ||
672 | } | ||
673 | free(p); | ||
674 | free(q); | ||
675 | return 0; //hmmm | ||
676 | } | ||
677 | |||
678 | int multiroot(int method, void f(int, double*, int, double*), | ||
679 | double epsabs, int maxit, | ||
680 | KRVEC(xi), RMAT(sol)) { | ||
681 | REQUIRES(solr == maxit && solc == 1+2*xin,BAD_SIZE); | ||
682 | DEBUGMSG("root_only_f"); | ||
683 | gsl_multiroot_function my_func; | ||
684 | // extract function from pars | ||
685 | my_func.f = only_f_aux_multiroot; | ||
686 | my_func.n = xin; | ||
687 | my_func.params = f; | ||
688 | size_t iter = 0; | ||
689 | int status; | ||
690 | const gsl_multiroot_fsolver_type *T; | ||
691 | gsl_multiroot_fsolver *s; | ||
692 | // Starting point | ||
693 | KDVVIEW(xi); | ||
694 | switch(method) { | ||
695 | case 0 : {T = gsl_multiroot_fsolver_hybrids;; break; } | ||
696 | case 1 : {T = gsl_multiroot_fsolver_hybrid; break; } | ||
697 | case 2 : {T = gsl_multiroot_fsolver_dnewton; break; } | ||
698 | case 3 : {T = gsl_multiroot_fsolver_broyden; break; } | ||
699 | default: ERROR(BAD_CODE); | ||
700 | } | ||
701 | s = gsl_multiroot_fsolver_alloc (T, my_func.n); | ||
702 | gsl_multiroot_fsolver_set (s, &my_func, V(xi)); | ||
703 | |||
704 | do { | ||
705 | status = gsl_multiroot_fsolver_iterate (s); | ||
706 | |||
707 | solp[iter*solc+0] = iter+1; | ||
708 | |||
709 | int k; | ||
710 | for(k=0;k<xin;k++) { | ||
711 | solp[iter*solc+k+1] = gsl_vector_get(s->x,k); | ||
712 | } | ||
713 | for(k=xin;k<2*xin;k++) { | ||
714 | solp[iter*solc+k+1] = gsl_vector_get(s->f,k-xin); | ||
715 | } | ||
716 | |||
717 | iter++; | ||
718 | if (status) /* check if solver is stuck */ | ||
719 | break; | ||
720 | |||
721 | status = | ||
722 | gsl_multiroot_test_residual (s->f, epsabs); | ||
723 | } | ||
724 | while (status == GSL_CONTINUE && iter < maxit); | ||
725 | |||
726 | int i,j; | ||
727 | for (i=iter; i<solr; i++) { | ||
728 | solp[i*solc+0] = iter; | ||
729 | for(j=1;j<solc;j++) { | ||
730 | solp[i*solc+j]=0.; | ||
731 | } | ||
732 | } | ||
733 | gsl_multiroot_fsolver_free(s); | ||
734 | OK | ||
735 | } | ||
736 | |||
737 | // working with the jacobian | ||
738 | |||
739 | typedef struct {int (*f)(int, double*, int, double *); | ||
740 | int (*jf)(int, double*, int, int, double*);} Tfjf; | ||
741 | |||
742 | int f_aux(const gsl_vector*x, void *pars, gsl_vector*y) { | ||
743 | Tfjf * fjf = ((Tfjf*) pars); | ||
744 | double* p = (double*)calloc(x->size,sizeof(double)); | ||
745 | double* q = (double*)calloc(y->size,sizeof(double)); | ||
746 | int k; | ||
747 | for(k=0;k<x->size;k++) { | ||
748 | p[k] = gsl_vector_get(x,k); | ||
749 | } | ||
750 | (fjf->f)(x->size,p,y->size,q); | ||
751 | for(k=0;k<y->size;k++) { | ||
752 | gsl_vector_set(y,k,q[k]); | ||
753 | } | ||
754 | free(p); | ||
755 | free(q); | ||
756 | return 0; | ||
757 | } | ||
758 | |||
759 | int jf_aux(const gsl_vector * x, void * pars, gsl_matrix * jac) { | ||
760 | Tfjf * fjf = ((Tfjf*) pars); | ||
761 | double* p = (double*)calloc(x->size,sizeof(double)); | ||
762 | double* q = (double*)calloc((jac->size1)*(jac->size2),sizeof(double)); | ||
763 | int i,j,k; | ||
764 | for(k=0;k<x->size;k++) { | ||
765 | p[k] = gsl_vector_get(x,k); | ||
766 | } | ||
767 | |||
768 | (fjf->jf)(x->size,p,jac->size1,jac->size2,q); | ||
769 | |||
770 | k=0; | ||
771 | for(i=0;i<jac->size1;i++) { | ||
772 | for(j=0;j<jac->size2;j++){ | ||
773 | gsl_matrix_set(jac,i,j,q[k++]); | ||
774 | } | ||
775 | } | ||
776 | free(p); | ||
777 | free(q); | ||
778 | return 0; | ||
779 | } | ||
780 | |||
781 | int fjf_aux(const gsl_vector * x, void * pars, gsl_vector * f, gsl_matrix * g) { | ||
782 | f_aux(x,pars,f); | ||
783 | jf_aux(x,pars,g); | ||
784 | return 0; | ||
785 | } | ||
786 | |||
787 | int multirootj(int method, int f(int, double*, int, double*), | ||
788 | int jac(int, double*, int, int, double*), | ||
789 | double epsabs, int maxit, | ||
790 | KRVEC(xi), RMAT(sol)) { | ||
791 | REQUIRES(solr == maxit && solc == 1+2*xin,BAD_SIZE); | ||
792 | DEBUGMSG("root_fjf"); | ||
793 | gsl_multiroot_function_fdf my_func; | ||
794 | // extract function from pars | ||
795 | my_func.f = f_aux; | ||
796 | my_func.df = jf_aux; | ||
797 | my_func.fdf = fjf_aux; | ||
798 | my_func.n = xin; | ||
799 | Tfjf stfjf; | ||
800 | stfjf.f = f; | ||
801 | stfjf.jf = jac; | ||
802 | my_func.params = &stfjf; | ||
803 | size_t iter = 0; | ||
804 | int status; | ||
805 | const gsl_multiroot_fdfsolver_type *T; | ||
806 | gsl_multiroot_fdfsolver *s; | ||
807 | // Starting point | ||
808 | KDVVIEW(xi); | ||
809 | switch(method) { | ||
810 | case 0 : {T = gsl_multiroot_fdfsolver_hybridsj;; break; } | ||
811 | case 1 : {T = gsl_multiroot_fdfsolver_hybridj; break; } | ||
812 | case 2 : {T = gsl_multiroot_fdfsolver_newton; break; } | ||
813 | case 3 : {T = gsl_multiroot_fdfsolver_gnewton; break; } | ||
814 | default: ERROR(BAD_CODE); | ||
815 | } | ||
816 | s = gsl_multiroot_fdfsolver_alloc (T, my_func.n); | ||
817 | |||
818 | gsl_multiroot_fdfsolver_set (s, &my_func, V(xi)); | ||
819 | |||
820 | do { | ||
821 | status = gsl_multiroot_fdfsolver_iterate (s); | ||
822 | |||
823 | solp[iter*solc+0] = iter+1; | ||
824 | |||
825 | int k; | ||
826 | for(k=0;k<xin;k++) { | ||
827 | solp[iter*solc+k+1] = gsl_vector_get(s->x,k); | ||
828 | } | ||
829 | for(k=xin;k<2*xin;k++) { | ||
830 | solp[iter*solc+k+1] = gsl_vector_get(s->f,k-xin); | ||
831 | } | ||
832 | |||
833 | iter++; | ||
834 | if (status) /* check if solver is stuck */ | ||
835 | break; | ||
836 | |||
837 | status = | ||
838 | gsl_multiroot_test_residual (s->f, epsabs); | ||
839 | } | ||
840 | while (status == GSL_CONTINUE && iter < maxit); | ||
841 | |||
842 | int i,j; | ||
843 | for (i=iter; i<solr; i++) { | ||
844 | solp[i*solc+0] = iter; | ||
845 | for(j=1;j<solc;j++) { | ||
846 | solp[i*solc+j]=0.; | ||
847 | } | ||
848 | } | ||
849 | gsl_multiroot_fdfsolver_free(s); | ||
850 | OK | ||
851 | } | ||
852 | |||
853 | //-------------- non linear least squares fitting ------------------- | ||
854 | |||
855 | int nlfit(int method, int f(int, double*, int, double*), | ||
856 | int jac(int, double*, int, int, double*), | ||
857 | double epsabs, double epsrel, int maxit, int p, | ||
858 | KRVEC(xi), RMAT(sol)) { | ||
859 | REQUIRES(solr == maxit && solc == 2+xin,BAD_SIZE); | ||
860 | DEBUGMSG("nlfit"); | ||
861 | const gsl_multifit_fdfsolver_type *T; | ||
862 | gsl_multifit_fdfsolver *s; | ||
863 | gsl_multifit_function_fdf my_f; | ||
864 | // extract function from pars | ||
865 | my_f.f = f_aux; | ||
866 | my_f.df = jf_aux; | ||
867 | my_f.fdf = fjf_aux; | ||
868 | my_f.n = p; | ||
869 | my_f.p = xin; // !!!! | ||
870 | Tfjf stfjf; | ||
871 | stfjf.f = f; | ||
872 | stfjf.jf = jac; | ||
873 | my_f.params = &stfjf; | ||
874 | size_t iter = 0; | ||
875 | int status; | ||
876 | |||
877 | KDVVIEW(xi); | ||
878 | //DMVIEW(cov); | ||
879 | |||
880 | switch(method) { | ||
881 | case 0 : { T = gsl_multifit_fdfsolver_lmsder; break; } | ||
882 | case 1 : { T = gsl_multifit_fdfsolver_lmder; break; } | ||
883 | default: ERROR(BAD_CODE); | ||
884 | } | ||
885 | |||
886 | s = gsl_multifit_fdfsolver_alloc (T, my_f.n, my_f.p); | ||
887 | gsl_multifit_fdfsolver_set (s, &my_f, V(xi)); | ||
888 | |||
889 | do { status = gsl_multifit_fdfsolver_iterate (s); | ||
890 | |||
891 | solp[iter*solc+0] = iter+1; | ||
892 | solp[iter*solc+1] = gsl_blas_dnrm2 (s->f); | ||
893 | |||
894 | int k; | ||
895 | for(k=0;k<xin;k++) { | ||
896 | solp[iter*solc+k+2] = gsl_vector_get(s->x,k); | ||
897 | } | ||
898 | |||
899 | iter++; | ||
900 | if (status) /* check if solver is stuck */ | ||
901 | break; | ||
902 | |||
903 | status = gsl_multifit_test_delta (s->dx, s->x, epsabs, epsrel); | ||
904 | } | ||
905 | while (status == GSL_CONTINUE && iter < maxit); | ||
906 | |||
907 | int i,j; | ||
908 | for (i=iter; i<solr; i++) { | ||
909 | solp[i*solc+0] = iter; | ||
910 | for(j=1;j<solc;j++) { | ||
911 | solp[i*solc+j]=0.; | ||
912 | } | ||
913 | } | ||
914 | |||
915 | //gsl_multifit_covar (s->J, 0.0, M(cov)); | ||
916 | |||
917 | gsl_multifit_fdfsolver_free (s); | ||
918 | OK | ||
919 | } | ||
920 | |||
921 | |||
922 | ////////////////////////////////////////////////////// | ||
923 | |||
924 | |||
925 | #define RAN(C,F) case C: { for(k=0;k<rn;k++) { rp[k]= F(gen); }; OK } | ||
926 | |||
927 | int random_vector_GSL(int seed, int code, RVEC(r)) { | ||
928 | DEBUGMSG("random_vector_GSL") | ||
929 | static gsl_rng * gen = NULL; | ||
930 | if (!gen) { gen = gsl_rng_alloc (gsl_rng_mt19937);} | ||
931 | gsl_rng_set (gen, seed); | ||
932 | int k; | ||
933 | switch (code) { | ||
934 | RAN(0,gsl_rng_uniform) | ||
935 | RAN(1,gsl_ran_ugaussian) | ||
936 | default: ERROR(BAD_CODE); | ||
937 | } | ||
938 | } | ||
939 | #undef RAN | ||
940 | |||
941 | ////////////////////////////////////////////////////// | ||
942 | |||
943 | #include "gsl-ode.c" | ||
944 | |||
945 | ////////////////////////////////////////////////////// | ||