ElementaryStepOptimizer.h

Go to the documentation of this file.

#ifndef READUCT_ELEMENTARYSTEPOPTIMIZATION_ELEMENTARYSTEPOPTIMIZER_H
#define READUCT_ELEMENTARYSTEPOPTIMIZATION_ELEMENTARYSTEPOPTIMIZER_H

#include "CostCalculators/ReaductDefault.h"
#include "ElementaryStepOptimizerSettings.h"
#include "EnergiesAndGradientsAlongSpline.h"
#include "ReactionProfile.h"
#include "RecurringProfileCalculator.h"
#include "TypeConverter.h"
#include "Utils/Geometry/AtomCollection.h"
#include "Utils/Optimizer/GradientBased/SteepestDescent.h"
#include "Utils/Settings.h"
#include "Utils/UniversalSettings/SettingsNames.h"
#include <Core/Interfaces/Calculator.h>
#include <Eigen/Core>
#include <utility>

namespace Scine {
namespace Readuct {
namespace ElementaryStepOptimization {

class ElementaryStepOptimizerBase {
 public:
  ElementaryStepOptimizerBase() = default;

  virtual ~ElementaryStepOptimizerBase() = default;

  virtual int optimize(Core::Log& log) = 0;

  virtual ReactionProfile& getReactionProfile() = 0;

  virtual void setSettings(const Utils::Settings& settings) = 0;

  virtual Utils::Settings getSettings() const = 0;
  virtual const std::shared_ptr<Utils::Settings> getCalculatorSettings() const = 0;
  virtual Utils::GradientBasedCheck getConvergenceCheck() const = 0;
};

template<class OptimizerType>
class ElementaryStepOptimizer : public ElementaryStepOptimizerBase {
 public:
  explicit ElementaryStepOptimizer(Core::Calculator& calculator, ReactionProfile initialReactionProfile)
    : _calculator(calculator), _profile(std::move(initialReactionProfile)){};

  virtual void setSettings(const Utils::Settings& settings) override {
    check.applySettings(settings);
    optimizer.applySettings(settings);
    numberEquidistantPoints = settings.getInt("num_integration_points");
  };

  virtual Utils::Settings getSettings() const override {
    auto settings = ElementaryStepOptimizerSettings<OptimizerType, Utils::GradientBasedCheck>(optimizer, check);
    settings.modifyInt("num_integration_points", numberEquidistantPoints);
    return settings;
  };

  const std::shared_ptr<Utils::Settings> getCalculatorSettings() const override {
    return std::make_shared<Utils::Settings>(_calculator.settings());
  };

  int optimize(Core::Log& log) override {
    // Preparations
    EnergiesAndGradientsAlongSpline valuesAlongSpline;
    const auto& elements = _profile.getMolecularSpline().getElements();
    const auto structure = _calculator.getStructure();
    if (!structure || structure->getElements() != elements) {
      Utils::PositionCollection positions = Utils::PositionCollection::Zero(elements.size(), 3); // dummy positions
      Utils::AtomCollection atoms(elements, positions);
      _calculator.setStructure(atoms);
    }

    // Get initial variable values and map to format required by update function
    auto& spline = _profile.getMolecularSpline().getBSpline();
    auto variables = TypeConverter::getInnerControlPointMatrix(spline);
    Eigen::VectorXd variablesVector;
    int nRows = variables.rows();
    int nColumns = variables.cols();
    variablesVector = Eigen::Map<const Eigen::VectorXd>(variables.data(), nRows * nColumns);

    // Create a profile calculator
    RecurringProfileCalculator profileCalculator(_calculator, numberEquidistantPoints);

    // Define update function
    auto const update = [&](const Eigen::VectorXd& parameters, double& value, Eigen::VectorXd& gradients) {
      Eigen::MatrixXd variables;
      variables = Eigen::Map<const Eigen::MatrixXd>(parameters.data(), nRows, nColumns);

      // Insert variables into Spline
      TypeConverter::setInnerControlPoints(spline, variables);

      // Evaluate cost value
      profileCalculator.calculateEnergiesAndGradients(spline);
      auto valuesAlongSpline = profileCalculator.valuesAlongSpline();
      CostBasedOptimization::ReaductDefault costCalculator;
      costCalculator.calculateCost(spline, valuesAlongSpline);
      value = costCalculator.getCost();

      // Evaluate gradients
      Eigen::MatrixXd fullDerivatives = costCalculator.getCostDerivatives();
      int numberInnerControlPoints = spline.controlPointCount() - 2;
      Eigen::MatrixXd derivatives = fullDerivatives.middleRows(1, numberInnerControlPoints);
      gradients = Eigen::Map<const Eigen::VectorXd>(derivatives.data(), nRows * nColumns);

      // Update reaction profile
      _profile.getProfileEnergies() = ProfileEnergies{profileCalculator.getCoordinates(), profileCalculator.getEnergies()};
    };

    // Optimize
    auto cycles = optimizer.optimize(variablesVector, update, check, log);

    return cycles;
  }

  ReactionProfile& getReactionProfile() override {
    return _profile;
  }

  Utils::GradientBasedCheck getConvergenceCheck() const override {
    return check;
  };

  OptimizerType optimizer;

  Utils::GradientBasedCheck check;

  int numberEquidistantPoints = numIntegrationPointsDefaultValue;

 private:
  Core::Calculator& _calculator;
  ReactionProfile _profile;
};
} // namespace ElementaryStepOptimization
} // namespace Readuct
} // namespace Scine

#endif // READUCT_ELEMENTARYSTEPOPTIMIZER_H