ScalarUtils.cpp

Go to the documentation of this file.

// Mantid Repository : https://github.com/mantidproject/mantid
//
// Copyright © 2018 ISIS Rutherford Appleton Laboratory UKRI,
//   NScD Oak Ridge National Laboratory, European Spallation Source,
//   Institut Laue - Langevin & CSNS, Institute of High Energy Physics, CAS
// SPDX - License - Identifier: GPL - 3.0 +
/* File: ScalarUtils.cpp */
 
#include "MantidGeometry/Crystal/ScalarUtils.h"
#include "MantidGeometry/Crystal/IndexingUtils.h"
#include "MantidGeometry/Crystal/OrientedLattice.h"
#include "MantidGeometry/Crystal/ReducedCell.h"
#include <algorithm>
#include <stdexcept>
 
using namespace Mantid::Geometry;
using Mantid::Kernel::DblMatrix;
using Mantid::Kernel::Matrix;
using Mantid::Kernel::V3D;
 
static const std::string BRAVAIS_TYPE[15] = {
    ReducedCell::CUBIC(), // F
    ReducedCell::CUBIC(), // I
    ReducedCell::CUBIC(), // P
 
    ReducedCell::HEXAGONAL(), // P
 
    ReducedCell::RHOMBOHEDRAL(), // R
 
    ReducedCell::TETRAGONAL(), // I
    ReducedCell::TETRAGONAL(), // P
 
    ReducedCell::ORTHORHOMBIC(), // F
    ReducedCell::ORTHORHOMBIC(), // I
    ReducedCell::ORTHORHOMBIC(), // C
    ReducedCell::ORTHORHOMBIC(), // P
 
    ReducedCell::MONOCLINIC(), // C
    ReducedCell::MONOCLINIC(), // I
    ReducedCell::MONOCLINIC(), // P
 
    ReducedCell::TRICLINIC() // P
};
 
static const std::string BRAVAIS_CENTERING[15] = {
    ReducedCell::F_CENTERED(), // cubic
    ReducedCell::I_CENTERED(), // cubic
    ReducedCell::P_CENTERED(), // cubic
 
    ReducedCell::P_CENTERED(), // hexagonal
 
    ReducedCell::R_CENTERED(), // rhombohedral
 
    ReducedCell::I_CENTERED(), // tetragonal
    ReducedCell::P_CENTERED(), // tetragonal
 
    ReducedCell::F_CENTERED(), // orthorhombic
    ReducedCell::I_CENTERED(), // orthorhombic
    ReducedCell::C_CENTERED(), // orthorhombic
    ReducedCell::P_CENTERED(), // orthorhombic
 
    ReducedCell::C_CENTERED(), // monoclinic
    ReducedCell::I_CENTERED(), // monoclinic
    ReducedCell::P_CENTERED(), // monoclinic
 
    ReducedCell::P_CENTERED() // triclinic
};
 
std::vector<ConventionalCell> ScalarUtils::GetCells(const DblMatrix &UB, bool best_only, bool allowPermutations) {
  std::vector<ConventionalCell> result;
 
  size_t num_lattices = 15;
  for (size_t i = 0; i < num_lattices; i++) {
    std::vector<ConventionalCell> temp = GetCells(UB, BRAVAIS_TYPE[i], BRAVAIS_CENTERING[i], allowPermutations);
    if (best_only) {
      ConventionalCell info = GetCellBestError(temp, true);
      temp.clear();
      temp.emplace_back(info);
    }
    for (auto &k : temp)
      AddIfBest(result, k);
  }
 
  return result;
}
 
std::vector<ConventionalCell> ScalarUtils::GetCells(const DblMatrix &UB, const std::string &cell_type,
                                                    const std::string &centering, bool allowPermutations) {
  std::vector<ConventionalCell> result;
 
  std::vector<DblMatrix> UB_list;
  if (allowPermutations) {
    double angle_tolerance = 2.0;
    double length_factor = 1.05;
    UB_list = GetRelatedUBs(UB, length_factor, angle_tolerance);
  } else {
    // Get exact form requested and not permutations
    UB_list.emplace_back(UB);
  }
 
  for (auto &k : UB_list) {
    std::vector<ConventionalCell> temp = GetCellsUBOnly(k, cell_type, centering, allowPermutations);
 
    for (auto &cell : temp)
      AddIfBest(result, cell);
  }
 
  return result;
}
 
std::vector<ConventionalCell> ScalarUtils::GetCellsUBOnly(const DblMatrix &UB, const std::string &cell_type,
                                                          const std::string &centering, bool allowPermutations) {
  std::vector<ConventionalCell> result;
 
  std::vector<double> lat_par;
  IndexingUtils::GetLatticeParameters(UB, lat_par);
 
  for (size_t i = 0; i <= ReducedCell::NUM_CELL_TYPES; i++) {
    ReducedCell rcell(i, lat_par[0], lat_par[1], lat_par[2], lat_par[3], lat_par[4], lat_par[5]);
 
    if (rcell.GetCentering() == centering && rcell.GetCellType() == cell_type) {
      ConventionalCell cell_info(UB, i, allowPermutations);
      result.emplace_back(cell_info);
    }
  }
 
  return result;
}
 
ConventionalCell ScalarUtils::GetCellForForm(const DblMatrix &UB, size_t form_num, bool allowPermutations) {
  ConventionalCell info(UB);
  ReducedCell form_0;
  ReducedCell form;
 
  double min_error = 1e20; // errors are usually < 10, so this is big enough
 
  std::vector<double> l_params;
  std::vector<DblMatrix> UB_list;
  if (allowPermutations) {
    double angle_tolerance = 2.0;
    double length_factor = 1.05;
    UB_list = GetRelatedUBs(UB, length_factor, angle_tolerance);
  } else {
    // Get exact form requested and not permutations
    UB_list.emplace_back(UB);
  }
  for (auto &ub : UB_list) {
    IndexingUtils::GetLatticeParameters(ub, l_params);
 
    form_0 = ReducedCell(0, l_params[0], l_params[1], l_params[2], l_params[3], l_params[4], l_params[5]);
 
    form = ReducedCell(form_num, l_params[0], l_params[1], l_params[2], l_params[3], l_params[4], l_params[5]);
 
    double error = form_0.WeightedDistance(form);
    if (error < min_error) {
      info = ConventionalCell(ub, form_num, allowPermutations);
      min_error = error;
    }
  }
 
  return info;
}
 
void ScalarUtils::RemoveHighErrorForms(std::vector<ConventionalCell> &list, double level) {
  if (list.empty()) // nothing to do
    return;
  const auto removeRangeBegin =
      std::remove_if(list.begin(), list.end(), [level](const auto &cell) { return cell.GetError() > level; });
  list.erase(removeRangeBegin, list.end());
}
 
ConventionalCell ScalarUtils::GetCellBestError(const std::vector<ConventionalCell> &list, bool use_triclinic) {
  if (list.empty()) {
    throw std::invalid_argument("GetCellBestError(): list is empty");
  }
 
  ConventionalCell info = list[0];
  double min_error = 1.0e20;
 
  bool min_found = false;
  for (const auto &cell : list) {
    std::string type = cell.GetCellType();
    double error = cell.GetError();
    if ((use_triclinic || type != ReducedCell::TRICLINIC()) && error < min_error) {
      info = cell;
      min_error = error;
      min_found = true;
    }
  }
 
  if (!min_found) {
    throw std::invalid_argument("GetCellBestError(): no allowed form with min error");
  }
 
  return info;
}
 
std::vector<DblMatrix> ScalarUtils::GetRelatedUBs(const DblMatrix &UB, double factor, double angle_tolerance) {
  std::vector<DblMatrix> result;
 
  V3D a, b, c, a_vec, b_vec, c_vec, // vectors for generating reflections
      m_a_vec, m_b_vec, m_c_vec;    // of pairs of sides
 
  V3D a_temp, b_temp, c_temp,       // vectors for generating handedness
      m_a_temp, m_b_temp, m_c_temp; // preserving permutations of sides
 
  OrientedLattice::GetABC(UB, a_vec, b_vec, c_vec);
 
  m_a_vec = a_vec * (-1.0);
  m_b_vec = b_vec * (-1.0);
  m_c_vec = c_vec * (-1.0);
  // make list of reflections of all pairs
  // of sides.  NOTE: These preserve the
  // ordering of magnitudes: |a|<=|b|<=|c|
  V3D reflections[4][3] = {
      {a_vec, b_vec, c_vec}, {m_a_vec, m_b_vec, c_vec}, {m_a_vec, b_vec, m_c_vec}, {a_vec, m_b_vec, m_c_vec}};
 
  // make list of the angles that are not
  // changed by each of the reflections.  IF
  // that angle is close to 90 degrees, then
  // we may need to switch between all angles
  // >= 90 and all angles < 90.  An angle
  // near 90 degrees may be mis-categorized
  // due to errors in the data.
  double alpha = b_vec.angle(c_vec) * 180.0 / M_PI;
  double beta = c_vec.angle(a_vec) * 180.0 / M_PI;
  double gamma = a_vec.angle(b_vec) * 180.0 / M_PI;
  double angles[4] = {90.0, gamma, beta, alpha};
 
  for (size_t row = 0; row < 4; row++) {
    if (fabs(angles[row] - 90.0) < angle_tolerance) // if nearly 90,
    {                                               // try related cell
      a_temp = reflections[row][0];                 // +cell <-> -cell
      b_temp = reflections[row][1];
      c_temp = reflections[row][2];
      // for each accepted reflection, try all
      // modified premutations that preserve the
      // handedness AND keep the cell edges
      // nearly ordered as a <= b <= c.
      m_a_temp = a_temp * (-1.0);
      m_b_temp = b_temp * (-1.0);
      m_c_temp = c_temp * (-1.0);
 
      V3D permutations[6][3] = {{a_temp, b_temp, c_temp},   {m_a_temp, c_temp, b_temp}, {b_temp, c_temp, a_temp},
                                {m_b_temp, a_temp, c_temp}, {c_temp, a_temp, b_temp},   {m_c_temp, b_temp, a_temp}};
 
      for (auto &permutation : permutations) {
        a = permutation[0];
        b = permutation[1];
        c = permutation[2];
        if (a.norm() <= factor * b.norm() && b.norm() <= factor * c.norm()) // could be Niggli within
        {                                                                   // experimental error
          Matrix<double> temp_UB(3, 3, false);
          OrientedLattice::GetUB(temp_UB, a, b, c);
          result.emplace_back(temp_UB);
        }
      }
    }
  }
 
  return result;
}
 
void ScalarUtils::AddIfBest(std::vector<ConventionalCell> &list, ConventionalCell &info) {
  size_t form_num = info.GetFormNum();
  double new_error = info.GetError();
  bool done = false;
  size_t i = 0;
  while (!done && i < list.size()) {
    if (list[i].GetFormNum() == form_num) // if found, replace if better
    {
      done = true;
      if (list[i].GetError() > new_error)
        list[i] = info;
    } else
      i++;
  }
 
  if (!done) // if never found, add to end of list
    list.emplace_back(info);
}