Merge pull request #70 from oberbichler/feature/steepest-decent

Steepest descent

include/eqlib/Armijo.h

@@ -51,10 +51,9 @@ class Armijo
 
         if (reset) {
             m_problem->set_f(f_init);
+            m_problem->set_x(m_x_init);
         }
 
-        m_problem->set_x(m_x_init);
-
         return alpha;
     }
 

include/eqlib/Problem.h

@@ -992,6 +992,38 @@ class Problem {
         set_x(Map<const Vector>(value, nb_variables()));
     }
 
+    void add_x(Ref<const Vector> delta) const
+    {
+        if (length(delta) != nb_variables()) {
+            throw std::runtime_error("Invalid size");
+        }
+
+        for (index i = 0; i < length(delta); i++) {
+            variable(i)->value() += delta[i];
+        }
+    }
+
+    void add_x(double* const delta) const
+    {
+        add_x(Map<const Vector>(delta, nb_variables()));
+    }
+
+    void sub_x(Ref<const Vector> delta) const
+    {
+        if (length(delta) != nb_variables()) {
+            throw std::runtime_error("Invalid size");
+        }
+
+        for (index i = 0; i < length(delta); i++) {
+            variable(i)->value() -= delta[i];
+        }
+    }
+
+    void sub_x(double* const delta) const
+    {
+        sub_x(Map<const Vector>(delta, nb_variables()));
+    }
+
     Vector variable_multipliers() const
     {
         Vector result(nb_variables());
@@ -1319,6 +1351,8 @@ class Problem {
             .def_property("variable_multipliers", py::overload_cast<>(&Type::variable_multipliers, py::const_), py::overload_cast<Ref<const Vector>>(&Type::set_variable_multipliers, py::const_))
             .def_property("equation_multipliers", py::overload_cast<>(&Type::equation_multipliers, py::const_), py::overload_cast<Ref<const Vector>>(&Type::set_equation_multipliers, py::const_))
             // methods
+            .def("add_x", py::overload_cast<Ref<const Vector>>(&Type::add_x, py::const_))
+            .def("sub_x", py::overload_cast<Ref<const Vector>>(&Type::sub_x, py::const_))
             .def("clone", &Type::clone)
             .def("remove_inactive_elements", &Type::remove_inactive_elements)
             .def("compute", &Type::compute<true>, "order"_a = 2, py::call_guard<py::gil_scoped_release>())

include/eqlib/SteepestDecent.h

@@ -0,0 +1,206 @@
+#pragma once
+
+#include "Define.h"
+
+#include "Armijo.h"
+#include "Log.h"
+#include "Problem.h"
+#include "Settings.h"
+#include "Timer.h"
+
+namespace eqlib {
+
+class SteepestDecent {
+private: // types
+    using Type = eqlib::SteepestDecent;
+
+private: // members
+    Pointer<Problem> m_problem;
+    index m_iterations;
+    index m_maxiter;
+    index m_fevals;
+    index m_gevals;
+    index m_hevals;
+    double m_rnorm;
+    double m_xnorm;
+    double m_rtol;
+    double m_xtol;
+    double m_damping;
+    index m_stopping_reason;
+    Armijo m_linesearch;
+
+private: // methods
+public: // constructor
+    SteepestDecent(Pointer<Problem> problem)
+        : m_problem(problem)
+        , m_maxiter(100)
+        , m_rtol(1e-6)
+        , m_xtol(1e-6)
+        , m_damping(0.0)
+        , m_linesearch(problem)
+    {
+        if (problem->is_constrained()) {
+            throw std::runtime_error("Constraints are not supported");
+        }
+    }
+
+public: // methods
+    index iterations() const noexcept
+    {
+        return m_iterations;
+    }
+
+    index fevals() const noexcept
+    {
+        return m_fevals;
+    }
+
+    index gevals() const noexcept
+    {
+        return m_gevals;
+    }
+
+    index hevals() const noexcept
+    {
+        return m_hevals;
+    }
+
+    index maxiter() const noexcept
+    {
+        return m_maxiter;
+    }
+
+    void set_maxiter(const index value) noexcept
+    {
+        m_maxiter = value;
+    }
+
+    double rnorm() const noexcept
+    {
+        return m_rnorm;
+    }
+
+    double rtol() const noexcept
+    {
+        return m_rtol;
+    }
+
+    void set_rtol(const double value) noexcept
+    {
+        m_rtol = value;
+    }
+
+    double xtol() const noexcept
+    {
+        return m_xtol;
+    }
+
+    void set_xtol(const double value) noexcept
+    {
+        m_xtol = value;
+    }
+
+    double damping() const noexcept
+    {
+        return m_damping;
+    }
+
+    void set_damping(const double value) noexcept
+    {
+        m_damping = value;
+    }
+
+    void run()
+    {
+        // setup
+
+        Log::task_begin("Solving nonlinear system using Steepest Decent...");
+
+        Timer timer;
+
+        const index n = m_problem->nb_variables();
+
+        Vector search_direction(n);
+
+        for (index iteration = 0;; iteration++) {
+            // check max iterations
+
+            if (iteration >= m_maxiter) {
+                m_stopping_reason = 2;
+                Log::warn(1, "Stopped because iteration >= {}", m_maxiter);
+                break;
+            }
+
+            Log::task_info("Iteration {}:", iteration + 1);
+
+            // compute g
+
+            Log::task_step("Computing system...");
+
+            m_problem->compute<false, 1>();
+            m_gevals += 1;
+
+            Log::task_info("The current value is {}", m_problem->f());
+
+            // check residual
+
+            Log::task_step("Computing residual...");
+
+            const double rnorm = m_problem->df().norm();
+
+            Log::task_info("The norm of the residual is {}", rnorm);
+
+            // check residual norm
+
+            if (rnorm < m_rtol) {
+                m_stopping_reason = 0;
+                Log::task_info("Stopped because rnorm < {}", m_rtol);
+                break;
+            }
+
+            // compute delta
+
+            search_direction = -m_problem->df();
+
+            const double alpha = m_linesearch.search(search_direction, 1.0, false);
+
+            // check x norm
+
+            const double xnorm = alpha * search_direction.norm();
+
+            Log::task_info("The norm of the step is {}", xnorm);
+
+            if (xnorm < m_xtol) {
+                Log::task_info("Stopped because xnorm < {}", m_xtol);
+                m_stopping_reason = 1;
+                break;
+            }
+        }
+
+        Log::task_end("System solved in {:.3f} sec", timer.ellapsed());
+    }
+
+public: // python
+    template <typename TModule>
+    static void register_python(TModule& m)
+    {
+        namespace py = pybind11;
+        using namespace pybind11::literals;
+
+        py::class_<Type>(m, "SteepestDecent")
+            .def(py::init<Pointer<eqlib::Problem>>(), "problem"_a)
+            .def("run", &Type::run, py::call_guard<py::gil_scoped_release>())
+            // properties
+            .def_property("damping", &Type::damping, &Type::set_damping)
+            .def_property("maxiter", &Type::maxiter, &Type::set_maxiter)
+            .def_property("rtol", &Type::rtol, &Type::set_rtol)
+            // read-only properties
+            .def_property_readonly("iterations", &Type::iterations)
+            .def_property_readonly("rnorm", &Type::rnorm)
+            .def_property_readonly("fevals", &Type::fevals)
+            .def_property_readonly("gevals", &Type::gevals)
+            .def_property_readonly("hevals", &Type::hevals);
+    }
+};
+
+} // namespace eqlib

include/eqlib/Variable.h

@@ -62,6 +62,11 @@ class Variable {
         return m_act_value;
     }
 
+    double& value() noexcept
+    {
+        return m_act_value;
+    }
+
     void set_value(const double value) noexcept
     {
         m_act_value = value;
@@ -107,6 +112,15 @@ class Variable {
         m_multiplier = value;
     }
 
+    void clamp() noexcept
+    {
+        if (value() < lower_bound()) {
+            set_value(lower_bound());
+        } else if (upper_bound() < value()) {
+            set_value(upper_bound());
+        }
+    }
+
     const std::string& name() const noexcept
     {
         return m_name;
@@ -159,19 +173,19 @@ class Variable {
         using Holder = Pointer<Type>;
 
         py::class_<Type, Holder>(m, "Variable")
-            .def(py::init<double, double, double, bool, double, std::string>(), "value"_a = 0.0, "lower_bound"_a = -infinity,
-                "upper_bound"_a = infinity, "is_active"_a = true, "multiplier"_a = 1.0, "name"_a = "")
+            // constructors
+            .def(py::init<double, double, double, bool, double, std::string>(), "value"_a = 0.0, "lower_bound"_a = -infinity, "upper_bound"_a = infinity, "is_active"_a = true, "multiplier"_a = 1.0, "name"_a = "")
             .def(py::init<>())
-            .def_property("value", &Type::value, &Type::set_value)
-            .def_property("lower_bound", &Type::lower_bound,
-                &Type::set_lower_bound)
-            .def_property("upper_bound", &Type::upper_bound,
-                &Type::set_upper_bound)
+            // methods
+            .def("__float__", &Type::operator double)
+            .def("clamp", &Type::clamp)
+            // properties
+            .def_property("value", py::overload_cast<>(&Type::value, py::const_), &Type::set_value)
+            .def_property("lower_bound", &Type::lower_bound, &Type::set_lower_bound)
+            .def_property("upper_bound", &Type::upper_bound, &Type::set_upper_bound)
             .def_property("is_active", &Type::is_active, &Type::set_active)
-            .def_property("multiplier", &Type::multiplier,
-                &Type::set_multiplier)
-            .def_property("name", &Type::name, &Type::set_name)
-            .def("__float__", &Type::operator double);
+            .def_property("multiplier", &Type::multiplier, &Type::set_multiplier)
+            .def_property("name", &Type::name, &Type::set_name);
     }
 };
 

src/Module.cpp

@@ -17,6 +17,7 @@
 #include <eqlib/Objective.h>
 #include <eqlib/Parameter.h>
 #include <eqlib/Problem.h>
+#include <eqlib/SteepestDecent.h>
 #include <eqlib/SparseLU.h>
 #include <eqlib/SparseStructure.h>
 #include <eqlib/Variable.h>
@@ -86,6 +87,9 @@ PYBIND11_MODULE(eqlib, m)
     // NewtonRaphson
     eqlib::NewtonRaphson::register_python(m);
 
+    // SteepestDecent
+    eqlib::SteepestDecent::register_python(m);
+
     // LinearSolver
     eqlib::LinearSolver::register_python(m);