diff --git a/docs/source/advancedcommands.rst b/docs/source/advancedcommands.rst index 5afa5e8f3..a01d3b529 100644 --- a/docs/source/advancedcommands.rst +++ b/docs/source/advancedcommands.rst @@ -1,6 +1,65 @@ Advanced Commands ***************** +.. _migp: + +Choice Variables +================ +If MOSEK 9 is installed, GPkit supports discretized free variables with the ``mosek_conif`` solver. +Choice variables are free in the sense of having multiple possible choices, but discretized in the +sense of having a limited set of possible solutions. + +.. literalinclude:: examples/migp.py + +If solved with the mosek_conif solver, the script above will print:: + + Optimal Cost + ------------ + [ 1.5 2.15 2.8 3.22 ... ] + + ~~~~~~~~ + WARNINGS + ~~~~~~~~ + No Dual Solution + ---------------- + This model has the discretized choice variables [x] 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"]`. + ~~~~~~~~ + + 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 : [ 1 + 1 + 1 + 2 + 2 + 2 + 2 + 2 + 2 + 3 + 3 ] + +Note that the optimal values for ``x`` are discretized, clicking from 1 +to 2 to 3 during the sweep, and that the solution has no dual variables. + + Derived Variables ================= diff --git a/docs/source/examples/migp.py b/docs/source/examples/migp.py new file mode 100644 index 000000000..372c7ab1f --- /dev/null +++ b/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()) diff --git a/docs/source/examples/migp_output.txt b/docs/source/examples/migp_output.txt new file mode 100644 index 000000000..8dfe175c2 --- /dev/null +++ b/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) + diff --git a/docs/source/installation.rst b/docs/source/installation.rst index 96dc5d84a..15fb59528 100644 --- a/docs/source/installation.rst +++ b/docs/source/installation.rst @@ -4,42 +4,44 @@ Installation ************ 1. If you are on Mac or Windows, we recommend installing `Anaconda `_. Alternatively, `install pip and create a virtual environment `_. -2. (optional) Install the MOSEK solver as directed below -3. Run ``pip install gpkit`` in the appropriate terminal or command prompt. -4. Open a Python prompt and run ``import gpkit`` to finish installation and run unit tests. +2. (optional) Install the MOSEK 9 solver with ``pip install Mosek``, then a license as described below +3. (optional) Install the MOSEK 8 solver as described below +4. Run ``pip install gpkit`` in the appropriate terminal or command prompt. +5. Open a Python prompt and run ``import gpkit`` to finish installation and run unit tests. If you encounter any bugs please email ``gpkit@mit.edu`` or `raise a GitHub issue `_. -Installing MOSEK -================ -GPkit interfaces with two off the shelf solvers: cvxopt, and MOSEK. +Installing MOSEK 8 +================== +GPkit interfaces with two off the shelf solvers: cvxopt, and MOSEK (versions 8 and 9). Cvxopt is open source and installed by default; MOSEK requires a commercial licence or (free) -academic license. +academic license. In MOSEK version 8 GPkit uses the command-line interface ``mskexpopt`` solver, while +in MOSEK 9 it uses the more active exponential-cone interface (and hence supports :ref:`migp`). Mac OS X - If ``which gcc`` does not return anything, install the `Apple Command Line Tools `_. - Download `MOSEK 8 `_, then: - Move the ``mosek`` folder to your home directory - Follow `these steps for Mac `_. - - Request an `academic license file `_ and put it in ``~/mosek/`` + - Request an `academic license file `_ and put it in ``~/mosek/`` Linux - Download `MOSEK 8 `_, then: - Move the ``mosek`` folder to your home directory - Follow `these steps for Linux `_. - - Request an `academic license file `_ and put it in ``~/mosek/`` + - Request an `academic license file `_ and put it in ``~/mosek/`` Windows + - Make sure ``gcc`` is on your system path. + - To do this, type ``gcc`` into a command prompt. + - If you get ``executable not found``, then install the 64-bit version (x86_64 installer architecture dropdown option) with GCC version 6.4.0 or older of `mingw `_. + - In an Anaconda command prompt (or equivalent), run ``cd C:\Program Files\mingw-w64\x86_64-6.4.0-posix-seh-rt_v5-rev0\`` (or whatever corresponds to the correct installation directory; note that if mingw is in ``Program Files (x86)`` instead of ``Program Files`` you've installed the 32-bit version by mistake) + - Run ``mingw-w64`` to add it to your executable path. For step 3 of the install process you'll need to run ``pip install gpkit`` from this prompt. - Download `MOSEK 8 `_, then: - Follow `these steps for Windows `_. - - Request an `academic license file `_ and put it in ``C:\Users\(your_username)\mosek\`` - - Make sure ``gcc`` is on your system path. - - To do this, type ``gcc`` into a command prompt. - - If you get ``executable not found``, then install the 64-bit version (x86_64 installer architecture dropdown option) with GCC version 6.4.0 or older of `mingw `_. - - In an Anaconda command prompt (or equivalent), run ``cd C:\Program Files\mingw-w64\x86_64-6.4.0-posix-seh-rt_v5-rev0\`` (or whatever corresponds to the correct installation directory; note that if mingw is in ``Program Files (x86)`` instead of ``Program Files`` you've installed the 32-bit version by mistake) - - Run ``mingw-w64`` to add it to your executable path. For step 3 of the install process you'll need to run ``pip install gpkit`` from this prompt. + - Request an `academic license file `_ and put it in ``C:\Users\(your_username)\mosek\`` Debugging your installation =========================== diff --git a/gpkit/constraints/gp.py b/gpkit/constraints/gp.py index 89a8fe294..b236596e7 100644 --- a/gpkit/constraints/gp.py +++ b/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): diff --git a/gpkit/small_classes.py b/gpkit/small_classes.py index b66d71464..81d5043fc 100644 --- a/gpkit/small_classes.py +++ b/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) diff --git a/gpkit/solution_array.py b/gpkit/solution_array.py index 7af996894..30fe3eda9 100644 --- a/gpkit/solution_array.py +++ b/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 diff --git a/gpkit/solvers/mosek_conif.py b/gpkit/solvers/mosek_conif.py index 75c4b123c..450acaef1 100644 --- a/gpkit/solvers/mosek_conif.py +++ b/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 diff --git a/gpkit/tests/helpers.py b/gpkit/tests/helpers.py index 206575032..dadad448b 100644 --- a/gpkit/tests/helpers.py +++ b/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 diff --git a/gpkit/tests/t_examples.py b/gpkit/tests/t_examples.py index 41db5f576..c84bb9cb6 100644 --- a/gpkit/tests/t_examples.py +++ b/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)