SpatialModel.h

Go to the documentation of this file.

#ifndef INCLUDE_MOLASSEMBLER_DISTANCE_GEOMETRY_SPATIAL_MODEL_H
#define INCLUDE_MOLASSEMBLER_DISTANCE_GEOMETRY_SPATIAL_MODEL_H

#include "Molassembler/DistanceGeometry/DistanceBoundsMatrix.h"
#include "Molassembler/Molecule.h"
#include "Molassembler/Graph.h"
#include "Molassembler/StereopermutatorList.h"
#include "Molassembler/AtomStereopermutator.h"

#include <unordered_map>

namespace Scine {
namespace Molassembler {
namespace DistanceGeometry {

class SpatialModel {
public:
  struct ModelGraphWriter;

  template<unsigned size>
  using BoundsMapType = std::unordered_map<
    std::array<AtomIndex, size>,
    ValueBounds,
    boost::hash<
      std::array<AtomIndex, size>
    >
  >;

  using FixedPositionsMapType = std::unordered_map<AtomIndex, Utils::Position>;

  using BoundsMatrix = Eigen::MatrixXd;

  struct BoundsMatrixHelper {
    inline explicit BoundsMatrixHelper(AtomIndex size) : matrix(size, size) { matrix.setZero(); }

    void add(AtomIndex i, AtomIndex j, const ValueBounds& bounds);
    void addMap(const BoundsMapType<2>& boundsMap);

    ValueBounds get(AtomIndex i, AtomIndex j) const;

    Eigen::MatrixXd matrix;
  };


  static constexpr double bondRelativeVariance = 0.01;
  static constexpr double angleRelativeVariance = 0.02;
  static constexpr double dihedralAbsoluteVariance = M_PI / 90; // ~ 2°

  static constexpr ValueBounds angleClampBounds {0.0, M_PI};
  static const ValueBounds defaultDihedralBounds;

/* Static functions */
  static double modelDistance(AtomIndex i, AtomIndex j, const PrivateGraph& graph);
  static double modelDistance(const BondIndex& bond, const PrivateGraph& graph);

  static std::vector<BondIndex> cycleConsistingOfExactly(
    const std::vector<AtomIndex>& atoms,
    const PrivateGraph& graph
  );

  static boost::optional<ValueBounds> coneAngle(
    const std::vector<AtomIndex>& baseConstituents,
    const ValueBounds& coneHeightBounds,
    const Graph& graph
  );

  static double spiroCrossAngle(double alpha, double beta);

  static ValueBounds siteDistanceFromCenter(
    const std::vector<AtomIndex>& siteAtomList,
    AtomIndex placement,
    const Graph& graph
  );

  static ValueBounds makeBoundsFromCentralValue(
    double centralValue,
    double absoluteVariance
  );

  static double smallestCycleDistortionMultiplier(
    const AtomIndex i,
    const Cycles& cycles
  );

  /*
   * @brief Generates a matrix of atom-pairwise distance bounds
   *
   * @complexity{@math{O(P_2 + P_3 + P_4)} where @math{P_i} is the number of
   * distinct paths of length @math{i} in the graph. That should scale at least
   * linearly in the number of vertices.}
   *
   * @param N size of the system
   * @param fixedPositionbounds Distances to enforce due to fixed positions
   * @param bondBounds Distances to enforce for bonds
   * @param angleBounds Angle bounds to enforce on atom index triples
   * @param dihedralBounds Dihedral bounds to enforce on atom index quadruplet
   *
   * @return A matrix of atom-pairwise distance bounds. Lower bounds are in the
   * strict lower triangle of the matrix, upper bounds in the strict upper
   * triangle.
   */
  static BoundsMatrix makePairwiseBounds(
    unsigned N,
    const BoundsMapType<2>& fixedPositionBounds,
    const BoundsMapType<2>& bondBounds,
    const BoundsMapType<3>& angleBounds,
    const BoundsMapType<4>& dihedralBounds
  );

  static double siteCentralAngle(
    AtomIndex placement,
    const Shapes::Shape& shape,
    const RankingInformation& ranking,
    const AtomStereopermutator::ShapeMap& shapeVertexMap,
    const std::pair<SiteIndex, SiteIndex>& sites,
    const PrivateGraph& inner
  );

  static double siteCentralAngle(
    const AtomStereopermutator& permutator,
    const std::pair<SiteIndex, SiteIndex>& sites,
    const PrivateGraph& inner
  );

  static ValueBounds modelSiteAngleBounds(
    const AtomStereopermutator& permutator,
    const std::pair<SiteIndex, SiteIndex>& sites,
    const double looseningMultiplier,
    const PrivateGraph& inner
  );

  static ChiralConstraint makeChiralConstraint(
    const AtomStereopermutator::MinimalChiralConstraint& minimalConstraint,
    const AtomStereopermutator& permutator,
    const double looseningMultiplier
  );

  static ValueBounds clamp(
    ValueBounds bounds,
    const ValueBounds& clampBounds
  );

  static void checkFixedPositionsPreconditions(
    const Molecule& molecule,
    const Configuration& configuration
  );


  SpatialModel(
    const Molecule& molecule,
    const Configuration& configuration
  );


  void setBondBoundsIfEmpty(
    std::array<AtomIndex, 2> bondIndices,
    ValueBounds bounds
  );

  void setAngleBoundsIfEmpty(
    std::array<AtomIndex, 3> angleIndices,
    ValueBounds bounds
  );

  void setDihedralBoundsIfEmpty(
    std::array<AtomIndex, 4> dihedralIndices,
    ValueBounds bounds
  );

  void addAtomStereopermutatorInformation(
    const AtomStereopermutator& permutator,
    const PrivateGraph& graph,
    double looseningMultiplier,
    const std::unordered_map<AtomIndex, Utils::Position>& fixedAngstromPositions,
    bool forceChiralConstraintEmission
  );

  void addBondStereopermutatorInformation(
    const BondStereopermutator& permutator,
    const AtomStereopermutator& stereopermutatorA,
    const AtomStereopermutator& stereopermutatorB,
    double looseningMultiplier,
    const std::unordered_map<AtomIndex, Utils::Position>& fixedAngstromPositions
  );

  bool modelPartiallyFixedBond(
    const BondStereopermutator& permutator,
    const std::pair<const AtomStereopermutator&, const AtomStereopermutator&>& atomStereopermutators,
    const std::unordered_map<AtomIndex, Utils::Position>& fixedAngstromPositions
  );


  std::vector<ChiralConstraint> getChiralConstraints() const;

  std::vector<DihedralConstraint> getDihedralConstraints() const;

  BoundsMatrix makePairwiseBounds() const;

  std::string dumpGraphviz() const;

  void writeGraphviz(const std::string& filename) const;

private:
  // Molecule closure
  const Molecule& molecule_;

  BoundsMapType<2> constraints_;
  BoundsMapType<2> bondBounds_;
  BoundsMapType<3> angleBounds_;
  BoundsMapType<4> dihedralBounds_;

  std::vector<ChiralConstraint> chiralConstraints_;
  std::vector<DihedralConstraint> dihedralConstraints_;

  void addDefaultAngles_();

  void addDefaultDihedrals_();

  void instantiateMissingAtomStereopermutators_(Random::Engine& engine);

  void modelBondDistances_(
    const FixedPositionsMapType& fixedAngstromPositions,
    double looseningFactor
  );

  void modelFlatCycles_(
    const FixedPositionsMapType& fixedAngstromPositions,
    double looseningFactor
  );

  void modelSpirocenters_(
    const FixedPositionsMapType& fixedAngstromPositions
  );
};

} // namespace DistanceGeometry
} // namespace Molassembler
} // namespace Scine

#endif