Merge pull request #3 from nunonmg/forced_budget

Add acyclic computation for SparseMAP

.github/workflows/wheels.yaml

@@ -0,0 +1,47 @@
+name: Build
+
+on: [push, pull_request, release]
+
+jobs:
+  build_wheels:
+    name: Build wheels
+    runs-on: ${{ matrix.os }}
+    strategy:
+      fail-fast: false
+      matrix:
+        os: [ubuntu-20.04, macos-10.15]
+        cibw_python: ["cp38-*", "cp39-*"]
+        cibw_manylinux: ["manylinux2014"]
+
+    steps:
+      - uses: actions/checkout@v2
+
+      - name: Build wheels
+        uses: pypa/cibuildwheel@v2.0.1
+        env:
+           CIBW_BUILD: ${{ matrix.cibw_python }}
+           CIBW_BUILD_VERBOSITY: 3
+           CIBW_BEFORE_BUILD_LINUX: |
+               pip install -U pip
+               yum install eigen3-devel -y
+           CIBW_BEFORE_BUILD_MACOS: |
+               pip install -U pip
+               brew install eigen
+           CIBW_ARCHS_MACOS: x86_64 arm64
+           CIBW_ARCHS_LINUX: auto64 
+           CIBW_MANYLINUX_X86_64_IMAGE: ${{ matrix.cibw_manylinux }}
+           MACOS_DEPLOYMENT_TARGET: 10.14
+           HOMEBREW_NO_INSTALL_CLEANUP: 1
+           CIBW_BEFORE_TEST_LINUX: |
+               pip install pytest
+               pip install torch==1.8.1+cpu -f https://download.pytorch.org/whl/lts/1.8/torch_lts.html
+           CIBW_BEFORE_TEST_MACOS: |
+               pip install pytest
+               pip install torch==1.8.1 -f https://download.pytorch.org/whl/lts/1.8/torch_lts.html
+           CIBW_TEST_COMMAND: |
+               pytest --pyargs lpsmap
+   
+      - uses: actions/upload-artifact@v2
+        with:
+          name: wheels
+          path: ./wheelhouse/*.whl

README.md

@@ -1,3 +1,6 @@
+![wheels](https://github.com/deep-spin/lp-sparsemap/actions/workflows/wheels.yaml/badge.svg)
+[![PyPI version](https://badge.fury.io/py/lp-sparsemap.svg)](https://badge.fury.io/py/lp-sparsemap)
+
 # LP-SparseMAP
 Differentiable sparse structured prediction in coarse factor graphs
 
@@ -28,26 +31,39 @@ Prediction. https://arxiv.org/abs/2001.04437
 
 *Requirements:*
  - Cython
- - [Eigen](https://gitlab.com/libeigen/eigen)
+ - [Eigen](https://gitlab.com/libeigen/eigen) (if it's a non-standard directory,
+   set `EIGEN_DIR=/path/to/eigen`.)
 
 For examples and tests: numpy, pytest.
 
 *Installation:*
 
 ```
-export EIGEN_DIR=/path/to/eigen
-python setup.py build_clib  # builds ad3 in-place
-pip install .               # builds lpsmap and installs
+pip install lp-sparsemap    # installs a wheel, if available.
 ```
 
-*In-place installation:*
+*In-place installation from source:*
 
 ```
+# export MACOS_DEPLOYMENT_TARGET=10.14  # on MacOS
 export EIGEN_DIR=/path/to/eigen
 python setup.py build_clib  # builds ad3 in-place
 pip install -e .            # builds lpsmap and creates a link
 ```
 
+*Using the Cython API from your own code.*
+
+You can add custom factors and other extensions by `cimport`ing the base classes
+provided. (See an example [in this
+project](https://github.com/deep-spin/sparse-marginalization-lvm/).)  The
+installed `lp-sparsemap` package provides a copy of `libad3` to statically link
+against. To get the path to it, use `lpsmap.config.get_libdir()`.  *Warning:*
+both `lp-sparsemap` as well as client libraries linking against it should be
+compiled with the same standard library implementation. On MacOS you may have
+issues unless `MACOS_DEPLOYMENT_TARGET >= 10.14`. If you get undefined symbol
+errors for AD3 symbols, try compiling your code with the same toolchain as the
+installed `lp-sparsemap`. (If in doubt, recompile both locally.)
+
 
 ## dysparsemap
 

ad3qp/ad3/Factor.h

@@ -190,7 +190,10 @@ class Factor
     void SetAdditionalLogPotentials(
       const vector<double>& additional_log_potentials)
     {
-        assert(additional_log_potentials.size() == GetNumAdditionals());
+        if (additional_log_potentials.size() != GetNumAdditionals()) {
+	    throw std::invalid_argument(
+	      "Invalid size of additional log potentials.");
+	}
         additional_log_potentials_ = additional_log_potentials;
     }
 

ad3qp/ad3/FactorGraph.cpp

@@ -875,6 +875,121 @@ FactorGraph::RunBranchAndBound(double cumulative_value,
     return status;
 }
 
+bool
+FactorGraph::CheckAcyclic()
+{
+    bool acyclic = true;
+
+    for (size_t i = 0; i < variables_.size(); ++i) {
+        if (variables_[i]->Degree() > 1) {
+            acyclic = false;
+            break;
+        }
+    }
+
+    return acyclic;
+}
+
+/* standalone MAP / SparseMAP function for acyclic factor graphs.
+ * should avoid the complication of the RunAD3 function.
+ * */
+int
+FactorGraph::SolveAcyclic(bool solve_qp,
+			  vector<double>* posteriors,
+			  vector<double>* additionals,
+			  double* value)
+{
+    *value = 0.0;    
+    posteriors->resize(variables_.size(), 0.0);
+
+    vector<int> add_offsets(factors_.size());
+    int offset = 0;
+    for (size_t j = 0; j < factors_.size(); ++j) {
+        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();
+	    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;
+	    }
+        }
+    }
+
+    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();
+        auto eta_uf = f->GetMutableCachedVariableLogPotentials();
+        f->ComputeCachedAdditionalLogPotentials(1.0);
+        eta_uf->resize(factor_degree);
+
+        for (size_t i = 0; i < factor_degree; ++i) {
+            (*eta_uf)[i] = f->GetVariable(i)->GetLogPotential();
+        }
+
+	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();
+	    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];
+            ++offset;
+	    if (solve_qp) {
+	        *value += mu_vf[i] * add_logp[i];
+	    }
+        }
+    }
+
+    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
 FactorGraph::RunAD3(double lower_bound,
                     vector<double>* posteriors,
@@ -895,7 +1010,6 @@ FactorGraph::RunAD3(double lower_bound,
     gettimeofday(&start, NULL);
 
     // initialize for QP factor solving
-
     if (solve_qp) {
         for (auto&& factor : factors_) {
             size_t factor_degree = factor->Degree();
@@ -907,6 +1021,15 @@ FactorGraph::RunAD3(double lower_bound,
         }
     }
 
+    // for acyclic graph, can avoid iterative algorithm.
+    if (autodetect_acyclic_ && CheckAcyclic()) {
+        if (verbosity_ > 1) {
+            std::cout << "Factor graph is acyclic; using this. \n"
+                      << std::endl;
+        }
+        return SolveAcyclic(solve_qp, posteriors, additional_posteriors, value);
+    }
+
     // Stopping criterion parameters.
     double residual_threshold = ad3_residual_threshold_; // 1e-6;
     // double gap_threshold = 1e-6;

ad3qp/ad3/FactorGraph.h

@@ -417,6 +417,7 @@ class FactorGraph
     }
     void SetEtaAD3(double eta) { ad3_eta_ = eta; }
     void AdaptEtaAD3(bool adapt) { ad3_adapt_eta_ = adapt; }
+    void SetAutodetectAcyclic(bool a) { autodetect_acyclic_ = a; }
     void SetResidualThresholdAD3(double threshold)
     {
         ad3_residual_threshold_ = threshold;
@@ -432,8 +433,12 @@ class FactorGraph
                           double* value)
     {
         double upper_bound;
-        return RunAD3(
-          -1e100, posteriors, additional_posteriors, value, &upper_bound);
+        return RunAD3(-1e100,
+                      posteriors,
+                      additional_posteriors,
+                      value,
+                      &upper_bound,
+                      /*solve_qp=*/false);
     }
 
     int SolveQP(vector<double>* posteriors,
@@ -526,6 +531,13 @@ class FactorGraph
                           double* best_lower_bound,
                           double* best_upper_bound);
 
+    bool CheckAcyclic();
+
+    int SolveAcyclic(bool solve_qp,
+		     vector<double>* posteriors,
+		     vector<double>* additionals,
+		     double* value);
+
 
   private:
     vector<BinaryVariable*> variables_;
@@ -546,6 +558,7 @@ class FactorGraph
     bool ad3_adapt_eta_;
     // Threshold for primal/dual residuals.
     double ad3_residual_threshold_;
+    bool autodetect_acyclic_;
 
     // Parameters for PSDD:
     int psdd_max_iterations_; // Maximum number of iterations.

lpsmap/__init__.py

@@ -1,8 +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, SequenceBudget
-
 from .api.api import FactorGraph
 
 try:

lpsmap/ad3ext/FactorSequenceBudget.h

@@ -30,11 +30,11 @@ using namespace std;
 namespace AD3 {
 
 class FactorSequenceBudget : public GenericFactor {
- protected:
-  virtual double GetNodeScore(int position,
-                                int state,
-                                const vector<double>& variable_log_potentials,
-                                const vector<double>& additional_log_potentials)
+  protected:
+    virtual double GetNodeScore(int position,
+				int state,
+				const vector<double>& variable_log_potentials,
+				const vector<double>& additional_log_potentials)
     {
         return variable_log_potentials[offset_states_[position] + state];
     }
@@ -221,8 +221,9 @@ class FactorSequenceBudget : public GenericFactor {
                        additional_posteriors);
       if (i > 0) {
         AddEdgePosterior(i, previous_state, state, weight,
-                       variable_posteriors,
-                       additional_posteriors);}
+                         variable_posteriors,
+                         additional_posteriors);
+      }
       previous_state = state;
     }
     AddEdgePosterior(sequence->size(), previous_state, 0, weight,