text: /*
   * Copyright 2009-2020 The VOTCA Development Team (http://www.votca.org)
   *
   * Licensed under the Apache License, Version 2.0 (the "License");
   * you may not use this file except in compliance with the License.
   *
   *     http://www.apache.org/licenses/LICENSE-2.0
   *
   * Unless required by applicable law or agreed to in writing, software
   * distributed under the License is distributed on an "AS IS" BASIS,
   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
   * See the License for the specific language governing permissions and
   * limitations under the License.
   *
   */
  #define BOOST_TEST_MAIN
  
  #define BOOST_TEST_MODULE bfgs_test
  
  // Standard includes
  #include <iostream>
  
  // Third party includes
  #include <boost/format.hpp>
  #include <boost/test/unit_test.hpp>
  
  // Local VOTCA includes
  #include "votca/xtp/adiis_costfunction.h"
  #include "votca/xtp/bfgs_trm.h"
  #include "votca/xtp/logger.h"
  #include "votca/xtp/optimiser_costfunction.h"
  
  using namespace votca::xtp;
  using namespace votca;
  using namespace std;
  
  BOOST_AUTO_TEST_SUITE(bfgs_test)
  
  BOOST_AUTO_TEST_CASE(parabola_test) {
    class parabola : public Optimiser_costfunction {
  
      double EvaluateCost(const Eigen::VectorXd& parameters) override {
        Eigen::VectorXd value = parameters;
        value(0) -= 2;
        double cost = value.cwiseAbs2().sum();
        return cost;
      }
  
      Eigen::VectorXd EvaluateGradient(
          const Eigen::VectorXd& parameters) override {
        Eigen::VectorXd gradient = 2 * parameters;
        gradient(0) -= 4;
        return gradient;
      }
  
      bool Converged(const Eigen::VectorXd&, double,
                     const Eigen::VectorXd& gradient) override {
        if (gradient.cwiseAbs().maxCoeff() < 1e-8) {
          return true;
        } else {
          return false;
        }
      }
  
      Index NumParameters() const override { return 5; }
    };
  
    parabola p5;
    BFGSTRM bfgstrm(p5);
    Logger log;
    bfgstrm.setLog(&log);
    bfgstrm.setNumofIterations(100);
    bfgstrm.setTrustRadius(0.1);
    bfgstrm.Optimize(5 * Eigen::VectorXd::Ones(5));
  
    Eigen::VectorXd ref = Eigen::VectorXd::Zero(5);
    ref(0) = 2;
    bool equal = bfgstrm.getParameters().isApprox(ref, 0.00001);
    if (!equal) {
      cout << "minimum found:" << endl;
      cout << bfgstrm.getParameters() << endl;
      cout << "minimum ref:" << endl;
      cout << ref << endl;
    } else {
      cout << bfgstrm.getIteration() << endl;
    }
    BOOST_CHECK_EQUAL(equal, 1);
  }
  
  BOOST_AUTO_TEST_CASE(booth_test) {
    class booth : public Optimiser_costfunction {
  
      double EvaluateCost(const Eigen::VectorXd& parameters) override {
        double x = parameters[0];
        double y = parameters[1];
  
        return (x + 2 * y - 7) * (x + 2 * y - 7) +
               (2 * x + y - 5) * (2 * x + y - 5);
      }
  
      Eigen::VectorXd EvaluateGradient(
          const Eigen::VectorXd& parameters) override {
        double x = parameters[0];
        double y = parameters[1];
        Eigen::VectorXd gradient = Eigen::VectorXd::Zero(2);
        gradient[0] = 2 * (5 * x + 4 * y - 17);
        gradient[1] = 2 * (4 * x + 5 * y - 19);
        return gradient;
      }
  
      bool Converged(const Eigen::VectorXd&, double,
                     const Eigen::VectorXd& gradient) override {
        if (gradient.cwiseAbs().maxCoeff() < 1e-8) {
          return true;
        } else {
          return false;
        }
      }
  
      Index NumParameters() const override { return 2; }
    };
  
    booth p2;
    BFGSTRM bfgstrm(p2);
    Logger log;
    bfgstrm.setLog(&log);
    bfgstrm.setNumofIterations(1000);
    bfgstrm.setTrustRadius(0.1);
    bfgstrm.Optimize(5 * Eigen::VectorXd::Ones(2));
  
    Eigen::VectorXd ref = Eigen::VectorXd::Zero(2);
    ref << 1, 3;
    bool equal = bfgstrm.getParameters().isApprox(ref, 0.00001);
    if (!equal) {
      cout << "minimum found:" << endl;
      cout << bfgstrm.getParameters() << endl;
      cout << "minimum ref:" << endl;
      cout << ref << endl;
    } else {
      cout << bfgstrm.getIteration() << endl;
    }
    BOOST_CHECK_EQUAL(equal, 1);
  }
  
  BOOST_AUTO_TEST_CASE(adiis_test) {
    Index size = 5;
  
    Eigen::VectorXd DiF = Eigen::VectorXd::Zero(size);
    DiF << 0.679243, 0.562675, 0.39399, -0.0258519, 0;
    Eigen::MatrixXd DiFj = Eigen::MatrixXd::Zero(size, size);
    DiFj << 0.613998, 0.192684, 0.0326914, -0.193661, 0, 0.192371, 0.0653754,
        0.0131825, -0.0631915, 0, 0.032739, 0.0131883, 0.0038493, -0.0110991, 0,
        -0.192873, -0.0631203, -0.0111063, 0.0633221, 0, 0, 0, 0, 0, 0;
  
    Logger log;
  
    ADIIS_costfunction a_cost = ADIIS_costfunction(DiF, DiFj);
    BFGSTRM optimizer = BFGSTRM(a_cost);
    optimizer.setNumofIterations(1000);
    optimizer.setTrustRadius(0.01);
    optimizer.setLog(&log);
    // Starting point: equal weights on all matrices
    Eigen::VectorXd coeffs = Eigen::VectorXd::Constant(size, 1.0 / (double)size);
    optimizer.Optimize(coeffs);
    bool success = optimizer.Success();
    coeffs = optimizer.getParameters().cwiseAbs2();
    double xnorm = coeffs.sum();
    coeffs /= xnorm;
  
    BOOST_CHECK_EQUAL(success, true);
    Eigen::VectorXd ref = Eigen::VectorXd::Zero(size);
    ref << 0, 0, 0, 0.40826075912352, 0.59173924087648;
    bool equal = coeffs.isApprox(ref, 0.00001);
    if (!equal) {
      cout << "minimum found:" << endl;
      cout << coeffs << endl;
      cout << "minimum ref:" << endl;
      cout << ref << endl;
    } else {
      cout << optimizer.getIteration() << endl;
    }
    BOOST_CHECK_EQUAL(equal, 1);
  }
  
  BOOST_AUTO_TEST_SUITE_END()
  