VMD.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 +
#include "MantidKernel/VMD.h"
#include "MantidKernel/Matrix.h"
#include "MantidKernel/StringTokenizer.h"
#include "MantidKernel/Strings.h"
#include "MantidKernel/Tolerance.h"
#include "MantidKernel/V3D.h"
 
#include <algorithm>
#include <cmath>
#include <cstddef>
#include <iterator>
#include <sstream>
#include <stdexcept>
 
using namespace Mantid::Kernel;
 
namespace Mantid::Kernel {
 
template <typename TYPE> VMDBase<TYPE>::VMDBase() : nd(1) {
  data = new TYPE[nd];
  for (size_t d = 0; d < nd; d++)
    data[d] = TYPE(0.0);
}
 
template <typename TYPE> VMDBase<TYPE>::VMDBase(size_t nd) : nd(nd) {
  if (nd == 0)
    throw std::invalid_argument("nd must be > 0");
  data = new TYPE[nd];
  for (size_t d = 0; d < nd; d++)
    data[d] = TYPE(0.0);
}
 
template <typename TYPE> VMDBase<TYPE>::VMDBase(double val0, double val1) : nd(2) {
  data = new TYPE[nd];
  data[0] = TYPE(val0);
  data[1] = TYPE(val1);
}
 
template <typename TYPE> VMDBase<TYPE>::VMDBase(double val0, double val1, double val2) : nd(3) {
  data = new TYPE[nd];
  data[0] = TYPE(val0);
  data[1] = TYPE(val1);
  data[2] = TYPE(val2);
}
 
template <typename TYPE> VMDBase<TYPE>::VMDBase(double val0, double val1, double val2, double val3) : nd(4) {
  data = new TYPE[nd];
  data[0] = TYPE(val0);
  data[1] = TYPE(val1);
  data[2] = TYPE(val2);
  data[3] = TYPE(val3);
}
 
template <typename TYPE>
VMDBase<TYPE>::VMDBase(double val0, double val1, double val2, double val3, double val4) : nd(5) {
  data = new TYPE[nd];
  data[0] = TYPE(val0);
  data[1] = TYPE(val1);
  data[2] = TYPE(val2);
  data[3] = TYPE(val3);
  data[4] = TYPE(val4);
}
 
template <typename TYPE>
VMDBase<TYPE>::VMDBase(double val0, double val1, double val2, double val3, double val4, double val5) : nd(6) {
  data = new TYPE[nd];
  data[0] = TYPE(val0);
  data[1] = TYPE(val1);
  data[2] = TYPE(val2);
  data[3] = TYPE(val3);
  data[4] = TYPE(val4);
  data[5] = TYPE(val5);
}
 
template <typename TYPE> VMDBase<TYPE>::VMDBase(const VMDBase &other) : nd(other.nd) {
  if (nd == 0)
    throw std::invalid_argument("nd must be > 0");
  data = new TYPE[nd];
  for (size_t d = 0; d < nd; d++)
    data[d] = other.data[d];
}
 
template <typename TYPE> VMDBase<TYPE> &VMDBase<TYPE>::operator=(const VMDBase &other) {
  if ((other.nd) != nd) {
    nd = other.nd;
    delete[] data;
    data = new TYPE[nd];
  }
  for (size_t d = 0; d < nd; d++)
    data[d] = other.data[d];
  return *this;
}
 
template <typename TYPE> VMDBase<TYPE>::VMDBase(VMDBase &&other) noexcept : nd(other.nd), data(other.data) {
  other.data = nullptr;
  other.nd = 0;
}
 
template <typename TYPE> VMDBase<TYPE> &VMDBase<TYPE>::operator=(VMDBase &&other) noexcept {
  if (this != &other) {
    this->nd = other.nd;
    other.nd = 0;
    delete[] this->data;
    this->data = other.data;
    other.data = nullptr;
  }
  return *this;
}
 
template <typename TYPE> VMDBase<TYPE>::VMDBase(size_t nd, const double *bareData) : nd(nd) {
  if (nd == 0)
    throw std::invalid_argument("nd must be > 0");
  data = new TYPE[nd];
  for (size_t d = 0; d < nd; d++)
    data[d] = TYPE(bareData[d]);
}
 
template <typename TYPE> VMDBase<TYPE>::VMDBase(size_t nd, const float *bareData) : nd(nd) {
  if (nd == 0)
    throw std::invalid_argument("nd must be > 0");
  data = new TYPE[nd];
  for (size_t d = 0; d < nd; d++)
    data[d] = TYPE(bareData[d]);
}
 
template <typename TYPE> VMDBase<TYPE>::VMDBase(const V3D &vector) : nd(3) {
  data = new TYPE[nd];
  for (size_t d = 0; d < nd; d++)
    data[d] = TYPE(vector[d]);
}
 
template <typename TYPE> VMDBase<TYPE>::VMDBase(const std::vector<double> &vector) : nd(vector.size()) {
  if (nd == 0)
    throw std::invalid_argument("nd must be > 0");
  data = new TYPE[nd];
  for (size_t d = 0; d < nd; d++)
    data[d] = TYPE(vector[d]);
}
 
template <typename TYPE> VMDBase<TYPE>::VMDBase(const std::vector<float> &vector) : nd(vector.size()) {
  if (nd == 0)
    throw std::invalid_argument("nd must be > 0");
  data = new TYPE[nd];
  for (size_t d = 0; d < nd; d++)
    data[d] = TYPE(vector[d]);
}
 
template <typename TYPE> VMDBase<TYPE>::VMDBase(const std::string &str) {
 
  StringTokenizer strs(str, ", ", StringTokenizer::TOK_IGNORE_EMPTY);
 
  std::vector<TYPE> vals;
  std::transform(strs.cbegin(), strs.cend(), std::back_inserter(vals), [](const std::string &token) {
    TYPE v;
    if (!Strings::convert(token, v))
      throw std::invalid_argument("VMDBase: Unable to convert the string '" + token + "' to a number.");
    return v;
  });
 
  nd = vals.size();
  if (nd == 0)
    throw std::invalid_argument("nd must be > 0");
  data = new TYPE[nd];
  std::copy(vals.cbegin(), vals.cend(), data);
}
 
template <typename TYPE> VMDBase<TYPE>::~VMDBase() { delete[] data; }
 
template <typename TYPE> size_t VMDBase<TYPE>::getNumDims() const { return nd; }
 
template <typename TYPE> size_t VMDBase<TYPE>::size() const { return nd; }
 
template <typename TYPE> const TYPE &VMDBase<TYPE>::operator[](const size_t index) const { return data[index]; }
 
template <typename TYPE> TYPE &VMDBase<TYPE>::operator[](const size_t index) { return data[index]; }
 
template <typename TYPE> const TYPE *VMDBase<TYPE>::getBareArray() const { return data; }
 
template <typename TYPE> std::string VMDBase<TYPE>::toString(const std::string &separator) const {
  std::ostringstream mess;
  for (size_t d = 0; d < nd; d++)
    mess << (d > 0 ? separator : "") << data[d];
  return mess.str();
}
 
template <typename TYPE> bool VMDBase<TYPE>::operator==(const VMDBase &v) const {
  if (v.nd != nd)
    return false;
  for (size_t d = 0; d < nd; d++)
    if ((std::fabs(data[d] - v.data[d]) > Tolerance))
      return false;
  return true;
}
 
template <typename TYPE> bool VMDBase<TYPE>::operator!=(const VMDBase &v) const { return !operator==(v); }
 
template <typename TYPE> VMDBase<TYPE> VMDBase<TYPE>::operator+(const VMDBase &v) const {
  VMDBase out(*this);
  out += v;
  return out;
}
 
template <typename TYPE> VMDBase<TYPE> &VMDBase<TYPE>::operator+=(const VMDBase &v) {
  if (v.nd != this->nd)
    throw std::runtime_error("Mismatch in number of dimensions in operation "
                             "between two VMDBase vectors.");
  for (size_t d = 0; d < nd; d++)
    data[d] += v.data[d];
  return *this;
}
 
template <typename TYPE> VMDBase<TYPE> VMDBase<TYPE>::operator-(const VMDBase &v) const {
  VMDBase out(*this);
  out -= v;
  return out;
}
 
template <typename TYPE> VMDBase<TYPE> &VMDBase<TYPE>::operator-=(const VMDBase &v) {
  if (v.nd != this->nd)
    throw std::runtime_error("Mismatch in number of dimensions in operation "
                             "between two VMDBase vectors.");
  for (size_t d = 0; d < nd; d++)
    data[d] -= v.data[d];
  return *this;
}
 
template <typename TYPE> VMDBase<TYPE> VMDBase<TYPE>::operator*(const VMDBase &v) const {
  VMDBase out(*this);
  out *= v;
  return out;
}
 
template <typename TYPE> VMDBase<TYPE> &VMDBase<TYPE>::operator*=(const VMDBase &v) {
  if (v.nd != this->nd)
    throw std::runtime_error("Mismatch in number of dimensions in operation "
                             "between two VMDBase vectors.");
  for (size_t d = 0; d < nd; d++)
    data[d] *= v.data[d];
  return *this;
}
 
template <typename TYPE> VMDBase<TYPE> VMDBase<TYPE>::operator/(const VMDBase &v) const {
  VMDBase out(*this);
  out /= v;
  return out;
}
 
template <typename TYPE> VMDBase<TYPE> &VMDBase<TYPE>::operator/=(const VMDBase &v) {
  if (v.nd != this->nd)
    throw std::runtime_error("Mismatch in number of dimensions in operation "
                             "between two VMDBase vectors.");
  for (size_t d = 0; d < nd; d++)
    data[d] /= v.data[d];
  return *this;
}
 
template <typename TYPE> VMDBase<TYPE> VMDBase<TYPE>::operator*(const double scalar) const {
  VMDBase out(*this);
  out *= scalar;
  return out;
}
 
template <typename TYPE> VMDBase<TYPE> &VMDBase<TYPE>::operator*=(const double scalar) {
  for (size_t d = 0; d < nd; d++)
    data[d] *= TYPE(scalar);
  return *this;
}
 
template <typename TYPE> VMDBase<TYPE> VMDBase<TYPE>::operator/(const double scalar) const {
  VMDBase out(*this);
  out /= scalar;
  return out;
}
 
template <typename TYPE> VMDBase<TYPE> &VMDBase<TYPE>::operator/=(const double scalar) {
  for (size_t d = 0; d < nd; d++)
    data[d] /= TYPE(scalar);
  return *this;
}
 
template <typename TYPE> TYPE VMDBase<TYPE>::scalar_prod(const VMDBase &v) const {
  TYPE out = 0;
  if (v.nd != this->nd)
    throw std::runtime_error("Mismatch in number of dimensions in operation "
                             "between two VMDBase vectors.");
  for (size_t d = 0; d < nd; d++)
    out += (data[d] * v.data[d]);
  return out;
}
 
template <typename TYPE> VMDBase<TYPE> VMDBase<TYPE>::cross_prod(const VMDBase &v) const {
  if (v.nd != this->nd)
    throw std::runtime_error("Mismatch in number of dimensions in operation "
                             "between two VMDBase vectors.");
  if (v.nd != 3)
    throw std::runtime_error("Cross product of vectors only works in 3 dimensions.");
  V3D a(data[0], data[1], data[2]);
  V3D b(v.data[0], v.data[1], v.data[2]);
  V3D c = a.cross_prod(b);
  VMDBase out(c);
  return out;
}
 
template <typename TYPE> TYPE VMDBase<TYPE>::length() const { return TYPE(std::sqrt(this->norm2())); }
 
template <typename TYPE> TYPE VMDBase<TYPE>::norm() const { return this->length(); }
 
template <typename TYPE> TYPE VMDBase<TYPE>::norm2() const { return this->scalar_prod(*this); }
 
template <typename TYPE> TYPE VMDBase<TYPE>::normalize() {
  TYPE length = this->length();
  for (size_t d = 0; d < nd; d++)
    data[d] /= length;
  return length;
}
 
template <typename TYPE> TYPE VMDBase<TYPE>::angle(const VMDBase &v) const {
  return TYPE(acos(this->scalar_prod(v) / (this->norm() * v.norm())));
}
 
//-------------------------------------------------------------------------------------------------
template <typename TYPE>
std::vector<VMDBase<TYPE>> VMDBase<TYPE>::makeVectorsOrthogonal(std::vector<VMDBase> &vectors) {
  if (vectors.size() != 2)
    throw std::runtime_error("VMDBase::makeVectorsOrthogonal(): Need 2 input vectors.");
  if (vectors[0].getNumDims() != 3 || vectors[1].getNumDims() != 3)
    throw std::runtime_error("VMDBase::makeVectorsOrthogonal(): Need 3D input vectors.");
  std::vector<V3D> in, out;
  for (size_t i = 0; i < vectors.size(); i++)
    in.emplace_back(vectors[i][0], vectors[i][1], vectors[i][2]);
  out = V3D::makeVectorsOrthogonal(in);
  std::vector<VMDBase> retVal;
  retVal.reserve(out.size());
  std::copy(std::make_move_iterator(out.begin()), std::make_move_iterator(out.end()), std::back_inserter(retVal));
  return retVal;
}
 
//-------------------------------------------------------------------------------------------------
template <typename TYPE> VMDBase<TYPE> VMDBase<TYPE>::getNormalVector(const std::vector<VMDBase<TYPE>> &vectors) {
  if (vectors.empty())
    throw std::invalid_argument("VMDBase::getNormalVector: Must give at least 1 vector");
  size_t nd = vectors[0].getNumDims();
  if (nd < 2)
    throw std::invalid_argument("VMDBase::getNormalVector: Must have at least 2 dimensions!");
  if (vectors.size() != nd - 1)
    throw std::invalid_argument("VMDBase::getNormalVector: Must have as many "
                                "N-1 vectors if there are N dimensions.");
  for (size_t i = 0; i < vectors.size(); i++)
    if (vectors[i].getNumDims() != nd)
      throw std::invalid_argument("VMDBase::getNormalVector: Inconsistent "
                                  "number of dimensions in the vectors given!");
 
  // Start the normal vector
  VMDBase normal = VMDBase(nd);
  TYPE sign = +1.0;
  for (size_t d = 0; d < nd; d++) {
    // Make the sub-determinant with the columns of every other dimension.
    Matrix<TYPE> mat(nd - 1, nd - 1);
    for (size_t row = 0; row < vectors.size(); row++) {
      VMDBase vec = vectors[row];
      size_t col = 0;
      for (size_t i = 0; i < nd; i++) {
        if (i != d) // Skip the column of this dimension
        {
          mat[row][col] = vec[i];
          col++;
        }
      }
    } // Building the matrix rows
 
    TYPE det = mat.determinant();
 
    // The determinant of the sub-matrix = the normal at that dimension
    normal[d] = sign * det;
 
    // Sign flips each time
    sign *= TYPE(-1.0);
  } // each dimension of the normal vector
 
  // Unity normal is better.
  double length = normal.normalize();
  if (length == 0)
    throw std::runtime_error("VMDBase::getNormalVector: 0-length normal found. "
                             "Are your vectors collinear?");
 
  return normal;
}
 
template class VMDBase<double>;
template class VMDBase<float>;
 
std::ostream &operator<<(std::ostream &os, const VMDBase<double> &v) {
  os << v.toString();
  return os;
}
 
std::ostream &operator<<(std::ostream &os, const VMDBase<float> &v) {
  os << v.toString();
  return os;
}
 
} // namespace Mantid::Kernel