From 8d97565c126916b6213ee4ccd300eabd26ed8093 Mon Sep 17 00:00:00 2001 From: rileyjmurray Date: Sun, 19 Jan 2020 10:49:11 -0800 Subject: [PATCH 01/21] gitignore --- .gitignore | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/.gitignore b/.gitignore index ba309decc..466bf4ff7 100644 --- a/.gitignore +++ b/.gitignore @@ -18,6 +18,7 @@ __pycache__/ .Python env/ bin/ +gpkit/env/* build/ develop-eggs/ dist/ @@ -30,6 +31,7 @@ var/ *.egg-info/ .installed.cfg *.egg +gpkit/_mosek/build/* # Installer logs pip-log.txt @@ -42,6 +44,7 @@ htmlcov/ .cache nosetests.xml coverage.xml +gplibrary/* # Translations *.mo @@ -76,3 +79,6 @@ test_reports_nounits/ # OSX .DS_Store + +# Development environment +.idea/ \ No newline at end of file From eacc720f6cf08e3284bb2c278baa816aeca51aaa Mon Sep 17 00:00:00 2001 From: rileyjmurray Date: Sun, 19 Jan 2020 13:16:28 -0800 Subject: [PATCH 02/21] fixed docstring error --- gpkit/_cvxopt.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/gpkit/_cvxopt.py b/gpkit/_cvxopt.py index 941a95026..3754c768b 100644 --- a/gpkit/_cvxopt.py +++ b/gpkit/_cvxopt.py @@ -18,9 +18,9 @@ def cvxoptimize(c, A, k, *args, **kwargs): --------- c : floats array of shape n Coefficients of each monomial - A : floats array of shape (m,n) + A : floats array of shape (n, m) Exponents of the various free variables for each monomial. - k : ints array of shape n + k : ints array of shape p+1 number of monomials (columns of F) present in each constraint Returns From 739a33c92df49bffeeeb7211f0f56214f925b2cb Mon Sep 17 00:00:00 2001 From: rileyjmurray Date: Sun, 19 Jan 2020 13:33:10 -0800 Subject: [PATCH 03/21] fixed inline comment error --- gpkit/constraints/gp.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/gpkit/constraints/gp.py b/gpkit/constraints/gp.py index 3bf5cc79e..0a4282ca0 100644 --- a/gpkit/constraints/gp.py +++ b/gpkit/constraints/gp.py @@ -93,7 +93,7 @@ def __init__(self, cost, constraints, substitutions, self.posynomials.extend(self.as_posyslt1(self.substitutions)) self.hmaps = [p.hmap for p in self.posynomials] ## Generate various maps into the posy- and monomials - # k [j]: number of monomials (columns of F) present in each constraint + # k [j]: number of monomials (rows of A) present in each constraint self.k = [len(hm) for hm in self.hmaps] p_idxs = [] # p_idxs [i]: posynomial index of each monomial self.m_idxs = [] # m_idxs [i]: monomial indices of each posynomial From 4c464d933e521485ee4561bd4a3c76cb3f3c8837 Mon Sep 17 00:00:00 2001 From: rileyjmurray Date: Sun, 19 Jan 2020 14:19:13 -0800 Subject: [PATCH 04/21] more docstring changes --- gpkit/_cvxopt.py | 5 +++-- gpkit/_mosek/mosek9.py | 0 2 files changed, 3 insertions(+), 2 deletions(-) create mode 100644 gpkit/_mosek/mosek9.py diff --git a/gpkit/_cvxopt.py b/gpkit/_cvxopt.py index 3754c768b..c223e5d5e 100644 --- a/gpkit/_cvxopt.py +++ b/gpkit/_cvxopt.py @@ -12,7 +12,7 @@ def cvxoptimize(c, A, k, *args, **kwargs): "[a,b] array of floats" indicates array-like data with shape [a,b] n is the number of monomials in the gp m is the number of variables in the gp - p is the number of posynomials in the gp + p is the number of posynomial constraints in the gp Arguments --------- @@ -21,7 +21,8 @@ def cvxoptimize(c, A, k, *args, **kwargs): A : floats array of shape (n, m) Exponents of the various free variables for each monomial. k : ints array of shape p+1 - number of monomials (columns of F) present in each constraint + k[0] is the number of monomials (rows of A) present in the objective + k[1:] is number of monomials (rows of A) present in each constraint Returns ------- diff --git a/gpkit/_mosek/mosek9.py b/gpkit/_mosek/mosek9.py new file mode 100644 index 000000000..e69de29bb From 1b3dd6c2dfaef3c73e0d6aa65f8b311b0d959745 Mon Sep 17 00:00:00 2001 From: rileyjmurray Date: Sun, 19 Jan 2020 20:35:04 -0800 Subject: [PATCH 05/21] basic, untested mosek 9+ interface --- gpkit/_cvxopt.py | 2 +- gpkit/_mosek/mosek9.py | 0 gpkit/_mosek/mosek_conif.py | 132 ++++++++++++++++++++++++++++++++++++ gplibrary | 1 + 4 files changed, 134 insertions(+), 1 deletion(-) delete mode 100644 gpkit/_mosek/mosek9.py create mode 100644 gpkit/_mosek/mosek_conif.py create mode 160000 gplibrary diff --git a/gpkit/_cvxopt.py b/gpkit/_cvxopt.py index c223e5d5e..7bd26620c 100644 --- a/gpkit/_cvxopt.py +++ b/gpkit/_cvxopt.py @@ -22,7 +22,7 @@ def cvxoptimize(c, A, k, *args, **kwargs): Exponents of the various free variables for each monomial. k : ints array of shape p+1 k[0] is the number of monomials (rows of A) present in the objective - k[1:] is number of monomials (rows of A) present in each constraint + k[1:] is the number of monomials (rows of A) present in each constraint Returns ------- diff --git a/gpkit/_mosek/mosek9.py b/gpkit/_mosek/mosek9.py deleted file mode 100644 index e69de29bb..000000000 diff --git a/gpkit/_mosek/mosek_conif.py b/gpkit/_mosek/mosek_conif.py new file mode 100644 index 000000000..fde8639cb --- /dev/null +++ b/gpkit/_mosek/mosek_conif.py @@ -0,0 +1,132 @@ +"Implements the GPkit interface to MOSEK (version >= 9) python-based Optimizer API" + +import numpy as np + + +def mskoptimize(c, A, k, p_idxs, *args, **kwargs): + """ + Definitions + ----------- + "[a,b] array of floats" indicates array-like data with shape [a,b] + n is the number of monomials in the gp + m is the number of variables in the gp + p is the number of posynomial constraints in the gp + + Arguments + --------- + c : floats array of shape n + Coefficients of each monomial + A : gpkit.small_classes.CootMatrix, of shape (n, m) + Exponents of the various free variables for each monomial. + k : ints array of shape p+1 + k[0] is the number of monomials (rows of A) present in the objective + k[1:] is the number of monomials (rows of A) present in each constraint + p_idxs : list of bool arrays, each of shape m + p_idxs[i] selects rows of A and entries of c of the i-th posynomial + fi(x) = c[p_idxs[i]] @ exp(A[p_idxs[i],:] @ x). The 0-th posynomial + gives the objective function, and the remaining posynomials should + be constrained to be <= 1. + + Returns + ------- + dict + Contains the following keys + "success": bool + "objective_sol" float + Optimal value of the objective + "primal_sol": floats array of size m + Optimal value of free variables. Note: not in logspace. + "dual_sol": floats array of size p + Optimal value of the dual variables, in logspace. + """ + import mosek + # + # Prepare some initial problem data + # + # Say that the optimization variable is y = [x;t], where + # x is of length m, and t is of length 3*n. Entry + # t[3*ell] is the epigraph variable for the ell-th monomial + # t[3*ell + 1] should be == 1. + # t[3*ell + 2] should be == A[ell, :] @ x. + # + c = np.array(c) + n, m = A.shape + p = len(p_idxs) - 1 + n_msk = m + 3*n + # + # Create MOSEK task. add variables and conic constraints. + # + env = mosek.Env() + task = env.Task(0, 0) + task.appendvars(n_msk) + task.putvarboundlist(np.arange(n_msk, dtype=int), + [mosek.boundkey.fr] * n_msk, + np.zeros(n_msk), + np.zeros(n_msk)) + for i in range(n): + idx = m + 3*i + task.appendcone(mosek.conetype.pexp, 0.0, np.arange(idx, idx + 3)) + # + # Calls to MOSEK's "putaijlist" + # + task.appendcons(2*n + p) + # Linear equations: t[3*ell + 1] == 1. + rows = [i for i in range(n)] + cols = (m + 3*np.arange(n) + 1).tolist() + vals = [1.0] * n + task.putaijlist(rows, cols, vals) + # Linear equations: A[ell,:] @ x - t[3*ell + 2] == 0 + cur_con_idx = n + rows = [cur_con_idx + r for r in A.rows] + task.putaijlist(rows, A.cols, A.vals) + rows = [cur_con_idx + i for i in range(n)] + cols = (m + 3*np.arange(n) + 2).tolist() + vals = [-1.0] * n + task.putaijlist(rows, cols, vals) + # Linear inequalities: c[sels] @ t[3 * sels] <= 1, for sels = np.nonzero(p_idxs[i]) + cur_con_idx = 2*n + rows, cols, vals = [], [], [] + for i in range(p): + sels = np.nonzero(p_idxs[i + 1])[0] + rows.extend([cur_con_idx] * sels.size) + cols.extend(m + 3 * sels) + vals.extend(c[sels]) + cur_con_idx += 1 + task.putaijlist(rows, cols, vals) + # + # Build the right-hand-sides of the [in]equality constraints + # + type_constraint = [mosek.boundkey.fx] * (2*n) + [mosek.boundkey.up] * p + h = np.concatenate([np.ones(n), np.zeros(n), np.ones(p)]) + task.putconboundlist(np.arange(h.size, dtype=int), type_constraint, h, h) + # + # Set the objective function + # + sels = np.nonzero(p_idxs[0])[0] + cols = (m + 3 * sels).tolist() + task.putclist(cols, c[sels].tolist()) + task.putobjsense(mosek.objsense.minimize) + # + # Set solver parameters, and call .solve(). + # + task.optimize() + # + # Recover the solution + # + msk_solsta = task.getsolsta(mosek.soltype.itr) + if msk_solsta == mosek.solsta.optimal: + str_status = 'optimal' + x = [0.] * m + task.getxxslice(mosek.soltype.itr, 0, len(x), x) + solution = {'status': str_status, 'primal': np.array(x)} + return solution + elif msk_solsta == mosek.solsta.prim_infeas_cer: + str_status = 'infeasible' + solution = {'status': str_status, 'primal': None} + return solution + elif msk_solsta == mosek.solsta.dual_infeas_cer: + str_status = 'unbounded' + solution = {'status': str_status, 'primal': None} + return solution + else: + raise RuntimeError('Unexpected solver status.') diff --git a/gplibrary b/gplibrary new file mode 160000 index 000000000..98a82984b --- /dev/null +++ b/gplibrary @@ -0,0 +1 @@ +Subproject commit 98a82984be73b23a766206de406f4ad857a73ce0 From c4a4bf8fbc8bcf1adceda15a4b4cbf42c8b3f57a Mon Sep 17 00:00:00 2001 From: rileyjmurray Date: Mon, 20 Jan 2020 11:05:34 -0800 Subject: [PATCH 06/21] working mosek interface, but dual variables dont match up with gpkit expectations. Need to reformulate the primal problem in order to align with gpkit expectations. This may cause things to break; committing current work so as to not lose it. --- gpkit/_mosek/mosek_conif.py | 35 ++++++++++++++++++++--------------- 1 file changed, 20 insertions(+), 15 deletions(-) diff --git a/gpkit/_mosek/mosek_conif.py b/gpkit/_mosek/mosek_conif.py index fde8639cb..cc9a9f41b 100644 --- a/gpkit/_mosek/mosek_conif.py +++ b/gpkit/_mosek/mosek_conif.py @@ -51,7 +51,7 @@ def mskoptimize(c, A, k, p_idxs, *args, **kwargs): # c = np.array(c) n, m = A.shape - p = len(p_idxs) - 1 + p = len(k) - 1 n_msk = m + 3*n # # Create MOSEK task. add variables and conic constraints. @@ -75,36 +75,36 @@ def mskoptimize(c, A, k, p_idxs, *args, **kwargs): cols = (m + 3*np.arange(n) + 1).tolist() vals = [1.0] * n task.putaijlist(rows, cols, vals) - # Linear equations: A[ell,:] @ x - t[3*ell + 2] == 0 + # Linear equations: A[ell,:] @ x - t[3*ell + 2] == -np.log(c[ell]) cur_con_idx = n - rows = [cur_con_idx + r for r in A.rows] - task.putaijlist(rows, A.cols, A.vals) + rows = [cur_con_idx + r for r in A.row] + task.putaijlist(rows, A.col, A.data) rows = [cur_con_idx + i for i in range(n)] cols = (m + 3*np.arange(n) + 2).tolist() vals = [-1.0] * n task.putaijlist(rows, cols, vals) - # Linear inequalities: c[sels] @ t[3 * sels] <= 1, for sels = np.nonzero(p_idxs[i]) + # Linear inequalities: 1 @ t[3 * sels] <= 1, for sels = np.nonzero(p_idxs[i]) cur_con_idx = 2*n rows, cols, vals = [], [], [] for i in range(p): - sels = np.nonzero(p_idxs[i + 1])[0] + sels = np.nonzero(p_idxs == (i + 1))[0] rows.extend([cur_con_idx] * sels.size) cols.extend(m + 3 * sels) - vals.extend(c[sels]) + vals.extend([1] * sels.size) cur_con_idx += 1 task.putaijlist(rows, cols, vals) # # Build the right-hand-sides of the [in]equality constraints # type_constraint = [mosek.boundkey.fx] * (2*n) + [mosek.boundkey.up] * p - h = np.concatenate([np.ones(n), np.zeros(n), np.ones(p)]) + h = np.concatenate([np.ones(n), -np.log(c), np.ones(p)]) task.putconboundlist(np.arange(h.size, dtype=int), type_constraint, h, h) # # Set the objective function # - sels = np.nonzero(p_idxs[0])[0] + sels = np.nonzero(p_idxs == 0)[0] cols = (m + 3 * sels).tolist() - task.putclist(cols, c[sels].tolist()) + task.putclist(cols, [1] * sels.size) task.putobjsense(mosek.objsense.minimize) # # Set solver parameters, and call .solve(). @@ -117,16 +117,21 @@ def mskoptimize(c, A, k, p_idxs, *args, **kwargs): if msk_solsta == mosek.solsta.optimal: str_status = 'optimal' x = [0.] * m - task.getxxslice(mosek.soltype.itr, 0, len(x), x) - solution = {'status': str_status, 'primal': np.array(x)} - return solution + task.getxxslice(mosek.soltype.itr, 0, m, x) + z = [0.] * p + task.getsucslice(mosek.soltype.itr, 2*n, 2*n + p, z) + z = np.array(z) + z[z < 0] = 0 + x = np.array(x) + solution = {'status': str_status, 'primal': x, 'la': z} elif msk_solsta == mosek.solsta.prim_infeas_cer: str_status = 'infeasible' solution = {'status': str_status, 'primal': None} - return solution elif msk_solsta == mosek.solsta.dual_infeas_cer: str_status = 'unbounded' solution = {'status': str_status, 'primal': None} - return solution else: raise RuntimeError('Unexpected solver status.') + task.__exit__(None, None, None) + env.__exit__(None, None, None) + return solution From b83c85d4f8111a52a7dfad779a173068ca84b79c Mon Sep 17 00:00:00 2001 From: rileyjmurray Date: Mon, 20 Jan 2020 11:17:04 -0800 Subject: [PATCH 07/21] make sure this class can see the mosek_conif solver --- gpkit/constraints/gp.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/gpkit/constraints/gp.py b/gpkit/constraints/gp.py index 0a4282ca0..a097b1c05 100644 --- a/gpkit/constraints/gp.py +++ b/gpkit/constraints/gp.py @@ -15,7 +15,7 @@ DEFAULT_SOLVER_KWARGS = {"cvxopt": {"kktsolver": "ldl"}} -SOLUTION_TOL = {"cvxopt": 1e-3, "mosek_cli": 1e-4, "mosek": 1e-5} +SOLUTION_TOL = {"cvxopt": 1e-3, "mosek_cli": 1e-4, "mosek": 1e-5, 'mosek_conif': 1e-5} def _get_solver(solver, kwargs): @@ -38,6 +38,9 @@ def _get_solver(solver, kwargs): elif solver == "mosek": from .._mosek import expopt solverfn = expopt.imize + elif solver =='mosek_conif': + from .._mosek import mosek_conif + solverfn = mosek_conif.mskoptimize elif hasattr(solver, "__call__"): solverfn = solver solver = solver.__name__ From 39e3b56f3110d73f6e132fd4269fb47f5e53a209 Mon Sep 17 00:00:00 2001 From: rileyjmurray Date: Mon, 20 Jan 2020 11:18:23 -0800 Subject: [PATCH 08/21] add finding mosek python interface to the build process --- gpkit/build.py | 17 ++++++++++++++++- 1 file changed, 16 insertions(+), 1 deletion(-) diff --git a/gpkit/build.py b/gpkit/build.py index be28a5249..b910c5dad 100644 --- a/gpkit/build.py +++ b/gpkit/build.py @@ -122,6 +122,21 @@ def look(self): pass +class MosekConif(SolverBackend): + + name = 'mosek_conif' + + def look(self): + "Attempts to import mosek." + try: + log("# Trying to import mosek...") + # Testing the import, so the variable is intentionally not used + import mosek # pylint: disable=unused-variable + return "in Python path" + except ImportError: + pass + + class Mosek(SolverBackend): "MOSEK finder and builder." name = "mosek" @@ -279,7 +294,7 @@ def build(): log("Started building gpkit...\n") log("Attempting to find and build solvers:\n") - solvers = [Mosek(), MosekCLI(), CVXopt()] + solvers = [MosekConif(), Mosek(), MosekCLI(), CVXopt()] installed_solvers = [solver.name for solver in solvers if solver.installed] From d29905d90650fc4deca8f9f4ef888f32f64caf71 Mon Sep 17 00:00:00 2001 From: rileyjmurray Date: Mon, 20 Jan 2020 11:23:36 -0800 Subject: [PATCH 09/21] specify solvers in all tests. minor additions for mosek_conif interface --- gpkit/tests/t_model.py | 52 +++++++++++++++++++++--------------------- 1 file changed, 26 insertions(+), 26 deletions(-) diff --git a/gpkit/tests/t_model.py b/gpkit/tests/t_model.py index d891c6a01..1b1e041e0 100644 --- a/gpkit/tests/t_model.py +++ b/gpkit/tests/t_model.py @@ -6,7 +6,7 @@ SignomialsEnabled, ArrayVariable, SignomialEquality) from gpkit.constraints.bounded import Bounded from gpkit.small_classes import CootMatrix -from gpkit.exceptions import InvalidGPConstraint +from gpkit.exceptions import InvalidGPConstraint, InvalidPosynomial from gpkit import NamedVariables, units, parse_variables from gpkit.constraints.relax import ConstraintsRelaxed from gpkit.constraints.relax import ConstraintsRelaxedEqually @@ -15,7 +15,7 @@ if sys.version_info >= (3, 0): unicode = str # pylint:disable=redefined-builtin,invalid-name -NDIGS = {"cvxopt": 4, "mosek": 5, "mosek_cli": 5} +NDIGS = {"cvxopt": 4, "mosek": 5, "mosek_cli": 5, 'mosek_conif': 4} # name: decimal places of accuracy # pylint: disable=invalid-name,attribute-defined-outside-init,unused-variable,undefined-variable,exec-used @@ -285,26 +285,26 @@ def test_sp_relaxation(self): with self.assertRaises(ValueError): mr1 = Model(x*r1.relaxvars, r1) # no 'prod' mr1 = Model(x*r1.relaxvars.prod()**10, r1) - cost1 = mr1.localsolve(verbosity=0)["cost"] + cost1 = mr1.localsolve(verbosity=0, solver=self.solver)["cost"] self.assertAlmostEqual(cost1/1024, 1, self.ndig) - m.debug(verbosity=0) + m.debug(verbosity=0, solver=self.solver) with SignomialsEnabled(): m = Model(x, [x+y >= z, x+y <= z/2, y <= x, y >= 1], {z: 3}) if self.solver != "cvxopt": - m.debug(verbosity=0) + m.debug(verbosity=0, solver=self.solver) r2 = ConstraintsRelaxed(m) self.assertEqual(len(r2.varkeys), 7) mr2 = Model(x*r2.relaxvars.prod()**10, r2) - cost2 = mr2.localsolve(verbosity=0)["cost"] + cost2 = mr2.localsolve(verbosity=0, solver=self.solver)["cost"] self.assertAlmostEqual(cost2/1024, 1, self.ndig) with SignomialsEnabled(): m = Model(x, [x+y >= z, x+y <= z/2, y <= x, y >= 1], {z: 3}) if self.solver != "cvxopt": - m.debug(verbosity=0) + m.debug(verbosity=0, solver=self.solver) r3 = ConstraintsRelaxedEqually(m) self.assertEqual(len(r3.varkeys), 4) mr3 = Model(x*r3.relaxvar**10, r3) - cost3 = mr3.localsolve(verbosity=0)["cost"] + cost3 = mr3.localsolve(verbosity=0, solver=self.solver)["cost"] self.assertAlmostEqual(cost3/(32*0.8786796585), 1, self.ndig) def test_sp_bounded(self): @@ -313,17 +313,17 @@ def test_sp_bounded(self): with SignomialsEnabled(): m = Model(x, [x + y >= 1, y <= 0.1]) # solves - cost = m.localsolve(verbosity=0)["cost"] + cost = m.localsolve(verbosity=0, solver=self.solver)["cost"] self.assertAlmostEqual(cost, 0.9, self.ndig) with SignomialsEnabled(): m = Model(x, [x + y >= 1]) # dual infeasible with self.assertRaises((RuntimeWarning, ValueError)): - m.localsolve(verbosity=0) + m.localsolve(verbosity=0, solver=self.solver), with SignomialsEnabled(): m = Model(x, Bounded([x + y >= 1], verbosity=0)) - sol = m.localsolve(verbosity=0) + sol = m.localsolve(verbosity=0, solver=self.solver) boundedness = sol["boundedness"] if "value near lower bound" in boundedness: self.assertIn(x.key, boundedness["value near lower bound"]) @@ -341,7 +341,7 @@ def test_values_vs_subs(self): y >= x - 1] m = Model(x + y*z, constraints) m.substitutions.update({"z": 5}) - sol = m.localsolve(verbosity=0) + sol = m.localsolve(verbosity=0, solver=self.solver) self.assertAlmostEqual(sol["cost"], 13, self.ndig) # Constant variable declaration method @@ -350,7 +350,7 @@ def test_values_vs_subs(self): constraints = [x + y >= z, y >= x - 1] m = Model(x + y*z, constraints) - sol = m.localsolve(verbosity=0) + sol = m.localsolve(verbosity=0, solver=self.solver) self.assertAlmostEqual(sol["cost"], 13, self.ndig) def test_initially_infeasible(self): @@ -363,7 +363,7 @@ def test_initially_infeasible(self): m = Model(1/x, [sigc, sigc2, y <= 0.5]) - sol = m.localsolve(verbosity=0) + sol = m.localsolve(verbosity=0, solver=self.solver) self.assertAlmostEqual(sol["cost"], 2**0.5, self.ndig) self.assertAlmostEqual(sol(y), 0.5, self.ndig) second_solve_key_names = [key.name[:5] @@ -383,14 +383,14 @@ def test_sp_substitutions(self): with SignomialsEnabled(): m = Model(x, [x + z >= y]) with self.assertRaises(ValueError): - m.localsolve(verbosity=0) + m.localsolve(verbosity=0, solver=self.solver) with SignomialsEnabled(): m = Model(x, [x + y >= z]) m.substitutions[y] = 1 m.substitutions[z] = 4 sol = m.solve(self.solver, verbosity=0) - self.assertAlmostEqual(sol["cost"], 3) + self.assertAlmostEqual(sol["cost"], 3, self.ndig) sys.stdout = old_stdout self.assertEqual(stringout.getvalue(), ( @@ -432,7 +432,7 @@ def test_impossible(self): m1 = Model(x, [x + y >= z, x >= y]) m1.substitutions.update({'x': 0, 'y': 0}) with self.assertRaises(ValueError): - _ = m1.localsolve() + _ = m1.localsolve(solver=self.solver) def test_trivial_sp(self): x = Variable('x') @@ -440,7 +440,7 @@ def test_trivial_sp(self): with SignomialsEnabled(): m = Model(x, [x >= 1-y, y <= 0.1]) with self.assertRaises(InvalidGPConstraint): - m.solve(verbosity=0) + m.solve(verbosity=0, solver=self.solver) sol = m.localsolve(self.solver, verbosity=0) self.assertAlmostEqual(sol["variables"]["x"], 0.9, self.ndig) with SignomialsEnabled(): @@ -455,7 +455,7 @@ def test_tautological_spconstraint(self): with SignomialsEnabled(): m = Model(x, [x >= 1-y, y <= 0.1, y >= z]) with self.assertRaises(InvalidGPConstraint): - m.solve(verbosity=0) + m.solve(verbosity=0, solver=self.solver) sol = m.localsolve(self.solver, verbosity=0) self.assertAlmostEqual(sol["variables"]["x"], 0.9, self.ndig) @@ -466,7 +466,7 @@ def test_relaxation(self): constraints = [y + x >= 2, y <= x] objective = x m = Model(objective, constraints) - m.localsolve(verbosity=0) + m.localsolve(verbosity=0, solver=self.solver) # issue #257 @@ -479,7 +479,7 @@ def test_relaxation(self): C <= 1] obj = 1/A[0] + 1/A[1] m = Model(obj, constraints) - m.localsolve(verbosity=0) + m.localsolve(verbosity=0, solver=self.solver) def test_issue180(self): L = Variable("L") @@ -574,7 +574,7 @@ def test_small_named_signomial(self): J = 0.01*(x - 1)**2 + nonzero_adder with NamedVariables("SmallSignomial"): m = Model(z, [z >= J]) - sol = m.localsolve(verbosity=0) + sol = m.localsolve(verbosity=0, solver=self.solver) self.assertAlmostEqual(sol['cost'], nonzero_adder, local_ndig) self.assertAlmostEqual(sol('x'), 0.98725425, self.ndig) @@ -585,7 +585,7 @@ def test_sigs_not_allowed_in_cost(self): J = 0.01*((x - 1)**2 + (y - 1)**2) + (x*y - 1)**2 m = Model(J) with self.assertRaises(TypeError): - m.localsolve(verbosity=0) + m.localsolve(verbosity=0, solver=self.solver) def test_partial_sub_signomial(self): "Test SP partial x0 initialization" @@ -607,7 +607,7 @@ def test_becomes_signomial(self): with SignomialsEnabled(): _ = m.gp() with self.assertRaises(RuntimeWarning): - _ = m.localsolve() + _ = m.localsolve(solver=self.solver) def test_reassigned_constant_cost(self): # for issue 1131 @@ -626,7 +626,7 @@ def test_unbounded_debugging(self): x = Variable("x") y = Variable("y") m = Model(x*y, [x*y**1.01 >= 100]) - with self.assertRaises((RuntimeWarning, ValueError)): + with self.assertRaises(InvalidPosynomial): m.solve(self.solver, verbosity=0) # test one-sided bound m = Model(x*y, Bounded(m, verbosity=0, lower=0.001)) @@ -645,7 +645,7 @@ def test_unbounded_debugging(self): def test_penalty_ccp_solve(self): "Test penalty convex-concave algorithm from [Lipp/Boyd 2016]." m = SPThing() - sol = m.localsolve(verbosity=0) + sol = m.localsolve(verbosity=0, solver=self.solver) sol_pccp = m.penalty_ccp_solve(verbosity=0) self.assertEqual(len(m.program.gps[-1].varkeys), 3) self.assertAlmostEqual(sol['cost'], sol_pccp['cost']) From 80a8fde793cfd2f6745488b9251d30318dbd0b7d Mon Sep 17 00:00:00 2001 From: Riley Murray Date: Mon, 20 Jan 2020 13:56:26 -0800 Subject: [PATCH 10/21] precision in some tests. line spacing issues. --- gpkit/tests/t_model.py | 19 ++++++++++++------- 1 file changed, 12 insertions(+), 7 deletions(-) diff --git a/gpkit/tests/t_model.py b/gpkit/tests/t_model.py index 1b1e041e0..fb15c8503 100644 --- a/gpkit/tests/t_model.py +++ b/gpkit/tests/t_model.py @@ -29,7 +29,7 @@ class TestGP(unittest.TestCase): name = "TestGP_" # solver and ndig get set in loop at bottom this file, a bit hacky solver = None - ndig = None + ndig = 4 def test_trivial_gp(self): """ @@ -64,7 +64,7 @@ def test_sigeq(self): SignomialEquality(x**2 + x, y)]) sol = m.localsolve(solver=self.solver, verbosity=0, mutategp=False) self.assertAlmostEqual(sol("x"), 0.1639472, self.ndig) - self.assertAlmostEqual(sol("y")[0], 0.1908254, self.ndig) + self.assertAlmostEqual(sol("y")[0], 0.1908254, 3) self.assertAlmostEqual(sol("c"), 0.2669448, self.ndig) # test right vector input to sigeq with SignomialsEnabled(): @@ -80,9 +80,9 @@ def test_sigeq(self): m = Model(c, [c >= (x + 0.25)**2 + (y - 0.5)**2, SignomialEquality(x**2 + x, y)]) sol = m.localsolve(solver=self.solver, verbosity=0) - self.assertAlmostEqual(sol("x"), 0.1639472, self.ndig) - self.assertAlmostEqual(sol("y"), 0.1908254, self.ndig) - self.assertAlmostEqual(sol("c"), 0.2669448, self.ndig) + self.assertAlmostEqual(sol("x"), 0.1639472, 3) + self.assertAlmostEqual(sol("y"), 0.1908254, 3) + self.assertAlmostEqual(sol("c"), 0.2669448, 3) def test_601(self): # tautological monomials should solve but not pass to the solver @@ -271,7 +271,7 @@ class TestSP(unittest.TestCase): """test case for SP class -- gets run for each installed solver""" name = "TestSP_" solver = None - ndig = None + ndig = 4 def test_sp_relaxation(self): w = Variable('w') @@ -626,7 +626,7 @@ def test_unbounded_debugging(self): x = Variable("x") y = Variable("y") m = Model(x*y, [x*y**1.01 >= 100]) - with self.assertRaises(InvalidPosynomial): + with self.assertRaises((RuntimeWarning, InvalidPosynomial)): m.solve(self.solver, verbosity=0) # test one-sided bound m = Model(x*y, Bounded(m, verbosity=0, lower=0.001)) @@ -650,6 +650,7 @@ def test_penalty_ccp_solve(self): self.assertEqual(len(m.program.gps[-1].varkeys), 3) self.assertAlmostEqual(sol['cost'], sol_pccp['cost']) + class TestModelSolverSpecific(unittest.TestCase): """test cases run only for specific solvers""" def test_cvxopt_kwargs(self): @@ -670,11 +671,13 @@ def setup(self, length): c = Variable("c", 17/4., "g") return [a >= c/b] + class Thing2(Model): "another thing for model testing" def setup(self): return [Thing(2), Model()] + class SPThing(Model): "a simple SP" def setup(self): @@ -686,6 +689,7 @@ def setup(self): self.cost = 1/z return constraints + class Box(Model): """simple box for model testing @@ -779,6 +783,7 @@ def test_duplicate_submodel_varnames(self): self.assertEqual(len(w.varkeys.keymap["m"]), 2) w2 = Widget() + TESTS = [TestModelSolverSpecific, TestModelNoSolve] MULTI_SOLVER_TESTS = [TestGP, TestSP] From 517a267aaa002f9fd1c6f310ca6acf7db8a9bee5 Mon Sep 17 00:00:00 2001 From: Riley Murray Date: Mon, 20 Jan 2020 13:57:01 -0800 Subject: [PATCH 11/21] switched to a log-sum-exp formulation; dual variables recovering correctly. --- gpkit/_mosek/mosek_conif.py | 70 ++++++++++++++++++++++--------------- 1 file changed, 42 insertions(+), 28 deletions(-) diff --git a/gpkit/_mosek/mosek_conif.py b/gpkit/_mosek/mosek_conif.py index cc9a9f41b..ef206df1b 100644 --- a/gpkit/_mosek/mosek_conif.py +++ b/gpkit/_mosek/mosek_conif.py @@ -22,10 +22,9 @@ def mskoptimize(c, A, k, p_idxs, *args, **kwargs): k[0] is the number of monomials (rows of A) present in the objective k[1:] is the number of monomials (rows of A) present in each constraint p_idxs : list of bool arrays, each of shape m - p_idxs[i] selects rows of A and entries of c of the i-th posynomial - fi(x) = c[p_idxs[i]] @ exp(A[p_idxs[i],:] @ x). The 0-th posynomial - gives the objective function, and the remaining posynomials should - be constrained to be <= 1. + sel = p_idxs == i selects rows of A and entries of c of the i-th posynomial + fi(x) = c[sel] @ exp(A[sel,:] @ x). The 0-th posynomial gives the objective + function, and the remaining posynomials should be constrained to be <= 1. Returns ------- @@ -43,39 +42,48 @@ def mskoptimize(c, A, k, p_idxs, *args, **kwargs): # # Prepare some initial problem data # - # Say that the optimization variable is y = [x;t], where - # x is of length m, and t is of length 3*n. Entry - # t[3*ell] is the epigraph variable for the ell-th monomial - # t[3*ell + 1] should be == 1. - # t[3*ell + 2] should be == A[ell, :] @ x. + # Say that the optimization variable is y = [x;t;z], where + # x is of length m, t is of length 3*n, and z is of length p+1. + # + # Triples + # t[3*ell] + # t[3*ell + 1] should be == 1. + # t[3*ell + 2] should be == A[ell, :] @ x + np.log(c[ell]) - z[posynom idx corresponding to ell]. + # should belong to the exponential cone. + # + # For each i from {0,...,p}, the "t" should also satisfy + # sum(t[3*ell] for ell correponsing to posynomial i) <= 1. + # + # The vector "z" should satisfy + # z[1:] <= 0 + # # c = np.array(c) n, m = A.shape p = len(k) - 1 - n_msk = m + 3*n + n_msk = m + 3*n + p + 1 # # Create MOSEK task. add variables and conic constraints. # env = mosek.Env() task = env.Task(0, 0) task.appendvars(n_msk) + bound_types = [mosek.boundkey.fr] * (m + 3*n + 1) + [mosek.boundkey.up] * p task.putvarboundlist(np.arange(n_msk, dtype=int), - [mosek.boundkey.fr] * n_msk, - np.zeros(n_msk), - np.zeros(n_msk)) + bound_types, np.zeros(n_msk), np.zeros(n_msk)) for i in range(n): idx = m + 3*i task.appendcone(mosek.conetype.pexp, 0.0, np.arange(idx, idx + 3)) # # Calls to MOSEK's "putaijlist" # - task.appendcons(2*n + p) + task.appendcons(2*n + p + 1) # Linear equations: t[3*ell + 1] == 1. rows = [i for i in range(n)] cols = (m + 3*np.arange(n) + 1).tolist() vals = [1.0] * n task.putaijlist(rows, cols, vals) - # Linear equations: A[ell,:] @ x - t[3*ell + 2] == -np.log(c[ell]) + # Linear equations: A[ell,:] @ x - t[3*ell + 2] - z[posynom idx corresponding to ell] == -np.log(c[ell]) cur_con_idx = n rows = [cur_con_idx + r for r in A.row] task.putaijlist(rows, A.col, A.data) @@ -83,11 +91,15 @@ def mskoptimize(c, A, k, p_idxs, *args, **kwargs): cols = (m + 3*np.arange(n) + 2).tolist() vals = [-1.0] * n task.putaijlist(rows, cols, vals) - # Linear inequalities: 1 @ t[3 * sels] <= 1, for sels = np.nonzero(p_idxs[i]) + rows = [cur_con_idx + i for i in range(n)] + cols = [m + 3*n + p_idxs[i] for i in range(n)] + vals = [-1.0] * n + task.putaijlist(rows, cols, vals) + # Linear inequalities: 1 @ t[3 * sels] <= 1, for sels = np.nonzero(p_idxs[i] == i) cur_con_idx = 2*n rows, cols, vals = [], [], [] - for i in range(p): - sels = np.nonzero(p_idxs == (i + 1))[0] + for i in range(p+1): + sels = np.nonzero(p_idxs == i)[0] rows.extend([cur_con_idx] * sels.size) cols.extend(m + 3 * sels) vals.extend([1] * sels.size) @@ -96,15 +108,13 @@ def mskoptimize(c, A, k, p_idxs, *args, **kwargs): # # Build the right-hand-sides of the [in]equality constraints # - type_constraint = [mosek.boundkey.fx] * (2*n) + [mosek.boundkey.up] * p - h = np.concatenate([np.ones(n), -np.log(c), np.ones(p)]) + type_constraint = [mosek.boundkey.fx] * (2*n) + [mosek.boundkey.up] * (p + 1) + h = np.concatenate([np.ones(n), -np.log(c), np.ones(p + 1)]) task.putconboundlist(np.arange(h.size, dtype=int), type_constraint, h, h) # # Set the objective function # - sels = np.nonzero(p_idxs == 0)[0] - cols = (m + 3 * sels).tolist() - task.putclist(cols, [1] * sels.size) + task.putclist([int(m + 3*n)], [1]) task.putobjsense(mosek.objsense.minimize) # # Set solver parameters, and call .solve(). @@ -116,14 +126,18 @@ def mskoptimize(c, A, k, p_idxs, *args, **kwargs): msk_solsta = task.getsolsta(mosek.soltype.itr) if msk_solsta == mosek.solsta.optimal: str_status = 'optimal' + # recover primal variables x = [0.] * m task.getxxslice(mosek.soltype.itr, 0, m, x) - z = [0.] * p - task.getsucslice(mosek.soltype.itr, 2*n, 2*n + p, z) - z = np.array(z) - z[z < 0] = 0 x = np.array(x) - solution = {'status': str_status, 'primal': x, 'la': z} + # recover dual variables for log-sum-exp epigraph constraints + # (skip epigraph of the objective function). + z_duals = [0.] * p + task.getsuxslice(mosek.soltype.itr, m + 3*n + 1, n_msk, z_duals) + z_duals = np.array(z_duals) + z_duals[z_duals < 0] = 0 + # wrap things up in a dictionary + solution = {'status': str_status, 'primal': x, 'la': z_duals} elif msk_solsta == mosek.solsta.prim_infeas_cer: str_status = 'infeasible' solution = {'status': str_status, 'primal': None} From 152511cfc4a3dea7cff608137cd809ce8ab9a640 Mon Sep 17 00:00:00 2001 From: Riley Murray Date: Mon, 20 Jan 2020 13:57:48 -0800 Subject: [PATCH 12/21] precision requirements --- gpkit/constraints/gp.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/gpkit/constraints/gp.py b/gpkit/constraints/gp.py index a097b1c05..1b103413c 100644 --- a/gpkit/constraints/gp.py +++ b/gpkit/constraints/gp.py @@ -15,7 +15,7 @@ DEFAULT_SOLVER_KWARGS = {"cvxopt": {"kktsolver": "ldl"}} -SOLUTION_TOL = {"cvxopt": 1e-3, "mosek_cli": 1e-4, "mosek": 1e-5, 'mosek_conif': 1e-5} +SOLUTION_TOL = {"cvxopt": 1e-3, "mosek_cli": 1e-4, "mosek": 1e-5, 'mosek_conif': 1e-3} def _get_solver(solver, kwargs): @@ -38,7 +38,7 @@ def _get_solver(solver, kwargs): elif solver == "mosek": from .._mosek import expopt solverfn = expopt.imize - elif solver =='mosek_conif': + elif solver == 'mosek_conif': from .._mosek import mosek_conif solverfn = mosek_conif.mskoptimize elif hasattr(solver, "__call__"): From 791d58ad7fe383fca6d92b2a598e419ab82f0c2c Mon Sep 17 00:00:00 2001 From: Riley Murray Date: Mon, 20 Jan 2020 14:19:18 -0800 Subject: [PATCH 13/21] precision and API changes to some tests. Remaining failing tests are because of more serious unit-conversion or string-printing issues. --- gpkit/tests/t_constraints.py | 2 +- gpkit/tests/t_examples.py | 2 +- gpkit/tests/t_model.py | 2 +- gpkit/tests/t_solution_array.py | 7 ++++--- 4 files changed, 7 insertions(+), 6 deletions(-) diff --git a/gpkit/tests/t_constraints.py b/gpkit/tests/t_constraints.py index 21c87d023..9dc90bd17 100644 --- a/gpkit/tests/t_constraints.py +++ b/gpkit/tests/t_constraints.py @@ -62,7 +62,7 @@ def test_equality_relaxation(self): m = Model(x, [x == 3, x == 4]) rc = ConstraintsRelaxed(m) m2 = Model(rc.relaxvars.prod() * x**0.01, rc) - self.assertAlmostEqual(m2.solve(verbosity=0)(x), 3, 5) + self.assertAlmostEqual(m2.solve(verbosity=0)(x), 3, places=3) def test_constraintget(self): x = Variable("x") diff --git a/gpkit/tests/t_examples.py b/gpkit/tests/t_examples.py index 46eaefa34..f1f69b3f0 100644 --- a/gpkit/tests/t_examples.py +++ b/gpkit/tests/t_examples.py @@ -113,7 +113,7 @@ def test_primal_infeasible_ex1(self, example): with self.assertRaises(RuntimeWarning) as cm: example.m.solve(verbosity=0) err = str(cm.exception) - if "mosek" in err: + if "mosek" in err and 'mosek_conif' not in err: self.assertIn("PRIM_INFEAS_CER", err) elif "cvxopt" in err: self.assertIn("unknown", err) diff --git a/gpkit/tests/t_model.py b/gpkit/tests/t_model.py index fb15c8503..efc567373 100644 --- a/gpkit/tests/t_model.py +++ b/gpkit/tests/t_model.py @@ -626,7 +626,7 @@ def test_unbounded_debugging(self): x = Variable("x") y = Variable("y") m = Model(x*y, [x*y**1.01 >= 100]) - with self.assertRaises((RuntimeWarning, InvalidPosynomial)): + with self.assertRaises((RuntimeWarning, ValueError, InvalidPosynomial)): m.solve(self.solver, verbosity=0) # test one-sided bound m = Model(x*y, Bounded(m, verbosity=0, lower=0.001)) diff --git a/gpkit/tests/t_solution_array.py b/gpkit/tests/t_solution_array.py index eb78e6a36..abbbdc4ae 100644 --- a/gpkit/tests/t_solution_array.py +++ b/gpkit/tests/t_solution_array.py @@ -15,7 +15,7 @@ def test_call(self): A = Variable('A', '-', 'Test Variable') prob = Model(A, [A >= 1]) sol = prob.solve(verbosity=0) - self.assertAlmostEqual(sol(A), 1.0, 10) + self.assertAlmostEqual(sol(A), 1.0, 8) def test_call_units(self): # test from issue541 @@ -35,7 +35,8 @@ def test_call_vector(self): self.assertEqual(type(solx), Quantity) self.assertEqual(type(sol["variables"][x]), np.ndarray) self.assertEqual(solx.shape, (n,)) - self.assertTrue((abs(solx - 2.5*np.ones(n)) < 1e-7).all()) + for i in range(n): + self.assertAlmostEqual(solx[i], 2.5, places=4) def test_subinto(self): Nsweep = 20 @@ -72,7 +73,7 @@ def test_units_sub(self): tminsub = 1000 * gpkit.ureg.lbf m.substitutions.update({Tmin: tminsub}) sol = m.solve(verbosity=0) - self.assertEqual(sol(Tmin), tminsub) + self.assertEqual(sol(Tmin) - tminsub, 0) self.assertFalse( "1000N" in sol.table().replace(" ", "").replace("[", "").replace("]", "")) From 36851575cac98caf69dda6789a735e12115e27d1 Mon Sep 17 00:00:00 2001 From: Riley Murray Date: Mon, 20 Jan 2020 14:56:17 -0800 Subject: [PATCH 14/21] gitignore --- .gitignore | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/.gitignore b/.gitignore index 466bf4ff7..ca54e26a5 100644 --- a/.gitignore +++ b/.gitignore @@ -45,6 +45,10 @@ htmlcov/ nosetests.xml coverage.xml gplibrary/* +gpkit/tests/*.csv +gpkit/tests/*.txt +gpkit/tests/*.mat +gpkit/tests/*.pkl # Translations *.mo From 6efe3e5984123c86c2857520c990c72b712b09df Mon Sep 17 00:00:00 2001 From: Riley Murray Date: Mon, 20 Jan 2020 15:01:18 -0800 Subject: [PATCH 15/21] status handling --- docs/source/examples/autosweep.py | 2 +- gpkit/_mosek/mosek_conif.py | 11 ++++------- gpkit/tests/t_examples.py | 4 +++- 3 files changed, 8 insertions(+), 9 deletions(-) diff --git a/docs/source/examples/autosweep.py b/docs/source/examples/autosweep.py index 27b2a64d3..11561284e 100644 --- a/docs/source/examples/autosweep.py +++ b/docs/source/examples/autosweep.py @@ -36,7 +36,7 @@ # this problem is two intersecting lines in logspace m2 = Model(A**2, [A >= (l/3)**2, A >= (l/3)**0.5 * units.m**1.5]) tol2 = {"mosek": 1e-12, "cvxopt": 1e-7, - "mosek_cli": 1e-6}[gpkit.settings["default_solver"]] + "mosek_cli": 1e-6, 'mosek_conif': 1e-6}[gpkit.settings["default_solver"]] # test Model method sol2 = m2.autosweep({l: [1, 10]}, tol2, verbosity=0) bst2 = sol2.bst diff --git a/gpkit/_mosek/mosek_conif.py b/gpkit/_mosek/mosek_conif.py index ef206df1b..156dc6961 100644 --- a/gpkit/_mosek/mosek_conif.py +++ b/gpkit/_mosek/mosek_conif.py @@ -125,7 +125,6 @@ def mskoptimize(c, A, k, p_idxs, *args, **kwargs): # msk_solsta = task.getsolsta(mosek.soltype.itr) if msk_solsta == mosek.solsta.optimal: - str_status = 'optimal' # recover primal variables x = [0.] * m task.getxxslice(mosek.soltype.itr, 0, m, x) @@ -137,15 +136,13 @@ def mskoptimize(c, A, k, p_idxs, *args, **kwargs): z_duals = np.array(z_duals) z_duals[z_duals < 0] = 0 # wrap things up in a dictionary - solution = {'status': str_status, 'primal': x, 'la': z_duals} + solution = {'status': 'optimal', 'primal': x, 'la': z_duals} elif msk_solsta == mosek.solsta.prim_infeas_cer: - str_status = 'infeasible' - solution = {'status': str_status, 'primal': None} + solution = {'status': 'infeasible', 'primal': None} elif msk_solsta == mosek.solsta.dual_infeas_cer: - str_status = 'unbounded' - solution = {'status': str_status, 'primal': None} + solution = {'status': 'unbounded', 'primal': None} else: - raise RuntimeError('Unexpected solver status.') + solution = {'status': 'unknown', 'primal': None} task.__exit__(None, None, None) env.__exit__(None, None, None) return solution diff --git a/gpkit/tests/t_examples.py b/gpkit/tests/t_examples.py index f1f69b3f0..0c410bd98 100644 --- a/gpkit/tests/t_examples.py +++ b/gpkit/tests/t_examples.py @@ -113,7 +113,9 @@ def test_primal_infeasible_ex1(self, example): with self.assertRaises(RuntimeWarning) as cm: example.m.solve(verbosity=0) err = str(cm.exception) - if "mosek" in err and 'mosek_conif' not in err: + if 'mosek_conif' in err: + self.assertIn('infeasible', err) + elif "mosek" in err: self.assertIn("PRIM_INFEAS_CER", err) elif "cvxopt" in err: self.assertIn("unknown", err) From 980727835dd3904c4091f22958b523efe2f6e265 Mon Sep 17 00:00:00 2001 From: Riley Murray Date: Mon, 20 Jan 2020 17:49:34 -0800 Subject: [PATCH 16/21] improved mosek implementation; monomial posynomial constraints are passed as linear inequalities. --- gpkit/_mosek/mosek_conif.py | 140 +++++++++++++++++++++++++++--------- 1 file changed, 107 insertions(+), 33 deletions(-) diff --git a/gpkit/_mosek/mosek_conif.py b/gpkit/_mosek/mosek_conif.py index 156dc6961..ab4863f41 100644 --- a/gpkit/_mosek/mosek_conif.py +++ b/gpkit/_mosek/mosek_conif.py @@ -59,67 +59,128 @@ def mskoptimize(c, A, k, p_idxs, *args, **kwargs): # # c = np.array(c) - n, m = A.shape - p = len(k) - 1 - n_msk = m + 3*n + p + 1 + # separate exponential from linear constraints. + # the objective is always kept in exponential form. + exp_posys = [0] + [i+1 for i, val in enumerate(k[1:]) if val > 1] + lin_posys = [i for i in range(len(k)) if i not in exp_posys] + if len(lin_posys) > 0: + A = A.tocsr() + lin_idxs = [] + for i in lin_posys: + temp = np.nonzero(p_idxs == i) + temp = temp[0] + lin_idxs.append(temp) + lin_idxs = np.concatenate(lin_idxs) + nonlin_idxs = np.ones(A.shape[0], dtype=bool) + nonlin_idxs[lin_idxs] = False + A_exp = A[nonlin_idxs, :].tocoo() + c_exp = c[nonlin_idxs] + A_lin = A[lin_idxs, :].tocoo() + c_lin = c[lin_idxs] + else: + c_lin = np.array([]) + A_exp = A + c_exp = c + + m = A.shape[1] + k_exp = [k[i] for i in exp_posys] + n_exp = sum(k_exp) + p_exp = len(k_exp) + n_msk = m + 3*n_exp + p_exp + exp_p_idx = [] + for i, ki in enumerate(k_exp): + exp_p_idx.extend([i] * ki) + exp_p_idx = np.array(exp_p_idx) # # Create MOSEK task. add variables and conic constraints. # env = mosek.Env() task = env.Task(0, 0) task.appendvars(n_msk) - bound_types = [mosek.boundkey.fr] * (m + 3*n + 1) + [mosek.boundkey.up] * p + bound_types = [mosek.boundkey.fr] * (m + 3*n_exp + 1) + [mosek.boundkey.up] * (p_exp - 1) task.putvarboundlist(np.arange(n_msk, dtype=int), bound_types, np.zeros(n_msk), np.zeros(n_msk)) - for i in range(n): + for i in range(n_exp): idx = m + 3*i task.appendcone(mosek.conetype.pexp, 0.0, np.arange(idx, idx + 3)) # - # Calls to MOSEK's "putaijlist" + # Affine constraints related to the exponential cone # - task.appendcons(2*n + p + 1) - # Linear equations: t[3*ell + 1] == 1. - rows = [i for i in range(n)] - cols = (m + 3*np.arange(n) + 1).tolist() - vals = [1.0] * n + task.appendcons(2*n_exp + p_exp) + # 1st n_exp: Linear equations: t[3*ell + 1] == 1 + rows = [i for i in range(n_exp)] + cols = (m + 3*np.arange(n_exp) + 1).tolist() + vals = [1.0] * n_exp task.putaijlist(rows, cols, vals) - # Linear equations: A[ell,:] @ x - t[3*ell + 2] - z[posynom idx corresponding to ell] == -np.log(c[ell]) - cur_con_idx = n - rows = [cur_con_idx + r for r in A.row] - task.putaijlist(rows, A.col, A.data) - rows = [cur_con_idx + i for i in range(n)] - cols = (m + 3*np.arange(n) + 2).tolist() - vals = [-1.0] * n + # 2nd n_exp: Linear equations: + # A[ell,:] @ x - t[3*ell + 2] - z[posynom idx corresponding to ell] == -np.log(c[ell]) + cur_con_idx = n_exp + rows = [cur_con_idx + r for r in A_exp.row] + task.putaijlist(rows, A_exp.col, A_exp.data) + rows = [cur_con_idx + i for i in range(n_exp)] + cols = (m + 3*np.arange(n_exp) + 2).tolist() + vals = [-1.0] * n_exp task.putaijlist(rows, cols, vals) - rows = [cur_con_idx + i for i in range(n)] - cols = [m + 3*n + p_idxs[i] for i in range(n)] - vals = [-1.0] * n + rows = [cur_con_idx + i for i in range(n_exp)] + cols = [m + 3*n_exp + exp_p_idx[i] for i in range(n_exp)] + vals = [-1.0] * n_exp task.putaijlist(rows, cols, vals) - # Linear inequalities: 1 @ t[3 * sels] <= 1, for sels = np.nonzero(p_idxs[i] == i) - cur_con_idx = 2*n + # last p_exp: Linear inequalities: + # 1 @ t[3 * sels] <= 1, for sels = np.nonzero(exp_p_idxs[i] == i) + cur_con_idx = 2*n_exp rows, cols, vals = [], [], [] - for i in range(p+1): - sels = np.nonzero(p_idxs == i)[0] + for i in range(p_exp): + sels = np.nonzero(exp_p_idx == i)[0] rows.extend([cur_con_idx] * sels.size) cols.extend(m + 3 * sels) vals.extend([1] * sels.size) cur_con_idx += 1 task.putaijlist(rows, cols, vals) + # Build the right-hand-sides of the [in]equality constraints + type_constraint = [mosek.boundkey.fx] * (2*n_exp) + [mosek.boundkey.up] * p_exp + h = np.concatenate([np.ones(n_exp), -np.log(c_exp), np.ones(p_exp)]) + task.putconboundlist(np.arange(h.size, dtype=int), type_constraint, h, h) # - # Build the right-hand-sides of the [in]equality constraints + # Affine constraints, not needing the exponential cone # - type_constraint = [mosek.boundkey.fx] * (2*n) + [mosek.boundkey.up] * (p + 1) - h = np.concatenate([np.ones(n), -np.log(c), np.ones(p + 1)]) - task.putconboundlist(np.arange(h.size, dtype=int), type_constraint, h, h) + cur_con_idx = 2*n_exp + p_exp + if c_lin.size > 0: + task.appendcons(c_lin.size) + rows = [cur_con_idx + r for r in A_lin.row] + task.putaijlist(rows, A_lin.col, A_lin.data) + type_constraint = [mosek.boundkey.up] * c_lin.size + con_indices = np.arange(cur_con_idx, cur_con_idx + c_lin.size, dtype=int) + h = -np.log(c_lin) + task.putconboundlist(con_indices, type_constraint, h, h) + cur_con_idx += c_lin.size # # Set the objective function # - task.putclist([int(m + 3*n)], [1]) + task.putclist([int(m + 3*n_exp)], [1]) task.putobjsense(mosek.objsense.minimize) # # Set solver parameters, and call .solve(). # + verbose = False + if 'verbose' in kwargs: + verbose = kwargs['verbose'] + if verbose: + import sys + + def streamprinter(text): + sys.stdout.write(text) + sys.stdout.flush() + + print('\n') + env.set_Stream(mosek.streamtype.log, streamprinter) + task.set_Stream(mosek.streamtype.log, streamprinter) + task.putintparam(mosek.iparam.infeas_report_auto, mosek.onoffkey.on) + task.putintparam(mosek.iparam.log_presolve, 0) + task.optimize() + + if verbose: + task.solutionsummary(mosek.streamtype.msg) # # Recover the solution # @@ -131,12 +192,25 @@ def mskoptimize(c, A, k, p_idxs, *args, **kwargs): x = np.array(x) # recover dual variables for log-sum-exp epigraph constraints # (skip epigraph of the objective function). - z_duals = [0.] * p - task.getsuxslice(mosek.soltype.itr, m + 3*n + 1, n_msk, z_duals) + z_duals = [0.] * (p_exp - 1) + task.getsuxslice(mosek.soltype.itr, m + 3*n_exp + 1, n_msk, z_duals) z_duals = np.array(z_duals) z_duals[z_duals < 0] = 0 + # recover dual variables for the remaining user-provided constraints + if c_lin.size > 0: + aff_duals = [0.] * c_lin.size + task.getsucslice(mosek.soltype.itr, 2*n_exp + p_exp, 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[exp_posys[1:]] = z_duals + merged_duals[lin_posys] = aff_duals + merged_duals = merged_duals[1:] + else: + merged_duals = z_duals # wrap things up in a dictionary - solution = {'status': 'optimal', 'primal': x, 'la': z_duals} + solution = {'status': 'optimal', 'primal': x, 'la': merged_duals} elif msk_solsta == mosek.solsta.prim_infeas_cer: solution = {'status': 'infeasible', 'primal': None} elif msk_solsta == mosek.solsta.dual_infeas_cer: From 8169e362549c372ea1dbd4729878bf1c22bba409 Mon Sep 17 00:00:00 2001 From: rileyjmurray Date: Mon, 20 Jan 2020 21:18:25 -0800 Subject: [PATCH 17/21] test precision --- gpkit/tests/t_model.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/gpkit/tests/t_model.py b/gpkit/tests/t_model.py index efc567373..da5035c21 100644 --- a/gpkit/tests/t_model.py +++ b/gpkit/tests/t_model.py @@ -64,7 +64,7 @@ def test_sigeq(self): SignomialEquality(x**2 + x, y)]) sol = m.localsolve(solver=self.solver, verbosity=0, mutategp=False) self.assertAlmostEqual(sol("x"), 0.1639472, self.ndig) - self.assertAlmostEqual(sol("y")[0], 0.1908254, 3) + self.assertAlmostEqual(sol("y")[0], 0.1908254, self.ndig) self.assertAlmostEqual(sol("c"), 0.2669448, self.ndig) # test right vector input to sigeq with SignomialsEnabled(): @@ -80,9 +80,9 @@ def test_sigeq(self): m = Model(c, [c >= (x + 0.25)**2 + (y - 0.5)**2, SignomialEquality(x**2 + x, y)]) sol = m.localsolve(solver=self.solver, verbosity=0) - self.assertAlmostEqual(sol("x"), 0.1639472, 3) - self.assertAlmostEqual(sol("y"), 0.1908254, 3) - self.assertAlmostEqual(sol("c"), 0.2669448, 3) + self.assertAlmostEqual(sol("x"), 0.1639472, self.ndig) + self.assertAlmostEqual(sol("y"), 0.1908254, self.ndig) + self.assertAlmostEqual(sol("c"), 0.2669448, self.ndig) def test_601(self): # tautological monomials should solve but not pass to the solver From a302f6ed06380d5755cf8b335c4346ff86a46f44 Mon Sep 17 00:00:00 2001 From: rileyjmurray Date: Mon, 20 Jan 2020 22:32:17 -0800 Subject: [PATCH 18/21] improved documentation for mskoptimize in mosek_conif.py --- gpkit/_mosek/mosek_conif.py | 211 +++++++++++++++++++----------------- 1 file changed, 111 insertions(+), 100 deletions(-) diff --git a/gpkit/_mosek/mosek_conif.py b/gpkit/_mosek/mosek_conif.py index ab4863f41..69b48ac61 100644 --- a/gpkit/_mosek/mosek_conif.py +++ b/gpkit/_mosek/mosek_conif.py @@ -21,7 +21,7 @@ def mskoptimize(c, A, k, p_idxs, *args, **kwargs): k : ints array of shape p+1 k[0] is the number of monomials (rows of A) present in the objective k[1:] is the number of monomials (rows of A) present in each constraint - p_idxs : list of bool arrays, each of shape m + p_idxs : ints array of shape n. sel = p_idxs == i selects rows of A and entries of c of the i-th posynomial fi(x) = c[sel] @ exp(A[sel,:] @ x). The 0-th posynomial gives the objective function, and the remaining posynomials should be constrained to be <= 1. @@ -30,133 +30,144 @@ def mskoptimize(c, A, k, p_idxs, *args, **kwargs): ------- dict Contains the following keys - "success": bool - "objective_sol" float - Optimal value of the objective - "primal_sol": floats array of size m - Optimal value of free variables. Note: not in logspace. - "dual_sol": floats array of size p - Optimal value of the dual variables, in logspace. + "status": string + "optimal", "infeasible", "unbounded", or "unknown". + "primal" np.ndarray or None + The values of the ``m`` primal variables. + "la": np.ndarray or None + The dual variables to the ``p`` posynomial constraints, when + those constraints are represented in log-sum-exp ("LSE") form. """ import mosek # - # Prepare some initial problem data + # Initial transformations of problem data. # - # Say that the optimization variable is y = [x;t;z], where - # x is of length m, t is of length 3*n, and z is of length p+1. + # separate monomial constraints (call them "lin"), from those which require + # an LSE representation (call those "lse"). # - # Triples - # t[3*ell] - # t[3*ell + 1] should be == 1. - # t[3*ell + 2] should be == A[ell, :] @ x + np.log(c[ell]) - z[posynom idx corresponding to ell]. - # should belong to the exponential cone. + # NOTE: the epigraph of the objective function always gets an "lse" + # representation, even if the objective is a monomial. # - # For each i from {0,...,p}, the "t" should also satisfy - # sum(t[3*ell] for ell correponsing to posynomial i) <= 1. - # - # The vector "z" should satisfy - # z[1:] <= 0 - # - # - c = np.array(c) - # separate exponential from linear constraints. - # the objective is always kept in exponential form. - exp_posys = [0] + [i+1 for i, val in enumerate(k[1:]) if val > 1] - lin_posys = [i for i in range(len(k)) if i not in exp_posys] + log_c = np.log(np.array(c)) + lse_posys = [0] + [i+1 for i, val in enumerate(k[1:]) if val > 1] + lin_posys = [i for i in range(len(k)) if i not in lse_posys] if len(lin_posys) > 0: A = A.tocsr() - lin_idxs = [] - for i in lin_posys: - temp = np.nonzero(p_idxs == i) - temp = temp[0] - lin_idxs.append(temp) - lin_idxs = np.concatenate(lin_idxs) - nonlin_idxs = np.ones(A.shape[0], dtype=bool) - nonlin_idxs[lin_idxs] = False - A_exp = A[nonlin_idxs, :].tocoo() - c_exp = c[nonlin_idxs] + lin_idxs = np.concatenate([np.nonzero(p_idxs == i)[0] for i in lin_posys]) + lse_idxs = np.ones(A.shape[0], dtype=bool) + lse_idxs[lin_idxs] = False + A_lse = A[lse_idxs, :].tocoo() + log_c_lse = log_c[lse_idxs] A_lin = A[lin_idxs, :].tocoo() - c_lin = c[lin_idxs] + log_c_lin = log_c[lin_idxs] else: - c_lin = np.array([]) - A_exp = A - c_exp = c - - m = A.shape[1] - k_exp = [k[i] for i in exp_posys] - n_exp = sum(k_exp) - p_exp = len(k_exp) - n_msk = m + 3*n_exp + p_exp - exp_p_idx = [] - for i, ki in enumerate(k_exp): - exp_p_idx.extend([i] * ki) - exp_p_idx = np.array(exp_p_idx) + log_c_lin = np.array([]) # A_lin won't be referenced later, so no need to define it. + A_lse = A + log_c_lse = log_c + k_lse = [k[i] for i in lse_posys] + n_lse = sum(k_lse) + p_lse = len(k_lse) + lse_p_idx = [] + for i, ki in enumerate(k_lse): + lse_p_idx.extend([i] * ki) + lse_p_idx = np.array(lse_p_idx) + # + # Create MOSEK task. Add variables, and conic constraints. + # + # Say that MOSEK's optimization variable is a block vector, [x;t;z], where ... + # x is the user-defined primal variable (length m), + # t is an auxiliary variable for exponential cones (length 3 * n_lse), and + # z is an epigraph variable for LSE terms (length p_lse). + # + # The variable z[0] is special, because it's the epigraph of the objective function + # in LSE form. The sign of this variable is not constrained. + # + # The variables z[1:] are epigraph terms for "log", in constraints that naturally + # write as LSE(Ai @ x + log_ci) <= 0. These variables need to be <= 0. # - # Create MOSEK task. add variables and conic constraints. + # The main constraints on (x, t, z) are described in next comment block. # env = mosek.Env() task = env.Task(0, 0) - task.appendvars(n_msk) - bound_types = [mosek.boundkey.fr] * (m + 3*n_exp + 1) + [mosek.boundkey.up] * (p_exp - 1) - task.putvarboundlist(np.arange(n_msk, dtype=int), - bound_types, np.zeros(n_msk), np.zeros(n_msk)) - for i in range(n_exp): + m = A.shape[1] + msk_nvars = m + 3 * n_lse + p_lse + task.appendvars(msk_nvars) + bound_types = [mosek.boundkey.fr] * (m + 3*n_lse + 1) + [mosek.boundkey.up] * (p_lse - 1) + task.putvarboundlist(np.arange(msk_nvars, dtype=int), + bound_types, np.zeros(msk_nvars), np.zeros(msk_nvars)) + for i in range(n_lse): idx = m + 3*i task.appendcone(mosek.conetype.pexp, 0.0, np.arange(idx, idx + 3)) # # Affine constraints related to the exponential cone # - task.appendcons(2*n_exp + p_exp) - # 1st n_exp: Linear equations: t[3*ell + 1] == 1 - rows = [i for i in range(n_exp)] - cols = (m + 3*np.arange(n_exp) + 1).tolist() - vals = [1.0] * n_exp - task.putaijlist(rows, cols, vals) - # 2nd n_exp: Linear equations: - # A[ell,:] @ x - t[3*ell + 2] - z[posynom idx corresponding to ell] == -np.log(c[ell]) - cur_con_idx = n_exp - rows = [cur_con_idx + r for r in A_exp.row] - task.putaijlist(rows, A_exp.col, A_exp.data) - rows = [cur_con_idx + i for i in range(n_exp)] - cols = (m + 3*np.arange(n_exp) + 2).tolist() - vals = [-1.0] * n_exp - task.putaijlist(rows, cols, vals) - rows = [cur_con_idx + i for i in range(n_exp)] - cols = [m + 3*n_exp + exp_p_idx[i] for i in range(n_exp)] - vals = [-1.0] * n_exp + # For each i in {0, ..., n_lse - 1}, we need + # t[3*i + 1] == 1, and + # t[3*i + 2] == A_lse[i, :] @ x + log_c_lse[i] - z[lse_p_idx[i]]. + # This contributes 2 * n_lse constraints. + # + # For each j from {0, ..., p_lse - 1}, the "t" should also satisfy + # sum(t[3*i] for i where i == lse_p_idx[j]) <= 1. + # This contributes another p_lse constraints. + # + # The above constraints imply that for ``sel = lse_p_idx == i``, we have + # LSE(A_lse[sel, :] @ x + log_c_lse[sel]) <= z[i]. + # + # We specify the necessary constraints to MOSEK in three phases. Over the + # course of these three phases, we make a total of five calls to "putaijlist" + # and a single call to "putconboundlist". + # + task.appendcons(2*n_lse + p_lse) + # 1st n_lse: Linear equations: t[3*i + 1] == 1 + rows = [i for i in range(n_lse)] + cols = (m + 3*np.arange(n_lse) + 1).tolist() + vals = [1.0] * n_lse task.putaijlist(rows, cols, vals) - # last p_exp: Linear inequalities: - # 1 @ t[3 * sels] <= 1, for sels = np.nonzero(exp_p_idxs[i] == i) - cur_con_idx = 2*n_exp + cur_con_idx = n_lse + # 2nd n_lse: Linear equations between (x,t,z). + rows = [cur_con_idx + r for r in A_lse.row] + task.putaijlist(rows, A_lse.col, A_lse.data) # coefficients on "x" + rows = [cur_con_idx + i for i in range(n_lse)] + cols = (m + 3*np.arange(n_lse) + 2).tolist() + vals = [-1.0] * n_lse + task.putaijlist(rows, cols, vals) # coefficients on "t" + rows = [cur_con_idx + i for i in range(n_lse)] + cols = [m + 3*n_lse + lse_p_idx[i] for i in range(n_lse)] + vals = [-1.0] * n_lse + task.putaijlist(rows, cols, vals) # coefficients on "z". + cur_con_idx = 2 * n_lse + # last p_lse: Linear inequalities on certain sums of "t". rows, cols, vals = [], [], [] - for i in range(p_exp): - sels = np.nonzero(exp_p_idx == i)[0] + for i in range(p_lse): + sels = np.nonzero(lse_p_idx == i)[0] rows.extend([cur_con_idx] * sels.size) cols.extend(m + 3 * sels) vals.extend([1] * sels.size) cur_con_idx += 1 task.putaijlist(rows, cols, vals) + cur_con_idx = 2 * n_lse + p_lse # Build the right-hand-sides of the [in]equality constraints - type_constraint = [mosek.boundkey.fx] * (2*n_exp) + [mosek.boundkey.up] * p_exp - h = np.concatenate([np.ones(n_exp), -np.log(c_exp), np.ones(p_exp)]) + type_constraint = [mosek.boundkey.fx] * (2*n_lse) + [mosek.boundkey.up] * p_lse + h = np.concatenate([np.ones(n_lse), -log_c_lse, np.ones(p_lse)]) task.putconboundlist(np.arange(h.size, dtype=int), type_constraint, h, h) # # Affine constraints, not needing the exponential cone # - cur_con_idx = 2*n_exp + p_exp - if c_lin.size > 0: - task.appendcons(c_lin.size) + # Require A_lin @ x <= -log_c_lin. + # + if log_c_lin.size > 0: + task.appendcons(log_c_lin.size) rows = [cur_con_idx + r for r in A_lin.row] task.putaijlist(rows, A_lin.col, A_lin.data) - type_constraint = [mosek.boundkey.up] * c_lin.size - con_indices = np.arange(cur_con_idx, cur_con_idx + c_lin.size, dtype=int) - h = -np.log(c_lin) + type_constraint = [mosek.boundkey.up] * log_c_lin.size + con_indices = np.arange(cur_con_idx, cur_con_idx + log_c_lin.size, dtype=int) + h = -log_c_lin task.putconboundlist(con_indices, type_constraint, h, h) - cur_con_idx += c_lin.size + cur_con_idx += log_c_lin.size # # Set the objective function # - task.putclist([int(m + 3*n_exp)], [1]) + task.putclist([int(m + 3*n_lse)], [1]) task.putobjsense(mosek.objsense.minimize) # # Set solver parameters, and call .solve(). @@ -192,19 +203,19 @@ def streamprinter(text): x = np.array(x) # recover dual variables for log-sum-exp epigraph constraints # (skip epigraph of the objective function). - z_duals = [0.] * (p_exp - 1) - task.getsuxslice(mosek.soltype.itr, m + 3*n_exp + 1, n_msk, z_duals) + 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 c_lin.size > 0: - aff_duals = [0.] * c_lin.size - task.getsucslice(mosek.soltype.itr, 2*n_exp + p_exp, cur_con_idx, aff_duals) + if log_c_lin.size > 0: + aff_duals = [0.] * log_c_lin.size + task.getsucslice(mosek.soltype.itr, 2*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[exp_posys[1:]] = z_duals + merged_duals[lse_posys[1:]] = z_duals merged_duals[lin_posys] = aff_duals merged_duals = merged_duals[1:] else: @@ -212,11 +223,11 @@ def streamprinter(text): # wrap things up in a dictionary solution = {'status': 'optimal', 'primal': x, 'la': merged_duals} elif msk_solsta == mosek.solsta.prim_infeas_cer: - solution = {'status': 'infeasible', 'primal': None} + solution = {'status': 'infeasible', 'primal': None, 'la': None} elif msk_solsta == mosek.solsta.dual_infeas_cer: - solution = {'status': 'unbounded', 'primal': None} + solution = {'status': 'unbounded', 'primal': None, 'la': None} else: - solution = {'status': 'unknown', 'primal': None} + solution = {'status': 'unknown', 'primal': None, 'la': None} task.__exit__(None, None, None) env.__exit__(None, None, None) return solution From 17be8784cdd9c0d95096a42de653be4d6e0fac77 Mon Sep 17 00:00:00 2001 From: Riley Murray Date: Tue, 21 Jan 2020 15:23:21 -0800 Subject: [PATCH 19/21] mosek_conif: check if an installed mosek python package is associated with MOSEK version >= 9. --- gpkit/_mosek/mosek_conif.py | 6 ++++++ gpkit/build.py | 7 +++++-- 2 files changed, 11 insertions(+), 2 deletions(-) diff --git a/gpkit/_mosek/mosek_conif.py b/gpkit/_mosek/mosek_conif.py index 69b48ac61..04dce1ab5 100644 --- a/gpkit/_mosek/mosek_conif.py +++ b/gpkit/_mosek/mosek_conif.py @@ -39,6 +39,12 @@ def mskoptimize(c, A, k, p_idxs, *args, **kwargs): those constraints are represented in log-sum-exp ("LSE") form. """ import mosek + if not hasattr(mosek.conetype, 'pexp'): + msg = """ + mosek_conif requires MOSEK version >= 9. The imported version of MOSEK does + not have attribute ``mosek.conetype.pexp``, which was introduced in MOSEK 9. + """ + raise RuntimeError(msg) # # Initial transformations of problem data. # diff --git a/gpkit/build.py b/gpkit/build.py index b910c5dad..4a7298fa7 100644 --- a/gpkit/build.py +++ b/gpkit/build.py @@ -127,12 +127,15 @@ class MosekConif(SolverBackend): name = 'mosek_conif' def look(self): - "Attempts to import mosek." + "Attempts to import mosek, version >= 9." try: log("# Trying to import mosek...") # Testing the import, so the variable is intentionally not used import mosek # pylint: disable=unused-variable - return "in Python path" + if hasattr(mosek.conetype, 'pexp'): + return "in Python path" + else: + pass except ImportError: pass From 1d79b07dc8a4669295f6eecf25fe89280f923ce8 Mon Sep 17 00:00:00 2001 From: Riley Murray Date: Tue, 21 Jan 2020 15:51:47 -0800 Subject: [PATCH 20/21] pylint --- docs/source/examples/autosweep.py | 3 +- gpkit/_mosek/mosek_conif.py | 73 +++++++++++++++++++------------ gpkit/build.py | 1 + gpkit/constraints/gp.py | 3 +- gpkit/tests/t_model.py | 2 +- 5 files changed, 50 insertions(+), 32 deletions(-) diff --git a/docs/source/examples/autosweep.py b/docs/source/examples/autosweep.py index 11561284e..2519db615 100644 --- a/docs/source/examples/autosweep.py +++ b/docs/source/examples/autosweep.py @@ -36,7 +36,8 @@ # this problem is two intersecting lines in logspace m2 = Model(A**2, [A >= (l/3)**2, A >= (l/3)**0.5 * units.m**1.5]) tol2 = {"mosek": 1e-12, "cvxopt": 1e-7, - "mosek_cli": 1e-6, 'mosek_conif': 1e-6}[gpkit.settings["default_solver"]] + "mosek_cli": 1e-6, + 'mosek_conif': 1e-6}[gpkit.settings["default_solver"]] # test Model method sol2 = m2.autosweep({l: [1, 10]}, tol2, verbosity=0) bst2 = sol2.bst diff --git a/gpkit/_mosek/mosek_conif.py b/gpkit/_mosek/mosek_conif.py index 04dce1ab5..9a3a75cf5 100644 --- a/gpkit/_mosek/mosek_conif.py +++ b/gpkit/_mosek/mosek_conif.py @@ -1,5 +1,6 @@ -"Implements the GPkit interface to MOSEK (version >= 9) python-based Optimizer API" - +"Implements the GPkit interface to MOSEK (v >= 9) python-based Optimizer API" +from __future__ import print_function +import sys import numpy as np @@ -22,9 +23,10 @@ def mskoptimize(c, A, k, p_idxs, *args, **kwargs): k[0] is the number of monomials (rows of A) present in the objective k[1:] is the number of monomials (rows of A) present in each constraint p_idxs : ints array of shape n. - sel = p_idxs == i selects rows of A and entries of c of the i-th posynomial - fi(x) = c[sel] @ exp(A[sel,:] @ x). The 0-th posynomial gives the objective - function, and the remaining posynomials should be constrained to be <= 1. + sel = p_idxs == i selects rows of A and entries of c of the i-th + posynomial fi(x) = c[sel] @ exp(A[sel,:] @ x). The 0-th posynomial gives + the objective function, and the remaining posynomials should be + constrained to be <= 1. Returns ------- @@ -41,15 +43,16 @@ def mskoptimize(c, A, k, p_idxs, *args, **kwargs): import mosek if not hasattr(mosek.conetype, 'pexp'): msg = """ - mosek_conif requires MOSEK version >= 9. The imported version of MOSEK does - not have attribute ``mosek.conetype.pexp``, which was introduced in MOSEK 9. + mosek_conif requires MOSEK version >= 9. The imported version of MOSEK + does not have attribute ``mosek.conetype.pexp``, which was introduced + in MOSEK 9. """ raise RuntimeError(msg) # # Initial transformations of problem data. # - # separate monomial constraints (call them "lin"), from those which require - # an LSE representation (call those "lse"). + # separate monomial constraints (call them "lin"), from those which + # require an LSE representation (call those "lse"). # # NOTE: the epigraph of the objective function always gets an "lse" # representation, even if the objective is a monomial. @@ -57,9 +60,10 @@ def mskoptimize(c, A, k, p_idxs, *args, **kwargs): log_c = np.log(np.array(c)) lse_posys = [0] + [i+1 for i, val in enumerate(k[1:]) if val > 1] lin_posys = [i for i in range(len(k)) if i not in lse_posys] - if len(lin_posys) > 0: + if lin_posys: A = A.tocsr() - lin_idxs = np.concatenate([np.nonzero(p_idxs == i)[0] for i in lin_posys]) + lin_idxs = np.concatenate( + [np.nonzero(p_idxs == i)[0] for i in lin_posys]) lse_idxs = np.ones(A.shape[0], dtype=bool) lse_idxs[lin_idxs] = False A_lse = A[lse_idxs, :].tocoo() @@ -67,7 +71,8 @@ def mskoptimize(c, A, k, p_idxs, *args, **kwargs): A_lin = A[lin_idxs, :].tocoo() log_c_lin = log_c[lin_idxs] else: - log_c_lin = np.array([]) # A_lin won't be referenced later, so no need to define it. + log_c_lin = np.array([]) + # A_lin won't be referenced later, so no need to define it. A_lse = A log_c_lse = log_c k_lse = [k[i] for i in lse_posys] @@ -80,25 +85,31 @@ def mskoptimize(c, A, k, p_idxs, *args, **kwargs): # # Create MOSEK task. Add variables, and conic constraints. # - # Say that MOSEK's optimization variable is a block vector, [x;t;z], where ... + # Say that MOSEK's optimization variable is a block vector, [x;t;z], + # where ... # x is the user-defined primal variable (length m), - # t is an auxiliary variable for exponential cones (length 3 * n_lse), and + # t is an auxiliary variable for exponential cones (length + # 3 * n_lse), and # z is an epigraph variable for LSE terms (length p_lse). # - # The variable z[0] is special, because it's the epigraph of the objective function - # in LSE form. The sign of this variable is not constrained. + # The variable z[0] is special, because it's the epigraph of the + # objective function in LSE form. The sign of this variable is + # not constrained. # - # The variables z[1:] are epigraph terms for "log", in constraints that naturally - # write as LSE(Ai @ x + log_ci) <= 0. These variables need to be <= 0. + # The variables z[1:] are epigraph terms for "log", in constraints + # that naturally write as LSE(Ai @ x + log_ci) <= 0. These variables + # need to be <= 0. # - # The main constraints on (x, t, z) are described in next comment block. + # The main constraints on (x, t, z) are described in next + # comment block. # env = mosek.Env() task = env.Task(0, 0) m = A.shape[1] msk_nvars = m + 3 * n_lse + p_lse task.appendvars(msk_nvars) - bound_types = [mosek.boundkey.fr] * (m + 3*n_lse + 1) + [mosek.boundkey.up] * (p_lse - 1) + bound_types = [mosek.boundkey.fr] * (m + 3*n_lse + 1) + bound_types.extend([mosek.boundkey.up] * (p_lse - 1)) task.putvarboundlist(np.arange(msk_nvars, dtype=int), bound_types, np.zeros(msk_nvars), np.zeros(msk_nvars)) for i in range(n_lse): @@ -116,16 +127,17 @@ def mskoptimize(c, A, k, p_idxs, *args, **kwargs): # sum(t[3*i] for i where i == lse_p_idx[j]) <= 1. # This contributes another p_lse constraints. # - # The above constraints imply that for ``sel = lse_p_idx == i``, we have + # The above constraints imply that for ``sel = lse_p_idx == i``, + # we have # LSE(A_lse[sel, :] @ x + log_c_lse[sel]) <= z[i]. # - # We specify the necessary constraints to MOSEK in three phases. Over the - # course of these three phases, we make a total of five calls to "putaijlist" - # and a single call to "putconboundlist". + # We specify the necessary constraints to MOSEK in three phases. + # Over the course of these three phases, we make a total of five + # calls to "putaijlist" and a single call to "putconboundlist". # task.appendcons(2*n_lse + p_lse) # 1st n_lse: Linear equations: t[3*i + 1] == 1 - rows = [i for i in range(n_lse)] + rows = list(range(n_lse)) cols = (m + 3*np.arange(n_lse) + 1).tolist() vals = [1.0] * n_lse task.putaijlist(rows, cols, vals) @@ -153,7 +165,8 @@ def mskoptimize(c, A, k, p_idxs, *args, **kwargs): task.putaijlist(rows, cols, vals) cur_con_idx = 2 * n_lse + p_lse # Build the right-hand-sides of the [in]equality constraints - type_constraint = [mosek.boundkey.fx] * (2*n_lse) + [mosek.boundkey.up] * p_lse + type_constraint = [mosek.boundkey.fx] * (2*n_lse) + type_constraint.extend([mosek.boundkey.up] * p_lse) h = np.concatenate([np.ones(n_lse), -log_c_lse, np.ones(p_lse)]) task.putconboundlist(np.arange(h.size, dtype=int), type_constraint, h, h) # @@ -166,7 +179,8 @@ def mskoptimize(c, A, k, p_idxs, *args, **kwargs): rows = [cur_con_idx + r for r in A_lin.row] task.putaijlist(rows, A_lin.col, A_lin.data) type_constraint = [mosek.boundkey.up] * log_c_lin.size - con_indices = np.arange(cur_con_idx, cur_con_idx + log_c_lin.size, dtype=int) + con_indices = np.arange(cur_con_idx, cur_con_idx + log_c_lin.size, + dtype=int) h = -log_c_lin task.putconboundlist(con_indices, type_constraint, h, h) cur_con_idx += log_c_lin.size @@ -182,9 +196,9 @@ def mskoptimize(c, A, k, p_idxs, *args, **kwargs): if 'verbose' in kwargs: verbose = kwargs['verbose'] if verbose: - import sys def streamprinter(text): + """A helper, for logging MOSEK output to sys.stdout.""" sys.stdout.write(text) sys.stdout.flush() @@ -216,7 +230,8 @@ def streamprinter(text): # recover dual variables for the remaining user-provided constraints if log_c_lin.size > 0: aff_duals = [0.] * log_c_lin.size - task.getsucslice(mosek.soltype.itr, 2*n_lse + p_lse, cur_con_idx, aff_duals) + task.getsucslice(mosek.soltype.itr, 2*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 diff --git a/gpkit/build.py b/gpkit/build.py index 4a7298fa7..8de5bf36c 100644 --- a/gpkit/build.py +++ b/gpkit/build.py @@ -123,6 +123,7 @@ def look(self): class MosekConif(SolverBackend): + "Find MOSEK version >= 9." name = 'mosek_conif' diff --git a/gpkit/constraints/gp.py b/gpkit/constraints/gp.py index 1b103413c..0b8eaa1c8 100644 --- a/gpkit/constraints/gp.py +++ b/gpkit/constraints/gp.py @@ -15,7 +15,8 @@ DEFAULT_SOLVER_KWARGS = {"cvxopt": {"kktsolver": "ldl"}} -SOLUTION_TOL = {"cvxopt": 1e-3, "mosek_cli": 1e-4, "mosek": 1e-5, 'mosek_conif': 1e-3} +SOLUTION_TOL = {"cvxopt": 1e-3, "mosek_cli": 1e-4, "mosek": 1e-5, + 'mosek_conif': 1e-3} def _get_solver(solver, kwargs): diff --git a/gpkit/tests/t_model.py b/gpkit/tests/t_model.py index da5035c21..0c9eca34c 100644 --- a/gpkit/tests/t_model.py +++ b/gpkit/tests/t_model.py @@ -319,7 +319,7 @@ def test_sp_bounded(self): with SignomialsEnabled(): m = Model(x, [x + y >= 1]) # dual infeasible with self.assertRaises((RuntimeWarning, ValueError)): - m.localsolve(verbosity=0, solver=self.solver), + m.localsolve(verbosity=0, solver=self.solver) with SignomialsEnabled(): m = Model(x, Bounded([x + y >= 1], verbosity=0)) From 3d7ce7179319aec601f87cbd78b0511ae4045a32 Mon Sep 17 00:00:00 2001 From: rileyjmurray Date: Tue, 21 Jan 2020 22:22:55 -0800 Subject: [PATCH 21/21] simplified implementation, by taking cues from bqpd's mosek9 branch (specifically commit 1e6791363a4fc2955ccc8b23085230b4d83352ec) --- gpkit/_mosek/mosek_conif.py | 88 ++++++++++++++++++------------------- 1 file changed, 42 insertions(+), 46 deletions(-) diff --git a/gpkit/_mosek/mosek_conif.py b/gpkit/_mosek/mosek_conif.py index 9a3a75cf5..2db6e717f 100644 --- a/gpkit/_mosek/mosek_conif.py +++ b/gpkit/_mosek/mosek_conif.py @@ -83,7 +83,8 @@ def mskoptimize(c, A, k, p_idxs, *args, **kwargs): lse_p_idx.extend([i] * ki) lse_p_idx = np.array(lse_p_idx) # - # Create MOSEK task. Add variables, and conic constraints. + # Create MOSEK task. Add variables, bound constraints, and conic + # constraints. # # Say that MOSEK's optimization variable is a block vector, [x;t;z], # where ... @@ -108,53 +109,51 @@ def mskoptimize(c, A, k, p_idxs, *args, **kwargs): m = A.shape[1] msk_nvars = m + 3 * n_lse + p_lse task.appendvars(msk_nvars) - bound_types = [mosek.boundkey.fr] * (m + 3*n_lse + 1) - bound_types.extend([mosek.boundkey.up] * (p_lse - 1)) - task.putvarboundlist(np.arange(msk_nvars, dtype=int), - bound_types, np.zeros(msk_nvars), np.zeros(msk_nvars)) - for i in range(n_lse): - idx = m + 3*i - task.appendcone(mosek.conetype.pexp, 0.0, np.arange(idx, idx + 3)) - # - # Affine constraints related to the exponential cone + # "x" is free + task.putvarboundlist(np.arange(m), [mosek.boundkey.fr] * m, + np.zeros(m), np.zeros(m)) + # t[3 * i + i] == 1, other components are free. + bound_types = [mosek.boundkey.fr, mosek.boundkey.fx, mosek.boundkey.fr] + task.putvarboundlist(np.arange(m, m + 3*n_lse), bound_types * n_lse, + np.ones(3*n_lse), np.ones(3*n_lse)) + # z[0] is free; z[1:] <= 0. + bound_types = [mosek.boundkey.fr] + [mosek.boundkey.up] * (p_lse - 1) + task.putvarboundlist(np.arange(m + 3*n_lse, msk_nvars), bound_types, + np.zeros(p_lse), np.zeros(p_lse)) + # t[3*i], t[3*i + 1], t[3*i + 2] belongs to the exponential cone + task.appendconesseq([mosek.conetype.pexp] * n_lse, [0.0] * n_lse, + [3] * n_lse, m) + # + # Exponential cone affine constraints (other than t[3*i + 1] == 1). # # For each i in {0, ..., n_lse - 1}, we need - # t[3*i + 1] == 1, and # t[3*i + 2] == A_lse[i, :] @ x + log_c_lse[i] - z[lse_p_idx[i]]. - # This contributes 2 * n_lse constraints. # # For each j from {0, ..., p_lse - 1}, the "t" should also satisfy # sum(t[3*i] for i where i == lse_p_idx[j]) <= 1. - # This contributes another p_lse constraints. - # - # The above constraints imply that for ``sel = lse_p_idx == i``, - # we have - # LSE(A_lse[sel, :] @ x + log_c_lse[sel]) <= z[i]. - # - # We specify the necessary constraints to MOSEK in three phases. - # Over the course of these three phases, we make a total of five - # calls to "putaijlist" and a single call to "putconboundlist". # - task.appendcons(2*n_lse + p_lse) - # 1st n_lse: Linear equations: t[3*i + 1] == 1 - rows = list(range(n_lse)) - cols = (m + 3*np.arange(n_lse) + 1).tolist() - vals = [1.0] * n_lse + # When combined with bound constraints ``t[3*i + 1] == 1``, the + # above constraints imply + # LSE(A_lse[sel, :] @ x + log_c_lse[sel]) <= z[i] + # for ``sel = lse_p_idx == i``. + # + task.appendcons(n_lse + p_lse) + # Linear equations between (x,t,z). + # start with coefficients on "x" + rows = [r for r in A_lse.row] + cols = [c for c in A_lse.col] + vals = [v for v in A_lse.data] + # add coefficients on "t" + rows += list(range(n_lse)) + cols += (m + 3*np.arange(n_lse) + 2).tolist() + vals += [-1.0] * n_lse + # add coefficients on "z" + rows += list(range(n_lse)) + cols += [m + 3*n_lse + lse_p_idx[i] for i in range(n_lse)] + vals += [-1.0] * n_lse task.putaijlist(rows, cols, vals) cur_con_idx = n_lse - # 2nd n_lse: Linear equations between (x,t,z). - rows = [cur_con_idx + r for r in A_lse.row] - task.putaijlist(rows, A_lse.col, A_lse.data) # coefficients on "x" - rows = [cur_con_idx + i for i in range(n_lse)] - cols = (m + 3*np.arange(n_lse) + 2).tolist() - vals = [-1.0] * n_lse - task.putaijlist(rows, cols, vals) # coefficients on "t" - rows = [cur_con_idx + i for i in range(n_lse)] - cols = [m + 3*n_lse + lse_p_idx[i] for i in range(n_lse)] - vals = [-1.0] * n_lse - task.putaijlist(rows, cols, vals) # coefficients on "z". - cur_con_idx = 2 * n_lse - # last p_lse: Linear inequalities on certain sums of "t". + # Linear inequalities on certain sums of "t". rows, cols, vals = [], [], [] for i in range(p_lse): sels = np.nonzero(lse_p_idx == i)[0] @@ -163,11 +162,9 @@ def mskoptimize(c, A, k, p_idxs, *args, **kwargs): vals.extend([1] * sels.size) cur_con_idx += 1 task.putaijlist(rows, cols, vals) - cur_con_idx = 2 * n_lse + p_lse # Build the right-hand-sides of the [in]equality constraints - type_constraint = [mosek.boundkey.fx] * (2*n_lse) - type_constraint.extend([mosek.boundkey.up] * p_lse) - h = np.concatenate([np.ones(n_lse), -log_c_lse, np.ones(p_lse)]) + type_constraint = [mosek.boundkey.fx] * n_lse + [mosek.boundkey.up] * p_lse + h = np.concatenate([-log_c_lse, np.ones(p_lse)]) task.putconboundlist(np.arange(h.size, dtype=int), type_constraint, h, h) # # Affine constraints, not needing the exponential cone @@ -179,8 +176,7 @@ def mskoptimize(c, A, k, p_idxs, *args, **kwargs): rows = [cur_con_idx + r for r in A_lin.row] task.putaijlist(rows, A_lin.col, A_lin.data) type_constraint = [mosek.boundkey.up] * log_c_lin.size - con_indices = np.arange(cur_con_idx, cur_con_idx + log_c_lin.size, - dtype=int) + con_indices = np.arange(cur_con_idx, cur_con_idx + log_c_lin.size) h = -log_c_lin task.putconboundlist(con_indices, type_constraint, h, h) cur_con_idx += log_c_lin.size @@ -230,7 +226,7 @@ def streamprinter(text): # recover dual variables for the remaining user-provided constraints if log_c_lin.size > 0: aff_duals = [0.] * log_c_lin.size - task.getsucslice(mosek.soltype.itr, 2*n_lse + p_lse, cur_con_idx, + 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