DistanceBoundsMatrix.h

Go to the documentation of this file.

#ifndef INCLUDE_MOLASSEMBLER_DISTANCE_GEOMETRY_DISTANCE_BOUNDS_MATRIX_H
#define INCLUDE_MOLASSEMBLER_DISTANCE_GEOMETRY_DISTANCE_BOUNDS_MATRIX_H

#include <Eigen/Core>

#include "Molassembler/DistanceGeometry/DistanceGeometry.h"
#include "Molassembler/Modeling/AtomInfo.h"
#include "Molassembler/Conformers.h"

namespace Scine {
namespace Molassembler {

namespace Random {
class Engine;
} // namespace Random

namespace outcome = OUTCOME_V2_NAMESPACE;

namespace DistanceGeometry {

class DistanceBoundsMatrix {
public:
  static constexpr double defaultLower = 0.0;
  static constexpr double defaultUpper = 100.0;

  static inline double& lowerBound(Eigen::Ref<Eigen::MatrixXd> matrix, const AtomIndex i, const AtomIndex j) {
    if(i < j) {
      return matrix(j, i);
    }

    return matrix(i, j);
  }

  static inline double& upperBound(Eigen::Ref<Eigen::MatrixXd> matrix, const AtomIndex i, const AtomIndex j) {
    if(i < j) {
      return matrix(i, j);
    }

    return matrix(j, i);
  }

  static void smooth(Eigen::Ref<Eigen::MatrixXd> matrix);

  using BoundsMatrix = Eigen::MatrixXd;

  DistanceBoundsMatrix();

  explicit DistanceBoundsMatrix(long unsigned N);

  explicit DistanceBoundsMatrix(Eigen::MatrixXd matrix);

  template<typename PrivateGraph>
  DistanceBoundsMatrix(
    const PrivateGraph& inner,
    BoundsMatrix bounds
  ) : matrix_(std::move(bounds)) {
    // Populate the lower bounds if no explicit information is present
    const AtomIndex N = inner.N();
    for(AtomIndex i = 0; i < N - 1; ++i) {
      for(AtomIndex j = i + 1; j < N; ++j) {
        // i < j in all cases -> lower bound is (j, i), upper bound is (i, j)
        if(matrix_(j, i) == 0) {
          matrix_(j, i) = (
            AtomInfo::vdwRadius(inner.elementType(i))
            + AtomInfo::vdwRadius(inner.elementType(j))
          );
        }

        // Ensure upper bound has a default value for Floyd's algorithm
        if(matrix_(i, j) == 0.0) {
          matrix_(i, j) = defaultUpper;
        }
      }
    }

    assert(boundInconsistencies() == 0);
  }

  bool setUpperBound(AtomIndex i, AtomIndex j, double newUpperBound);
  bool setLowerBound(AtomIndex i, AtomIndex j, double newLowerBound);

  void smooth();


  inline double upperBound(AtomIndex i, AtomIndex j) const {
    if(i < j) {
      return matrix_(i, j);
    }

    return matrix_(j, i);
  }

  inline double lowerBound(AtomIndex i, AtomIndex j) const {
    if(i < j) {
      return matrix_(j, i);
    }

    return matrix_(i, j);
  }

  unsigned boundInconsistencies() const;

  const Eigen::MatrixXd& access() const;

  outcome::result<Eigen::MatrixXd> makeDistanceMatrix(
    Random::Engine& engine,
    Partiality partiality = Partiality::All
  ) const noexcept;

  Eigen::MatrixXd makeSquaredBoundsMatrix() const;

  unsigned N() const;

private:
  Eigen::MatrixXd matrix_;

  inline double& lowerBound_(const AtomIndex i, const AtomIndex j) {
    return matrix_(
      std::max(i, j),
      std::min(i, j)
    );
  }

  inline double& upperBound_(const AtomIndex i, const AtomIndex j) {
    return matrix_(
      std::min(i, j),
      std::max(i, j)
    );
  }

};

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

#endif