Solution postsolve: dual for the main constraint only #184

Also: fix storage of sparsity patterns in MP2NL #237

include/mp/arrayref.h

@@ -112,7 +112,7 @@ namespace mp {
     }
 
     template <class AR>
-    void init_from_rvalue(AR&& other) {
+    void init_from_rvalue(AR&& other) noexcept {
       if (other.save_.size()) {
         save_ = std::move(other.save_);
         data_ = save_.data();

include/mp/valcvt-base.h

@@ -375,6 +375,17 @@ class VCString {
 // etc...
 
 
+/// How to resolve the conflict
+/// when several values are gathered in one
+/// source or dest node
+enum class ValueResolution {
+  ValResAll,       // All values transferred
+                   // and SetVal() decides what to do
+                   // (can max out non-0 values.)
+  ValResDefault = ValResAll,
+  ValResLast       // just the last conection is transferred
+};
+
 /// ValuePresolver interface.
 /// Addresses value pre- / postsolve (solutions, basis, etc)
 class BasicValuePresolver : public EnvKeeper {
@@ -390,10 +401,12 @@ class BasicValuePresolver : public EnvKeeper {
 #define PRESOLVE_KIND(name, ValType) \
   virtual MVOverEl<ValType> \
     Presolve ## name ( \
-      const MVOverEl<ValType> & mv) = 0; \
+      const MVOverEl<ValType> & mv, \
+      ValueResolution = ValueResolution::ValResDefault) = 0; \
   virtual MVOverEl<ValType> \
     Postsolve ## name ( \
-      const MVOverEl<ValType> & mv) = 0;
+      const MVOverEl<ValType> & mv, \
+      ValueResolution = ValueResolution::ValResDefault) = 0;
 
   LIST_PRESOLVE_METHODS
 

include/mp/valcvt-link.h

@@ -45,8 +45,10 @@ class BasicLink {
   /// Postsolves should usually loop the range backwards
 #undef PRESOLVE_KIND
 #define PRESOLVE_KIND(name, ValType) \
-  virtual void Presolve ## name (LinkIndexRange ) = 0; \
-  virtual void Postsolve ## name(LinkIndexRange ) = 0;
+  virtual void Presolve ## name (LinkIndexRange, \
+    ValueResolution ) = 0; \
+  virtual void Postsolve ## name(LinkIndexRange, \
+    ValueResolution ) = 0;
 
   LIST_PRESOLVE_METHODS
 
@@ -159,9 +161,9 @@ class CopyLink : public BasicLink {
   /// Copy everything, MaxAmongNon0 should not apply
 #undef PRESOLVE_KIND
 #define PRESOLVE_KIND(name, ValType) \
-  void Presolve ## name (LinkIndexRange ir) override \
+  void Presolve ## name (LinkIndexRange ir, ValueResolution ) override \
   { CopySrcDest <ValType>(ir); } \
-  void Postsolve ## name(LinkIndexRange ir) override \
+  void Postsolve ## name(LinkIndexRange ir, ValueResolution ) override \
   { CopyDestSrc <ValType>(ir); }
 
   LIST_PRESOLVE_METHODS
@@ -228,10 +230,24 @@ class Many2ManyLink : public BasicLink {
            && entries_.back().first.TryExtendBy(be.first))
           )) {
       entries_.push_back(be);             // Add new entry
+      if_consider_for_last_link_.push_back(false);
       RegisterLinkIndex(entries_.size()-1);
     }
   }
 
+  /// Mark entry as contining the last link of a group.
+  /// needed for solution pre-postsolve where
+  /// we actually need the last link of a group,
+  /// corr. to the main constraint.
+  /// Example: a nonlinear con is redefined into a bunch
+  /// of smaller constraints, but only the main
+  /// resulting con receives/sends the dual value.
+  void MarkLastEntryAsLastInGroup() {
+    assert(if_consider_for_last_link_.size() == entries_.size());
+    assert(entries_.size());
+    if_consider_for_last_link_.back() = true;
+  }
+
   /// Get source/target nodes for a given link entry.
   /// This is used for graph export
   void ExportEntryItems(EntryItems& ei, int i) const override {
@@ -245,62 +261,80 @@ class Many2ManyLink : public BasicLink {
   /// All pre- / postsolves just take max from non-0
 #undef PRESOLVE_KIND
 #define PRESOLVE_KIND(name, ValType) \
-  void Presolve ## name (LinkIndexRange ir) override \
-  { DistributeFromSrc2Dest <ValType>(ir); } \
-  void Postsolve ## name(LinkIndexRange ir) override \
-  { CollectFromDest2Src <ValType>(ir); }
+  void Presolve ## name (LinkIndexRange ir, ValueResolution vr) override \
+  { DistributeFromSrc2Dest <ValType>(ir, vr); } \
+  void Postsolve ## name(LinkIndexRange ir, ValueResolution vr) override \
+  { CollectFromDest2Src <ValType>(ir, vr); }
 
   LIST_PRESOLVE_METHODS
 
 protected:
   /// Distribute values of type T from nr1 to nr2
   template <class T>
-  void Distr(NodeRange nr1, NodeRange nr2) {
+  void Distr(NodeRange nr1, NodeRange nr2, ValueResolution vr) {
     auto ir1 = nr1.GetIndexRange();
     auto ir2 = nr2.GetIndexRange();
     for (auto i0=ir1.beg_; i0!=ir1.end_; ++i0) {
       // Need reference here for reference counting
       // in PresolveNames():
       const auto& val = nr1.GetValueNode()->
           GetVal<T>(i0);
-      for (auto i=ir2.beg_; i!=ir2.end_; ++i)
-        nr2.GetValueNode()->SetVal(i, val);
+      if (ValueResolution::ValResAll == vr) {
+        for (auto i=ir2.beg_; i!=ir2.end_; ++i)
+          nr2.GetValueNode()->SetVal(i, val);
+      } else {
+        assert(ValueResolution::ValResLast == vr);
+        nr2.GetValueNode()->SetVal(ir2.end_ - 1, val);
+      }
     }
   }
 
   /// Collect values of type T from nr2 to nr1
   template <class T>
-  void Collect(NodeRange nr1, NodeRange nr2) {
+  void Collect(NodeRange nr1, NodeRange nr2, ValueResolution vr) {
     auto ir1 = nr1.GetIndexRange();
     auto ir2 = nr2.GetIndexRange();
     auto& vec2 = nr2.GetValueNode()->GetValVec<T>();
     for (auto i0=ir1.beg_; i0!=ir1.end_; ++i0) {
-      for (auto i=ir2.beg_; i!=ir2.end_; ++i)
-        nr1.GetValueNode()->SetVal(i0, vec2.at(i));
+      if (ValueResolution::ValResAll == vr) {
+        for (auto i=ir2.beg_; i!=ir2.end_; ++i) {
+          nr1.GetValueNode()->SetVal(i0, vec2.at(i));
+        }
+      } else {
+        assert(ValueResolution::ValResLast == vr);
+        nr1.GetValueNode()->SetVal(i0, vec2.at(ir2.end_ - 1));
+      }
     }
   }
 
   /// Src -> dest for the entries index range ir
   template <class T>
-  void DistributeFromSrc2Dest(LinkIndexRange ir) {
+  void DistributeFromSrc2Dest(LinkIndexRange ir, ValueResolution vr) {
     for (int i=ir.beg_; i!=ir.end_; ++i) {
-      const auto& br = entries_[i];
-      Distr<T>(br.first, br.second);
+      if (ValueResolution::ValResAll == vr
+          || if_consider_for_last_link_[i]) {
+        const auto& br = entries_[i];
+        Distr<T>(br.first, br.second, vr);
+      }
     }
   }
 
   /// Collect src <- dest for index range ir. Loop backwards
   template <class T>
-  void CollectFromDest2Src(LinkIndexRange ir) {
+  void CollectFromDest2Src(LinkIndexRange ir, ValueResolution vr) {
     for (int i=ir.end_; (i--)!=ir.beg_; ) {
-      const auto& br = entries_[i];
-      Collect<T>(br.first, br.second);
+      if (ValueResolution::ValResAll == vr
+          || if_consider_for_last_link_[i]) {
+        const auto& br = entries_[i];
+        Collect<T>(br.first, br.second, vr);
+      }
     }
   }
 
 
 private:
   CollectionOfEntries entries_;
+  std::vector<bool> if_consider_for_last_link_;
 };
 
 
@@ -387,10 +421,10 @@ class BasicIndivEntryLink : public BasicLink {
   /// and calls the derived class' method for each.
 #undef PRESOLVE_KIND
 #define PRESOLVE_KIND(name, ValType) \
-  void Presolve ## name (LinkIndexRange ir) override { \
+  void Presolve ## name (LinkIndexRange ir, ValueResolution) override { \
     for (int i=ir.beg_; i!=ir.end_; ++i) \
       MPD( Presolve ## name ## Entry(entries_.at(i)) ); } \
-  void Postsolve ## name(LinkIndexRange ir) override { \
+  void Postsolve ## name(LinkIndexRange ir, ValueResolution) override { \
     for (int i=ir.end_; i--!=ir.beg_; ) \
       MPD( Postsolve ## name ## Entry(entries_.at(i)) ); }
 
@@ -514,6 +548,7 @@ class AutoLinkScope {
           cvt_.GetOne2ManyLink().AddEntry(   // use One2ManyLink
                 { cvt_.GetAutoLinkSource(), t } );
         }
+        cvt_.GetOne2ManyLink().MarkLastEntryAsLastInGroup();
       }
     }
     cvt_.TurnOffAutoLinking();

include/mp/valcvt-node.h

@@ -156,7 +156,7 @@ class ValueNode {
 
   /////////////////////// Access value vectors ///////////////////////
 
-  /// Assign from ArrayRef<int>. Always copy the values.
+  /// Assign from vector<int>.
   /// Some (target) node assignments can be longer vectors:
   /// e.g., Gurobi adds variables for FeasRelax.
   ValueNode& operator=(std::vector<int> ai)
@@ -382,7 +382,7 @@ ValueNode CreateArray(BasicValuePresolver& vp) { return ValueNode{vp}; }
 /// @return always true currently
 template <class Vec>
 inline
-bool CopyRange(Vec& src, Vec& dest, NodeIndexRange ir, int i1) {
+bool CopyRange(const Vec& src, Vec& dest, NodeIndexRange ir, int i1) {
   assert(ir.end_ <= (int)src.size());
   assert(i1 + ir.Size() <= (int)dest.size());
   std::copy(src.begin()+ir.beg_, src.begin()+ir.end_,

include/mp/valcvt.h

@@ -104,13 +104,15 @@ class ValuePresolverImpl : public BasicValuePresolver {
 #define PRESOLVE_KIND(name, ValType) \
   MVOverEl<ValType> \
     Presolve ## name ( \
-      const MVOverEl<ValType> & mv) override { \
-    return RunPresolve(&BasicLink::Presolve ## name, mv); \
+      const MVOverEl<ValType> & mv, \
+      ValueResolution vr = ValueResolution::ValResDefault) override { \
+    return RunPresolve(&BasicLink::Presolve ## name, mv, vr); \
   } \
   MVOverEl<ValType> \
     Postsolve ## name ( \
-      const MVOverEl<ValType> & mv) override { \
-    return RunPostsolve(&BasicLink::Postsolve ## name, mv); \
+      const MVOverEl<ValType> & mv, \
+      ValueResolution vr = ValueResolution::ValResDefault) override { \
+    return RunPostsolve(&BasicLink::Postsolve ## name, mv, vr); \
   }
 
   LIST_PRESOLVE_METHODS
@@ -131,25 +133,27 @@ class ValuePresolverImpl : public BasicValuePresolver {
 
 protected:
   /// Helper type: virtual member function pointer
-  using LinkFn = void (BasicLink::*)(LinkIndexRange);
+  using LinkFn = void (BasicLink::*)(LinkIndexRange, ValueResolution);
 
   /// Generic value presolve loop: forward
   template <class ModelValues>
-  ModelValues RunPresolve(LinkFn fn, const ModelValues& mv) const {
+  ModelValues RunPresolve(LinkFn fn, const ModelValues& mv,
+                          ValueResolution vr) const {
     CleanUpValueNodes();
     src_ = mv;
     for (const auto& br: brl_)
-      (br.b_.*fn)(br.ir_);
+      (br.b_.*fn)(br.ir_, vr);
     return dest_;
   }
 
   /// Generic value postsolve loop: backward
   template <class ModelValues>
-  ModelValues RunPostsolve(LinkFn fn, const ModelValues& mv) const {
+  ModelValues RunPostsolve(LinkFn fn, const ModelValues& mv,
+                           ValueResolution vr) const {
     CleanUpValueNodes();
     dest_ = mv;
     for (auto rit=brl_.rbegin(); rit!=brl_.rend(); ++rit) {
-      (rit->b_.*fn)(rit->ir_);
+      (rit->b_.*fn)(rit->ir_, vr);
     }
     return src_;
   }
@@ -261,9 +265,12 @@ class ValuePresolver : public ValuePresolverImpl {
 
   /// Override PresolveSolution().
   /// Update warm start values to fit into their bounds.
+  /// @param vr: ignored for simplicity
   MVOverEl<double> PresolveSolution (
-      const MVOverEl<double>& mv) override {
-    auto result = ValuePresolverImpl::PresolveSolution(mv);
+      const MVOverEl<double>& mv,
+      ValueResolution  = ValueResolution::ValResLast) override {
+    auto result = ValuePresolverImpl::PresolveSolution(
+        mv, ValueResolution::ValResLast);
     auto& x = result.GetVarValues()();
     const auto& lbs = model_.GetVarLBs();
     const auto& ubs = model_.GetVarUBs();
@@ -283,8 +290,10 @@ class ValuePresolver : public ValuePresolverImpl {
   /// 1. Call preposts_() if provided.
   /// 2. Check solution if checker provided.
   /// 3. mv's ExtraData() is passed to the checker.
+  /// @param vr: ignored
   MVOverEl<double> PostsolveSolution (
-      const MVOverEl<double>& mv) override {
+      const MVOverEl<double>& mv,
+      ValueResolution  = ValueResolution::ValResLast) override {
     auto mv_copy = mv;
     if (preposts_)
       preposts_(mv_copy);
@@ -299,7 +308,8 @@ class ValuePresolver : public ValuePresolverImpl {
         solchk_(mx(), mv_copy.GetConValues(), objs, mv_copy.ExtraData());
       }
     }
-    return ValuePresolverImpl::PostsolveSolution(mv_copy);
+    return ValuePresolverImpl::PostsolveSolution(
+        mv_copy, ValueResolution::ValResLast);
   }
 
 

solvers/mp2nl/mp2nlmodelapi.h

@@ -1335,6 +1335,16 @@ class MP2NLModelAPI
     LinPartRefOrStorage(ArrayRef<double> c, ArrayRef<int> v)
         : coefs_(std::move(c)), vars_(std::move(v))
     { assert(Check()); }
+    /// Copy construct: delete, too risky with ArrayRef
+    LinPartRefOrStorage(const LinPartRefOrStorage& lp) noexcept = delete;
+    /// operator=(const &): delete, too risky
+    LinPartRefOrStorage& operator=(const LinPartRefOrStorage& lp) noexcept
+        = delete;
+    /// Move construct
+    LinPartRefOrStorage(LinPartRefOrStorage&& ) noexcept = default;
+    /// operator=(&&)
+    LinPartRefOrStorage& operator=(LinPartRefOrStorage&& ) noexcept
+        = default;
     /// Check
     bool Check() const { return coefs_.size()==vars_.size(); }
     /// Size
@@ -1359,35 +1369,28 @@ class MP2NLModelAPI
              bool fLogical,
              bool f_m = 0, LinPartRefOrStorage lpe = {}
              //, StaticItemTypeID iid, ExpressionTypeID eid
-             )
+             ) noexcept
         : disp_(disp), p_item_(pitem), f_logical_(fLogical),
-        f_marked_(f_m), lin_part_ext_(lpe)
+        f_marked_(f_m), lin_part_ext_(std::move(lpe))
         // no storing, itemID_(iid), exprID_(eid)
     { }
-    /// Copy construct
-    ItemInfo(const ItemInfo& ii)
-        : ItemInfo(ii.disp_, ii.p_item_, ii.f_logical_,
-                   ii.f_marked_, ii.lin_part_ext_) { }
-    /// operator=
-    ItemInfo& operator=(const ItemInfo& ii) {
-      assert(&disp_==&ii.disp_);
-      p_item_ = ii.p_item_;
-      assert(f_logical_==ii.f_logical_);
-      f_marked_ = ii.f_marked_;
-      lin_part_ext_ = ii.lin_part_ext_;
-      return *this;
-    }
+    /// Copy construct: delete, too risky with \a lin_part_ext_
+    ItemInfo(const ItemInfo& ii) = delete;
+    /// operator=(const &): delete
+    ItemInfo& operator=(const ItemInfo& ii) = delete;
     /// Move construct
-    ItemInfo(ItemInfo&& ii)
+    ItemInfo(ItemInfo&& ii) noexcept
         : ItemInfo(ii.disp_, ii.p_item_, ii.f_logical_,
                    ii.f_marked_, std::move(ii.lin_part_ext_)) { }
     /// operator=(&&)
-    ItemInfo& operator=(ItemInfo&& ii) {
+    ItemInfo& operator=(ItemInfo&& ii) noexcept {
       assert(&disp_==&ii.disp_);
       p_item_ = ii.p_item_;
       assert(f_logical_==ii.f_logical_);
       f_marked_ = ii.f_marked_;
+      auto sz = ii.lin_part_ext_.size();
       lin_part_ext_ = std::move(ii.lin_part_ext_);
+      assert(lin_part_ext_.size() == sz);
       return *this;
     }
     /// Get dispatcher