Merge pull request #74 from oberbichler/fix/gil-lock

Add membrane

external_libraries/hyperjet/detail/fwd.h

@@ -13,6 +13,7 @@
 
 #include <assert.h>
 #include <ostream>
+#include <sstream>
 #include <string>
 #include <type_traits>
 
@@ -32,6 +33,8 @@ HYPERJET_INLINE index length(const T& container)
     return static_cast<index>(container.size());
 }
 
+constexpr index Dynamic = -1;
+
 constexpr index init_size(const int size)
 {
     return size != -1 ? size : 0;

external_libraries/hyperjet/detail/hyperjet.h

@@ -256,10 +256,102 @@ class HyperJet {
         return result;
     }
 
+    auto backward(const std::vector<HyperJet<TScalar, Dynamic>>& xs) const
+    {
+        const index size = xs[0].size();
+
+        auto result = HyperJet<TScalar, Dynamic>::zero(size);
+
+        result.f() = f();
+
+        backward_to(xs, result.g(), result.h(), true);
+
+        return result;
+    }
+
+    void backward_to(const std::vector<HyperJet<TScalar, Dynamic>>& xs, Eigen::Ref<Eigen::VectorXd> g, Eigen::Ref<Eigen::MatrixXd> h, const bool full) const
+    {
+        const index size = g.size();
+
+        for (index i = 0; i < size; i++) {
+            for (index r = 0; r < this->size(); r++) {
+                assert(xs[r].size() == size);
+                g(i) += this->g(r) * xs[r].g(i);
+            }
+
+            for (index j = i; j < size; j++) {
+                for (index r = 0; r < this->size(); r++) {
+                    h(i, j) += this->g(r) * xs[r].h(i, j);
+
+                    for (index s = 0; s < this->size(); s++) {
+                        h(i, j) += this->h(r, s) * xs[r].g(i) * xs[s].g(j);
+                    }
+                }
+            }
+        }
+
+        if (!full) {
+            return;
+        }
+
+        for (index i = 0; i < size; i++) {
+            for (index j = 0; j < i; j++) {
+                h(i, j) = h(j, i);
+            }
+        }
+    }
+
+    auto forward(const std::vector<HyperJet<TScalar, Dynamic>>& xs) const
+    {
+        const index size = length(xs);
+
+        auto result = HyperJet<TScalar, Dynamic>::zero(size);
+
+        result.f() = f();
+
+        forward_to(xs, result.g(), result.h(), true);
+
+        return result;
+    }
+
+    void forward_to(const std::vector<HyperJet<TScalar, Dynamic>>& xs, Eigen::Ref<Eigen::VectorXd> g, Eigen::Ref<Eigen::MatrixXd> h, const bool full) const
+    {
+        const index size = g.size();
+
+        for (index i = 0; i < size; i++) {
+            g(i) += xs[i].g().dot(this->g());
+
+            for (index j = i; j < size; j++) {
+                for (index r = 0; r < this->size(); r++) {
+                    for (index s = 0; s < this->size(); s++) {
+                        h(i, j) += this->g(r) * xs[i].h(r, s) + this->h(r, s) * xs[i].g(r) * xs[j].g(r);
+                    }
+                }
+            }
+        }
+
+        if (!full) {
+            return;
+        }
+
+        for (index i = 0; i < size; i++) {
+            for (index j = 0; j < i; j++) {
+                h(i, j) = h(j, i);
+            }
+        }
+    }
+
+    std::string to_string() const
+    {
+        std::stringstream ss;
+        ss << *this;
+        return ss.str();
+    }
+
 public: // operators
     friend std::ostream& operator<<(std::ostream& out, const HyperJet<TScalar, TSize>& value)
     {
-        out << "HyperJet<" << value.m_f << ">";
+        out << value.m_f << "hj";
         return out;
     }
 
@@ -778,7 +870,7 @@ class HyperJet {
             .def("__len__", &Type::size)
             .def("__pow__", &Type::pow<double>)
             .def("__pow__", &Type::pow<index>)
-            // .def("__repr__", &Type::to_string)
+            .def("__repr__", &Type::to_string)
             .def("abs", &Type::abs)
             .def("acos", &Type::acos)
             .def("arccos", &Type::acos)
@@ -787,8 +879,12 @@ class HyperJet {
             .def("arctan2", &Type::atan2)
             .def("asin", &Type::asin)
             .def("atan", &Type::atan)
+            .def("backward", &Type::backward, "xs"_a)
+            .def("backward_to", &Type::backward_to, "xs"_a, "g"_a, "h"_a, "full"_a=true)
             .def("cos", &Type::cos)
             .def("enlarge", py::overload_cast<index, index>(&Type::enlarge, py::const_), "left"_a = 0, "right"_a = 0)
+            .def("forward", &Type::forward, "xs"_a)
+            .def("forward_to", &Type::forward_to, "xs"_a, "g"_a, "h"_a, "full"_a=true)
             .def("sin", &Type::sin)
             .def("sqrt", &Type::sqrt)
             .def("tan", &Type::tan)

external_libraries/hyperjet/detail/jet.h

@@ -50,7 +50,7 @@ class Jet {
     HYPERJET_INLINE static Type empty()
     {
         assert(TSize != -1);
-        return Type(Scalar(), Vector(TSize), Matrix(TSize, TSize));
+        return Type(Scalar(), Vector(TSize));
     }
 
     HYPERJET_INLINE static Type empty(const index size)
@@ -189,10 +189,64 @@ class Jet {
         return result;
     }
 
+    auto backward(const std::vector<Jet<TScalar, Dynamic>>& xs) const
+    {
+        const index size = xs[0].size();
+
+        auto result = Jet<TScalar, Dynamic>::zero(size);
+
+        result.f() = f();
+
+        backward_to(xs, result.g());
+
+        return result;
+    }
+
+    void backward_to(const std::vector<Jet<TScalar, Dynamic>>& xs, Eigen::Ref<Eigen::VectorXd> g) const
+    {
+        const index size = g.size();
+
+        for (index i = 0; i < size; i++) {
+            for (index r = 0; r < this->size(); r++) {
+                assert(xs[r].size() == size);
+                g(i) += this->g(r) * xs[r].g(i);
+            }
+        }
+    }
+
+    auto forward(const std::vector<Jet<TScalar, Dynamic>>& xs) const
+    {
+        const index size = length(xs);
+
+        auto result = Jet<TScalar, Dynamic>::zero(size);
+
+        result.f() = f();
+
+        forward_to(xs, result.g());
+
+        return result;
+    }
+
+    void forward_to(const std::vector<Jet<TScalar, Dynamic>>& xs, Eigen::Ref<Eigen::VectorXd> g) const
+    {
+        const index size = g.size();
+
+        for (index i = 0; i < size; i++) {
+            g(i) = xs[i].g().dot(this->g());
+        }
+    }
+
+    std::string to_string() const
+    {
+        std::stringstream ss;
+        ss << *this;
+        return ss.str();
+    }
+
 public: // operators
     friend std::ostream& operator<<(std::ostream& out, const Jet<TScalar, TSize>& value)
     {
-        out << "Jet<" << value.m_f << ">";
+        out << value.m_f << "j";
         return out;
     }
 
@@ -677,7 +731,7 @@ class Jet {
             .def("__len__", &Type::size)
             .def("__pow__", &Type::pow<double>)
             .def("__pow__", &Type::pow<index>)
-            // .def("__repr__", &Type::to_string)
+            .def("__repr__", &Type::to_string)
             .def("abs", &Type::abs)
             .def("acos", &Type::acos)
             .def("arccos", &Type::acos)
@@ -686,8 +740,12 @@ class Jet {
             .def("arctan2", &Type::atan2)
             .def("asin", &Type::asin)
             .def("atan", &Type::atan)
+            .def("backward", &Type::backward, "xs"_a)
+            .def("backward_to", &Type::backward_to, "xs"_a, "g"_a)
             .def("cos", &Type::cos)
             .def("enlarge", py::overload_cast<index, index>(&Type::enlarge, py::const_), "left"_a = 0, "right"_a = 0)
+            .def("forward", &Type::forward, "xs"_a)
+            .def("forward_to", &Type::forward_to, "xs"_a, "g"_a)
             .def("sin", &Type::sin)
             .def("sqrt", &Type::sqrt)
             .def("tan", &Type::tan)

external_libraries/hyperjet/detail/space.h

@@ -26,6 +26,13 @@ struct Space<0, TScalar, TSize> {
     template <index TRows, index TCols>
     using Matrix = Eigen::Matrix<Scalar, TRows, TCols>;
 
+    static Scalar empty()
+    {
+        static_assert(-1 <= TSize);
+
+        return Scalar();
+    }
+
     static Scalar constant(TScalar value)
     {
         static_assert(-1 <= TSize);
@@ -40,6 +47,20 @@ struct Space<0, TScalar, TSize> {
         return value;
     }
 
+    template <index TOffset, index TDerivedSize>
+    HYPERJET_INLINE static Eigen::Matrix<Scalar, 1, TDerivedSize> variables(Eigen::Matrix<TScalar, 1, TDerivedSize> value)
+    {
+        static_assert(-1 <= TSize);
+
+        Eigen::Matrix<Scalar, 1, TDerivedSize> result;
+
+        for (index i = 0; i < TDerivedSize; i++) {
+            result(i) = variable(TOffset + i, value(i));
+        }
+
+        return result;
+    }
+
     static TScalar f(const Scalar& variable)
     {
         return variable;
@@ -84,6 +105,13 @@ struct Space<1, TScalar, TSize> {
     template <index TRows, index TCols>
     using Matrix = Eigen::Matrix<Scalar, TRows, TCols>;
 
+    static Scalar empty()
+    {
+        static_assert(-1 <= TSize);
+
+        return Scalar::empty();
+    }
+
     static Scalar constant(TScalar value)
     {
         static_assert(-1 <= TSize);
@@ -100,6 +128,20 @@ struct Space<1, TScalar, TSize> {
         return result;
     }
 
+    template <index TOffset, index TDerivedSize>
+    HYPERJET_INLINE static Eigen::Matrix<Scalar, 1, TDerivedSize> variables(Eigen::Matrix<TScalar, 1, TDerivedSize> value)
+    {
+        static_assert(-1 <= TSize);
+
+        Eigen::Matrix<Scalar, 1, TDerivedSize> result;
+
+        for (index i = 0; i < TDerivedSize; i++) {
+            result(i) = variable(TOffset + i, value(i));
+        }
+
+        return result;
+    }
+
     static TScalar f(const Scalar& variable)
     {
         return variable.f();
@@ -150,6 +192,13 @@ struct Space<2, TScalar, TSize> {
     template <index TRows, index TCols>
     using Matrix = Eigen::Matrix<Scalar, TRows, TCols>;
 
+    static Scalar empty()
+    {
+        static_assert(-1 <= TSize);
+
+        return Scalar::empty();
+    }
+
     static Scalar constant(TScalar value)
     {
         static_assert(-1 <= TSize);
@@ -166,6 +215,20 @@ struct Space<2, TScalar, TSize> {
         return result;
     }
 
+    template <index TOffset, index TDerivedSize>
+    HYPERJET_INLINE static Eigen::Matrix<Scalar, 1, TDerivedSize> variables(Eigen::Matrix<TScalar, 1, TDerivedSize> value)
+    {
+        static_assert(-1 <= TSize);
+
+        Eigen::Matrix<Scalar, 1, TDerivedSize> result;
+
+        for (index i = 0; i < TDerivedSize; i++) {
+            result(i) = variable(TOffset + i, value(i));
+        }
+
+        return result;
+    }
+
     static TScalar f(const Scalar& variable)
     {
         return variable.f();

include/eqlib/LambdaConstraint.h

@@ -33,6 +33,7 @@ class LambdaConstraint : public Constraint {
 public: // methods
     void compute(Ref<Vector> rs, const std::vector<Ref<Vector>>& gs, const std::vector<Ref<Matrix>>& hs) const override
     {
+        pybind11::gil_scoped_acquire acquire;
         m_compute(m_equations, m_variables, rs, gs, hs);
     }
 

include/eqlib/LambdaObjective.h

@@ -30,6 +30,7 @@ class LambdaObjective : public Objective {
 public: // methods
     double compute(Ref<Vector> g, Ref<Matrix> h) const override
     {
+        pybind11::gil_scoped_acquire acquire;
         return m_compute(m_variables, g, h);
     }
 

include/eqlib/Problem.h

@@ -438,7 +438,11 @@ class Problem {
 
                 index index = m_structure_hm.get_index(row.global, col.global);
 
-                data.hm_value(index) += h(row.local, col.local);
+                if (row.local < col.local) {
+                    data.hm_value(index) += h(row.local, col.local);
+                } else {
+                    data.hm_value(index) += h(col.local, row.local);
+                }
             }
         }
 

include/eqlib/SparseStructure.h

@@ -78,13 +78,11 @@ class SparseStructure {
 
     double density() const noexcept
     {
-        const index size = rows() * cols();
-
-        if (size == 0) {
+        if (rows() == 0 || cols() == 0) {
             return 0;
         }
 
-        return (double)nb_nonzeros() / size;
+        return (double)nb_nonzeros() / rows() / cols();
     }
 
     const std::vector<TIndex>& ia() const noexcept