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>
26 #include <gsl/gsl_sf_trig.h>
27 #include <hkl/hkl-pseudoaxis-auto.h>
31 /*********************************************/
32 /* methods use to solve numerical pseudoAxes */
33 /*********************************************/
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
48 static void find_degenerated_axes(HklPseudoAxisEngine
*self
,
49 gsl_multiroot_function
*func
,
50 gsl_vector
const *x
, gsl_vector
const *f
,
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
)
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");
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
,
97 gsl_multiroot_fsolver_type
const *T
;
98 gsl_multiroot_fsolver
*s
;
100 size_t len
= HKL_LIST_LEN(self
->axes
);
102 double *x_data0
= alloca(len
* sizeof(*x_data0
));
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
;
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 */
126 status
= gsl_multiroot_fsolver_iterate(s
);
127 if (status
|| iter
% 100 == 0) {
128 /* Restart from another point. */
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);
138 fprintf(stdout
, "\nstatus : %d iter : %d", status
, iter
);
140 fprintf(stdout
, " %.7f", s
->f
->data
[i
]);
141 fprintf(stdout
, "\n");
144 if (status
!= GSL_CONTINUE
) {
145 find_degenerated_axes(self
, f
, s
->x
, s
->f
, degenerated
);
148 /* print the test header */
149 fprintf(stdout
, "\n");
151 fprintf(stdout
, "\t f(%d)", i
);
153 fprintf(stdout
, "\t \"%s\"", ((HklParameter
*)self
->axes
[i
])->name
);
155 /* set the geometry from the gsl_vector */
156 /* in a futur version the geometry must contain a gsl_vector */
158 x_data
= (double *)s
->x
->data
;
161 hkl_axis_set_value(self
->axes
[i
], x_data0
[i
]);
163 hkl_axis_set_value(self
->axes
[i
], x_data
[i
]);
165 hkl_geometry_update(self
->geometry
);
171 gsl_multiroot_fsolver_free(s
);
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.
189 static void change_sector(double x
[], double const x0
[],
190 int const sector
[], size_t n
)
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
,
223 int res
= HKL_SUCCESS
;
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
){
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
]);
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
);
246 fprintf(stdout
, "\t FAIL");
248 fprintf(stdout
, "\t SUCCESS");
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
)
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
);
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
)
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
));
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
);
319 res
= find_first_geometry(self
, &f
, degenerated
);
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
]);
330 for (i
=0; i
<op_len
[0]; ++i
)
331 perm_r(len
, op_len
, p
, 0, i
, &f
, x0
, _x
, _f
);
339 int hkl_pseudo_axis_engine_mode_set_real(HklPseudoAxisEngineMode
*self
,
340 HklPseudoAxisEngine
*engine
,
341 HklGeometry
*geometry
,
342 HklDetector
*detector
,
349 if(!self
|| !engine
|| !geometry
|| !detector
|| !sample
)
352 for(i
=0;i
<HKL_LIST_LEN(self
->functions
);++i
)
353 res
&= solve_function(engine
, self
->functions
[i
]);
356 hkl_pseudo_axis_engine_fprintf(stdout
, engine
);