src/hkl-pseudoaxis-auto.c

   1 #include <string.h>
   2 #include <gsl/gsl_sf_trig.h>
   3 #include <hkl/hkl-pseudoaxis-auto.h>
   4
   5 /*********************************************/
   6 /* methods use to solve numerical pseudoAxes */
   7 /*********************************************/
   8
   9 /**
  10  * @brief This private method find the degenerated axes.
  11  *
  12  * @param func the gsl_multiroopt_function to test
  13  * @param x the starting point
  14  * @param f the result of the function evaluation.
  15  *
  16  * with this method we can see if an axis is degenerated or not.
  17  * A degenerated axis is an axis with no effect on the function evaluation.
  18  * In the Jacobian matrix all elements of a columnn is null.
  19  * Once we know this the axis is mark as degenerated and we do not need to
  20  * change is sector.
  21  */
  22 static void find_degenerated_axes(HklPseudoAxisEngine *self,
  23                                   gsl_multiroot_function *func,
  24                                   gsl_vector const *x, gsl_vector const *f,
  25                                   int degenerated[])
  26 {
  27         gsl_matrix *J;
  28         size_t i, j;
  29
  30         memset(degenerated, 0, x->size * sizeof(int));
  31         J = gsl_matrix_alloc(x->size, f->size);
  32
  33         gsl_multiroot_fdjacobian(func, x, f, GSL_SQRT_DBL_EPSILON, J);
  34         for(j=0; j<x->size && !degenerated[j]; ++j) {
  35                 for(i=0; i<f->size; ++i)
  36                         if (fabs(gsl_matrix_get(J, i, j)) > HKL_EPSILON)
  37                                 break;
  38                 if (i == f->size)
  39                         degenerated[j] = 1;
  40         }
  41
  42         /*
  43           hkl_pseudoAxisEngine_fprintf(func->params, stdout);
  44           fprintf(stdout, "\n");
  45           for(i=0; i<x->size; ++i)
  46           fprintf(stdout, " %d", degenerated[i]);
  47           for(i=0;i<x->size;++i) {
  48           fprintf(stdout, "\n   ");
  49           for(j=0;j<f->size;++j)
  50           fprintf(stdout, " %f", gsl_matrix_get(J, i, j));
  51           }
  52           fprintf(stdout, "\n");
  53         */
  54         gsl_matrix_free(J);
  55 }
  56
  57 /**
  58  * @brief this private method try to find the first solution
  59  *
  60  * @param self the current HklPseudoAxeEngine.
  61  * @param f The function to use for the computation.
  62  *
  63  * If a solution was found it also check for degenerated axes.
  64  * A degenerated axes is an Axes with no effect on the function.
  65  * @see find_degenerated
  66  * @return HKL_SUCCESS (0) or HKL_FAIL (-1).
  67  */
  68 static int find_first_geometry(HklPseudoAxisEngine *self,
  69                                gsl_multiroot_function *f,
  70                                int degenerated[])
  71 {
  72         gsl_multiroot_fsolver_type const *T;
  73         gsl_multiroot_fsolver *s;
  74         gsl_vector *x;
  75         double *x_data;
  76         double x_data0[self->axes_len];
  77         size_t iter = 0;
  78         int status;
  79         int res = HKL_FAIL;
  80         size_t i;
  81
  82         // get the starting point from the geometry
  83         // must be put in the auto_set method
  84         x = gsl_vector_alloc(self->axes_len);
  85         x_data = (double *)x->data;
  86         for(i=0; i<self->axes_len; ++i)
  87                 x_data[i] = self->axes[i]->config.value;
  88
  89         // keep a copy of the first axes positions to deal with degenerated axes
  90         memcpy(x_data0, x_data, self->axes_len * sizeof(double));
  91
  92         // Initialize method
  93         T = gsl_multiroot_fsolver_hybrid;
  94         s = gsl_multiroot_fsolver_alloc (T, self->axes_len);
  95         gsl_multiroot_fsolver_set (s, f, x);
  96
  97         // iterate to find the solution
  98         do {
  99                 ++iter;
 100                 status = gsl_multiroot_fsolver_iterate(s);
 101                 if (status || iter % 1000 == 0) {
 102                         // Restart from another point.
 103                         for(i=0; i<self->axes_len; ++i)
 104                                 x_data[i] = (double)rand() / RAND_MAX * 180. / M_PI;
 105                         gsl_multiroot_fsolver_set(s, f, x);
 106                         status = gsl_multiroot_fsolver_iterate(s);
 107                 }
 108                 status = gsl_multiroot_test_residual (s->f, HKL_EPSILON);
 109         } while (status == GSL_CONTINUE && iter < 1000);
 110
 111         if (status != GSL_CONTINUE) {
 112                 find_degenerated_axes(self, f, s->x, s->f, degenerated);
 113
 114                 // set the geometry from the gsl_vector
 115                 // in a futur version the geometry must contain a gsl_vector
 116                 // to avoid this.
 117                 x_data = (double *)s->x->data;
 118                 for(i=0; i<self->axes_len; ++i) {
 119                         HklAxis *axis = self->axes[i];
 120                         if (degenerated[i])
 121                                 axis->config.value = x_data0[i];
 122                         else
 123                                 axis->config.value = x_data[i];
 124
 125                         axis->config.dirty = 1;
 126                 }
 127                 hkl_geometry_update(self->geometry);
 128                 res = HKL_SUCCESS;
 129         }
 130
 131         // release memory
 132         gsl_vector_free(x);
 133         gsl_multiroot_fsolver_free(s);
 134
 135         return res;
 136 }
 137
 138 /**
 139  * @brief This private method change the sector of angles.
 140  *
 141  * @param x The vector of changed angles.
 142  * @param x0 The vector of angles to change.
 143  * @param sector the sector vector operation.
 144  * @param n the size of all vectors.
 145  *
 146  * 0 -> no change
 147  * 1 -> pi - angle
 148  * 2 -> pi + angle
 149  * 3 -> -angle
 150  */
 151 static void change_sector(double x[], double const x0[],
 152                           int const sector[], size_t n)
 153 {
 154         size_t i;
 155
 156         for(i=0; i<n; ++i) {
 157                 switch (sector[i]) {
 158                 case 0:
 159                         x[i] = x0[i];
 160                         break;
 161                 case 1:
 162                         x[i] = M_PI - x0[i];
 163                         break;
 164                 case 2:
 165                         x[i] = M_PI + x0[i];
 166                         break;
 167                 case 3:
 168                         x[i] = -x0[i];
 169                         break;
 170                 }
 171         }
 172 }
 173
 174 /**
 175  * @brief Test if an angle combination is compatible with q function.
 176  *
 177  * @param x The vector of angles to test.
 178  * @param function The gsl_multiroot_function used for the test.
 179  * @param f a gsl_vector use to compute the result (optimization)
 180  */
 181 static int test_sector(gsl_vector const *x,
 182                        gsl_multiroot_function *function,
 183                        gsl_vector *f)
 184 {
 185         size_t i;
 186         double *f_data = f->data;
 187
 188         function->f(x, function->params, f);
 189
 190         for(i=0; i<f->size; ++i)
 191                 if (fabs(f_data[i]) > HKL_EPSILON)
 192                         return HKL_FAIL;
 193
 194         return HKL_SUCCESS;
 195 }
 196
 197 /**
 198  * @brief compute the permutation and test its validity.
 199  *
 200  * @param axes_len number of axes
 201  * @param op_len number of operation per axes. (4 for now)
 202  * @param p The vector containing the current permutation.
 203  * @param axes_idx The index of the axes we are permution.
 204  * @param op the current operation to set.
 205  * @param f The function for the validity test.
 206  * @param x0 The starting point of all geometry permutations.
 207  * @param _x a gsl_vector use to compute the sectors (optimization)
 208  * @param _f a gsl_vector use during the sector test (optimization)
 209  */
 210 static void perm_r(size_t axes_len, int op_len[], int p[], int axes_idx,
 211                    int op, gsl_multiroot_function *f, double x0[],
 212                    gsl_vector *_x, gsl_vector *_f)
 213 {
 214         int i;
 215
 216         p[axes_idx++] = op;
 217         if (axes_idx == axes_len) {
 218                 double *x_data = _x->data;
 219                 change_sector(x_data, x0, p, axes_len);
 220                 if (HKL_SUCCESS == test_sector(_x, f, _f))
 221                         hkl_pseudoAxisEngine_add_geometry(f->params, x_data);
 222         } else
 223                 for (i=0; i<op_len[axes_idx]; ++i)
 224                         perm_r(axes_len, op_len, p, axes_idx, i, f, x0, _x, _f);
 225 }
 226
 227 /**
 228  * @brief Find all numerical solutions of a mode.
 229  *
 230  * @param self the current HklPseudoAxisEngine
 231  * @param function The mode function
 232  *
 233  * @return HKL_SUCCESS (0) or HKL_FAIL (-1)
 234  *
 235  * This method find a first solution with a numerical method from the
 236  * GSL library (the multi root solver hybrid). Then it multiplicates the
 237  * solutions from this starting point using cosinus/sinus properties.
 238  * It addes all valid solutions to the self->geometries.
 239  */
 240 int hkl_pseudoAxeEngine_solve_function(HklPseudoAxisEngine *self,
 241                                        HklPseudoAxisEngineFunction function)
 242 {
 243
 244         size_t i;
 245         size_t n = self->axes_len;
 246         int p[n];
 247         double x0[n];
 248         int degenerated[n];
 249         int op_len[n];
 250         int res;
 251         gsl_vector *_x; /* use to compute sectors in perm_r (avoid copy) */
 252         gsl_vector *_f; /* use to test sectors in perm_r (avoid copy) */
 253
 254         _x = gsl_vector_alloc(n);
 255         _f = gsl_vector_alloc(n);
 256         gsl_multiroot_function f = {function, self->axes_len, self};
 257         res = find_first_geometry(self, &f, degenerated);
 258         if (res == HKL_SUCCESS) {
 259                 memset(p, 0, sizeof(p));
 260                 /* use first solution as starting point for permutations */
 261                 for(i=0; i<n; ++i){
 262                         x0[i] = self->axes[i]->config.value;
 263                         if (degenerated[i])
 264                                 op_len[i] = 1;
 265                         else
 266                                 op_len[i] = 4;
 267                 }
 268                 for (i=0; i<op_len[0]; ++i)
 269                         perm_r(n, op_len, p, 0, i, &f, x0, _x, _f);
 270         }
 271         gsl_vector_free(_f);
 272         gsl_vector_free(_x);
 273         return res;
 274 }