Added min set cover ILP model

Author: alexandrutomescu

examples/min_set_cover.py

@@ -0,0 +1,25 @@
+import flowpaths as fp
+
+def main():
+    universe = [1,2,3,4,5]
+    subsets = [{1,2,3}, {2,4}, {3,5}, {1,4,5}]
+    subset_weights = [1,2,1,1]
+
+    msc_solver = fp.MinSetCover(
+        universe=universe,
+        subsets=subsets,
+        subset_weights=subset_weights
+        )
+    
+    msc_solver.solve()
+    process_solution(msc_solver)
+
+def process_solution(model: fp.MinGenSet):
+    if model.is_solved():
+        print(model.get_solution(as_subsets=True))
+    else:
+        print("Model could not be solved.")
+    print(model.solve_statistics)
+
+if __name__ == "__main__":
+    main()

flowpaths/__init__.py

@@ -8,6 +8,7 @@
 from .nodeexpandeddigraph import NodeExpandedDiGraph
 from .utils import graphutils as graphutils
 from .mingenset import MinGenSet
+from .minsetcover import MinSetCover
 
 __all__ = [
     "AbstractPathModelDAG",
@@ -20,4 +21,5 @@
     "NodeExpandedDiGraph",
     "graphutils",
     "MinGenSet",
+    "MinSetCover",
 ]

flowpaths/minsetcover.py

@@ -0,0 +1,142 @@
+import flowpaths.utils.solverwrapper as sw
+import time
+
+class MinSetCover():
+    def __init__(
+        self,
+        universe: list,
+        subsets: list,
+        subset_weights: list = None
+        ):
+        """
+        This class solves the minimum set cover problem. Given a universe `universe` and a list of subsets `subsets`,
+        the goal is to find the smallest list of subsets `set_cover` such that: 
+        
+        - every element in `universe` is in at least one subset in `set_cover`.
+        - the sum of the weights of the subsets in `set_cover` is minimized.
+
+        Parameters
+        ----------
+
+        - `universe` : list`
+
+            The universe of elements that must be covered.
+
+        - `subsets : list`
+
+            A list of subsets that can be used to cover the universe.
+
+        - `subset_weights : list`
+
+            The weight of each subset, as a list in the same order that the subsets appear in the list `subsets`. 
+            If not provided, each subset is assumed to have a weight of 1.
+
+        """
+        
+        self.universe = universe
+        self.subsets = subsets
+        self.subset_weights = subset_weights
+        self.set_cover = []
+        self.set_cover_indices = []
+        self.set_cover_weights = []
+
+        self.__is_solved = None
+        self.__solution = None
+
+        self.__encode_set_cover()
+        
+    def __encode_set_cover(self):
+        """
+        This function encodes the set cover problem as an integer linear program.
+        """
+        self.solver = sw.SolverWrapper()
+
+        self.subset_indexes = [(i)   for i in range(len(self.subsets))]
+
+        self.subset_vars = self.solver.add_variables(
+            self.subset_indexes, 
+            name_prefix="subset",
+            lb=0, 
+            ub=1,
+            var_type="integer"
+        )
+
+        # Every element of the universe must be in at least one subset
+        for element in self.universe:            
+            self.solver.add_constraint(
+                self.solver.quicksum(
+                    self.subset_vars[i]
+                    for i in range(len(self.subsets)) if element in self.subsets[i]
+                )
+                >= 1,
+                name=f"total",
+            )
+
+        # Objective function
+        self.solver.set_objective(
+            self.solver.quicksum(
+                self.subset_weights[i] * self.subset_vars[i]
+                for i in range(len(self.subsets))
+            )
+        )
+        
+    def solve(self) -> bool:
+        """
+        Solves the minimum set cover problem. 
+
+        Returns
+        -------
+        - bool
+            
+            `True` if the model was solved, `False` otherwise.
+        """
+        start_time = time.time()
+
+        self.solver.optimize()
+        if self.solver.get_model_status() == "kOptimal":
+            subset_cover_sol = self.solver.get_variable_values("subset", [int])
+            self.__solution = [i for i in range(len(self.subsets)) if subset_cover_sol[i] == 1]
+            self.__is_solved = True
+            self.solve_statistics = {
+                "solve_time": time.time() - start_time,
+                "num_elements": len(self.__solution),
+                "status": self.solver.get_model_status(),
+            }
+            return True
+        else:
+            self.solve_statistics = {
+                "solve_time": time.time() - start_time,
+                "status": self.solver.get_model_status(),
+            }
+            return False
+
+    def is_solved(self):
+        """
+        Returns `True` if the model was solved, `False` otherwise.
+        """
+        if self.__is_solved is None:
+            raise Exception("Model not yet solved. If you want to solve it, call the `solve` method first.")
+        
+        return self.__is_solved
+    
+    def check_is_solved(self):
+        if not self.is_solved():
+            raise Exception(
+                "Model not solved. If you want to solve it, call the solve method first. \
+                  If you already ran the solve method, then the model is infeasible, or you need to increase parameter time_limit."
+            )
+
+    def get_solution(self, as_subsets: bool = False):
+        """
+        Returns the solution to the minimum generating set problem, if the model was solved. 
+        
+        !!! warning "Warning"
+            Call the `solve` method first.
+        """
+        if self.__solution is not None:
+            if not as_subsets:
+                return self.__solution
+            else:
+                return [self.subsets[i] for i in self.__solution]
+        
+        self.check_is_solved()