Multi-obj: HiGHS .objweight suffix 'intuitive way' #240 #239

include/mp/backend-std.h

@@ -287,7 +287,7 @@ class StdBackend :
 
   /// Standard extras
   virtual void InputStdExtras() {
-    if (multiobj()) {
+    if (multiobj() && multiobj_has_native()) {
       if (auto suf = ReadSuffix(suf_objpriority))
         ObjPriorities( suf );
       if (auto suf = ReadSuffix(suf_objweight))

include/mp/flat/converter.h

@@ -10,6 +10,7 @@
 #include "mp/env.h"
 #include "mp/format.h"
 #include "mp/solver-base.h"
+#include "mp/suffix.h"
 #include "mp/flat/converter_model.h"
 #include "mp/flat/convert_functional.h"
 #include "mp/flat/constr_keeper.h"
@@ -1294,13 +1295,33 @@ class FlatConverter :
           std::move(key), std::move(msg), replace);
   }
 
+  /// Provide suffix getters
+  void SetSuffixGetters(
+      std::function<ArrayRef<int>(const SuffixDef<int>& )> sgi,
+      std::function<ArrayRef<double>(const SuffixDef<double>& )> sgd)
+  { suf_get_int_=sgi; suf_get_dbl_=sgd; }
+
+
+public:
+  /// Read int suffix
+  ArrayRef<int> ReadIntSuffix(const SuffixDef<int>& sd)
+  { assert(suf_get_int_); return suf_get_int_(sd); }
+
+  /// Read double suffix
+  ArrayRef<double> ReadDblSuffix(const SuffixDef<double>& sd)
+  { assert(suf_get_dbl_); return suf_get_dbl_(sd); }
+
 
 private:
   /// We store ModelApi in the converter for speed.
   /// Should be before constraints
   ModelAPIType modelapi_;
   /// solve iteration
   int n_solve_iter_ {0};
+  /// Suffix getter int
+  std::function<ArrayRef<int>(const SuffixDef<int>& )> suf_get_int_;
+  /// Suffix getter double
+  std::function<ArrayRef<double>(const SuffixDef<double>& )> suf_get_dbl_;
   /// ValuePresolver: should be init before constraint keepers
   /// and links
   pre::ValuePresolver value_presolver_

include/mp/flat/converter_multiobj.h

@@ -104,15 +104,56 @@ class MOManager {
     MPD(set_skip_pushing_objs());  // could have a cleaner system of linking
       // via custom link restoring original objective values,
       // instead of current manual postsolving in ValuePresolver::PostsolveSolution().
-    obj_new_ = MPD( get_objectives() );   // no linking
+    const auto& obj_orig = MPD( get_objectives() );   // no linking
+    ///////////////// Read / set default suffixes ///////////////////
+    std::vector<double> objpr = MPD( ReadDblSuffix( {"objpriority", suf::OBJ} ) );
+    objpr.resize(obj_orig.size(), 0.0);               // blend objectives by default
+    std::vector<double> objwgt = MPD( ReadDblSuffix( {"objweight", suf::OBJ} ) );
+    objwgt.resize(obj_orig.size(), 1.0);
+    std::vector<double> objtola = MPD( ReadDblSuffix( {"objabstol", suf::OBJ} ) );
+    objtola.resize(obj_orig.size(), 0.0);
+    std::vector<double> objtolr = MPD( ReadDblSuffix( {"objreltol", suf::OBJ} ) );
+    objtolr.resize(obj_orig.size(), 0.0);
+    std::map<double, std::vector<int>, std::greater<double> > pr_map;      // Decreasing order
+    for (int i=0; i<objpr.size(); ++i)
+      pr_map[objpr[i]].push_back(i);
+    obj_new_ = {};         ////////////////// Aggregate new objectives ///////////////////
+    obj_new_.reserve(pr_map.size());
+    obj_new_tola_.reserve(pr_map.size());
+    obj_new_tolr_.reserve(pr_map.size());
+    for (const auto& pr_level: pr_map) {
+      const auto& i0_vec = pr_level.second;
+      obj_new_.push_back(obj_orig.at(i0_vec.front()));
+      obj_new_.back().GetLinTerms() *= objwgt.at(i0_vec.front());   // Use weight
+      obj_new_.back().GetQPTerms() *= objwgt.at(i0_vec.front());
+      obj_new_tola_.push_back(objtola.at(i0_vec.front()));
+      obj_new_tolr_.push_back(objtolr.at(i0_vec.front()));
+      for (auto i0i=i0_vec.size(); --i0i; ) {
+        // Add next objective with weight and sense factor
+        double sensef
+            = (obj_orig.at(i0_vec.front()).obj_sense() == obj_orig.at(i0_vec[i0i]).obj_sense())
+                  ? 1.0 : -1.0;
+        auto lt1 = obj_orig.at(i0_vec[i0i]).GetLinTerms();
+        lt1 *= sensef * objwgt.at(i0_vec[i0i]);
+        auto qt1 = obj_orig.at(i0_vec[i0i]).GetQPTerms();
+        qt1 *= sensef * objwgt.at(i0_vec[i0i]);
+        obj_new_.back().GetLinTerms().add(lt1);
+        obj_new_.back().GetQPTerms().add(qt1);
+        // Max the tolerances
+        obj_new_tola_.back() = std::max(obj_new_tola_.back(), objtola.at(i0_vec[i0i]));
+        obj_new_tolr_.back() = std::max(obj_new_tolr_.back(), objtolr.at(i0_vec[i0i]));
+      }
+      obj_new_.back().GetLinTerms().sort_terms();
+      obj_new_.back().GetQPTerms().sort_terms();
+    }
     if (MPD( GetEnv() ).verbose_mode())
       MPD( GetEnv() ).Print(
           "\n\n"
           "==============================================================================\n"
-          "MULTI-OBJECTIVE MODE: starting with {} objectives ...\n"
+          "MULTI-OBJECTIVE MODE: starting with {} objectives ({} combined) ...\n"
           "==============================================================================\n"
           "==============================================================================\n\n"
-          , obj_new_.size());
+          , obj_orig.size(), obj_new_.size());
   }
 
   /// Do prepare next solve
@@ -143,6 +184,10 @@ class MOManager {
   void RestrictLastObjVal() {
     assert(i_current_obj_ && i_current_obj_<obj_new_.size());
     const auto& obj_last = obj_new_[i_current_obj_-1];
+    auto lim = objval_last_;
+    auto diff = std::max(                 // Apply degradation tolerance
+        obj_new_tola_[i_current_obj_-1], std::fabs(lim) * obj_new_tolr_[i_current_obj_-1]);
+    lim += diff * (obj::MAX==obj_last.obj_sense() ? -1.0 : 1.0);
     if (obj_last.GetQPTerms().size()) {
       if (obj::MAX == obj_last.obj_sense())
         MPD( GetModelAPI() ).AddConstraint(
@@ -169,6 +214,8 @@ class MOManager {
 private:
   MOManagerStatus status_ {MOManagerStatus::NOT_SET};
   std::vector<QuadraticObjective> obj_new_;     // ranked aggregated objectives
+  std::vector<double> obj_new_tola_;
+  std::vector<double> obj_new_tolr_;
   int i_current_obj_ {-1};
   double objval_last_ {};
 };

include/mp/flat/problem_flattener.h

@@ -118,6 +118,9 @@ class ProblemFlattener :
 public:
   /// Convert the whole model, e.g., after reading from NL
   void ConvertModel() override {
+    GetFlatCvt().SetSuffixGetters(
+        [this](const SuffixDef<int>& sd) { return GetModel().ReadIntSuffix(sd); },
+        [this](const SuffixDef<double>& sd) { return GetModel().ReadDblSuffix(sd); });
     GetFlatCvt().StartModelInput();
     MP_DISPATCH( ConvertStandardItems() );
     GetFlatCvt().FinishModelInput();      // Chance to flush to the Backend

include/mp/problem.h

@@ -1220,8 +1220,10 @@ class BasicProblem : public ExprFactory, public SuffixManager {
   /// take that
   ArrayRef<double> ReadDblSuffix(const SuffixDef<double>& sufdef) {
     auto suf_dbl = ReadSuffix_OneTypeOnly(sufdef);
+    std::printf("   SUFD %s: %ld elements\n", sufdef.name(), suf_dbl.size());
     if (!suf_dbl) {
       auto suf_int = ReadSuffix_OneTypeOnly(sufdef.to_type<int>());
+      std::printf("        SUFI %s: %ld elements\n", sufdef.name(), suf_int.size());
       if (suf_int)
         return std::vector<double>(suf_int.begin(), suf_int.end());
     }

include/mp/solver-opt.h

@@ -146,6 +146,7 @@ class SolverOption {
   { return inline_synonyms_; }
   /// Add additional "inline" synonyms
   void add_synonyms_front(const char* names_list);
+  /// Add "inline" synonyms to the back
   void add_synonyms_back(const char* names_list);
 
   /// Is hidden?
@@ -155,19 +156,25 @@ class SolverOption {
   bool is_wildcard() const { return wc_headtails_.size(); }
   /// Checks if matches, then saves key & body
   bool wc_match(const std::string& key);
+  /// Wildcard head
   const std::string& wc_head() const
   { assert(is_wildcard()); return wc_headtails_[0].first; }
+  /// Wildcard tail
   const std::string& wc_tail() const
   { assert(is_wildcard()); return wc_headtails_[0].second; }
+  /// Last wildcard key
   const std::string& wc_key_last() const { return wc_key_last_; }
+  /// Last wildcard keybody
   const std::string& wc_keybody_last() const { return wc_body_last_; }
   /// Printing last parsed wc key in std form
   std::string wc_key_last__std_form() const
   { return wc_head() + wc_body_last_ + wc_tail(); }
 
-  /// Return/set the option description.
+  /// Option description.
   const char *description() const { return description_.c_str(); }
+  /// Set option description
   void set_description(const char* d) { description_=d; }
+  /// Add to option description
   void add_to_description(const char* d) { description_ += d; }
 
   /// Append the formatted description to the writer

solvers/gurobi/gurobibackend.cc

@@ -2101,16 +2101,7 @@ void GurobiBackend::InitCustomOptions() {
   /// Option "multiobj" is created internally if
   /// std feature MULTIOBJ is set.
   /// Change the help text
-  ReplaceOptionDescription("obj:multi",
-                           "0*/1: Whether to do multi-objective optimization.\n"
-                           "When obj:multi = 1 and several objectives are present, suffixes "
-                           ".objpriority, .objweight, .objreltol, and .objabstol on the "
-                           "objectives are relevant.  Objectives with greater .objpriority "
-                           "values (integer values) have higher priority.  Objectives with "
-                           "the same .objpriority are weighted by .objweight.  Objectives "
-                           "with positive .objabstol or .objreltol are allowed to be "
-                           "degraded by lower priority objectives by amounts not exceeding "
-                           "the .objabstol (absolute) and .objreltol (relative) limits. "
+  AddToOptionDescription("obj:multi",
                            "The objectives must all be linear.  Objective-specific "
                            "convergence tolerances and method values may be assigned via "
                            "keywords of the form obj_n_<name>, such as obj_1_method for the "

solvers/highsmp/highsmpmodelapi.cc

@@ -36,6 +36,7 @@ void HighsModelAPI::SetLinearObjective( int iobj, const LinearObjective& lo ) {
     for (auto i=lo.vars().size(); i--; )
       objc_new[lo.vars()[i]] = lo.coefs()[i];
     HIGHS_CCALL(Highs_changeColsCostByRange(lp(), 0, NumVars()-1, objc_new.data()));
+    std::printf("     HIGHS OPBJ SENSE:  %d \n", lo.obj_sense());
     HIGHS_CCALL(Highs_changeObjectiveSense(lp(), 
                                            obj::Type::MAX==lo.obj_sense() ?
                                                kHighsObjSenseMaximize : kHighsObjSenseMinimize) );

solvers/visitor/visitormodelapi.h

@@ -31,12 +31,14 @@ class VisitorModelAPI :
   /// After
   void FinishProblemModificationPhase();
 
-  /// TODO Implement the following functions using the solver's API
+  /// TODO Implement adding variables
   void AddVariables(const VarArrayDef& );
+  /// TODO Implement setting (also changing) a linear (part of the) objective
   void SetLinearObjective( int iobj, const LinearObjective& lo );
   /// Whether accepting quadratic objectives:
   /// 0 - no, 1 - convex, 2 - nonconvex
   static int AcceptsQuadObj() { return 0; }
+  /// TODO Implement setting (also changing) a quadratic objective
   void SetQuadraticObjective(int iobj, const QuadraticObjective& qo);
 
   //////////////////////////// GENERAL CONSTRAINTS ////////////////////////////

src/solver.cc

@@ -646,10 +646,16 @@ void BasicSolver::InitMetaInfoAndOptions(
 
   if ((flags & MULTIPLE_OBJ) != 0) {
     AddOption(OptionPtr(new BoolOption(multiobj_, "obj:multi multiobj",
-      "0*/1:  Whether to use multi-objective optimization. "
-      "If set to 1 multi-objective optimization is performed using "
-      "lexicographic method with the first objective treated as the most "
-      "important, then the second objective and so on.")));
+      "0*/1:  Whether to use multi-objective optimization.\n"
+                                       "\n"
+                                       "When obj:multi = 1 and several objectives are present, suffixes "
+                                       ".objpriority, .objweight, .objreltol, and .objabstol on the "
+                                       "objectives are relevant.  Objectives with greater .objpriority "
+                                       "values (integer values) have higher priority.  Objectives with "
+                                       "the same .objpriority are weighted by .objweight.  Objectives "
+                                       "with positive .objabstol or .objreltol are allowed to be "
+                                       "degraded by lower priority objectives by amounts not exceeding "
+                                       "the .objabstol (absolute) and .objreltol (relative) limits. "                                       )));
   }
 
   AddIntOption("tech:timing timing tech:report_times report_times",

test/end2end/cases/categorized/fast/multi_obj/modellist.json

@@ -45,8 +45,23 @@
       "_sobj[4]": 53.19255547602
     }
   },
+  {
+    "name" : "dietobj_1000 multiobj=1 obj_pr_1",
+    "Comment": "Modify obj priorities and weights via suffixes in dietobj_pr_1",
+    "tags" : ["linear", "continuous", "multiobj"],
+    "files" : ["dietobj_1000.mod", "dietobj.dat", "dietobj_pr_1.inc"],
+    "options": { "ANYSOLVER_options": "multiobj=1" },
+    "values": {
+      "total_cost[\"A&P\"]": 924.1546153846151,
+      "total_cost[\"JEWEL\"]": 925.6467032967034,
+      "total_cost[\"VONS\"]": 918.8280219780219,
+      "total_number": 32.604395604395734,
+      "_sobj[4]": 32.604395604395734
+    }
+  },
   {
     "name" : "dietobj multiobj=1 obj:2:priority=10",
+    "Comment": "Modify obj priority via an option",
     "tags" : ["linear", "continuous", "multiobj", "obj_priority"],
     "files" : ["dietobj.mod", "dietobj.dat"],
     "options": { "ANYSOLVER_options": "multiobj=1 obj:2:priority=10" },

test/end2end/scripts/python/Model.py

@@ -51,7 +51,7 @@ class ModelTags(enum.Enum):
     mipstart=20008
 
     multiobj=30009
-    obj_priority=30010
+    obj_priority=30010           # understands obj:2:priority etc.
 
     multisol = 40011
     sstatus = 40012              # Basis I/O

test/end2end/scripts/python/Solver.py

@@ -474,7 +474,11 @@ def __init__(self, exeName, timeout=None, nthreads=None, otherOptions=None):
         stags = {ModelTags.continuous, ModelTags.integer, ModelTags.binary,
                   ModelTags.linear, ModelTags.quadratic, ModelTags.sos,
                  ModelTags.return_mipgap,
-                 ModelTags.sstatus}
+                 ModelTags.sstatus,
+                 ModelTags.multisol,
+                 ModelTags.multiobj,
+                 ModelTags.iis,
+                 ModelTags.feasrelax}
         super().__init__(exeName, timeout, nthreads, otherOptions, stags)
 
     def _doParseSolution(self, st, stdout=None):
@@ -828,7 +832,8 @@ def __init__(self, exeName, timeout=None, nthreads=None,
                     ModelTags.logical,
                     ModelTags.plinear,
                     ModelTags.nonlinear,
-                    ModelTags.log
+                    ModelTags.log,
+                    ModelTags.multiobj
                 }
             # Direct/FlatConverter drivers with non-convex quadratics:
             if ModelTags.quadraticnonconvex in stags:
@@ -980,7 +985,7 @@ def __init__(self, exeName, timeout=None, nthreads=None,
                  ModelTags.return_mipgap,
 
                  ModelTags.warmstart, ModelTags.mipstart,
-                 ModelTags.multiobj, ModelTags.obj_priority,
+                 ModelTags.multiobj,
                  ModelTags.multisol, ModelTags.sstatus,
                  ModelTags.iis, ModelTags.fixmodel,