SupercellEnumerator.hh

Go to the documentation of this file.

#ifndef SupercellEnumerator_HH
#define SupercellEnumerator_HH

#include "casm/external/Eigen/Dense"

#include "casm/misc/cloneable_ptr.hh"
#include "casm/symmetry/SymGroup.hh"
#include "casm/crystallography/Lattice.hh"
#include "casm/crystallography/BasicStructure.hh"
#include "casm/crystallography/Site.hh"
#include "casm/casm_io/Log.hh"

#include "casm/container/Counter.hh"

namespace CASM {

  class PrimClex;

  class ScelEnumProps {

  public:

    typedef long size_type;

    ScelEnumProps(size_type begin_volume,
                  size_type end_volume,
                  std::string dirs = "abc",
                  Eigen::Matrix3i generating_matrix = Eigen::Matrix3i::Identity()) :
      m_begin_volume(begin_volume),
      m_end_volume(end_volume),
      m_dims(dirs.size()),
      m_dirs(dirs) {

      for(int i = 0; i < m_dirs.size(); i++) {
        if(m_dirs[i] != 'a' && m_dirs[i] != 'b' && m_dirs[i] != 'c') {
          std::string msg = "Error constructing ScelEnumProps: an element of dirs != 'a', 'b', or 'c'";
          default_err_log().error("Constructing ScelEnumProps");
          default_err_log() << msg << "\n" << std::endl;
          throw std::invalid_argument(msg);
        }
      }

      // add missing directions to 'dirs',
      // then permute 'G' so that the specified 'dirs' are first
      while(m_dirs.size() != 3) {
        if(std::find(m_dirs.begin(), m_dirs.end(), 'a') == m_dirs.end()) {
          m_dirs.push_back('a');
        }
        if(std::find(m_dirs.begin(), m_dirs.end(), 'b') == m_dirs.end()) {
          m_dirs.push_back('b');
        }
        if(std::find(m_dirs.begin(), m_dirs.end(), 'c') == m_dirs.end()) {
          m_dirs.push_back('c');
        }
      }
      Eigen::Matrix3i P = Eigen::Matrix3i::Zero();
      P(m_dirs[0] - 'a', 0) = 1;
      P(m_dirs[1] - 'a', 1) = 1;
      P(m_dirs[2] - 'a', 2) = 1;

      m_gen_mat = generating_matrix * P;
    }


    size_type begin_volume() const {
      return m_begin_volume;
    }

    size_type end_volume() const {
      return m_end_volume;
    }

    int dims() const {
      return m_dims;
    }

    std::string dirs() const {
      return m_dirs;
    }

    Eigen::Matrix3i generating_matrix() const {
      return m_gen_mat;
    }

  private:

    size_type m_begin_volume;
    size_type m_end_volume;
    int m_dims;
    std::string m_dirs;
    Eigen::Matrix3i m_gen_mat;

  };

  Eigen::Matrix3i make_unit_cell(PrimClex &primclex, const jsonParser &json);

  ScelEnumProps make_scel_enum_props(PrimClex &primclex, const jsonParser &input);

  class HermiteCounter {
  public:
    typedef Eigen::VectorXi::Scalar value_type;
    typedef CASM::Index Index;

    //HermiteCounter() {};

    HermiteCounter(int init_determinant, int init_dim);

    //You probably will never need these. They're just here for testing more than anything.
    //Either way, they're safe to call.
    Index position() const;
    Eigen::VectorXi diagonal() const;
    //value_type low() const;
    //value_type high() const;

    Eigen::MatrixXi current() const;

    value_type determinant() const;

    Index dim() const;

    void reset_current();

    void next_determinant();

    void jump_to_determinant(value_type new_det);

    HermiteCounter &operator++();

    Eigen::MatrixXi operator()() const;


  private:

    Index m_pos;

    Eigen::VectorXi m_diagonal;

    EigenVectorXiCounter m_upper_tri;

    Index _increment_diagonal();

  };

  namespace HermiteCounter_impl {
    HermiteCounter::Index _spill_factor(Eigen::VectorXi &diag, HermiteCounter::Index position, HermiteCounter::value_type attempt);

    HermiteCounter::Index next_spill_position(Eigen::VectorXi &diag, HermiteCounter::Index position);

    HermiteCounter::Index upper_size(HermiteCounter::Index init_dim);

    EigenVectorXiCounter _upper_tri_counter(const Eigen::VectorXi &current_diag);

    Eigen::MatrixXi _zip_matrix(const Eigen::VectorXi &current_diag, const Eigen::VectorXi &current_upper_tri);

    Eigen::MatrixXi _expand_dims_old(const Eigen::MatrixXi &hermit_mat, const Eigen::VectorXi &active_dims);

    Eigen::MatrixXi _expand_dims(const Eigen::MatrixXi &H, const Eigen::MatrixXi &G);

    Eigen::VectorXi _canonical_unroll(const Eigen::MatrixXi &hermit_mat);

    bool _canonical_compare(const Eigen::MatrixXi &H0, const Eigen::MatrixXi &H1);
  }

  //******************************************************************************************************************//

  template <typename UnitType> class SupercellEnumerator;

  template <typename UnitType>
  class SupercellIterator {

  public:

    EIGEN_MAKE_ALIGNED_OPERATOR_NEW

    typedef std::forward_iterator_tag iterator_category;
    typedef int difference_type;
    typedef UnitType value_type;
    typedef const UnitType &reference;
    typedef const UnitType *pointer;

    SupercellIterator<UnitType>() {}

    SupercellIterator<UnitType>(const SupercellEnumerator<UnitType> &enumerator,
                                int volume,
                                int dims);

    //required for all iterators
    //SupercellIterator<UnitType>(const SupercellIterator<UnitType> &B);

    //required for all iterators?
    SupercellIterator<UnitType> &operator=(const SupercellIterator<UnitType> &B);

    bool operator==(const SupercellIterator<UnitType> &B) const;

    bool operator!=(const SupercellIterator<UnitType> &B) const;

    reference operator*() const;

    pointer operator->() const;

    Eigen::Matrix3i matrix() const;

    HermiteCounter::value_type volume() const;

    const SupercellEnumerator<UnitType> &enumerator() const;

    SupercellIterator<UnitType> &operator++();

    //SupercellIterator<UnitType> operator++(int);  TODO


  private:

    void _increment();

    //bool _is_canonical() const;

    void _try_increment();

    void _update_super();


    mutable bool m_super_updated;

    const SupercellEnumerator<UnitType> *m_enum;

    notstd::cloneable_ptr<HermiteCounter> m_current;

    mutable UnitType m_super;

    std::vector<Eigen::Matrix3i> m_canon_hist;

  };


  template<typename UnitType>
  class SupercellEnumerator {

  public:

    typedef long int size_type;

    typedef SupercellIterator<UnitType> const_iterator;

    SupercellEnumerator(UnitType unit,
                        const ScelEnumProps &enum_props,
                        double tol);

    SupercellEnumerator(UnitType unit,
                        const SymGroup &point_grp,
                        const ScelEnumProps &enum_props);


    const UnitType &unit() const;

    const Lattice &lattice() const;

    const SymGroup &point_group() const;

    void begin_volume(size_type _begin_volume);

    size_type begin_volume() const;

    void end_volume(size_type _end_volume);

    size_type end_volume() const;

    const Eigen::Matrix3i &gen_mat() const;

    int dimension() const;

    const_iterator begin() const;

    const_iterator end() const;

    const_iterator cbegin() const;

    const_iterator cend() const;

    const_iterator citerator(size_type volume) const;

  private:

    const UnitType m_unit;

    Lattice m_lat;

    SymGroup m_point_group;         //factor group...?

    const int m_begin_volume;

    const int m_end_volume;

    const Eigen::Matrix3i m_gen_mat;

    const int m_dims;

  };


  template<typename UnitType>
  Eigen::Matrix3i enforce_min_volume(
    const UnitType &unit,
    const Eigen::Matrix3i &T,
    const SymGroup &point_grp,
    Index volume,
    bool fix_shape = false);



  template<typename UnitType>
  SupercellIterator<UnitType>::SupercellIterator(const SupercellEnumerator<UnitType> &enumerator,
                                                 int volume,
                                                 int dims):
    m_super_updated(false),
    m_enum(&enumerator),
    m_current(notstd::make_cloneable<HermiteCounter>(volume, dims)) {
    if(enumerator.begin_volume() > enumerator.end_volume()) {
      throw std::runtime_error("The beginning volume of the SupercellEnumerator cannot be greater than the end volume!");
    }

    if(dims < 1) {
      throw std::runtime_error("Dimensions to count over must be greater than 0!");
    }
  }

  /*
  template<typename UnitType>
  SupercellIterator<UnitType>::SupercellIterator(const SupercellIterator &B)
  {
      *this = B;
  };
  */

  template<typename UnitType>
  SupercellIterator<UnitType> &SupercellIterator<UnitType>::operator=(const SupercellIterator &B) {
    m_enum = B.m_enum;
    m_current = B.m_current;
    m_super_updated = false;

    m_canon_hist = B.m_canon_hist;
    return *this;
  }

  template<typename UnitType>
  bool SupercellIterator<UnitType>::operator==(const SupercellIterator &B) const {
    return (m_enum == B.m_enum) && (matrix() - B.matrix()).isZero();
  }

  template<typename UnitType>
  bool SupercellIterator<UnitType>::operator!=(const SupercellIterator &B) const {
    return !(*this == B);
  }

  template<typename UnitType>
  typename SupercellIterator<UnitType>::reference SupercellIterator<UnitType>::operator*() const {
    if(!m_super_updated) {
      m_super = make_supercell(m_enum->unit(), matrix());
      m_super_updated = true;
    }
    return m_super;
  }

  template<typename UnitType>
  typename SupercellIterator<UnitType>::pointer SupercellIterator<UnitType>::operator->() const {
    if(!m_super_updated) {
      m_super = make_supercell(m_enum->unit(), matrix());
      m_super_updated = true;
    }
    return &m_super;
  }

  template<typename UnitType>
  HermiteCounter::value_type SupercellIterator<UnitType>::volume() const {
    return m_current->determinant();
  }

  template<typename UnitType>
  Eigen::Matrix3i SupercellIterator<UnitType>::matrix() const {
    Eigen::Matrix3i expanded = HermiteCounter_impl::_expand_dims((*m_current)(), m_enum->gen_mat());
    return canonical_hnf(expanded, m_enum->point_group(), m_enum->lattice());
  }

  template<typename UnitType>
  const SupercellEnumerator<UnitType> &SupercellIterator<UnitType>::enumerator() const {
    return *m_enum;
  }

  // prefix
  template<typename UnitType>
  SupercellIterator<UnitType> &SupercellIterator<UnitType>::operator++() {
    _increment();
    return *this;
  }

  /*
  // postfix
  template<typename UnitType>
  SupercellIterator<UnitType> SupercellIterator<UnitType>::operator++(int)
  {
      SupercellIterator result(*this);
      _increment();
      return result;
  }
  */


  template<typename UnitType>
  void SupercellIterator<UnitType>::_increment() {
    m_canon_hist.push_back(matrix());
    HermiteCounter::value_type last_determinant = m_current->determinant();
    ++(*m_current);

    if(last_determinant != m_current->determinant()) {
      m_canon_hist.clear();
    }

    while(std::find(m_canon_hist.begin(), m_canon_hist.end(), matrix()) != m_canon_hist.end()) {
      ++(*m_current);
    }

    m_super_updated = false;
  }

  //********************************************************************************************************//

  template<typename UnitType>
  const UnitType &SupercellEnumerator<UnitType>::unit() const {
    return m_unit;
  }

  template<typename UnitType>
  const Lattice &SupercellEnumerator<UnitType>::lattice() const {
    return m_lat;
  }

  template<typename UnitType>
  const SymGroup &SupercellEnumerator<UnitType>::point_group() const {
    return m_point_group;
  }

  template<typename UnitType>
  const Eigen::Matrix3i &SupercellEnumerator<UnitType>::gen_mat() const {
    return m_gen_mat;
  }

  template<typename UnitType>
  int SupercellEnumerator<UnitType>::dimension() const {
    return m_dims;
  }

  /*
  template<typename UnitType>
  void SupercellEnumerator<UnitType>::begin_volume(size_type _begin_volume)
  {
      m_begin_volume = _begin_volume;
  }
  */

  template<typename UnitType>
  typename SupercellEnumerator<UnitType>::size_type SupercellEnumerator<UnitType>::begin_volume() const {
    return m_begin_volume;
  }

  /*
  template<typename UnitType>
  void SupercellEnumerator<UnitType>::end_volume(size_type _end_volume)
  {
      m_end_volume = _end_volume;
  }
  */

  template<typename UnitType>
  typename SupercellEnumerator<UnitType>::size_type SupercellEnumerator<UnitType>::end_volume() const {
    return m_end_volume;
  }

  template<typename UnitType>
  typename SupercellEnumerator<UnitType>::const_iterator SupercellEnumerator<UnitType>::begin() const {
    return const_iterator(*this, m_begin_volume, dimension());
  }

  template<typename UnitType>
  typename SupercellEnumerator<UnitType>::const_iterator SupercellEnumerator<UnitType>::end() const {
    return const_iterator(*this, m_end_volume, dimension());
  }

  template<typename UnitType>
  typename SupercellEnumerator<UnitType>::const_iterator SupercellEnumerator<UnitType>::cbegin() const {
    return const_iterator(*this, m_begin_volume, dimension());
  }

  template<typename UnitType>
  typename SupercellEnumerator<UnitType>::const_iterator SupercellEnumerator<UnitType>::cend() const {
    return const_iterator(*this, m_end_volume, dimension());
  }

  template<typename UnitType>
  typename SupercellEnumerator<UnitType>::const_iterator SupercellEnumerator<UnitType>::citerator(size_type volume) const {
    return SupercellIterator<UnitType>(*this, volume, dimension());
  }


  // declare specializations for Lattice

  template<>
  SupercellEnumerator<Lattice>::SupercellEnumerator(Lattice unit,
                                                    const ScelEnumProps &enum_props,
                                                    double tol);

  template<>
  SupercellEnumerator<Lattice>::SupercellEnumerator(Lattice unit,
                                                    const SymGroup &point_grp,
                                                    const ScelEnumProps &enum_props);

  template<>
  Eigen::Matrix3i enforce_min_volume<Lattice>(
    const Lattice &unit,
    const Eigen::Matrix3i &T,
    const SymGroup &point_grp,
    Index volume,
    bool fix_shape);


  Eigen::Matrix3i canonical_hnf(const Eigen::Matrix3i &T, const SymGroup &effective_pg, const Lattice &ref_lattice);

}

#endif