add a few scripts to clean the files and indent using emacs.
[hkl.git] / hkl / hkl-pseudoaxis-auto.c
blobde75bcf52bc98bb96d4c97628d801e3d3265292a
1 /* This file is part of the hkl library.
3 * The hkl library is free software: you can redistribute it and/or modify
4 * it under the terms of the GNU General Public License as published by
5 * the Free Software Foundation, either version 3 of the License, or
6 * (at your option) any later version.
8 * The hkl library is distributed in the hope that it will be useful,
9 * but WITHOUT ANY WARRANTY; without even the implied warranty of
10 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
11 * GNU General Public License for more details.
13 * You should have received a copy of the GNU General Public License
14 * along with the hkl library. If not, see <http://www.gnu.org/licenses/>.
16 * Copyright (C) 2003-2010 Synchrotron SOLEIL
17 * L'Orme des Merisiers Saint-Aubin
18 * BP 48 91192 GIF-sur-YVETTE CEDEX
20 * Authors: Picca Frédéric-Emmanuel <picca@synchrotron-soleil.fr>
22 #include <string.h>
23 #ifndef _MSC_VER
24 # include <alloca.h>
25 #endif
26 #include <gsl/gsl_sf_trig.h>
27 #include <hkl/hkl-pseudoaxis-auto.h>
29 /* #define DEBUG */
31 /*********************************************/
32 /* methods use to solve numerical pseudoAxes */
33 /*********************************************/
35 /**
36 * @brief This private method find the degenerated axes.
38 * @param func the gsl_multiroopt_function to test
39 * @param x the starting point
40 * @param f the result of the function evaluation.
42 * with this method we can see if an axis is degenerated or not.
43 * A degenerated axis is an axis with no effect on the function evaluation.
44 * In the Jacobian matrix all elements of a columnn is null.
45 * Once we know this the axis is mark as degenerated and we do not need to
46 * change is sector.
48 static void find_degenerated_axes(HklPseudoAxisEngine *self,
49 gsl_multiroot_function *func,
50 gsl_vector const *x, gsl_vector const *f,
51 int degenerated[])
53 gsl_matrix *J;
54 size_t i, j;
56 memset(degenerated, 0, x->size * sizeof(int));
57 J = gsl_matrix_alloc(x->size, f->size);
59 gsl_multiroot_fdjacobian(func, x, f, GSL_SQRT_DBL_EPSILON, J);
60 for(j=0; j<x->size && !degenerated[j]; ++j) {
61 for(i=0; i<f->size; ++i)
62 if (fabs(gsl_matrix_get(J, i, j)) > HKL_EPSILON)
63 break;
64 if (i == f->size)
65 degenerated[j] = 1;
68 #ifdef DEBUG
69 fprintf(stdout, "\nLooks for degenerated axes\n");
70 for(i=0; i<x->size; ++i)
71 fprintf(stdout, " %d", degenerated[i]);
72 for(i=0;i<x->size;++i) {
73 fprintf(stdout, "\n ");
74 for(j=0;j<f->size;++j)
75 fprintf(stdout, " %f", gsl_matrix_get(J, i, j));
77 fprintf(stdout, "\n");
78 #endif
79 gsl_matrix_free(J);
82 /**
83 * @brief this private method try to find the first solution
85 * @param self the current HklPseudoAxeEngine.
86 * @param f The function to use for the computation.
88 * If a solution was found it also check for degenerated axes.
89 * A degenerated axes is an Axes with no effect on the function.
90 * @see find_degenerated
91 * @return HKL_SUCCESS (0) or HKL_FAIL (-1).
93 static int find_first_geometry(HklPseudoAxisEngine *self,
94 gsl_multiroot_function *f,
95 int degenerated[])
97 gsl_multiroot_fsolver_type const *T;
98 gsl_multiroot_fsolver *s;
99 gsl_vector *x;
100 size_t len = HKL_LIST_LEN(self->axes);
101 double *x_data;
102 double *x_data0 = alloca(len * sizeof(*x_data0));
103 size_t iter = 0;
104 int status;
105 int res = HKL_FAIL;
106 size_t i;
108 /* get the starting point from the geometry */
109 /* must be put in the auto_set method */
110 x = gsl_vector_alloc(len);
111 x_data = (double *)x->data;
112 for(i=0; i<len; ++i)
113 x_data[i] = hkl_axis_get_value(self->axes[i]);
115 /* keep a copy of the first axes positions to deal with degenerated axes */
116 memcpy(x_data0, x_data, len * sizeof(double));
118 /* Initialize method */
119 T = gsl_multiroot_fsolver_hybrid;
120 s = gsl_multiroot_fsolver_alloc (T, len);
121 gsl_multiroot_fsolver_set (s, f, x);
123 /* iterate to find the solution */
124 do {
125 ++iter;
126 status = gsl_multiroot_fsolver_iterate(s);
127 if (status || iter % 100 == 0) {
128 /* Restart from another point. */
129 for(i=0; i<len; ++i)
130 x_data[i] = (double)rand() / RAND_MAX * 180. / M_PI;
131 gsl_multiroot_fsolver_set(s, f, x);
132 gsl_multiroot_fsolver_iterate(s);
134 status = gsl_multiroot_test_residual (s->f, HKL_EPSILON);
135 } while (status == GSL_CONTINUE && iter < 1000);
137 #ifdef DEBUG
138 fprintf(stdout, "\nstatus : %d iter : %d", status, iter);
139 for(i=0; i<len; ++i)
140 fprintf(stdout, " %.7f", s->f->data[i]);
141 fprintf(stdout, "\n");
142 #endif
144 if (status != GSL_CONTINUE) {
145 find_degenerated_axes(self, f, s->x, s->f, degenerated);
147 #ifdef DEBUG
148 /* print the test header */
149 fprintf(stdout, "\n");
150 for(i=0; i<len; ++i)
151 fprintf(stdout, "\t f(%d)", i);
152 for(i=0; i<len; ++i)
153 fprintf(stdout, "\t \"%s\"", ((HklParameter *)self->axes[i])->name);
154 #endif
155 /* set the geometry from the gsl_vector */
156 /* in a futur version the geometry must contain a gsl_vector */
157 /* to avoid this. */
158 x_data = (double *)s->x->data;
159 for(i=0; i<len; ++i)
160 if (degenerated[i])
161 hkl_axis_set_value(self->axes[i], x_data0[i]);
162 else
163 hkl_axis_set_value(self->axes[i], x_data[i]);
165 hkl_geometry_update(self->geometry);
166 res = HKL_SUCCESS;
169 /* release memory */
170 gsl_vector_free(x);
171 gsl_multiroot_fsolver_free(s);
173 return res;
177 * @brief This private method change the sector of angles.
179 * @param x The vector of changed angles.
180 * @param x0 The vector of angles to change.
181 * @param sector the sector vector operation.
182 * @param n the size of all vectors.
184 * 0 -> no change
185 * 1 -> pi - angle
186 * 2 -> pi + angle
187 * 3 -> -angle
189 static void change_sector(double x[], double const x0[],
190 int const sector[], size_t n)
192 size_t i;
194 for(i=0; i<n; ++i) {
195 switch (sector[i]) {
196 case 0:
197 x[i] = x0[i];
198 break;
199 case 1:
200 x[i] = M_PI - x0[i];
201 break;
202 case 2:
203 x[i] = M_PI + x0[i];
204 break;
205 case 3:
206 x[i] = -x0[i];
207 break;
213 * @brief Test if an angle combination is compatible with q function.
215 * @param x The vector of angles to test.
216 * @param function The gsl_multiroot_function used for the test.
217 * @param f a gsl_vector use to compute the result (optimization)
219 static int test_sector(gsl_vector const *x,
220 gsl_multiroot_function *function,
221 gsl_vector *f)
223 int res = HKL_SUCCESS;
224 size_t i;
225 double *f_data = f->data;
227 function->f(x, function->params, f);
229 for(i=0; i<f->size; ++i)
230 if (fabs(f_data[i]) > HKL_EPSILON){
231 res = HKL_FAIL;
232 break;
235 #ifdef DEBUG
236 fprintf(stdout, "\n");
237 for(i=0; i<f->size; ++i)
238 if(fabs(f_data[i]) < HKL_EPSILON)
239 fprintf(stdout, "\t%f *", f_data[i]);
240 else
241 fprintf(stdout, "\t%f", f_data[i]);
242 for(i=0; i<f->size; ++i)
243 fprintf(stdout, "\t%f", gsl_sf_angle_restrict_symm(x->data[i]) * HKL_RADTODEG);
245 if(res == HKL_FAIL)
246 fprintf(stdout, "\t FAIL");
247 else
248 fprintf(stdout, "\t SUCCESS");
249 #endif
251 return res;
255 * @brief compute the permutation and test its validity.
257 * @param axes_len number of axes
258 * @param op_len number of operation per axes. (4 for now)
259 * @param p The vector containing the current permutation.
260 * @param axes_idx The index of the axes we are permution.
261 * @param op the current operation to set.
262 * @param f The function for the validity test.
263 * @param x0 The starting point of all geometry permutations.
264 * @param _x a gsl_vector use to compute the sectors (optimization)
265 * @param _f a gsl_vector use during the sector test (optimization)
267 static void perm_r(size_t axes_len, size_t op_len[], int p[], size_t axes_idx,
268 int op, gsl_multiroot_function *f, double x0[],
269 gsl_vector *_x, gsl_vector *_f)
271 size_t i;
273 p[axes_idx++] = op;
274 if (axes_idx == axes_len) {
275 double *x_data = _x->data;
276 change_sector(x_data, x0, p, axes_len);
277 if (HKL_SUCCESS == test_sector(_x, f, _f))
278 hkl_pseudo_axis_engine_add_geometry(f->params, x_data);
279 } else
280 for (i=0; i<op_len[axes_idx]; ++i)
281 perm_r(axes_len, op_len, p, axes_idx, i, f, x0, _x, _f);
285 * @brief Find all numerical solutions of a mode.
287 * @param self the current HklPseudoAxisEngine
288 * @param function The mode function
290 * @return HKL_SUCCESS (0) or HKL_FAIL (-1)
292 * This method find a first solution with a numerical method from the
293 * GSL library (the multi root solver hybrid). Then it multiplicates the
294 * solutions from this starting point using cosinus/sinus properties.
295 * It addes all valid solutions to the self->geometries.
297 static int solve_function(HklPseudoAxisEngine *self,
298 HklFunction function)
301 size_t i;
302 size_t len = HKL_LIST_LEN(self->axes);
303 int *p = alloca(len * sizeof(*p));
304 double *x0 = alloca(len * sizeof(*x0));
305 int *degenerated = alloca(len * sizeof(*degenerated));
306 size_t *op_len = alloca(len * sizeof(*op_len));
307 int res;
308 gsl_vector *_x; /* use to compute sectors in perm_r (avoid copy) */
309 gsl_vector *_f; /* use to test sectors in perm_r (avoid copy) */
310 gsl_multiroot_function f;
312 _x = gsl_vector_alloc(len);
313 _f = gsl_vector_alloc(len);
315 f.f = function;
316 f.n = len;
317 f.params = self;
319 res = find_first_geometry(self, &f, degenerated);
320 if (!res) {
321 memset(p, 0, sizeof(p));
322 /* use first solution as starting point for permutations */
323 for(i=0; i<len; ++i){
324 x0[i] = hkl_axis_get_value(self->axes[i]);
325 if (degenerated[i])
326 op_len[i] = 1;
327 else
328 op_len[i] = 4;
330 for (i=0; i<op_len[0]; ++i)
331 perm_r(len, op_len, p, 0, i, &f, x0, _x, _f);
334 gsl_vector_free(_f);
335 gsl_vector_free(_x);
336 return res;
339 int hkl_pseudo_axis_engine_mode_set_real(HklPseudoAxisEngineMode *self,
340 HklPseudoAxisEngine *engine,
341 HklGeometry *geometry,
342 HklDetector *detector,
343 HklSample *sample,
344 HklError **error)
346 size_t i;
347 int res = HKL_FAIL;
349 if(!self || !engine || !geometry || !detector || !sample)
350 return res;
352 for(i=0;i<HKL_LIST_LEN(self->functions);++i)
353 res &= solve_function(engine, self->functions[i]);
355 #ifdef DEBUG
356 hkl_pseudo_axis_engine_fprintf(stdout, engine);
357 #endif
359 return res;