Implement force budget option and acyclic SparseMAP implementation

ad3qp/ad3/Factor.cpp

@@ -1013,7 +1013,7 @@ FactorBUDGET::SolveMAP(const vector<double>& variable_log_potentials,
         sum = 0.0;
         for (size_t k = 0; k < GetBudget(); ++k) {
             valaux = -scores[k].first;
-            if (valaux < 0.0)
+            if (valaux < 0.0 && !ForcedBudget())
                 break;
             int f = scores[k].second;
             (*variable_posteriors)[f] = negated_[f] ? 0.0 : 1.0;
@@ -1045,7 +1045,7 @@ FactorBUDGET::SolveQP(const vector<double>& variable_log_potentials,
     double budget = GetBudget();
     tight_ = false;
 
-    if (s > budget) {
+    if (s > budget || ForcedBudget()) {
         tight_ = true;
 
         FlipNegatedForQP(variable_log_potentials, variable_posteriors);

ad3qp/ad3/Factor.h

@@ -25,6 +25,7 @@
 #include <iostream>
 #include <stdexcept>
 #include <vector>
+#include <sstream>
 
 using namespace std;
 
@@ -190,7 +191,13 @@ class Factor
     void SetAdditionalLogPotentials(
       const vector<double>& additional_log_potentials)
     {
-        assert(additional_log_potentials.size() == GetNumAdditionals());
+        if (additional_log_potentials.size() != GetNumAdditionals()) {
+	    std::stringstream ss;
+	    ss << "Invalid size of additional log potentials: "
+	       << additional_log_potentials.size() << " vs "
+	       << GetNumAdditionals() << std::endl;
+	    throw std::invalid_argument(ss.str());
+	}
         additional_log_potentials_ = additional_log_potentials;
     }
 
@@ -574,6 +581,10 @@ class FactorBUDGET : public Factor
     size_t GetBudget() { return budget_; }
     void SetBudget(size_t budget) { budget_ = budget; }
 
+    // Get/set forced budget (false means can be <=, true means must be =).
+    bool ForcedBudget() { return forced_budget_; }
+    void ForceBudget(bool forced) { forced_budget_ = forced; }
+
     // Add evidence information to the factor.
     int AddEvidence(vector<bool>* active_links,
                     vector<int>* evidence,
@@ -602,6 +613,9 @@ class FactorBUDGET : public Factor
                      vector<double>& out_add);
 
   private:
+    // Forced budget (false means can be <=, true means must be =).
+    bool forced_budget_;
+
     // Budget value.
     size_t budget_;
 

ad3qp/ad3/FactorGraph.cpp

@@ -890,15 +890,16 @@ FactorGraph::CheckAcyclic()
     return acyclic;
 }
 
-/* standalone argmax function for acyclic factor graphs.
- * should avoid the cmplication of the RunAD3 function.
+/* standalone MAP / SparseMAP function for acyclic factor graphs.
+ * should avoid the complication of the RunAD3 function.
  * */
-double
-FactorGraph::MaximizeAcyclic(vector<double>* posteriors,
-                             vector<double>* additionals)
+int
+FactorGraph::SolveAcyclic(bool solve_qp,
+			  vector<double>* posteriors,
+			  vector<double>* additionals,
+			  double* value)
 {
-
-    double value = 0;
+    *value = 0.0;    
     posteriors->resize(variables_.size(), 0.0);
 
     vector<int> add_offsets(factors_.size());
@@ -907,20 +908,27 @@ FactorGraph::MaximizeAcyclic(vector<double>* posteriors,
         add_offsets[j] = offset;
         offset += factors_[j]->GetNumAdditionals();
     }
-
     additionals->resize(offset);
 
     // deal with zero-degree variables if any
     for (size_t i = 0; i < variables_.size(); ++i) {
         auto v = variables_[i];
         if (v->Degree() == 0) {
             auto eta_ui = v->GetLogPotential();
-            (*posteriors)[i] = (eta_ui > 0) ? 1.0 : 0.0;
-            value += (eta_ui > 0) ? eta_ui : 0;
+	    if (solve_qp) {
+                double u = max(0.0, min(eta_ui, 1.0));
+		(*posteriors)[i] = u;
+		*value += u * eta_ui - 0.5 * u * u;
+	    } else {
+	        (*posteriors)[i] = (eta_ui > 0) ? 1.0 : 0.0;
+		*value += (eta_ui > 0) ? eta_ui : 0;
+	    }
         }
     }
 
-    // run maximize for each factor
+    bool fractional = false;
+
+    // run QP for each factor
     for (size_t j = 0; j < factors_.size(); ++j) {
         Factor* f = factors_[j];
         size_t factor_degree = f->Degree();
@@ -932,25 +940,54 @@ FactorGraph::MaximizeAcyclic(vector<double>* posteriors,
             (*eta_uf)[i] = f->GetVariable(i)->GetLogPotential();
         }
 
-        double value_f;
-        f->SolveMAPCached(&value_f);
-        value += value_f;
+	if (solve_qp) {
+	    // Solve the QP.
+	    f->SolveQPCached();
+	} else {
+            // Solve MAP.
+	    double value_f;
+	    f->SolveMAPCached(&value_f);
+	    *value += value_f;
+	}
 
         const vector<double>& mu_uf = f->GetCachedVariablePosteriors();
-
         for (size_t i = 0; i < factor_degree; ++i) {
             auto k = f->GetVariable(i)->GetId();
-            (*posteriors)[k] = mu_uf[i];
+	    double logp = f->GetVariable(i)->GetLogPotential();
+	    if (solve_qp) {
+	        double u_i = mu_uf[i];
+		(*posteriors)[k] = u_i;
+		if (!NEARLY_BINARY((*posteriors)[i], 1e-12))
+		    fractional = true;
+	        *value += u_i * logp - .5 * u_i * u_i;
+	    } else {
+	        (*posteriors)[k] = mu_uf[i];
+	    }
         }
 
+	auto add_logp = f->GetAdditionalLogPotentials();
         const auto& mu_vf = f->GetCachedAdditionalPosteriors();
         offset = add_offsets[j];
         for (size_t i = 0; i < mu_vf.size(); ++i) {
-            (*additionals)[offset] = mu_vf[i];
+	    (*additionals)[offset] = mu_vf[i];
             ++offset;
+	    if (solve_qp) {
+	        *value += mu_vf[i] * add_logp[i];
+	    }
         }
     }
-    return value;
+
+    if (!fractional) {
+      if (verbosity_ > 1) {
+	std::cout << "Solution is integer." << std::endl;
+      }
+      return STATUS_OPTIMAL_INTEGER;
+    } else {
+      if (verbosity_ > 1) {
+	std::cout << "Solution is fractional." << std::endl;
+      }
+      return STATUS_OPTIMAL_FRACTIONAL;
+    }
 }
 
 int
@@ -985,13 +1022,12 @@ FactorGraph::RunAD3(double lower_bound,
     }
 
     // for acyclic graph, can avoid iterative algorithm.
-    if (!solve_qp && autodetect_acyclic_ && CheckAcyclic()) {
+    if (autodetect_acyclic_ && CheckAcyclic()) {
         if (verbosity_ > 1) {
             std::cout << "Factor graph is acyclic; using this. \n"
                       << std::endl;
         }
-        *value = MaximizeAcyclic(posteriors, additional_posteriors);
-        return STATUS_OPTIMAL_INTEGER;
+        return SolveAcyclic(solve_qp, posteriors, additional_posteriors, value);
     }
 
     // Stopping criterion parameters.

ad3qp/ad3/FactorGraph.h

@@ -272,20 +272,24 @@ class FactorGraph
     // Create a new BUDGET factor.
     Factor* CreateFactorBUDGET(const vector<BinaryVariable*>& variables,
                                int budget,
+			       bool force_budget = false,
                                bool owned_by_graph = true)
     {
         vector<bool> negated;
-        return CreateFactorBUDGET(variables, negated, budget, owned_by_graph);
+        return CreateFactorBUDGET(variables, negated, budget, force_budget,
+				  owned_by_graph);
     }
 
     Factor* CreateFactorBUDGET(const vector<BinaryVariable*>& variables,
                                const vector<bool>& negated,
                                int budget,
+			       bool force_budget = false,
                                bool owned_by_graph = true)
     {
         Factor* factor = new FactorBUDGET;
         DeclareFactor(factor, variables, negated, owned_by_graph);
         static_cast<FactorBUDGET*>(factor)->SetBudget(budget);
+        static_cast<FactorBUDGET*>(factor)->ForceBudget(force_budget);
         return factor;
     }
 
@@ -529,8 +533,10 @@ class FactorGraph
 
     bool CheckAcyclic();
 
-    double MaximizeAcyclic(vector<double>* posteriors,
-                           vector<double>* additionals);
+    int SolveAcyclic(bool solve_qp,
+		     vector<double>* posteriors,
+		     vector<double>* additionals,
+		     double* value);
 
 
   private:

lpsmap/__init__.py

@@ -1,7 +1,7 @@
 from .version import __version__
 from .api.factors import (Xor, Or, AtMostOne, Imply, XorOut, OrOut,
                          AndOut, Budget, Pair)
-from .api.factors_extension import Sequence, DepTree
+from .api.factors_extension import Sequence, DepTree, SequenceBudget
 from .api.api import FactorGraph
 
 try: