Merge ae241c0 into 40a9466

docs/source/examples/migp.py

@@ -0,0 +1,10 @@
+import numpy as np
+from gpkit import *
+
+x = Variable("x", choices=range(1,4))
+num = Variable("numerator", np.linspace(0.5, 7, 11))
+
+m = Model(x + num/x)
+sol = m.solve(verbosity=0)
+
+print(sol.table())

docs/source/examples/migp_output.txt

@@ -0,0 +1,59 @@
+
+Optimal Cost
+------------
+ [ 1.41      2.14      2.68      3.13     ... ]
+
+~~~~~~~~
+WARNINGS
+~~~~~~~~
+Freed Choice Variables
+----------------------
+This model has the discretized choice variables [x], but since the 'cvxopt' solver doesn't support discretization they were treated as continuous variables.
+~~~~~~~~
+
+Swept Variables
+---------------
+numerator : [ 0.5
+              1.15
+              1.8
+              2.45
+              3.1
+              3.75
+              4.4
+              5.05
+              5.7
+              6.35
+              7         ]
+
+Free Variables
+--------------
+x : [ 0.707
+      1.07
+      1.34
+      1.57
+      1.76
+      1.94
+      2.1
+      2.25
+      2.39
+      2.52
+      2.65      ]
+
+Variable Sensitivities
+----------------------
+numerator : [ +0.5
+              +0.5
+              +0.5
+              +0.5
+              +0.5
+              +0.5
+              +0.5
+              +0.5
+              +0.5
+              +0.5
+              +0.5      ]
+
+Most Sensitive Constraints (in last sweep)
+------------------------------------------
+(none)
+

gpkit/constraints/gp.py

@@ -69,6 +69,7 @@ class GeometricProgram:
     >>> gp.solve()
     """
     _result = solve_log = solver_out = model = v_ss = nu_by_posy = None
+    choicevaridxs = integersolve = None
 
     def __init__(self, cost, constraints, substitutions, *, checkbounds=True):
         self.cost, self.substitutions = cost, substitutions
@@ -155,6 +156,8 @@ def gen(self):
                 m_idx += 1
         self.p_idxs = np.array(self.p_idxs, "int32")  # to use array equalities
         self.varidxs = {vk: i for i, vk in enumerate(self.varlocs)}
+        self.choicevaridxs = {vk: i for i, vk in enumerate(self.varlocs)
+                              if vk.choices}
         for j, (var, locs) in enumerate(self.varlocs.items()):
             row.extend(locs)
             col.extend([j]*len(locs))
@@ -191,6 +194,9 @@ def solve(self, solver=None, *, verbosity=1, gen_result=True, **kwargs):
 
         solverargs = DEFAULT_SOLVER_KWARGS.get(solvername, {})
         solverargs.update(kwargs)
+        if self.choicevaridxs and solvername == "mosek_conif":
+            solverargs["choicevaridxs"] = self.choicevaridxs
+            self.integersolve = True
         starttime = time()
         solver_out, infeasibility, original_stdout = {}, None, sys.stdout
         try:
@@ -304,6 +310,27 @@ def _compile_result(self, solver_out):
         result["constants"] = KeyDict(self.substitutions)
         result["variables"] = KeyDict(result["freevariables"])
         result["variables"].update(result["constants"])
+        result["soltime"] = solver_out["soltime"]
+        if self.integersolve:
+            result["choicevariables"] = KeyDict( \
+                {k: v for k, v in result["freevariables"].items()
+                 if k in self.choicevaridxs})
+            result["warnings"] = {"No Dual Solution": [(\
+                "This model has the discretized choice variables"
+                " %s and hence no dual solution. You can fix those variables"
+                " to their optimal value and get sensitivities to the resulting"
+                " continuous problem by updating your model's substitions with"
+                " `sol[\"choicevariables\"]`."
+                % sorted(self.choicevaridxs.keys()), self.choicevaridxs)]}
+            return SolutionArray(result)
+        elif self.choicevaridxs:
+            result["warnings"] = {"Freed Choice Variables": [(\
+                "This model has the discretized choice variables"
+                " %s, but since the '%s' solver doesn't support discretization"
+                " they were treated as continuous variables."
+                % (sorted(self.choicevaridxs.keys()), solver_out["solver"]),
+                self.choicevaridxs)]}
+
         result["sensitivities"] = {"constraints": {}}
         la, self.nu_by_posy = self._generate_nula(solver_out)
         cost_senss = sum(nu_i*exp for (nu_i, exp) in zip(self.nu_by_posy[0],
@@ -342,7 +369,6 @@ def _compile_result(self, solver_out):
         result["sensitivities"]["variables"] = KeyDict(gpv_ss)
         result["sensitivities"]["constants"] = \
             result["sensitivities"]["variables"]  # NOTE: backwards compat.
-        result["soltime"] = solver_out["soltime"]
         return SolutionArray(result)
 
     def check_solution(self, cost, primal, nu, la, tol, abstol=1e-20):
@@ -379,6 +405,8 @@ def almost_equal(num1, num2):
                 raise Infeasible("Primal solution violates constraint: %s is "
                                  "greater than 1" % primal_exp_vals[mi].sum())
         # check dual sol #
+        if self.integersolve:
+            return
         # note: follows dual formulation in section 3.1 of
         # http://web.mit.edu/~whoburg/www/papers/hoburg_phd_thesis.pdf
         if not almost_equal(self.nu_by_posy[0].sum(), 1):

gpkit/small_classes.py

@@ -74,9 +74,10 @@ def __init__(self, output=None, *, verbosity=0):
 
     def write(self, writ):
         "Append and potentially write the new line."
+        if writ[:2] == "b'":
+            writ = writ[2:-1]
         if writ != "\n":
-            writ = writ.rstrip("\n")
-            self.append(str(writ))
+            self.append(writ.rstrip("\n"))
         if self.verbosity > 0:  # pragma: no cover
             self.output.write(writ)
 

gpkit/solution_array.py

@@ -243,7 +243,9 @@ def warnings_table(self, _, **kwargs):
         if len(data_vec) == 0:
             continue
         if not hasattr(data_vec, "shape"):
-            data_vec = [data_vec]
+            data_vec = [data_vec]  # not a sweep
+        if all((data == data_vec[0]).all() for data in data_vec[1:]):
+            data_vec = [data_vec[0]]  # warnings identical across all sweeps
         for i, data in enumerate(data_vec):
             if len(data) == 0:
                 continue
@@ -340,7 +342,8 @@ class SolutionArray(DictOfLists):
     """
     modelstr = ""
     _name_collision_varkeys = None
-    table_titles = {"sweepvariables": "Swept Variables",
+    table_titles = {"choicevariables": "Choice Variables",
+                    "sweepvariables": "Swept Variables",
                     "freevariables": "Free Variables",
                     "constants": "Fixed Variables",  # TODO: change everywhere
                     "variables": "Variables"}
@@ -683,7 +686,7 @@ def table(self, showvars=(),
         -------
         str
         """
-        if sortmodelsbysenss:
+        if sortmodelsbysenss and "sensitivities" in self:
             kwargs["sortmodelsbysenss"] = self["sensitivities"]["models"]
         else:
             kwargs["sortmodelsbysenss"] = False
@@ -700,7 +703,10 @@ def table(self, showvars=(),
         showvars = self._parse_showvars(showvars)
         strs = []
         for table in tables:
-            if table == "cost":
+            if "sensitivities" not in self and ("sensitivities" in table or
+                                                "constraints" in table):
+                continue
+            elif table == "cost":
                 cost = self["cost"]  # pylint: disable=unsubscriptable-object
                 if kwargs.get("latex", None):  # cost is not printed for latex
                     continue

gpkit/solvers/mosek_conif.py

@@ -180,6 +180,35 @@ def optimize(*, c, A, k, p_idxs, **kwargs):
         task.putconboundlist(con_indices, type_constraint, h, h)
         cur_con_idx += log_c_lin.size
     #
+    #   Constrain choice variables to discrete choices
+    #
+    choicevaridxs = kwargs.get("choicevaridxs", {})
+    if choicevaridxs:
+        n_choicevars = 0
+        for var, idx in choicevaridxs.items():
+            choices = sorted(var.choices)
+            m_choices = len(choices) - 1  # the first option is the default
+            choiceidxs = np.arange(msk_nvars + n_choicevars,
+                                   msk_nvars + n_choicevars + m_choices)
+            n_choicevars += m_choices
+            task.appendvars(m_choices)
+            task.putvartypelist(choiceidxs,
+                                [mosek.variabletype.type_int]*m_choices)
+            task.putvarboundlist(choiceidxs, [mosek.boundkey.ra]*m_choices,
+                                 np.zeros(m_choices), np.ones(m_choices))
+            task.appendcons(m_choices)
+            for i in range(m_choices - 1):
+                # each larger choice requires those before it
+                task.putarow(cur_con_idx + i, choiceidxs[i:i+2], [1.0, -1.0])
+                task.putconbound(cur_con_idx + i, mosek.boundkey.lo, 0.0, 0.0)
+                cur_con_idx += 1
+            base = np.log(choices[0])
+            logdiffs = np.diff(np.log(choices)).tolist()
+            task.putarow(cur_con_idx, choiceidxs.tolist() + [idx],
+                         logdiffs + [-1])  # choices are summed
+            task.putconbound(cur_con_idx, mosek.boundkey.fx, -base, -base)
+            cur_con_idx += 1
+    #
     #   Set the objective function
     #
     task.putclist([int(m + 3*n_lse)], [1])
@@ -212,40 +241,54 @@ def streamprinter(text):
     #
     #   Recover the solution
     #
-    msk_solsta = task.getsolsta(mosek.soltype.itr)
+    if choicevaridxs:
+        sol = mosek.soltype.itg
+        optimal = mosek.solsta.integer_optimal
+    else:
+        sol = mosek.soltype.itr
+        optimal = mosek.solsta.optimal
+    msk_solsta = task.getsolsta(sol)
     if msk_solsta == mosek.solsta.prim_infeas_cer:
         raise PrimalInfeasible()
     if msk_solsta == mosek.solsta.dual_infeas_cer:
         raise DualInfeasible()
-    if msk_solsta != mosek.solsta.optimal:  # pragma: no cover
+    if msk_solsta != optimal:  # pragma: no cover
         raise UnknownInfeasible("solution status: ", msk_solsta)
 
     # recover primal variables
     x = [0.] * m
-    task.getxxslice(mosek.soltype.itr, 0, m, x)
+    task.getxxslice(sol, 0, m, x)
+    # recover binary variables
+    # xbin = [0.] * (n_choicevars)
+    # task.getxxslice(sol, msk_nvars, msk_nvars + n_choicevars, xbin)
+    # wrap things up in a dictionary
+    solution = {"status": "optimal", "primal": np.array(x),
+                "objective": np.exp(task.getprimalobj(sol))}
     # recover dual variables for log-sum-exp epigraph constraints
     # (skip epigraph of the objective function).
-    z_duals = [0.] * (p_lse - 1)
-    task.getsuxslice(mosek.soltype.itr, m + 3*n_lse + 1, msk_nvars, z_duals)
-    z_duals = np.array(z_duals)
-    z_duals[z_duals < 0] = 0
-    # recover dual variables for the remaining user-provided constraints
-    if log_c_lin is not None:
-        aff_duals = [0.] * log_c_lin.size
-        task.getsucslice(mosek.soltype.itr, n_lse + p_lse, cur_con_idx,
-                         aff_duals)
-        aff_duals = np.array(aff_duals)
-        aff_duals[aff_duals < 0] = 0
-        # merge z_duals with aff_duals
-        merged_duals = np.zeros(len(k))
-        merged_duals[lse_posys[1:]] = z_duals  # skipping the cost
-        merged_duals[lin_posys] = aff_duals
-        merged_duals = merged_duals[1:]
+    if choicevaridxs:  # no dual solution
+        solution["la"] = []
+        solution["nu"] = []
     else:
-        merged_duals = z_duals
-    # wrap things up in a dictionary
-    solution = {"status": "optimal", "primal": np.array(x), "la": merged_duals,
-                "objective": np.exp(task.getprimalobj(mosek.soltype.itr))}
+        z_duals = [0.] * (p_lse - 1)
+        task.getsuxslice(mosek.soltype.itr, m + 3*n_lse + 1, msk_nvars, z_duals)
+        z_duals = np.array(z_duals)
+        z_duals[z_duals < 0] = 0
+        # recover dual variables for the remaining user-provided constraints
+        if log_c_lin is None:
+            solution["la"] = z_duals
+        else:
+            aff_duals = [0.] * log_c_lin.size
+            task.getsucslice(mosek.soltype.itr, n_lse + p_lse, cur_con_idx,
+                             aff_duals)
+            aff_duals = np.array(aff_duals)
+            aff_duals[aff_duals < 0] = 0
+            # merge z_duals with aff_duals
+            merged_duals = np.zeros(len(k))
+            merged_duals[lse_posys[1:]] = z_duals  # skipping the cost
+            merged_duals[lin_posys] = aff_duals
+            solution["la"] = merged_duals[1:]
+
     task.__exit__(None, None, None)
     env.__exit__(None, None, None)
     return solution

gpkit/tests/helpers.py

@@ -84,7 +84,10 @@ def test(self):
         filepath = ("".join([path, os.sep, "%s_output.txt" % name])
                     if name not in imported else None)
         with StdoutCaptured(logfilepath=filepath):
-            imported[name] = importlib.import_module(name)
+            if name in imported:
+                importlib.reload(imported[name])
+            else:
+                imported[name] = importlib.import_module(name)
         getattr(self, name)(imported[name])
     return test
 

gpkit/tests/t_examples.py

@@ -92,6 +92,13 @@ def test_evaluated_free_variables(self, example):
     def test_external_constraint(self, example):
         pass
 
+    def test_migp(self, example):
+        if settings["default_solver"] == "mosek_conif":
+            assert_logtol(example.sol(example.x), [1]*3 + [2]*6 + [3]*2)
+        else:
+            assert_logtol(example.sol(example.x),
+                          np.sqrt(example.sol(example.num)))
+
     def test_external_function(self, example):
         external_code = example.external_code
         self.assertEqual(external_code(0), 0)