Start implementing branch and bound

LoopModels · Sep 12, 2023 · 76041b7 · 76041b7
1 parent dd2f0ef
commit 76041b7
Show file tree

Hide file tree

Showing 4 changed files with 59 additions and 25 deletions.
diff --git a/include/Math/Array.hpp b/include/Math/Array.hpp
@@ -89,7 +89,7 @@ template <class T, class S> struct POLY_MATH_GSL_POINTER Array {
   constexpr auto operator=(Array &&) noexcept -> Array & = default;
   constexpr Array(const T *p, S s) : ptr(p), sz(s) {}
   constexpr Array(NotNull<const T> p, S s) : ptr(p), sz(s) {}
-  constexpr Array(const T *p, Row r, Col c) : ptr(p), sz(dimension<S>(r, c)) {}
+  constexpr Array(const T *p, Row r, Col c) : ptr(p), sz(S{r, c}) {}
   constexpr Array(NotNull<const T> p, Row r, Col c)
     : ptr(p), sz(dimension<S>(r, c)) {}
   template <std::convertible_to<S> V>

diff --git a/include/Math/BoxOpt.hpp b/include/Math/BoxOpt.hpp
@@ -32,6 +32,12 @@ class BoxTransformView {
   [[nodiscard]] constexpr auto view() const -> BoxTransformView {
     return *this;
   }
+  [[nodiscard]] constexpr auto getLowerBounds() -> MutPtrVector<int32_t> {
+    return {i32 + Ntotal, Ntotal};
+  }
+  [[nodiscard]] constexpr auto getUpperBounds() -> MutPtrVector<int32_t> {
+    return {i32 + ptrdiff_t(2) * Ntotal, Ntotal};
+  }
 
 protected:
   /// Ntotal offsets
@@ -62,12 +68,6 @@ class BoxTransformView {
   [[nodiscard]] constexpr auto getInds() -> MutPtrVector<int32_t> {
     return {i32, Ntotal};
   }
-  [[nodiscard]] constexpr auto getLowerBounds() -> MutPtrVector<int32_t> {
-    return {i32 + Ntotal, Ntotal};
-  }
-  [[nodiscard]] constexpr auto getUpperBounds() -> MutPtrVector<int32_t> {
-    return {i32 + ptrdiff_t(2) * Ntotal, Ntotal};
-  }
   [[nodiscard]] constexpr auto offs() -> MutPtrVector<double> {
     return {f64, Ntotal};
   }
@@ -153,6 +153,10 @@ class BoxTransform : public BoxTransformView {
     offs()[idx] = newOff;
   }
   constexpr ~BoxTransform() {
+    if (f64 == nullptr) {
+      invariant(i32 == nullptr);
+      return;
+    }
     delete[] f64;
     delete[] i32;
   }

diff --git a/include/Math/BoxOptInt.hpp b/include/Math/BoxOptInt.hpp
@@ -0,0 +1,32 @@
+#pragma once
+#include "Math/Array.hpp"
+#include "Math/BoxOpt.hpp"
+#include "Math/Constructors.hpp"
+#include <cstdint>
+
+namespace poly::math {
+
+// recursive branch and bound
+constexpr auto branchandbound(utils::Arena<> *alloc, MutPtrVector<double> x,
+                              BoxTransformView trf, double lower, double upper,
+                              const auto &f) -> double {}
+// set floor and ceil
+// Assumes that integer floor of values is optimal
+constexpr auto minimizeIntSol(utils::Arena<> *alloc, MutPtrVector<int32_t> r,
+                              int32_t lb, int32_t ub, const auto &f) {
+  // goal is to shrink all bounds such that lb==ub, i.e. we have all
+  // integer solutions.
+  unsigned N = r.size();
+  auto s = alloc->scope();
+  MutPtrVector<double> x{vector<double>(alloc, N)};
+  // MutDensePtrMatrix<double> X{matrix<double>(alloc, Row{2}, Col{N})};
+  poly::math::BoxTransform box{N, lb, ub};
+  double lower = minimize(alloc, x, box, f);
+  double upper = f(Elementwise{[](auto &x) { return std::floor(x); },
+                               BoxTransformVector{x, box}});
+  double opt = branchandbound(alloc, x, box, lower, upper, f);
+  r << box.getLowerBounds();
+  return opt;
+}
+
+} // namespace poly::math
diff --git a/include/Math/Math.hpp b/include/Math/Math.hpp
@@ -88,7 +88,7 @@ concept BinaryFuncOfElts =
   return true;
 }
 
-template <typename Op, typename A> struct ElementwiseUnaryOp {
+template <typename Op, typename A> struct Elementwise {
   using value_type =
     decltype(std::declval<Op>()(std::declval<utils::eltype_t<A>>()));
   [[no_unique_address]] Op op;
@@ -382,9 +382,8 @@ static_assert(!AbstractMatrix<StridedVector<int64_t>>);
 
 // static_assert(std::is_trivially_copyable_v<MutStridedVector<int64_t>>);
 static_assert(std::is_trivially_copyable_v<
-              ElementwiseUnaryOp<std::negate<>, StridedVector<int64_t>>>);
-static_assert(
-  Trivial<ElementwiseUnaryOp<std::negate<>, StridedVector<int64_t>>>);
+              Elementwise<std::negate<>, StridedVector<int64_t>>>);
+static_assert(Trivial<Elementwise<std::negate<>, StridedVector<int64_t>>>);
 
 constexpr void swap(MutPtrMatrix<int64_t> A, Row i, Row j) {
   if (i == j) return;
@@ -513,38 +512,37 @@ inline auto operator<<(std::ostream &os, const T &A) -> std::ostream & {
 
 constexpr auto operator-(const AbstractVector auto &a) {
   auto AA{a.view()};
-  return ElementwiseUnaryOp<std::negate<>, decltype(AA)>{.op = std::negate<>{},
-                                                         .a = AA};
+  return Elementwise<std::negate<>, decltype(AA)>{.op = std::negate<>{},
+                                                  .a = AA};
 }
 constexpr auto operator-(const AbstractMatrix auto &a) {
   auto AA{a.view()};
-  return ElementwiseUnaryOp<std::negate<>, decltype(AA)>{.op = std::negate<>{},
-                                                         .a = AA};
+  return Elementwise<std::negate<>, decltype(AA)>{.op = std::negate<>{},
+                                                  .a = AA};
 }
 
 constexpr auto operator!(const AbstractVector auto &a) {
   auto AA{a.view()};
-  return ElementwiseUnaryOp<std::logical_not<>, decltype(AA)>{
-    .op = std::negate<>{}, .a = AA};
+  return Elementwise<std::logical_not<>, decltype(AA)>{.op = std::negate<>{},
+                                                       .a = AA};
 }
 constexpr auto operator!(const AbstractMatrix auto &a) {
   auto AA{a.view()};
-  return ElementwiseUnaryOp<std::logical_not<>, decltype(AA)>{
-    .op = std::negate<>{}, .a = AA};
+  return Elementwise<std::logical_not<>, decltype(AA)>{.op = std::negate<>{},
+                                                       .a = AA};
 }
 
 constexpr auto operator~(const AbstractVector auto &a) {
   auto AA{a.view()};
-  return ElementwiseUnaryOp<std::bit_not<>, decltype(AA)>{.op = std::negate<>{},
-                                                          .a = AA};
+  return Elementwise<std::bit_not<>, decltype(AA)>{.op = std::negate<>{},
+                                                   .a = AA};
 }
 constexpr auto operator~(const AbstractMatrix auto &a) {
   auto AA{a.view()};
-  return ElementwiseUnaryOp<std::bit_not<>, decltype(AA)>{.op = std::negate<>{},
-                                                          .a = AA};
+  return Elementwise<std::bit_not<>, decltype(AA)>{.op = std::negate<>{},
+                                                   .a = AA};
 }
-static_assert(
-  AbstractMatrix<ElementwiseUnaryOp<std::negate<>, PtrMatrix<int64_t>>>);
+static_assert(AbstractMatrix<Elementwise<std::negate<>, PtrMatrix<int64_t>>>);
 static_assert(AbstractMatrix<Array<int64_t, SquareDims>>);
 static_assert(AbstractMatrix<ManagedArray<int64_t, SquareDims>>);