src/hkl-sample.c

   1 #include <gsl/gsl_multimin.h>
   2
   3 #include <hkl/hkl-sample.h>
   4 #include <hkl/hkl-matrix.h>
   5
   6 /* private */
   7
   8 static void *copy_ref(void const *item)
   9 {
  10         HklSampleReflection *copy = NULL;
  11         HklSampleReflection const *src = item;
  12
  13         copy = malloc(sizeof(*copy));
  14         if (!copy)
  15                 die("Can not allocate memory for a HklSampleReflection");
  16
  17         copy->geometry = hkl_geometry_new_copy(src->geometry);
  18         copy->detector = src->detector;
  19         copy->hkl = src->hkl;
  20         copy->_hkl = src->_hkl;
  21
  22         return copy;
  23 }
  24
  25 static void free_ref(void *item)
  26 {
  27         HklSampleReflection *ref = item;
  28         hkl_geometry_free(ref->geometry);
  29         free(ref);
  30 }
  31
  32 static int hkl_sample_compute_UB(HklSample *sample)
  33 {
  34         HklMatrix B;
  35
  36         if (hkl_lattice_get_B(sample->lattice, &B))
  37                 return HKL_FAIL;
  38
  39         sample->UB = sample->U;
  40         hkl_matrix_times_smatrix(&sample->UB, &B);
  41
  42         return HKL_SUCCESS;
  43 }
  44
  45 static double mono_crystal_fitness(gsl_vector const *x, void *params)
  46 {
  47         size_t i, j;
  48         double fitness;
  49         double euler_x;
  50         double euler_y;
  51         double euler_z;
  52         HklSample *sample = params;
  53
  54         euler_x = gsl_vector_get(x, 0);
  55         euler_y = gsl_vector_get(x, 1);
  56         euler_z = gsl_vector_get(x, 2);
  57         sample->lattice->a->value = gsl_vector_get(x, 3);
  58         sample->lattice->b->value = gsl_vector_get(x, 4);
  59         sample->lattice->c->value = gsl_vector_get(x, 5);
  60         sample->lattice->alpha->value = gsl_vector_get(x, 6);
  61         sample->lattice->beta->value = gsl_vector_get(x, 7);
  62         sample->lattice->gamma->value = gsl_vector_get(x, 8);
  63         hkl_matrix_from_euler(&sample->U, euler_x, euler_y, euler_z);
  64         if (hkl_sample_compute_UB(sample))
  65                 return GSL_NAN;
  66
  67         fitness = 0.;
  68         for(i=0; i<sample->reflections->len; ++i) {
  69                 HklVector UBh;
  70                 HklSampleReflection *reflection;
  71
  72                 reflection = hkl_list_get_by_idx(sample->reflections, i);
  73                 UBh = reflection->hkl;
  74                 hkl_matrix_times_vector(&sample->UB, &UBh);
  75
  76                 for(j=0; j<3; ++j) {
  77                         double tmp = UBh.data[j] - reflection->_hkl.data[j];
  78                         fitness += tmp * tmp;
  79                 }
  80         }
  81         return fitness;
  82 }
  83
  84 /* public */
  85
  86 HklSample* hkl_sample_new(char const *name, HklSampleType type)
  87 {
  88         HklSample *sample = NULL;
  89         sample = malloc(sizeof(*sample));
  90         if (!sample)
  91                 die("Cannot allocate memory for a Sample");
  92         sample->name = name;
  93         sample->type = type;
  94         sample->lattice = hkl_lattice_new_default();
  95         hkl_matrix_init(&sample->U,1, 0, 0, 0, 1, 0, 0, 0, 1);
  96         hkl_matrix_init(&sample->UB,1, 0, 0, 0, 1, 0, 0, 0, 1);
  97         hkl_sample_compute_UB(sample);
  98         sample->reflections = hkl_list_new_managed(&copy_ref, &free_ref);
  99
 100         return sample;
 101 }
 102
 103 HklSample *hkl_sample_new_copy(HklSample const *src)
 104 {
 105         HklSample *copy;
 106
 107         copy = malloc(sizeof(*copy));
 108         if (!copy)
 109                 die("Cannot allocate memory for a Sample");
 110         copy->name = src->name;
 111         copy->type = src->type;
 112         copy->lattice = hkl_lattice_new_copy(src->lattice);
 113         copy->U = src->U;
 114         copy->UB = src->UB;
 115         copy->reflections = hkl_list_new_copy(src->reflections);
 116
 117         return copy;
 118 }
 119
 120 void hkl_sample_free(HklSample *sample)
 121 {
 122         hkl_lattice_free(sample->lattice);
 123         hkl_list_free(sample->reflections);
 124         free(sample);
 125 }
 126
 127 HklSampleReflection *hkl_sample_add_reflection(HklSample *sample,
 128                 HklGeometry *g, HklDetector const *det,
 129                 double h, double k, double l)
 130 {
 131         HklSampleReflection *ref;
 132         HklHolder *holder_d;
 133         HklHolder *holder_s;
 134         HklVector ki;
 135         HklQuaternion q;
 136
 137         if (fabs(h) < HKL_EPSILON
 138                         && fabs(k) < HKL_EPSILON
 139                         && fabs(l) < HKL_EPSILON)
 140                 return NULL;
 141
 142         ref = malloc(sizeof(*ref));
 143         if (!ref)
 144                 die("Cannot allocate memory for an HklSampleReflection");
 145
 146         hkl_geometry_update(g);
 147
 148         ref->geometry = hkl_geometry_new_copy(g);
 149         ref->detector = *det;
 150         ref->hkl.data[0] = h;
 151         ref->hkl.data[1] = k;
 152         ref->hkl.data[2] = l;
 153
 154         // compute the _hkl using only the axes of the geometry
 155         holder_d = &ref->geometry->holders[det->idx];
 156         holder_s = &ref->geometry->holders[0];
 157
 158         // compute Q from angles
 159         hkl_source_compute_ki(&ref->geometry->source, &ki);
 160         ref->_hkl = ki;
 161         hkl_vector_rotated_quaternion(&ref->_hkl, &holder_d->q);
 162         hkl_vector_minus_vector(&ref->_hkl, &ki);
 163
 164         q = holder_s->q;
 165         hkl_quaternion_conjugate(&q);
 166         hkl_vector_rotated_quaternion(&ref->_hkl, &q);
 167
 168         hkl_list_append(sample->reflections, ref);
 169
 170         return ref;
 171 }
 172
 173 HklSampleReflection* hkl_sample_get_reflection(HklSample *sample,
 174                 size_t idx)
 175 {
 176         return hkl_list_get_by_idx(sample->reflections, idx);
 177 }
 178
 179 int hkl_sample_del_reflection(HklSample *sample, size_t idx)
 180 {
 181         return hkl_list_del_by_idx(sample->reflections, idx);
 182 }
 183
 184 int hkl_sample_compute_UB_busing_levy(HklSample *sample,
 185                 size_t idx1, size_t idx2)
 186 {
 187         if (idx1 < sample->reflections->len
 188                         && idx2 < sample->reflections->len) {
 189
 190                 HklSampleReflection *r1;
 191                 HklSampleReflection *r2;
 192
 193                 r1 = hkl_list_get_by_idx(sample->reflections, idx1);
 194                 r2 = hkl_list_get_by_idx(sample->reflections, idx2);
 195
 196                 if (!hkl_vector_is_colinear(&r1->hkl, &r2->hkl)) {
 197                         HklVector h1c;
 198                         HklVector h2c;
 199                         HklMatrix B;
 200                         HklMatrix Tc;
 201
 202                         // Compute matrix Tc from r1 and r2.
 203                         h1c = r1->hkl;
 204                         h2c = r2->hkl;
 205                         hkl_lattice_get_B(sample->lattice, &B);
 206                         hkl_matrix_times_vector(&B, &h1c);
 207                         hkl_matrix_times_vector(&B, &h2c);
 208                         hkl_matrix_from_two_vector(&Tc, &h1c, &h2c);
 209                         hkl_matrix_transpose(&Tc);
 210
 211                         // compute U
 212                         hkl_matrix_from_two_vector(&sample->U,
 213                                         &r1->_hkl, &r2->_hkl);
 214                         hkl_matrix_times_smatrix(&sample->U, &Tc);
 215                 } else
 216                         return HKL_FAIL;
 217         } else
 218                 return HKL_FAIL;
 219
 220         return HKL_SUCCESS;
 221 }
 222
 223 void hkl_sample_affine(HklSample *sample)
 224 {
 225         gsl_multimin_fminimizer_type const *T = gsl_multimin_fminimizer_nmsimplex;
 226         gsl_multimin_fminimizer *s = NULL;
 227         gsl_vector *ss, *x;
 228         gsl_multimin_function minex_func;
 229         size_t iter = 0;
 230         int status;
 231         double size;
 232
 233         // Starting point
 234         x = gsl_vector_alloc (9);
 235         gsl_vector_set (x, 0, 10 * HKL_DEGTORAD);
 236         gsl_vector_set (x, 1, 10 * HKL_DEGTORAD);
 237         gsl_vector_set (x, 2, 10 * HKL_DEGTORAD);
 238         gsl_vector_set (x, 3, sample->lattice->a->value);
 239         gsl_vector_set (x, 4, sample->lattice->b->value);
 240         gsl_vector_set (x, 5, sample->lattice->c->value);
 241         gsl_vector_set (x, 6, sample->lattice->alpha->value);
 242         gsl_vector_set (x, 7, sample->lattice->beta->value);
 243         gsl_vector_set (x, 8, sample->lattice->gamma->value);
 244
 245         // Set initial step sizes to 1
 246         ss = gsl_vector_alloc (9);
 247         gsl_vector_set (ss, 0, 1 * HKL_DEGTORAD);
 248         gsl_vector_set (ss, 1, 1 * HKL_DEGTORAD);
 249         gsl_vector_set (ss, 2, 1 * HKL_DEGTORAD);
 250         gsl_vector_set (ss, 3, !sample->lattice->a->not_to_fit);
 251         gsl_vector_set (ss, 4, !sample->lattice->b->not_to_fit);
 252         gsl_vector_set (ss, 5, !sample->lattice->c->not_to_fit);
 253         gsl_vector_set (ss, 6, !sample->lattice->alpha->not_to_fit);
 254         gsl_vector_set (ss, 7, !sample->lattice->beta->not_to_fit);
 255         gsl_vector_set (ss, 8, !sample->lattice->gamma->not_to_fit);
 256
 257         // Initialize method and iterate
 258         minex_func.n = 9;
 259         minex_func.f = &mono_crystal_fitness;
 260         minex_func.params = sample;
 261         s = gsl_multimin_fminimizer_alloc (T, 9);
 262         gsl_set_error_handler_off();
 263         gsl_multimin_fminimizer_set (s, &minex_func, x, ss);
 264         do {
 265                 ++iter;
 266                 status = gsl_multimin_fminimizer_iterate(s);
 267                 if (status)
 268                         break;
 269                 size = gsl_multimin_fminimizer_size (s);
 270                 status = gsl_multimin_test_size (size, HKL_EPSILON / 2.);
 271         } while (status == GSL_CONTINUE && iter < 10000);
 272         gsl_vector_free(x);
 273         gsl_vector_free(ss);
 274         gsl_multimin_fminimizer_free(s);
 275         gsl_set_error_handler (NULL);
 276 }
 277
 278 /*
 279 void hkl_sample_fprintf(HklSample *sample, FILE *f)
 280 {
 281         fprintf(f, "\"%s\"\n", sample->name);
 282 }
 283 */