Hull.cc

Go to the documentation of this file.

#include "casm/hull/Hull.hh"
 
#include <iterator>
 
#include "casm/casm_io/dataformatter/DataFormatter_impl.hh"
#include "casm/clex/PrimClex.hh"
#include "casm/database/ConfigDatabase.hh"
#include "casm/external/Eigen/Dense"
#include "casm/external/qhull/libqhullcpp/QhullFacetList.h"
#include "casm/external/qhull/libqhullcpp/QhullVertexSet.h"
#include "casm/misc/PCA.hh"
 
namespace CASM {
 
template class DataFormatter<std::pair<std::string, std::pair<bool, bool>>>;
 
namespace Hull_impl {
 
void _validate_input(const DB::Selection<Configuration> &selection,
                     const Hull::CompCalculator &comp_calculator,
                     const Hull::EnergyCalculator &energy_calculator);
}  // namespace Hull_impl
 
Hull::Hull(const DB::Selection<Configuration> &_selection,
           const CompCalculator &_comp_calculator,
           const EnergyCalculator &_energy_calculator,
           double _singular_value_tol, double _bottom_facet_tol)
    : m_selection(_selection),
      m_comp_calculator(_comp_calculator.clone()),
      m_energy_calculator(_energy_calculator.clone()) {
  Hull_impl::_validate_input(m_selection, *m_comp_calculator,
                             *m_energy_calculator);
 
  // get the number of configurations and compositions (full composition space)
  Index Nselected = m_selection.selected_size();
  Index Ncomp = composition(*m_selection.selected().begin()).size();
 
  // generate initial set of points (col vector matrix)
  Eigen::MatrixXd mat(Ncomp + 1, Nselected);
 
  Index i = 0;
  for (const auto &config : m_selection.selected()) {
    mat.block(0, i, Ncomp, 1) = composition(config);
    mat(Ncomp, i) = energy(config);
    ++i;
  }
 
  // principal component analysis to get rotation matrix
  PCA pca(mat.topRows(Ncomp), _singular_value_tol);
  m_reduce = pad(pca.reduce(), 1);
 
  Eigen::MatrixXd reduced_mat = m_reduce * mat;
 
  // Construct Qhull PointCoordiantes, with correct dimension
  orgQhull::PointCoordinates points(reduced_mat.rows(), "");
 
  // use matrix data to set point coordinates
  points.append(reduced_mat.size(), reduced_mat.data());
 
  // calculate hull
  std::string qh_command = "";
  m_hull.runQhull(points.comment().c_str(), points.dimension(), points.count(),
                  &*points.coordinates(), qh_command.c_str());
 
  // check for errors
  if (m_hull.hasQhullMessage()) {
    std::cerr << "\nQhull message:\n" << m_hull.qhullMessage();
    m_hull.clearQhullMessage();
    throw std::runtime_error("Qhull Error (see message)");
  }
 
  // want to find minimum distance from point to a bottom facet
 
  // V: any vector from point to facet
  // N: outward unit normal of facet
  // D: unit 'down' direction (by convention always [0, 0, ..., -1])
  //
  // distance to facet along normal: V.dot(N) = a (Qhull reports this as a
  // negative number) proj of D onto N: D.dot(N) = b (if b is positive, then
  // this is 'bottom' facet) distance along D to facet: a/b
 
  // collect bottom facets (along with 'b')
  {
    auto begin = m_hull.facetList().begin();
    auto end = m_hull.facetList().end();
    int dim = m_hull.dimension();
    double b;
    Eigen::Map<const Eigen::VectorXd> outnorm((*begin).hyperplane().begin(),
                                              dim);
 
    for (auto facet_it = begin; facet_it != end; ++facet_it) {
      // reset location of outnorm data
      new (&outnorm) Eigen::Map<const Eigen::VectorXd>(
          (*facet_it).hyperplane().begin(), dim);
      b = -outnorm(dim - 1);
 
      if (b > _bottom_facet_tol) {
        m_bottom_facets.push_back(std::make_pair(*facet_it, b));
      }
    }
  }
 
  // find set of bottom vertices using bottom facets
  // use pointers for the set comparison
  for (auto facet_it = m_bottom_facets.begin();
       facet_it != m_bottom_facets.end(); ++facet_it) {
    orgQhull::QhullVertexSet vertices = (*facet_it).first.vertices();
    for (auto vertex_it = vertices.begin(); vertex_it != vertices.end();
         ++vertex_it) {
      m_bottom_vertices.insert(*vertex_it);
    }
  }
}
 
const orgQhull::Qhull &Hull::data() const { return m_hull; }
 
const Eigen::MatrixXd &Hull::reduce() const { return m_reduce; }
 
const Eigen::Transpose<const Eigen::MatrixXd> Hull::expand() const {
  return m_reduce.transpose();
}
 
const Configuration &Hull::configuration(
    const orgQhull::QhullPoint &point) const {
  auto it = m_selection.selected().begin();
  std::advance(it, point.id());
  return *it;
}
 
const Configuration &Hull::configuration(
    const orgQhull::QhullVertex &vertex) const {
  return configuration(vertex.point());
}
 
Eigen::VectorXd Hull::mu(const orgQhull::QhullFacet facet) const {
  // mu_i = dG/dx_i
  // hyperplane normal = [dx_i, dx_j, ... dG] (normalized)
 
  int dim = m_hull.dimension();
  Eigen::Map<const Eigen::VectorXd> reduced_hyperplane(
      facet.hyperplane().begin(), dim);
 
  Eigen::VectorXd _mu =
      Eigen::VectorXd::Constant(dim - 1, reduced_hyperplane(dim - 1));
 
  return _mu.cwiseQuotient(expand().leftCols(dim - 1) *
                           reduced_hyperplane.head(dim - 1));
}
 
const Configuration &Hull::groundstate(const Eigen::VectorXd &mu) const {
  // mu_i = dG/dx_i
  // hyperplane normal = [dx_i, dx_j, ... dG] (normalized)
 
  // normal vector to hyperplane
  Eigen::VectorXd hyperplane(mu.size() + 1);
  hyperplane.head(mu.size()) = mu.cwiseInverse();
  hyperplane(mu.size()) = 1.0;
  hyperplane.normalize();
 
  Eigen::VectorXd reduced_hyperplane = reduce() * hyperplane;
 
  double max_offset = std::numeric_limits<double>::min();
  orgQhull::QhullVertex _groundstate = *m_bottom_vertices.begin();
 
  // m_bottom_vertices: a std::set<orgQhull::QhullVertex*>
  for (auto it = m_bottom_vertices.begin(); it != m_bottom_vertices.end();
       ++it) {
    Eigen::Map<const Eigen::VectorXd> reduced_vertex_point(it->point().begin(),
                                                           it->point().size());
 
    double offset = reduced_hyperplane.dot(reduced_vertex_point);
 
    if (offset > max_offset) {
      max_offset = offset;
      _groundstate = *it;
    }
  }
 
  return configuration(_groundstate);
}
 
double Hull::energy(const Configuration &config) const {
  return (*m_energy_calculator)(config);
}
 
Eigen::VectorXd Hull::composition(const Configuration &config) const {
  return (*m_comp_calculator)(config);
}
 
Eigen::VectorXd Hull::point(const Configuration &config) const {
  Eigen::VectorXd _point(m_reduce.cols());
  _point.head(m_reduce.cols() - 1) = composition(config);
  _point.tail(1)(0) = energy(config);
  return _point;
}
 
Eigen::VectorXd Hull::reduced_point(const Configuration &config) const {
  return reduce() * point(config);
}
 
double Hull::dist_to_hull(const Configuration &config) const {
  return dist_to_hull(reduced_point(config));
}
 
double Hull::dist_to_hull(Eigen::VectorXd _reduced_point) const {
  // want to find minimum distance from point to a bottom facet
 
  // V: any vector from point to facet
  // N: outward unit normal of facet
  // D: unit 'down' direction (by convention always [0, 0, ..., -1])
  //
  // distance to facet along normal: V.dot(N) = a (Qhull reports this as a
  // negative number) proj of D onto N: D.dot(N) = b (if b is positive, then
  // this is 'bottom' facet) distance along D to facet: a/b
 
  orgQhull::QhullPoint qpoint(_reduced_point.size(), _reduced_point.data());
 
  double a, b;
  double dist_to_hull = std::numeric_limits<double>::max();
 
  for (auto facet_it = m_bottom_facets.begin();
       facet_it != m_bottom_facets.end(); ++facet_it) {
    // Qhull distance is negative for internal points
    a = -(facet_it->first).distance(qpoint);
    b = facet_it->second;
    if (a / b < dist_to_hull) {
      dist_to_hull = a / b;
    }
  }
 
  return dist_to_hull;
}
 
namespace Hull_impl {
 
void _validate_input(const DB::Selection<Configuration> &selection,
                     const Hull::CompCalculator &comp_calculator,
                     const Hull::EnergyCalculator &energy_calculator) {
  typedef std::map<std::string, std::pair<bool, bool>> CheckMap;
 
  CheckMap invalid_data;
  for (const auto &config : selection.selected()) {
    if (!comp_calculator.validate(config) ||
        !energy_calculator.validate(config)) {
      invalid_data[config.name()] = std::make_pair(
          comp_calculator.validate(config), energy_calculator.validate(config));
    }
  }
 
  if (invalid_data.size()) {
    std::cerr << "Invalid data for hull construction:\n";
    typedef GenericDatumFormatter<std::string, CheckMap::value_type> Formatter;
 
    DataFormatter<CheckMap::value_type> f(
        Formatter("configname", "",
                  [](CheckMap::value_type v) { return v.first; }),
        Formatter("composition", "",
                  [](CheckMap::value_type v) {
                    return v.second.first ? "OK" : "invalid";
                  }),
        Formatter("energy", "", [](CheckMap::value_type v) {
          return v.second.second ? "OK" : "invalid";
        }));
 
    std::cerr << f(invalid_data.begin(), invalid_data.end()) << "\n";
 
    std::stringstream ss;
    ss << "Error in Hull(): Invalid composition or energy data. \n"
          "Make sure you have set composition axes, all selected "
          "configurations\n"
          "have calculation results, and you have set your chemical reference.";
    throw std::runtime_error(ss.str());
  }
}
 
}  // namespace Hull_impl
 
}  // namespace CASM