Minimum Feedback Arc Set (MFAS) for 1dsfm

gtsam/sfm/mfas.cpp

@@ -0,0 +1,101 @@
+/*
+  This file defines functions used to solve a Minimum feedback arc set (MFAS)
+  problem. This code was forked and modified from Kyle Wilson's repository at
+  https://github.com/wilsonkl/SfM_Init.
+  Copyright (c) 2014, Kyle Wilson
+  All rights reserved.
+*/
+
+#include <gtsam/sfm/mfas.h>
+
+#include <map>
+#include <set>
+#include <vector>
+
+using std::map;
+using std::pair;
+using std::set;
+using std::vector;
+
+namespace gtsam {
+namespace mfas {
+
+void flipNegEdges(vector<KeyPair> &edges, vector<double> &weights) {
+  // now renumber the edges
+  for (unsigned int i = 0; i < edges.size(); ++i) {
+    if (weights[i] < 0.0) {
+      Key tmp = edges[i].second;
+      edges[i].second = edges[i].first;
+      edges[i].first = tmp;
+      weights[i] *= -1;
+    }
+  }
+}
+
+void mfasRatio(const std::vector<KeyPair> &edges,
+               const std::vector<double> &weights, const KeyVector &nodes,
+               FastMap<Key, int> &ordered_positions) {
+  // initialize data structures
+  FastMap<Key, double> win_deg;
+  FastMap<Key, double> wout_deg;
+  FastMap<Key, vector<pair<int, double> > > inbrs;
+  FastMap<Key, vector<pair<int, double> > > onbrs;
+
+  // stuff data structures
+  for (unsigned int ii = 0; ii < edges.size(); ++ii) {
+    int i = edges[ii].first;
+    int j = edges[ii].second;
+    double w = weights[ii];
+
+    win_deg[j] += w;
+    wout_deg[i] += w;
+    inbrs[j].push_back(pair<int, double>(i, w));
+    onbrs[i].push_back(pair<int, double>(j, w));
+  }
+
+  unsigned int ordered_count = 0;
+  while (ordered_count < nodes.size()) {
+    // choose an unchosen node
+    Key choice;
+    double max_score = 0.0;
+    for (auto node : nodes) {
+      if (ordered_positions.find(node) == ordered_positions.end()) {
+        // is this a source
+        if (win_deg[node] < 1e-8) {
+          choice = node;
+          break;
+        } else {
+          double score = (wout_deg[node] + 1) / (win_deg[node] + 1);
+          if (score > max_score) {
+            max_score = score;
+            choice = node;
+          }
+        }
+      }
+    }
+    // find its inbrs, adjust their wout_deg
+    for (auto it = inbrs[choice].begin(); it != inbrs[choice].end(); ++it)
+      wout_deg[it->first] -= it->second;
+    // find its onbrs, adjust their win_deg
+    for (auto it = onbrs[choice].begin(); it != onbrs[choice].end(); ++it)
+      win_deg[it->first] -= it->second;
+
+    ordered_positions[choice] = ordered_count++;
+  }
+}
+
+void outlierWeights(const std::vector<KeyPair> &edges,
+                    const std::vector<double> &weight,
+                    const FastMap<Key, int> &ordered_positions,
+                    FastMap<KeyPair, double> &outlier_weights) {
+  // find the outlier edges
+  for (unsigned int i = 0; i < edges.size(); ++i) {
+    int x0 = ordered_positions.at(edges[i].first);
+    int x1 = ordered_positions.at(edges[i].second);
+    if ((x1 - x0) * weight[i] < 0)
+      outlier_weights[edges[i]] += weight[i] > 0 ? weight[i] : -weight[i];
+  }
+}
+
+}  // namespace mfas
+}  // namespace gtsam

gtsam/sfm/mfas.h

@@ -0,0 +1,65 @@
+/*
+  This file defines functions used to solve a Minimum feedback arc set (MFAS)
+  problem. This code was forked and modified from Kyle Wilson's repository at
+  https://github.com/wilsonkl/SfM_Init.
+  Copyright (c) 2014, Kyle Wilson
+  All rights reserved.
+
+  Given a weighted directed graph, the objective in a Minimum feedback arc set
+  problem is to obtain a graph that does not contain any cycles by removing
+  edges such that the total weight of removed edges is minimum.
+*/
+#ifndef __MFAS_H__
+#define __MFAS_H__
+
+#include <gtsam/inference/Key.h>
+
+#include <map>
+#include <vector>
+
+namespace gtsam {
+
+using KeyPair = std::pair<Key, Key>;
+
+namespace mfas {
+
+/*
+ * Given a vector of KeyPairs that constitutes edges in a graph and the weights
+ * corresponding to these edges, this function changes all the weights to
+ * positive numbers by flipping the direction of the edges that have a negative
+ * weight. The changes are made in place.
+ * @param edges reference to vector of KeyPairs
+ * @param weights weights corresponding to edges
+ */
+void flipNegEdges(std::vector<KeyPair> &edges, std::vector<double> &weights);
+
+/*
+ * Computes the MFAS ordering, ie an ordering of the nodes in the graph such
+ * that the source of any edge appears before its destination in the ordering.
+ * The weight of edges that are removed to obtain this ordering is minimized.
+ * @param edges: edges in the graph
+ * @param weights: weights corresponding to the edges (have to be positive)
+ * @param nodes: nodes in the graph
+ * @param ordered_positions: map from node to position in the ordering (0
+ * indexed)
+ */
+void mfasRatio(const std::vector<KeyPair> &edges,
+               const std::vector<double> &weights, const KeyVector &nodes,
+               FastMap<Key, int> &ordered_positions);
+
+/*
+ * Returns the weights of edges that are not consistent with the input ordering.
+ * @param edges in the graph
+ * @param weights of the edges in the graph
+ * @param ordered_positions: ordering (obtained from MFAS solution)
+ * @param outlier_weights: reference to a map from edges to their "outlier
+ * weights"
+ */
+void outlierWeights(const std::vector<KeyPair> &edges,
+                    const std::vector<double> &weight,
+                    const FastMap<Key, int> &ordered_positions,
+                    FastMap<KeyPair, double> &outlier_weights);
+
+}  // namespace mfas
+}  // namespace gtsam
+#endif  // __MFAS_H__

gtsam/sfm/tests/testMFAS.cpp

@@ -0,0 +1,127 @@
+#include <gtsam/sfm/mfas.h>
+
+#include <CppUnitLite/TestHarness.h>
+
+using namespace std;
+using namespace gtsam;
+
+/* We (partially) use the example from the paper on 1dsfm
+ * (https://research.cs.cornell.edu/1dsfm/docs/1DSfM_ECCV14.pdf, Fig 1, Page 5)
+ * for the unit tests here. The only change is that we leave out node 4 and use
+ * only nodes 0-3. This not only makes the test easier to understand but also
+ * avoids an ambiguity in the ground truth ordering that arises due to
+ * insufficient edges in the geaph. */
+
+// edges in the graph - last edge from node 3 to 0 is an outlier
+vector<KeyPair> graph = {make_pair(3, 2), make_pair(0, 1), make_pair(3, 1),
+                         make_pair(1, 2), make_pair(0, 2), make_pair(3, 0)};
+// nodes in the graph
+KeyVector nodes = {Key(0), Key(1), Key(2), Key(3)};
+// weights from projecting in direction-1 (bad direction, outlier accepted)
+vector<double> weights1 = {2, 1.5, 0.5, 0.25, 1, 0.75};
+// weights from projecting in direction-2 (good direction, outlier rejected)
+vector<double> weights2 = {0.5, 0.75, -0.25, 0.75, 1, 0.5};
+
+// Testing the flipNegEdges function for weights1
+TEST(MFAS, FlipNegEdges) {
+  vector<KeyPair> graph_copy = graph;
+  vector<double> weights1_positive = weights1;
+  mfas::flipNegEdges(graph_copy, weights1_positive);
+
+  // resulting graph and edges must be of same size
+  EXPECT_LONGS_EQUAL(graph_copy.size(), graph.size());
+  EXPECT_LONGS_EQUAL(weights1_positive.size(), weights1.size());
+
+  for (unsigned int i = 0; i < weights1.size(); i++) {
+    if (weights1[i] < 0) {
+      // if original weight was negative, edges must be flipped and new weight
+      // must be positive
+      EXPECT_DOUBLES_EQUAL(weights1_positive[i], -weights1[i], 1e-6);
+      EXPECT(graph_copy[i].first == graph[i].second &&
+             graph_copy[i].second == graph[i].first);
+    } else {
+      // unchanged if original weight was positive
+      EXPECT_DOUBLES_EQUAL(weights1_positive[i], weights1[i], 1e-6);
+      EXPECT(graph_copy[i].first == graph[i].first &&
+             graph_copy[i].second == graph[i].second);
+    }
+  }
+}
+
+// test the ordering and the outlierWeights function using weights2 - outlier
+// edge is rejected when projected in a direction that gives weights2
+TEST(MFAS, OrderingWeights2) {
+  vector<KeyPair> graph_copy = graph;
+  vector<double> weights2_positive = weights2;
+  mfas::flipNegEdges(graph_copy, weights2_positive);
+  FastMap<Key, int> ordered_positions;
+  // compute ordering from positive edge weights
+  mfas::mfasRatio(graph_copy, weights2_positive, nodes, ordered_positions);
+
+  // expected ordering in this example
+  FastMap<Key, int> gt_ordered_positions;
+  gt_ordered_positions[0] = 0;
+  gt_ordered_positions[1] = 1;
+  gt_ordered_positions[3] = 2;
+  gt_ordered_positions[2] = 3;
+
+  // check if the expected ordering is obtained
+  for (auto node : nodes) {
+    EXPECT_LONGS_EQUAL(gt_ordered_positions[node], ordered_positions[node]);
+  }
+
+  // testing the outlierWeights method
+  FastMap<KeyPair, double> outlier_weights;
+  mfas::outlierWeights(graph_copy, weights2_positive, gt_ordered_positions,
+                       outlier_weights);
+  // since edge between 3 and 0 is inconsistent with the ordering, it must have
+  // positive outlier weight, other outlier weights must be zero
+  for (auto &edge : graph_copy) {
+    if (edge == make_pair(Key(3), Key(0)) ||
+        edge == make_pair(Key(0), Key(3))) {
+      EXPECT_DOUBLES_EQUAL(outlier_weights[edge], 0.5, 1e-6);
+    } else {
+      EXPECT_DOUBLES_EQUAL(outlier_weights[edge], 0, 1e-6);
+    }
+  }
+}
+
+// test the ordering function and the outlierWeights method using
+// weights2 (outlier edge is accepted when projected in a direction that
+// produces weights2)
+TEST(MFAS, OrderingWeights1) {
+  vector<KeyPair> graph_copy = graph;
+  vector<double> weights1_positive = weights1;
+  mfas::flipNegEdges(graph_copy, weights1_positive);
+  FastMap<Key, int> ordered_positions;
+  // compute ordering from positive edge weights
+  mfas::mfasRatio(graph_copy, weights1_positive, nodes, ordered_positions);
+
+  // expected "ground truth" ordering in this example
+  FastMap<Key, int> gt_ordered_positions;
+  gt_ordered_positions[3] = 0;
+  gt_ordered_positions[0] = 1;
+  gt_ordered_positions[1] = 2;
+  gt_ordered_positions[2] = 3;
+
+  // check if the expected ordering is obtained
+  for (auto node : nodes) {
+    EXPECT_LONGS_EQUAL(gt_ordered_positions[node], ordered_positions[node]);
+  }
+
+  // since all edges (including the outlier) are consistent with this ordering,
+  // all outlier_weights must be zero
+  FastMap<KeyPair, double> outlier_weights;
+  mfas::outlierWeights(graph, weights1_positive, gt_ordered_positions,
+                       outlier_weights);
+  for (auto &edge : graph) {
+    EXPECT_DOUBLES_EQUAL(outlier_weights[edge], 0, 1e-6);
+  }
+}
+
+/* ************************************************************************* */
+int main() {
+  TestResult tr;
+  return TestRegistry::runAllTests(tr);
+}
+/* ************************************************************************* */