IndexingUtils.cpp 109 KB
Newer Older
1
2
3
// Mantid Repository : https://github.com/mantidproject/mantid
//
// Copyright © 2018 ISIS Rutherford Appleton Laboratory UKRI,
4
5
//   NScD Oak Ridge National Laboratory, European Spallation Source,
//   Institut Laue - Langevin & CSNS, Institute of High Energy Physics, CAS
6
// SPDX - License - Identifier: GPL - 3.0 +
7
#include "MantidGeometry/Crystal/IndexingUtils.h"
8
#include "MantidGeometry/Crystal/NiggliCell.h"
9
#include "MantidKernel/EigenConversionHelpers.h"
Lynch, Vickie's avatar
Lynch, Vickie committed
10
#include "MantidKernel/Quat.h"
11

12
13
#include <boost/math/special_functions/round.hpp>
#include <boost/numeric/conversion/cast.hpp>
14

15
#include <gsl/gsl_errno.h>
Lynch, Vickie's avatar
Lynch, Vickie committed
16
17
18
#include <gsl/gsl_fft_real.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_matrix.h>
19
20
#include <gsl/gsl_sys.h>
#include <gsl/gsl_vector.h>
21
22
23

#include <algorithm>
#include <cmath>
24
25

using namespace Mantid::Geometry;
26
using Mantid::Kernel::DblMatrix;
Lynch, Vickie's avatar
Lynch, Vickie committed
27
using Mantid::Kernel::Matrix;
28
using Mantid::Kernel::Quat;
Lynch, Vickie's avatar
Lynch, Vickie committed
29
using Mantid::Kernel::V3D;
30

31
namespace {
32
33
const constexpr double DEG_TO_RAD = M_PI / 180.;
const constexpr double RAD_TO_DEG = 180. / M_PI;
Lynch, Vickie's avatar
Lynch, Vickie committed
34
} // namespace
35

36
37
38
/**
  STATIC method Find_UB: Calculates the matrix that most nearly indexes
  the specified q_vectors, given the lattice parameters.  The sum of the
39
  squares of the residual errors is returned.  This method first sorts the
40
  specified q_vectors in order of increasing magnitude.  It then searches
41
  through all possible orientations to find an initial UB that indexes the
42
  lowest magnitude peaks.
43
44
     The resolution of the search through possible orientations is specified
  by the degrees_per_step parameter.  Approximately 2-4 degrees_per_step is
45
46
47
48
  usually adequate.  NOTE: This is an expensive calculation which takes
  approximately 1 second using 2 degrees_per_step.  However, the execution
  time is O(n^3) so decreasing the resolution to 1 degree per step will take
  about 8 seconds, etc.  It should not be necessary to decrease this value
49
50
51
  below 1 degree per step, and users will have to be VERY patient, if it is
  decreased much below 1 degree per step.
    The number of peaks used to obtain an initial indexing is specified by
52
53
  the "NumInitial" parameter.  Good values for this are typically around
  10-15, though with accurate peak positions, and good values for the lattice
54
55
  paramters, as few as 2 can be used.  Using substantially more than 15 peaks
  initially typically has no benefit and increases execution time.
56

57
  @param  UB                  3x3 matrix that will be set to the UB matrix
58
  @param  q_vectors           std::vector of V3D objects that contains the
59
                              list of q_vectors that are to be indexed
60
                              NOTE: There must be at least 2 q_vectors.
61
  @param  lattice             The orientated lattice with the lattice
62
63
                              parameters a,b,c and alpha, beta, gamma. The found
                              UB and errors will be set on this lattice.
64
65
  @param  required_tolerance  The maximum allowed deviation of Miller indices
                              from integer values for a peak to be indexed.
66
  @param  base_index          The sequence number of the peak that should
67
68
                              be used as the central peak.  On the first
                              scan for a UB matrix that fits the data,
69
                              the remaining peaks in the list of q_vectors
70
71
72
73
74
75
76
77
                              will be shifted by -base_peak, where base_peak
                              is the q_vector with the specified base index.
                              If fewer than 5 peaks are specified in the
                              q_vectors list, this parameter is ignored.
                              If this parameter is -1, and there are at least
                              four peaks in the q_vector list, then a base
                              index will be calculated internally.  In most
                              cases, it should suffice to set this to -1.
78
79
80
81
  @param  num_initial         The number of low |Q| peaks that should be
                              used to scan for an initial orientation matrix.
  @param  degrees_per_step    The number of degrees between different
                              orientations used during the initial scan.
82
83
  @param  fixAll              Fix the lattice parameters and do not optimise
                              the UB matrix.
84
  @param  iterations          Number of refinements of UB
85
86

  @return  This will return the sum of the squares of the residual errors.
87
88
89
90

  @throws  std::invalid_argument exception if UB is not a 3X3 matrix,
                                 if there are not at least 2 q vectors,
                                 if num_initial is < 2, or
91
                                 if the required_tolerance or degrees_per_step
92
                                 is <= 0.
93
*/
94
95
double IndexingUtils::Find_UB(DblMatrix &UB, const std::vector<V3D> &q_vectors, OrientedLattice &lattice,
                              double required_tolerance, int base_index, size_t num_initial, double degrees_per_step,
96
                              bool fixAll, int iterations) {
97
98
  if (UB.numRows() != 3 || UB.numCols() != 3) {
    throw std::invalid_argument("Find_UB(): UB matrix NULL or not 3X3");
99
100
  }

101
  if (q_vectors.size() < 2) {
102
    throw std::invalid_argument("Find_UB(): Two or more indexed peaks needed to find UB");
103
104
  }

105
  if (required_tolerance <= 0) {
106
    throw std::invalid_argument("Find_UB(): required_tolerance must be positive");
107
108
  }

109
  if (num_initial < 2) {
110
    throw std::invalid_argument("Find_UB(): number of peaks for inital scan must be > 2");
111
112
  }

113
114
  if (degrees_per_step <= 0) {
    throw std::invalid_argument("Find_UB(): degrees_per_step must be positive");
115
116
  }

117
  // get list of peaks, sorted on |Q|
118
  std::vector<V3D> sorted_qs;
119
120
121
122
123
124
125
126
127
128
  sorted_qs.reserve(q_vectors.size());

  if (q_vectors.size() > 5) // shift to be centered on peak (we lose
                            // one peak that way, so require > 5)
  {
    std::vector<V3D> shifted_qs(q_vectors);
    size_t mid_ind = q_vectors.size() / 3;
    // either do an initial sort and use
    // default mid index, or use the index
    // specified by the base_peak parameter
129
    if (base_index < 0 || base_index >= static_cast<int>(q_vectors.size())) {
130
      std::sort(shifted_qs.begin(), shifted_qs.end(), V3D::compareMagnitude);
131
    } else {
132
133
      mid_ind = base_index;
    }
134
    V3D mid_vec(shifted_qs[mid_ind]);
135

136
137
138
    for (size_t i = 0; i < shifted_qs.size(); i++) {
      if (i != mid_ind) {
        V3D shifted_vec(shifted_qs[i]);
139
        shifted_vec -= mid_vec;
140
        sorted_qs.emplace_back(shifted_vec);
141
142
      }
    }
143
  } else {
144
    std::copy(q_vectors.cbegin(), q_vectors.cend(), std::back_inserter(sorted_qs));
145
146
  }

147
  std::sort(sorted_qs.begin(), sorted_qs.end(), V3D::compareMagnitude);
148

149
  if (num_initial > sorted_qs.size())
150
151
152
    num_initial = sorted_qs.size();

  std::vector<V3D> some_qs;
153
  some_qs.reserve(q_vectors.size());
154

155
  for (size_t i = 0; i < num_initial; i++)
156
    some_qs.emplace_back(sorted_qs[i]);
157

158
  ScanFor_UB(UB, some_qs, lattice, degrees_per_step, required_tolerance);
159

160
  double fit_error = 0;
161
162
  std::vector<V3D> miller_ind;
  std::vector<V3D> indexed_qs;
163
164
  miller_ind.reserve(q_vectors.size());
  indexed_qs.reserve(q_vectors.size());
165

166
167
168
  // now gradually bring in the remaining
  // peaks and re-optimize the UB to index
  // them as well
169
  while (!fixAll && num_initial < sorted_qs.size()) {
170
    num_initial = std::lround(1.5 * static_cast<double>(num_initial + 3));
171
172
173
    // add 3, in case we started with
    // a very small number of peaks!
    if (num_initial >= sorted_qs.size())
174
      num_initial = sorted_qs.size();
175

176
    for (size_t i = some_qs.size(); i < num_initial; i++)
177
      some_qs.emplace_back(sorted_qs[i]);
178
179
    for (int counter = 0; counter < iterations; counter++) {
      try {
180
        GetIndexedPeaks(UB, some_qs, required_tolerance, miller_ind, indexed_qs, fit_error);
181
182
183
        Matrix<double> temp_UB(3, 3, false);
        fit_error = Optimize_UB(temp_UB, miller_ind, indexed_qs);
        UB = temp_UB;
Lynch, Vickie's avatar
Lynch, Vickie committed
184
      } catch (...) {
185
186
        // failed to fit using these peaks, so add some more and try again
      }
187
188
189
    }
  }

190
  std::vector<double> sigabc(7);
191
  if (!fixAll && q_vectors.size() >= 3) // try one last refinement using all peaks
192
  {
193
194
    for (int counter = 0; counter < iterations; counter++) {
      try {
195
        GetIndexedPeaks(UB, q_vectors, required_tolerance, miller_ind, indexed_qs, fit_error);
196
197
198
        Matrix<double> temp_UB = UB;
        fit_error = Optimize_UB(temp_UB, miller_ind, indexed_qs, sigabc);
        UB = temp_UB;
Lynch, Vickie's avatar
Lynch, Vickie committed
199
      } catch (...) {
200
201
        // failed to improve UB using these peaks, so just return the current UB
      }
202
203
    }
  }
204
205
206
  // Regardless of how we got the UB, find the
  // sum-squared errors for the indexing in
  // HKL space.
207
  GetIndexedPeaks(UB, q_vectors, required_tolerance, miller_ind, indexed_qs, fit_error);
208
209
  // set the error on the lattice parameters
  lattice.setUB(UB);
210
  lattice.setError(sigabc[0], sigabc[1], sigabc[2], sigabc[3], sigabc[4], sigabc[5]);
211
212
213
  return fit_error;
}

214
215
216
217
/**
    STATIC method Find_UB: This method will attempt to calculate the matrix
  that most nearly indexes the specified q_vectors, given only a range of
  possible unit cell edge lengths.  If successful, the matrix should
218
  correspond to the Niggli reduced cell.
219
220
     The resolution of the search through possible orientations is specified
  by the degrees_per_step parameter.  Approximately 1-3 degrees_per_step is
221
222
223
224
  usually adequate.  NOTE: This is an expensive calculation which takes
  approximately 1 second using 1 degree_per_step.  However, the execution
  time is O(n^3) so decreasing the resolution to 0.5 degree per step will take
  about 8 seconds, etc.  It should not be necessary to decrease this value
225
226
227
  below 1 degree per step, and users will have to be VERY patient, if it is
  decreased much below 1 degree per step.
    The number of peaks used to obtain an initial indexing is specified by
228
  the "NumInitial" parameter.  Good values for this are typically around
229
  15-25.  The specified q_vectors must correspond to a single crystal.  If
230
  several crystallites are present or there are other sources of "noise"
231
232
233
234
235
  leading to invalid peaks, this method will not work well.  The method that
  uses lattice parameters may be better in such cases.  Alternatively, adjust
  the list of specified q_vectors so it does not include noise peaks or peaks
  from more than one crystal, by increasing the threshold for what counts
  as a peak, or by other methods.
236

237
  @param  UB                  3x3 matrix that will be set to the UB matrix
238
  @param  q_vectors           std::vector of V3D objects that contains the
239
                              list of q_vectors that are to be indexed
240
                              NOTE: There must be at least 3 q_vectors.
241
242
243
244
245
246
247
248
249
250
  @param  min_d               Lower bound on shortest unit cell edge length.
                              This does not have to be specified exactly but
                              must be strictly less than the smallest edge
                              length, in Angstroms.
  @param  max_d               Upper bound on longest unit cell edge length.
                              This does not have to be specified exactly but
                              must be strictly more than the longest edge
                              length in angstroms.
  @param  required_tolerance  The maximum allowed deviation of Miller indices
                              from integer values for a peak to be indexed.
251
  @param  base_index          The sequence number of the peak that should
252
253
                              be used as the central peak.  On the first
                              scan for a UB matrix that fits the data,
254
                              the remaining peaks in the list of q_vectors
255
256
257
258
259
260
261
262
263
264
265
266
267
268
                              will be shifted by -base_peak, where base_peak
                              is the q_vector with the specified base index.
                              If fewer than 6 peaks are specified in the
                              q_vectors list, this parameter is ignored.
                              If this parameter is -1, and there are at least
                              five peaks in the q_vector list, then a base
                              index will be calculated internally.  In most
                              cases, it should suffice to set this to -1.
  @param  num_initial         The number of low |Q| peaks that should be
                              used to scan for an initial orientation matrix.
  @param  degrees_per_step    The number of degrees between different
                              orientations used during the initial scan.

  @return  This will return the sum of the squares of the residual errors.
269
270
271

  @throws  std::invalid_argument exception if UB is not a 3X3 matrix,
                                 if there are not at least 3 q vectors,
272
                                 if min_d >= max_d or min_d <= 0
273
                                 if num_initial is < 3, or
274
                                 if the required_tolerance or degrees_per_step
275
                                 is <= 0.
276
*/
277
278
double IndexingUtils::Find_UB(DblMatrix &UB, const std::vector<V3D> &q_vectors, double min_d, double max_d,
                              double required_tolerance, int base_index, size_t num_initial, double degrees_per_step) {
279
280
  if (UB.numRows() != 3 || UB.numCols() != 3) {
    throw std::invalid_argument("Find_UB(): UB matrix NULL or not 3X3");
281
282
  }

283
  if (q_vectors.size() < 3) {
284
    throw std::invalid_argument("Find_UB(): Three or more indexed peaks needed to find UB");
285
286
  }

287
288
  if (min_d >= max_d || min_d <= 0) {
    throw std::invalid_argument("Find_UB(): Need 0 < min_d < max_d");
289
290
  }

291
  if (required_tolerance <= 0) {
292
    throw std::invalid_argument("Find_UB(): required_tolerance must be positive");
293
294
  }

295
  if (num_initial < 3) {
296
    throw std::invalid_argument("Find_UB(): number of peaks for inital scan must be > 2");
297
298
  }

299
300
  if (degrees_per_step <= 0) {
    throw std::invalid_argument("Find_UB(): degrees_per_step must be positive");
301
302
  }

303
  // get list of peaks, sorted on |Q|
304
  std::vector<V3D> sorted_qs;
305
306
307
308
309
310
311
312
313
314
  sorted_qs.reserve(q_vectors.size());

  if (q_vectors.size() > 5) // shift to be centered on peak (we lose
                            // one peak that way, so require > 5)
  {
    std::vector<V3D> shifted_qs(q_vectors);
    size_t mid_ind = q_vectors.size() / 2;
    // either do an initial sort and use
    // default mid index, or use the index
    // specified by the base_peak parameter
315
    if (base_index < 0 || base_index >= static_cast<int>(q_vectors.size())) {
316
      std::sort(shifted_qs.begin(), shifted_qs.end(), V3D::compareMagnitude);
317
    } else {
318
319
      mid_ind = base_index;
    }
320
    V3D mid_vec(shifted_qs[mid_ind]);
321

322
323
324
    for (size_t i = 0; i < shifted_qs.size(); i++) {
      if (i != mid_ind) {
        V3D shifted_vec(shifted_qs[i]);
325
        shifted_vec -= mid_vec;
326
        sorted_qs.emplace_back(shifted_vec);
327
328
      }
    }
329
  } else {
330
    std::copy(q_vectors.cbegin(), q_vectors.cend(), std::back_inserter(sorted_qs));
331
332
  }

333
  std::sort(sorted_qs.begin(), sorted_qs.end(), V3D::compareMagnitude);
334

335
  if (num_initial > sorted_qs.size())
336
337
338
    num_initial = sorted_qs.size();

  std::vector<V3D> some_qs;
339
  some_qs.reserve(q_vectors.size());
340

341
  for (size_t i = 0; i < num_initial; i++)
342
    some_qs.emplace_back(sorted_qs[i]);
343
  std::vector<V3D> directions;
344
  ScanFor_Directions(directions, some_qs, min_d, max_d, required_tolerance, degrees_per_step);
345

346
  if (directions.size() < 3) {
347
    throw std::runtime_error("Find_UB(): Could not find at least three possible lattice directions");
348
  }
349

350
  std::sort(directions.begin(), directions.end(), V3D::compareMagnitude);
351

352
  if (!FormUB_From_abc_Vectors(UB, directions, 0, min_d, max_d)) {
353
    throw std::runtime_error("Find_UB(): Could not find independent a, b, c directions");
354
  }
355
356
357
  // now gradually bring in the remaining
  // peaks and re-optimize the UB to index
  // them as well
358
359
  std::vector<V3D> miller_ind;
  std::vector<V3D> indexed_qs;
360
361
  miller_ind.reserve(q_vectors.size());
  indexed_qs.reserve(q_vectors.size());
362

363
364
365
  Matrix<double> temp_UB(3, 3, false);
  double fit_error = 0;
  while (num_initial < sorted_qs.size()) {
366
    num_initial = std::lround(1.5 * static_cast<double>(num_initial + 3));
367
368
369
    // add 3, in case we started with
    // a very small number of peaks!
    if (num_initial >= sorted_qs.size())
370
371
      num_initial = sorted_qs.size();

372
    for (size_t i = some_qs.size(); i < num_initial; i++)
373
      some_qs.emplace_back(sorted_qs[i]);
374

375
    GetIndexedPeaks(UB, some_qs, required_tolerance, miller_ind, indexed_qs, fit_error);
376
377
    try {
      fit_error = Optimize_UB(temp_UB, miller_ind, indexed_qs);
378
      UB = temp_UB;
Lynch, Vickie's avatar
Lynch, Vickie committed
379
    } catch (...) {
380
381
382
383
384
      // failed to improve with these peaks, so continue with more peaks
      // if possible
    }
  }

385
  if (q_vectors.size() >= 5) // try one last refinement using all peaks
386
  {
387
    GetIndexedPeaks(UB, q_vectors, required_tolerance, miller_ind, indexed_qs, fit_error);
388
389
    try {
      fit_error = Optimize_UB(temp_UB, miller_ind, indexed_qs);
390
      UB = temp_UB;
Lynch, Vickie's avatar
Lynch, Vickie committed
391
    } catch (...) {
392
      // failed to improve with all peaks, so return the UB we had
393
394
395
    }
  }

396
  if (NiggliCell::MakeNiggliUB(UB, temp_UB))
397
398
    UB = temp_UB;

399
400
401
  return fit_error;
}

402
403
404
/**
    STATIC method Find_UB: This method will attempt to calculate the matrix
  that most nearly indexes the specified q_vectors, using FFTs to find
405
  patterns in projected Q-vectors, given only a range of possible unit cell
406
  edge lengths. If successful, the resulting matrix should correspond to
407
  the Niggli reduced cell.
408
    The resolution of the search through possible orientations is specified
409
  by the degrees_per_step parameter.  One to two degrees per step is usually
410
411
  adequate.
  NOTE: The execution time is O(n^3) where n is the number of degrees per
412
413
414
  step, so decreasing the resolution to 0.5 degree per step will take
  about 8 times longer than using 1 degree per step.  It should not be
  necessary to decrease this value below 1 degree per step, and users will
415
416
417
  have to be VERY patient, if it is decreased much below 1 degree per step.
    The specified q_vectors should correspond to a single crystal, for this
  to work reliably.
418

419
  @param  UB                  3x3 matrix that will be set to the UB matrix
420
  @param  q_vectors           std::vector of V3D objects that contains the
421
422
                              list of q_vectors that are to be indexed
                              NOTE: There must be at least 4 q_vectors and it
423
                              really should have at least 10 or more peaks
424
425
426
427
428
429
430
431
432
433
434
435
436
                              for this to work quite consistently.
  @param  min_d               Lower bound on shortest unit cell edge length.
                              This does not have to be specified exactly but
                              must be strictly less than the smallest edge
                              length, in Angstroms.
  @param  max_d               Upper bound on longest unit cell edge length.
                              This does not have to be specified exactly but
                              must be strictly more than the longest edge
                              length in angstroms.
  @param  required_tolerance  The maximum allowed deviation of Miller indices
                              from integer values for a peak to be indexed.
  @param  degrees_per_step    The number of degrees between different
                              orientations used during the initial scan.
437
  @param  iterations          Number of refinements of UB
438
439

  @return  This will return the sum of the squares of the residual errors.
440
441
442

  @throws  std::invalid_argument exception if UB is not a 3X3 matrix,
                                 if there are not at least 3 q vectors,
443
444
445
446
447
                                 if min_d >= max_d or min_d <= 0,
                                 if the required_tolerance or degrees_per_step
                                 is <= 0,
                                 if at least three possible a,b,c directions
                                 were not found,
448
                                 or if a valid UB matrix could not be formed
449
450
                                 from the a,b,c directions that were found.
*/
451
452
double IndexingUtils::Find_UB(DblMatrix &UB, const std::vector<V3D> &q_vectors, double min_d, double max_d,
                              double required_tolerance, double degrees_per_step, int iterations) {
453
454
  if (UB.numRows() != 3 || UB.numCols() != 3) {
    throw std::invalid_argument("Find_UB(): UB matrix NULL or not 3X3");
455
456
  }

457
  if (q_vectors.size() < 4) {
458
    throw std::invalid_argument("Find_UB(): Four or more indexed peaks needed to find UB");
459
460
  }

461
462
  if (min_d >= max_d || min_d <= 0) {
    throw std::invalid_argument("Find_UB(): Need 0 < min_d < max_d");
463
464
  }

465
  if (required_tolerance <= 0) {
466
    throw std::invalid_argument("Find_UB(): required_tolerance must be positive");
467
468
  }

469
470
  if (degrees_per_step <= 0) {
    throw std::invalid_argument("Find_UB(): degrees_per_step must be positive");
471
472
473
474
  }

  std::vector<V3D> directions;

475
476
477
478
479
  // NOTE: we use a somewhat higher tolerance when
  // finding individual directions since it is easier
  // to index one direction individually compared to
  // indexing three directions simultaneously.
  size_t max_indexed =
480
      FFTScanFor_Directions(directions, q_vectors, min_d, max_d, 0.75f * required_tolerance, degrees_per_step);
481

482
  if (max_indexed == 0) {
483
    throw std::invalid_argument("Find_UB(): Could not find any a,b,c vectors to index Qs");
484
485
  }

486
  if (directions.size() < 3) {
487
    throw std::invalid_argument("Find_UB(): Could not find enough a,b,c vectors");
488
489
  }

490
  std::sort(directions.begin(), directions.end(), V3D::compareMagnitude);
491

492
  double min_vol = min_d * min_d * min_d / 4.0;
493

494
495
  if (!FormUB_From_abc_Vectors(UB, directions, q_vectors, required_tolerance, min_vol)) {
    throw std::invalid_argument("Find_UB(): Could not form UB matrix from a,b,c vectors");
496
497
  }

498
  Matrix<double> temp_UB(3, 3, false);
499
  double fit_error = 0;
500
  if (q_vectors.size() >= 5) // repeatedly refine UB
501
502
503
  {
    std::vector<V3D> miller_ind;
    std::vector<V3D> indexed_qs;
504
505
    miller_ind.reserve(q_vectors.size());
    indexed_qs.reserve(q_vectors.size());
506

507
    GetIndexedPeaks(UB, q_vectors, required_tolerance, miller_ind, indexed_qs, fit_error);
508

509
    for (int counter = 0; counter < iterations; counter++) {
510
511
      try {
        fit_error = Optimize_UB(temp_UB, miller_ind, indexed_qs);
512
        UB = temp_UB;
513
        GetIndexedPeaks(UB, q_vectors, required_tolerance, miller_ind, indexed_qs, fit_error);
Lynch, Vickie's avatar
Lynch, Vickie committed
514
      } catch (std::exception &) {
515
        // failed to improve with all peaks, so just keep the UB we had
516
517
518
519
      }
    }
  }

520
  if (NiggliCell::MakeNiggliUB(UB, temp_UB))
521
522
    UB = temp_UB;

523
524
525
  return fit_error;
}

526
/**
527
  STATIC method Optimize_UB: Calculates the matrix that most nearly maps
528
529
530
  the specified hkl_vectors to the specified q_vectors.  The calculated
  UB minimizes the sum squared differences between UB*(h,k,l) and the
  corresponding (qx,qy,qz) for all of the specified hkl and Q vectors.
531
532
  The sum of the squares of the residual errors is returned.  This method is
  used to optimize the UB matrix once an initial indexing has been found.
533

534
  @param  UB           3x3 matrix that will be set to the UB matrix
535
  @param  hkl_vectors  std::vector of V3D objects that contains the
536
                       list of hkl values
537
  @param  q_vectors    std::vector of V3D objects that contains the list of
538
539
                       q_vectors that are indexed by the corresponding hkl
                       vectors.
540
  @param  sigabc      error in the crystal lattice parameter values if length
541
542
                      is at least 6. NOTE: Calculation of these errors is based
  on
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
                      SCD FORTRAN code base at IPNS. Contributors to the least
                      squares application(1979) are J.Marc Overhage, G.Anderson,
                      P. C. W. Leung, R. G. Teller, and  A. J. Schultz
  NOTE: The number of hkl_vectors and q_vectors must be the same, and must
        be at least 3.

  @return  This will return the sum of the squares of the residual differences
           between the Q vectors provided and the UB*hkl values, in
           reciprocal space.

  @throws  std::invalid_argument exception if there are not at least 3
                                 hkl and q vectors, or if the numbers of
                                 hkl and q vectors are not the same, or if
                                 the UB matrix is not a 3x3 matrix.

  @throws  std::runtime_error    exception if the QR factorization fails or
                                 the UB matrix can't be calculated or if
                                 UB is a singular matrix.
*/
562
double IndexingUtils::Optimize_UB(DblMatrix &UB, const std::vector<V3D> &hkl_vectors, const std::vector<V3D> &q_vectors,
563
                                  std::vector<double> &sigabc) {
564
565
566
  double result = 0;
  result = Optimize_UB(UB, hkl_vectors, q_vectors);

567
  if (sigabc.size() < 6) {
568
569
    sigabc.clear();
    return result;
570
571
572
573
574
575
576
577
  } else
    for (int i = 0; i < 6; i++)
      sigabc[i] = 0.0;

  size_t nDOF = 3 * (hkl_vectors.size() - 3);
  DblMatrix HKLTHKL(3, 3);
  for (int r = 0; r < 3; r++)
    for (int c = 0; c < 3; c++)
578
      for (const auto &hkl_vector : hkl_vectors) {
579
        HKLTHKL[r][c] += hkl_vector[r] * hkl_vector[c]; // rounded??? to nearest integer
580
      }
581
582
583

  HKLTHKL.Invert();

584
585
  double SMALL = 1.525878906E-5;
  Matrix<double> derivs(3, 7);
586
  std::vector<double> latOrig, latNew;
587
  GetLatticeParameters(UB, latOrig);
588

589
590
  for (int r = 0; r < 3; r++) {
    for (int c = 0; c < 3; c++) {
591
592

      UB[r][c] += SMALL;
593
      GetLatticeParameters(UB, latNew);
594
595
      UB[r][c] -= SMALL;

596
      for (size_t l = 0; l < 7; l++)
597
        derivs[c][l] = (latNew[l] - latOrig[l]) / SMALL;
598
599
    }

600
    for (size_t l = 0; l < std::min<size_t>(static_cast<size_t>(7), sigabc.size()); l++)
601
602
      for (int m = 0; m < 3; m++)
        for (int n = 0; n < 3; n++)
603
604
605
          sigabc[l] += (derivs[m][l] * HKLTHKL[m][n] * derivs[n][l]);
  }

606
  double delta = result / static_cast<double>(nDOF);
607

608
609
  for (size_t i = 0; i < std::min<size_t>(7, sigabc.size()); i++)
    sigabc[i] = sqrt(delta * sigabc[i]);
610
611
612
613

  return result;
}

614
615
616
617
618
619
620
621
622
/**
  STATIC method Optimize_UB: Calculates the matrix that most nearly maps
  the specified hkl_vectors to the specified q_vectors.  The calculated
  UB minimizes the sum squared differences between UB*(h,k,l) and the
  corresponding (qx,qy,qz) for all of the specified hkl and Q vectors.
  The sum of the squares of the residual errors is returned.  This method is
  used to optimize the UB matrix once an initial indexing has been found.

  @param  UB           3x3 matrix that will be set to the UB matrix
623
  @param  hkl_vectors  std::vector of V3D objects that contains the
624
                       list of hkl values
625
  @param  q_vectors    std::vector of V3D objects that contains the list of
626
627
                       q_vectors that are indexed by the corresponding hkl
                       vectors.
628
629
  NOTE: The number of hkl_vectors and q_vectors must be the same, and must
        be at least 3.
630
631
632
633
634

  @return  This will return the sum of the squares of the residual differences
           between the Q vectors provided and the UB*hkl values, in
           reciprocal space.

635
636
637
638
  @throws  std::invalid_argument exception if there are not at least 3
                                 hkl and q vectors, or if the numbers of
                                 hkl and q vectors are not the same, or if
                                 the UB matrix is not a 3x3 matrix.
639

640
  @throws  std::runtime_error    exception if the QR factorization fails or
641
                                 the UB matrix can't be calculated or if
642
                                 UB is a singular matrix.
643
*/
644
double IndexingUtils::Optimize_UB(DblMatrix &UB, const std::vector<V3D> &hkl_vectors,
645
646
647
                                  const std::vector<V3D> &q_vectors) {
  if (UB.numRows() != 3 || UB.numCols() != 3) {
    throw std::invalid_argument("Optimize_UB(): UB matrix NULL or not 3X3");
648
649
  }

650
  if (hkl_vectors.size() < 3) {
651
    throw std::invalid_argument("Optimize_UB(): Three or more indexed peaks needed to find UB");
652
  }
653

654
  if (hkl_vectors.size() != q_vectors.size()) {
655
    throw std::invalid_argument("Optimize_UB(): Number of hkl_vectors != number of q_vectors");
656
  }
657

658
659
  gsl_matrix *H_transpose = gsl_matrix_alloc(hkl_vectors.size(), 3);
  gsl_vector *tau = gsl_vector_alloc(3);
660
661

  double sum_sq_error = 0;
662
663
664
665
666
  // Make the H-transpose matrix from the
  // hkl vectors and form QR factorization
  for (size_t row = 0; row < hkl_vectors.size(); row++) {
    for (size_t col = 0; col < 3; col++) {
      gsl_matrix_set(H_transpose, row, col, (hkl_vectors[row])[col]);
667
668
    }
  }
669

670
  int returned_flag = gsl_linalg_QR_decomp(H_transpose, tau);
671

672
673
674
  if (returned_flag != 0) {
    gsl_matrix_free(H_transpose);
    gsl_vector_free(tau);
675
    throw std::runtime_error("Optimize_UB(): gsl QR_decomp failed, invalid hkl values");
676
  }
677
678
679
680
681
682
  // solve for each row of UB, using the
  // QR factorization of and accumulate the
  // sum of the squares of the residuals
  gsl_vector *UB_row = gsl_vector_alloc(3);
  gsl_vector *q = gsl_vector_alloc(q_vectors.size());
  gsl_vector *residual = gsl_vector_alloc(q_vectors.size());
683
684
685

  bool found_UB = true;

686
687
688
  for (size_t row = 0; row < 3; row++) {
    for (size_t i = 0; i < q_vectors.size(); i++) {
      gsl_vector_set(q, i, (q_vectors[i])[row] / (2.0 * M_PI));
689
    }
690

691
    returned_flag = gsl_linalg_QR_lssolve(H_transpose, tau, q, UB_row, residual);
692
    if (returned_flag != 0) {
693
      found_UB = false;
694
    }
695

696
697
    for (size_t i = 0; i < 3; i++) {
      double value = gsl_vector_get(UB_row, i);
698
      if (!std::isfinite(value))
699
700
        found_UB = false;
    }
701

702
    V3D row_values(gsl_vector_get(UB_row, 0), gsl_vector_get(UB_row, 1), gsl_vector_get(UB_row, 2));
703
    UB.setRow(row, row_values);
704

705
    for (size_t i = 0; i < q_vectors.size(); i++) {
706
      sum_sq_error += gsl_vector_get(residual, i) * gsl_vector_get(residual, i);
707
    }
708
709
  }

710
711
712
713
714
  gsl_matrix_free(H_transpose);
  gsl_vector_free(tau);
  gsl_vector_free(UB_row);
  gsl_vector_free(q);
  gsl_vector_free(residual);
715

716
  if (!found_UB) {
717
    throw std::runtime_error("Optimize_UB(): Failed to find UB, invalid hkl or Q values");
718
  }
719

720
721
  if (!CheckUB(UB)) {
    throw std::runtime_error("Optimize_UB(): The optimized UB is not valid");
722
  }
723

724
725
  return sum_sq_error;
}
726

727
728
729
730
731
732
733
/**
 STATIC method Optimize_UB: Calculates the matrix that most nearly maps
 the specified hkl_vectors to the specified q_vectors.  The calculated
 UB minimizes the sum squared differences between UB*(h,k,l) and the
 corresponding (qx,qy,qz) for all of the specified hkl and Q vectors.
 The sum of the squares of the residual errors is returned.  This method is
 used to optimize the UB matrix once an initial indexing has been found.
Lynch, Vickie's avatar
Lynch, Vickie committed
734

735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
 @param  UB           3x3 matrix that will be set to the UB matrix
 @param  ModUB        3x3 matrix that will be set to the ModUB matrix
 @param  hkl_vectors  std::vector of V3D objects that contains the
 list of hkl values
 @param  mnp_vectors  std::vector of V3D objects that contains the
 @param  ModDim       int value from 1 to 3, defines the number of dimensions
 of modulation.
 @param  q_vectors    std::vector of V3D objects that contains the list of
 q_vectors that are indexed by the corresponding hkl
 vectors.
 @param  sigabc       error in the crystal lattice parameter values if length
 is at least 6. NOTE: Calculation of these errors is based
 on
 SCD FORTRAN code base at IPNS. Contributors to the least
 squares application(1979) are J.Marc Overhage, G.Anderson,
 P. C. W. Leung, R. G. Teller, and  A. J. Schultz
 NOTE: The number of hkl_vectors and q_vectors must be the same, and must
 be at least 3.
 @param  sigq         error in the modulation vectors.
Lynch, Vickie's avatar
Lynch, Vickie committed
754
755


756
757
758
 @return  This will return the sum of the squares of the residual differences
 between the Q vectors provided and the UB*hkl values, in
 reciprocal space.
Lynch, Vickie's avatar
Lynch, Vickie committed
759

760
761
762
763
 @throws  std::invalid_argument exception if there are not at least 3
 hkl and q vectors, or if the numbers of
 hkl and q vectors are not the same, or if
 the UB matrix is not a 3x3 matrix.
Lynch, Vickie's avatar
Lynch, Vickie committed
764

765
766
767
768
769
 @throws  std::runtime_error    exception if the QR factorization fails or
 the UB matrix can't be calculated or if
 UB is a singular matrix.
 */

770
771
772
double IndexingUtils::Optimize_6dUB(DblMatrix &UB, DblMatrix &ModUB, const std::vector<V3D> &hkl_vectors,
                                    const std::vector<V3D> &mnp_vectors, const int &ModDim,
                                    const std::vector<V3D> &q_vectors, std::vector<double> &sigabc,
Lynch, Vickie's avatar
Lynch, Vickie committed
773
774
775
776
777
778
                                    std::vector<double> &sigq) {

  if (ModDim != 0 && ModDim != 1 && ModDim != 2 && ModDim != 3)
    throw std::invalid_argument("invalid Value for Modulation Dimension");

  double result = 0;
779
  result = Optimize_6dUB(UB, ModUB, hkl_vectors, mnp_vectors, ModDim, q_vectors);
Lynch, Vickie's avatar
Lynch, Vickie committed
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811

  if (sigabc.size() < static_cast<size_t>(6)) {
    sigabc.clear();
    return result;
  } else
    for (int i = 0; i < 6; i++)
      sigabc[i] = 0.0;

  size_t nDOF = 3 * (hkl_vectors.size() - 3);
  DblMatrix HKLTHKL(3, 3);
  for (int r = 0; r < 3; r++)
    for (int c = 0; c < 3; c++)
      for (const auto &hkl_vector : hkl_vectors) {
        HKLTHKL[r][c] += hkl_vector[r] * hkl_vector[c];
      }

  HKLTHKL.Invert();

  double SMALL = 1.525878906E-5;
  Matrix<double> derivs(3, 7);
  std::vector<double> latOrig, latNew;
  GetLatticeParameters(UB, latOrig);

  for (int r = 0; r < 3; r++) {
    for (int c = 0; c < 3; c++) {

      UB[r][c] += SMALL;
      GetLatticeParameters(UB, latNew);
      UB[r][c] -= SMALL;

      for (size_t l = 0; l < 7; l++)
        derivs[c][l] = (latNew[l] - latOrig[l]) / SMALL;
Lynch, Vickie's avatar
Lynch, Vickie committed
812
    }
Lynch, Vickie's avatar
Lynch, Vickie committed
813

814
    for (size_t l = 0; l < std::min<size_t>(static_cast<size_t>(7), sigabc.size()); l++)
Lynch, Vickie's avatar
Lynch, Vickie committed
815
816
817
818
      for (int m = 0; m < 3; m++)
        for (int n = 0; n < 3; n++)
          sigabc[l] += (derivs[m][l] * HKLTHKL[m][n] * derivs[n][l]);
  }
819

Lynch, Vickie's avatar
Lynch, Vickie committed
820
  double delta = result / static_cast<double>(nDOF);
821

Lynch, Vickie's avatar
Lynch, Vickie committed
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
  for (size_t i = 0; i < std::min<size_t>(7, sigabc.size()); i++)
    sigabc[i] = sqrt(delta * sigabc[i]);

  DblMatrix UBinv = UB;
  UBinv.Invert();

  if (sigq.size() < static_cast<size_t>(3)) {
    sigq.clear();
    return result;
  }

  else {
    sigq[0] = sqrt(delta) * latOrig[0];
    sigq[1] = sqrt(delta) * latOrig[1];
    sigq[2] = sqrt(delta) * latOrig[2];
  }

  return delta;
}
841

842
/**
Lynch, Vickie's avatar
Lynch, Vickie committed
843
844
845
846
847
848
849
850
851
852
853
854
 STATIC method Optimize_6dUB: Calculates the 6-dimensional matrix that most
 nearly maps the specified hkl_vectors and mnp_vectors to the specified
 q_vectors.  The calculated UB minimizes the sum squared differences between
 UB|ModUB*(h,k,l,m,n,p) and the corresponding (qx,qy,qz) for all of the
 specified hklmnp and Q vectors. The sum of the squares of the residual errors
 is returned.  This method is used to optimize the UB matrix and ModUB matrix
 once an initial indexing has been found.

 ModUB matrix is a 3x3 defines the modulation vectors in Q-sample. each colomn
 corresponds to a modulation vector.


855
856
857
858
859
860
861
862
863
864
865
 @param  UB           3x3 matrix that will be set to the UB matrix
 @param  ModUB        3x3 matrix that will be set to the ModUB matrix
 @param  hkl_vectors  std::vector of V3D objects that contains the
 list of hkl values
 @param  mnp_vectors  std::vector of V3D objects that contains the
 list of mnp values
 @param  ModDim       int value from 1 to 3, defines the number of dimensions
 of modulation.
 @param  q_vectors    std::vector of V3D objects that contains the list of
 q_vectors that are indexed by the corresponding hkl
 vectors.
Lynch, Vickie's avatar
Lynch, Vickie committed
866
867
868
 NOTE: The number of hkl_vectors and mnp_vectors and q_vectors must be the same,
 and must be at least 4.

869
870
871
 @return  This will return the sum of the squares of the residual differences
 between the Q vectors provided and the UB*hkl values, in
 reciprocal space.
Lynch, Vickie's avatar
Lynch, Vickie committed
872

873
874
875
876
877
 @throws  std::invalid_argument exception if there are not at least 3
 hkl and q vectors, or if the numbers of
 hkl and q vectors are not the same, or if
 the UB matrix is not a 3x3 matrix.
 or if the ModDim is not 1 or 2 or 3.
Lynch, Vickie's avatar
Lynch, Vickie committed
878

879
880
881
 @throws  std::runtime_error    exception if the QR factorization fails or
 the UB matrix can't be calculated or if
 UB is a singular matrix.
Lynch, Vickie's avatar
Lynch, Vickie committed
882

883
884
885
 Created by Shiyun Jin on 7/16/18.
 */

886
887
double IndexingUtils::Optimize_6dUB(DblMatrix &UB, DblMatrix</