Properties.h

Go to the documentation of this file.

#ifndef INCLUDE_SHAPE_CONSTEXPR_PROPERTIES_H
#define INCLUDE_SHAPE_CONSTEXPR_PROPERTIES_H

#include "Molassembler/Shapes/Data.h"
#include "Molassembler/Shapes/constexpr/Data.h"

#include "Molassembler/Temple/constexpr/Bitset.h"
#include "Molassembler/Temple/constexpr/DynamicSet.h"
#include "Molassembler/Temple/constexpr/FloatingPointComparison.h"
#include "Molassembler/Temple/constexpr/Numeric.h"
#include "Molassembler/Temple/Cache.h"

namespace Scine {
namespace Molassembler {
namespace Shapes {

namespace ConstexprProperties {

constexpr double floatingPointEqualityTolerance = 1e-4;

template<typename ShapeClass>
constexpr double calculateSmallestAngle() {
  double smallestAngle = M_PI;

  for(unsigned i = 0; i < ShapeClass::size; ++i) {
    for(unsigned j = i + 1; j < ShapeClass::size; ++j) {
      double returnedAngle = ShapeClass::angleFunction(i, j);
      if(returnedAngle < smallestAngle) {
        smallestAngle = returnedAngle;
      }
    }
  }

  return smallestAngle;
}

template<typename ShapeClass>
struct AngleBoundsFunctor {
  static constexpr std::pair<double, double> value() {
    double smallestAngle = M_PI;
    double largestAngle = 0;

    for(unsigned i = 0; i < ShapeClass::size; ++i) {
      for(unsigned j = i + 1; j < ShapeClass::size; ++j) {
        double returnedAngle = ShapeClass::angleFunction(i, j);
        if(returnedAngle < smallestAngle) {
          smallestAngle = returnedAngle;
        }

        if(returnedAngle > largestAngle) {
          largestAngle = returnedAngle;
        }
      }
    }

    return {smallestAngle, largestAngle};
  }
};

template<typename ... ShapeClass>
struct minAngleFunctor {
  static constexpr double value() {
    const std::array<double, sizeof...(ShapeClass)> smallestAngles {{
      calculateSmallestAngle<ShapeClass>()...
    }};

    // C++17 min_element (isn't constexpr before)
    double minElement = smallestAngles.at(0);

    for(unsigned i = 1; i < sizeof...(ShapeClass); ++i) {
      if(smallestAngles.at(i) < minElement) {
        minElement = smallestAngles.at(i);
      }
    }

    return minElement;
  }
};

template<typename T, size_t size>
using ArrayType = Temple::Array<T, size>;

template<typename ShapeClass>
constexpr auto startingIndexSequence() {
  return Temple::iota<ArrayType, unsigned, ShapeClass::size>();
}

template<typename ShapeClass, size_t... Indices>
constexpr ArrayType<unsigned, ShapeClass::size> applyRotationImpl(
  const ArrayType<unsigned, ShapeClass::size>& indices,
  const unsigned rotationFunctionIndex,
  std::index_sequence<Indices...> /* inds */
) {
  return {
    indices.at(
      ShapeClass::rotations.at(rotationFunctionIndex).at(Indices)
    )...
  };
}

template<typename ShapeClass>
constexpr ArrayType<unsigned, ShapeClass::size> applyRotation(
  const ArrayType<unsigned, ShapeClass::size>& indices,
  const unsigned rotationFunctionIndex
) {
  return applyRotationImpl<ShapeClass>(
    indices,
    rotationFunctionIndex,
    std::make_index_sequence<ShapeClass::size>{}
  );
}

template<typename ShapeClass, unsigned rotationFunctionIndex>
constexpr unsigned rotationPeriodicityImpl(
  const ArrayType<unsigned, ShapeClass::size>& runningIndices,
  const unsigned count
) {
  if(
    Temple::arraysEqual(
      runningIndices,
      startingIndexSequence<ShapeClass>()
    )
  ) {
    return count;
  }

  return rotationPeriodicityImpl<ShapeClass, rotationFunctionIndex>(
    applyRotation<ShapeClass>(
      runningIndices,
      rotationFunctionIndex
    ),
    count + 1
  );
}

template<typename ShapeClass, unsigned rotationFunctionIndex>
constexpr unsigned rotationPeriodicity() {
  return rotationPeriodicityImpl<ShapeClass, rotationFunctionIndex>(
    applyRotation<ShapeClass>(
      startingIndexSequence<ShapeClass>(),
      rotationFunctionIndex
    ),
    1
  );
}

template<typename ShapeClass, size_t ...Inds>
constexpr ArrayType<unsigned, ShapeClass::rotations.size()>
rotationPeriodicitiesImpl(std::index_sequence<Inds...> /* inds */) {
  return { rotationPeriodicity<ShapeClass, Inds>()... };
}

template<typename ShapeClass>
constexpr ArrayType<unsigned, ShapeClass::rotations.size()> rotationPeriodicities() {
  return rotationPeriodicitiesImpl<ShapeClass>(
    std::make_index_sequence<ShapeClass::rotations.size()>{}
  );
}

template<typename ShapeClass>
struct allRotationPeriodicities {
  static constexpr ArrayType<
    unsigned,
    ShapeClass::rotations
  > value = rotationPeriodicitiesImpl<ShapeClass>(
    std::make_index_sequence<ShapeClass::rotations.size()>{}
  );
};

template<typename ... ShapeClasses>
struct maxShapeSizeFunctor {
  static constexpr unsigned value() {
    ArrayType<unsigned, sizeof...(ShapeClasses)> sizes {
      ShapeClasses::size...
    };

    return Temple::max(sizes);
  }
};

constexpr unsigned maxShapeSize = Temple::Tuples::unpackToFunction<
  Shapes::Data::allShapeDataTypes,
  maxShapeSizeFunctor
>();

template<typename ShapeClass>
constexpr Temple::Vector getCoordinates(const unsigned indexInSymmetry) {
  if(indexInSymmetry != ORIGIN_PLACEHOLDER) {
    return ShapeClass::coordinates.at(indexInSymmetry);
  }

  return {0, 0, 0};
}

constexpr double getTetrahedronVolume(
  const Temple::Vector& i,
  const Temple::Vector& j,
  const Temple::Vector& k,
  const Temple::Vector& l
) {
  return (i - l).dot(
    (j - l).cross(k - l)
  );
}

// Typedef to avoid reusing C-Style function ptr type
using AngleFunctionPtr = double(*)(const unsigned, const unsigned);

template<size_t size>
constexpr double calculateAngularDistortion(
  const ArrayType<unsigned, size>& indexMapping,
  const size_t sourceSymmetrySize,
  const AngleFunctionPtr sourceAngleFunction,
  const AngleFunctionPtr targetAngleFunction
) {
  double distortionSum = 0;

  for(unsigned i = 0; i < sourceSymmetrySize; ++i) {
    for(unsigned j = i + 1; j < sourceSymmetrySize; ++j) {
      distortionSum += Temple::Math::abs(
        sourceAngleFunction(i, j)
        - targetAngleFunction(
          indexMapping.at(i),
          indexMapping.at(j)
        )
      );
    }
  }

  return distortionSum;
}

template<size_t size>
constexpr unsigned propagateSymmetryPosition(
  const unsigned symmetryPosition,
  const ArrayType<unsigned, size>& indexMapping
) {
  if(symmetryPosition != ORIGIN_PLACEHOLDER) {
    return indexMapping.at(symmetryPosition);
  }

  return ORIGIN_PLACEHOLDER;
}

template<typename ShapeClassFrom, typename ShapeClassTo>
constexpr double calculateChiralDistortion(
  const ArrayType<
    unsigned,
    Temple::Math::max(
      ShapeClassFrom::size,
      ShapeClassTo::size
    )
  >& indexMapping
) {
  double chiralDistortion = 0;

  // C++17:
  // for(const auto& tetrahedron : ShapeClassFrom::tetrahedra) {
  for(unsigned i = 0; i < ShapeClassFrom::tetrahedra.size(); ++i) {
    const auto& tetrahedron = ShapeClassFrom::tetrahedra.at(i);

    chiralDistortion += Temple::Math::abs(
      getTetrahedronVolume(
        getCoordinates<ShapeClassFrom>(tetrahedron.at(0)),
        getCoordinates<ShapeClassFrom>(tetrahedron.at(1)),
        getCoordinates<ShapeClassFrom>(tetrahedron.at(2)),
        getCoordinates<ShapeClassFrom>(tetrahedron.at(3))
      ) - getTetrahedronVolume(
        getCoordinates<ShapeClassTo>(
          propagateSymmetryPosition(tetrahedron.at(0), indexMapping)
        ),
        getCoordinates<ShapeClassTo>(
          propagateSymmetryPosition(tetrahedron.at(1), indexMapping)
        ),
        getCoordinates<ShapeClassTo>(
          propagateSymmetryPosition(tetrahedron.at(2), indexMapping)
        ),
        getCoordinates<ShapeClassTo>(
          propagateSymmetryPosition(tetrahedron.at(3), indexMapping)
        )
      )
    );
  }

  return chiralDistortion;
}

template<size_t size>
constexpr ArrayType<unsigned, size> symPosMapping(
  const ArrayType<unsigned, size>& mapping
) {
  /* Creates the list of indices in the target shape. Why is this necessary?
   *
   * E.g. An index mapping from linear to T-shaped. The individual
   * shape-internal numbering schemes are shown for the shape positions.
   *
   *  1  –▶  0
   *  |      |
   * (_)    (_) – 1 (new)
   *  |      |                Line pos. 0 to
   *  0  –▶  2                pos. 2 in Tshaped
   *                                 |  ┌– Line pos. 1 to pos. 0 in Tshaped
   *                                 |  |  ┌– This position is new
   * This mapping is represented as {2, 0, 1}.
   *
   * This function writes the indices of original mapping into the target
   * shape's indexing scheme.
   *
   * For this example, this returns {1, 2, 0}:
   *
   *  1 (at pos 0 in internal indexing scheme)
   *  |
   * (_) – 2 (etc.)
   *  |
   *  0
   *
   * The closely related mapping {0, 2, 1} yields target indices {0, 2, 1}.
   *
   * Which of these properties are related by target shape rotations?
   *
   *
   *     mapping       target indices
   * -----------------------------------
   *    {2, 0, 1}   =>   {1, 2, 0}
   *        ▲                ▲
   *        |                |
   *        X                | (C2 rotation)
   *        |                |
   *        ▼                ▼
   *    {0, 2, 1}   =>   {0, 2, 1}
   *
   */
  ArrayType<unsigned, size> symmetryPositions;

  for(unsigned i = 0; i < size; ++i) {
    symmetryPositions.at(
      mapping.at(i)
    ) = i;
  }

  return symmetryPositions;
}

template<typename ShapeClass>
constexpr unsigned maxRotations() {
  return Temple::reduce(rotationPeriodicities<ShapeClass>(), 1u, std::multiplies<>());
}

// Some helper types for use in generateAllRotations
template<typename ShapeClass>
using IndicesList = ArrayType<unsigned, ShapeClass::size>;

using IndexListStorageType = unsigned;

template<typename ShapeClass>
using RotationsSetType = Temple::DynamicSet<
  IndicesList<ShapeClass>,
  maxRotations<ShapeClass>() * 2
>;

template<typename ShapeClass>
using ChainStructuresArrayType = Temple::DynamicArray<
  IndicesList<ShapeClass>,
  maxRotations<ShapeClass>() * 2 // factor is entirely arbitrary
>;

template<typename ShapeClass>
using ChainArrayType = Temple::DynamicArray<
  unsigned,
  maxRotations<ShapeClass>() * 2 // factor is entirely arbitrary
>;

template<typename ShapeClass>
constexpr auto generateAllRotations(const IndicesList<ShapeClass>& indices) {
  RotationsSetType<ShapeClass> rotations;

  rotations.insert(indices);

  ChainStructuresArrayType<ShapeClass> chainStructures;
  chainStructures.push_back(indices);

  ChainArrayType<ShapeClass> chain {0u};

  while(chain.front() < ShapeClass::rotations.size()) {

    auto generated = applyRotation<ShapeClass>(
      chainStructures.back(),
      chain.back()
    );

    if(!rotations.contains(generated)) {
      rotations.insert(generated);
      chainStructures.push_back(generated);
      chain.push_back(0);
    } else {
      // collapse the chain until we are at an incrementable position (if need be)
      while(
        chain.size() > 1
        && chain.back() == ShapeClass::rotations.size() - 1
      ) {
        chain.pop_back();
        chainStructures.pop_back();
      }

      // increment
      ++chain.back();
    }
  }

  return rotations;
}

struct MappingsReturnType {
  static constexpr size_t maxMappingsSize = 50;

  using MappingsList = Temple::DynamicSet<
    Temple::DynamicArray<unsigned, maxShapeSize>,
    maxMappingsSize
  >;

  MappingsList mappings;
  double angularDistortion, chiralDistortion;

  constexpr MappingsReturnType()
    : angularDistortion(std::numeric_limits<double>::max()),
      chiralDistortion(std::numeric_limits<double>::max())
  {}

  constexpr MappingsReturnType(
    MappingsList&& passMappings,
    double&& passAngularDistortion,
    double&& passChiralDistortion
  ) : mappings(passMappings),
      angularDistortion(passAngularDistortion),
      chiralDistortion(passChiralDistortion)
  {}

  constexpr bool operator == (const MappingsReturnType& other) const {
    return (
      mappings == other.mappings
      && angularDistortion == other.angularDistortion
      && chiralDistortion == other.chiralDistortion
    );
  }

  constexpr bool operator < (const MappingsReturnType& other) const {
    return (
      std::tie(mappings, angularDistortion, chiralDistortion)
      < std::tie(other.mappings, other.angularDistortion, other.chiralDistortion)
    );
  }
};

template<class ShapeClassFrom, class ShapeClassTo>
constexpr auto symmetryTransitionMappings() {
  static_assert(
    (
      ShapeClassTo::size == ShapeClassFrom::size + 1
      || ShapeClassTo::size == ShapeClassFrom::size
    ),
    "This function can handle only cases of equal or increasing shape size"
  );

  typename MappingsReturnType::MappingsList bestMappings;

  double lowestAngleDistortion = 100;
  double lowestChiralDistortion = 100;

  auto indexMapping = startingIndexSequence<ShapeClassTo>();

  Temple::Bitset<Temple::Math::factorial(ShapeClassTo::size)> encounteredMappings;

  Temple::Floating::ExpandedRelativeEqualityComparator<double> comparator {
    floatingPointEqualityTolerance
  };

  do {
    auto mapped = symPosMapping(indexMapping);
    auto storageVersion = Temple::permutationIndex(mapped);

    if(!encounteredMappings.test(storageVersion)) {
      double angularDistortion = calculateAngularDistortion(
        indexMapping,
        ShapeClassFrom::size,
        ShapeClassFrom::angleFunction,
        ShapeClassTo::angleFunction
      );

      double chiralDistortion = calculateChiralDistortion<
        ShapeClassFrom,
        ShapeClassTo
      >(indexMapping);

      /* Summary of cases:
       * - If any improvement is made on angular distortion, clear all and add
       * - If angular distortion is equal but chiral distortion is improved,
       *   clear all and add
       * - If angular distortion and chiral distortion are equal, just add
       *
       * The boolean cases below are, AFAICT, the fewest comparisons.
       * Feel free to see if you can find some way with fewer.
       */

      bool addMapping = (
        comparator.isLessThan(angularDistortion, lowestAngleDistortion)
        || (
          comparator.isEqual(angularDistortion, lowestAngleDistortion)
          && comparator.isLessOrEqual(chiralDistortion, lowestChiralDistortion)
        )
      );

      bool clearExisting = (
        addMapping
        && !(
          comparator.isEqual(angularDistortion, lowestAngleDistortion)
          && comparator.isEqual(chiralDistortion, lowestChiralDistortion)
        )
      );

      if(clearExisting) {
        bestMappings.clear();
        lowestAngleDistortion = angularDistortion;
        lowestChiralDistortion = chiralDistortion;
      }

      if(addMapping) {
        bestMappings.insert(indexMapping);
      }

      // Add all rotations to the encountered mappings
      auto allRotations = generateAllRotations<ShapeClassTo>(mapped);

      for(const auto& rotation : allRotations) {
        encounteredMappings.set(
          Temple::permutationIndex(rotation)
        );
      }
    }
  } while(Temple::inPlaceNextPermutation(indexMapping));

  return MappingsReturnType(
    std::move(bestMappings),
    std::move(lowestAngleDistortion),
    std::move(lowestChiralDistortion)
  );
}

template<class ShapeClassFrom, class ShapeClassTo>
constexpr auto ligandLossMappings(const unsigned deletedSymmetryPosition) {

  static_assert(
    ShapeClassFrom::size == ShapeClassTo::size + 1,
    "Ligand loss pathway calculation must involve shape size decrease"
  );

  assert(deletedSymmetryPosition < ShapeClassFrom::size);


  typename MappingsReturnType::MappingsList bestMappings;

  double lowestAngleDistortion = 100;
  double lowestChiralDistortion = 100;

  // Construct the initial index mapping
  ArrayType<unsigned, ShapeClassTo::size> indexMapping;
  for(unsigned i = 0; i < deletedSymmetryPosition; ++i) {
    indexMapping.at(i) = i;
  }
  for(unsigned i = deletedSymmetryPosition; i < ShapeClassFrom::size - 1; ++i) {
    indexMapping.at(i) = i + 1;
  }

  /* NOTE: From here the algorithm is identical to symmetryTransitionMappings
   * save that symmetryTo and symmetryFrom are swapped in all occasions
   */
  Temple::Bitset<Temple::Math::factorial(ShapeClassFrom::size)> encounteredMappings;

  Temple::Floating::ExpandedRelativeEqualityComparator<double> comparator {
    floatingPointEqualityTolerance
  };

  do {
    auto mapped = symPosMapping(indexMapping);
    auto storageVersion = Temple::permutationIndex(mapped);

    if(!encounteredMappings.test(storageVersion)) {
      double angularDistortion = calculateAngularDistortion(
        indexMapping,
        ShapeClassTo::size,
        ShapeClassTo::angleFunction,
        ShapeClassFrom::angleFunction
      );

      double chiralDistortion = calculateChiralDistortion<
        ShapeClassTo,
        ShapeClassFrom
      >(indexMapping);

      bool addMapping = (
        comparator.isLessThan(angularDistortion, lowestAngleDistortion)
        || (
          comparator.isEqual(angularDistortion, lowestAngleDistortion)
          && comparator.isLessOrEqual(chiralDistortion, lowestChiralDistortion)
        )
      );

      bool clearExisting = (
        addMapping
        && !(
          comparator.isEqual(angularDistortion, lowestAngleDistortion)
          && comparator.isEqual(chiralDistortion, lowestChiralDistortion)
        )
      );

      if(clearExisting) {
        bestMappings.clear();
        lowestAngleDistortion = angularDistortion;
        lowestChiralDistortion = chiralDistortion;
      }

      if(addMapping) {
        bestMappings.insert(indexMapping);
      }

      // Add all rotations to the encountered mappings
      auto allRotations = generateAllRotations<ShapeClassTo>(indexMapping);

      for(const auto& rotation : allRotations) {
        encounteredMappings.set(
          Temple::permutationIndex(rotation)
        );
      }
    }
  } while(Temple::inPlaceNextPermutation(indexMapping));

  return MappingsReturnType(
    std::move(bestMappings),
    std::move(lowestAngleDistortion),
    std::move(lowestChiralDistortion)
  );
}

/* Pre-compute all ligand gain situations */
template<typename SymmetrySource, typename SymmetryTarget>
constexpr
std::enable_if_t<
  (
    SymmetryTarget::size <= 7
    && (
      SymmetrySource::size == SymmetryTarget::size
      || SymmetrySource::size + 1 == SymmetryTarget::size
    )
  ),
  Temple::Optional<MappingsReturnType>
> calculateMapping() {
  return Temple::Optional<MappingsReturnType> {
    symmetryTransitionMappings<SymmetrySource, SymmetryTarget>()
  };
}

template<typename SymmetrySource, typename SymmetryTarget>
constexpr
std::enable_if_t<
  !(
    SymmetryTarget::size <= 7
    && (
      SymmetrySource::size == SymmetryTarget::size
      || SymmetrySource::size + 1 == SymmetryTarget::size
    )
  ),
  Temple::Optional<MappingsReturnType>
> calculateMapping() {
  return Temple::Optional<MappingsReturnType> {};
}

template<typename ShapeClass>
constexpr unsigned numUnlinkedStereopermutations(
  const unsigned nIdenticalLigands
) {
  unsigned count = 1;

  auto indices = startingIndexSequence<ShapeClass>();

  for(unsigned i = 0; i < nIdenticalLigands; ++i) {
    indices.at(i) = 0;
  }

  Temple::Bitset<Temple::Math::factorial(ShapeClass::size)> rotations;

  auto initialRotations = generateAllRotations<ShapeClass>(indices);

  for(const auto& rotation : initialRotations) {
    rotations.set(Temple::permutationIndex(rotation));
  }

  while(Temple::inPlaceNextPermutation(indices)) {
    if(!rotations.test(Temple::permutationIndex(indices))) {
      auto allRotations = generateAllRotations<ShapeClass>(indices);
      for(const auto& rotation : allRotations) {
        rotations.set(Temple::permutationIndex(rotation));
      }

      ++count;
    }
  }

  return count;
}

template<typename ShapeClass>
constexpr bool hasMultipleUnlinkedStereopermutations(const unsigned nIdenticalLigands) {
  if(nIdenticalLigands == ShapeClass::size) {
    return false;
  }

  auto indices = startingIndexSequence<ShapeClass>();

  for(unsigned i = 0; i < nIdenticalLigands; ++i) {
    indices.at(i) = 0;
  }

  Temple::Bitset<Temple::Math::factorial(ShapeClass::size)> rotations;

  auto initialRotations = generateAllRotations<ShapeClass>(indices);

  for(const auto& rotation : initialRotations) {
    rotations.set(Temple::permutationIndex(rotation));
  }

  while(Temple::inPlaceNextPermutation(indices)) {
    if(!rotations.test(Temple::permutationIndex(indices))) {
      return true;
    }
  }

  return false;
}

} // namespace ConstexprProperties
} // namespace Shapes
} // namespace Molassembler
} // namespace Scine

#endif