Hot-keys on this page
r m x p toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
1"""Implement the GeometricProgram class"""
2import sys
3import warnings
4from time import time
5from collections import defaultdict
6import numpy as np
7from ..small_classes import CootMatrix, SolverLog, Numbers, FixedScalar
8from ..keydict import KeyDict
9from ..solution_array import SolutionArray
10from .set import ConstraintSet
11from ..exceptions import (InvalidPosynomial, Infeasible, UnknownInfeasible,
12 PrimalInfeasible, DualInfeasible, UnboundedGP,
13 InvalidLicense)
16DEFAULT_SOLVER_KWARGS = {"cvxopt": {"kktsolver": "ldl"}}
17SOLUTION_TOL = {"cvxopt": 1e-3, "mosek_cli": 1e-4, "mosek_conif": 1e-3}
20class MonoEqualityIndexes:
21 "Class to hold MonoEqualityIndexes"
23 def __init__(self):
24 self.all = set()
25 self.first_half = set()
28def _get_solver(solver, kwargs):
29 """Get the solverfn and solvername associated with solver"""
30 if solver is None:
31 from .. import settings
32 try:
33 solver = settings["default_solver"]
34 except KeyError:
35 raise ValueError("No default solver was set during build, so"
36 " solvers must be manually specified.")
37 if solver == "cvxopt":
38 from ..solvers.cvxopt import optimize
39 elif solver == "mosek_cli":
40 from ..solvers.mosek_cli import optimize_generator
41 optimize = optimize_generator(**kwargs)
42 elif solver == "mosek_conif":
43 from ..solvers.mosek_conif import optimize
44 elif hasattr(solver, "__call__"):
45 solver, optimize = solver.__name__, solver
46 else:
47 raise ValueError("Unknown solver '%s'." % solver)
48 return solver, optimize
51class GeometricProgram:
52 # pylint: disable=too-many-instance-attributes
53 """Standard mathematical representation of a GP.
55 Attributes with side effects
56 ----------------------------
57 `solver_out` and `solve_log` are set during a solve
58 `result` is set at the end of a solve if solution status is optimal
60 Examples
61 --------
62 >>> gp = gpkit.geometric_program.GeometricProgram(
63 # minimize
64 x,
65 [ # subject to
66 1/x # <= 1, implicitly
67 ], {})
68 >>> gp.solve()
69 """
70 _result = solve_log = solver_out = model = v_ss = nu_by_posy = None
72 def __init__(self, cost, constraints, substitutions, *, checkbounds=True):
73 self.cost, self.substitutions = cost, substitutions
74 for key, sub in self.substitutions.items():
75 if isinstance(sub, FixedScalar):
76 sub = sub.value
77 if hasattr(sub, "units"):
78 sub = sub.to(key.units or "dimensionless").magnitude
79 self.substitutions[key] = sub
80 if not isinstance(sub, (Numbers, np.ndarray)):
81 raise ValueError("substitution {%s: %s} has invalid value type"
82 " %s." % (key, sub, type(sub)))
83 cost_hmap = cost.hmap.sub(self.substitutions, cost.varkeys)
84 if any(c <= 0 for c in cost_hmap.values()):
85 raise InvalidPosynomial("a GP's cost must be Posynomial")
86 hmapgen = ConstraintSet.as_hmapslt1(constraints, self.substitutions)
87 self.hmaps = [cost_hmap] + list(hmapgen)
88 self.gen() # Generate various maps into the posy- and monomials
89 if checkbounds:
90 self.check_bounds(err_on_missing_bounds=True)
92 def check_bounds(self, *, err_on_missing_bounds=False):
93 "Checks if any variables are unbounded, through equality constraints."
94 missingbounds = {}
95 for var, locs in self.varlocs.items():
96 upperbound, lowerbound = False, False
97 for i in locs:
98 if i not in self.meq_idxs.all:
99 if self.exps[i][var] > 0: # pylint:disable=simplifiable-if-statement
100 upperbound = True
101 else:
102 lowerbound = True
103 if upperbound and lowerbound:
104 break
105 if not upperbound:
106 missingbounds[(var, "upper")] = "."
107 if not lowerbound:
108 missingbounds[(var, "lower")] = "."
109 if not missingbounds:
110 return {} # all bounds found in inequalities
111 meq_bounds = gen_meq_bounds(missingbounds, self.exps, self.meq_idxs)
112 fulfill_meq_bounds(missingbounds, meq_bounds)
113 if missingbounds and err_on_missing_bounds:
114 raise UnboundedGP(
115 "\n\n".join("%s has no %s bound%s" % (v, b, x)
116 for (v, b), x in missingbounds.items()))
117 return missingbounds
119 def gen(self):
120 "Generates nomial and solve data (A, p_idxs) from posynomials"
121 # k [posys]: number of monomials (rows of A) present in each constraint
122 self.k = [len(hmap) for hmap in self.hmaps]
123 # m_idxs [mons]: monomial indices of each posynomial
124 self.m_idxs = []
125 # p_idxs [mons]: posynomial index of each monomial
126 self.p_idxs = []
127 # cs, exps [mons]: coefficient and exponents of each monomial
128 self.cs, self.exps = [], []
129 # varlocs: {vk: monomial indices of each variables' location}
130 self.varkeys = self.varlocs = defaultdict(list)
131 # meq_idxs: {all indices of equality mons} and {just the first halves}
132 self.meq_idxs = MonoEqualityIndexes()
133 m_idx = 0
134 row, col, data = [], [], []
135 for p_idx, (N_mons, hmap) in enumerate(zip(self.k, self.hmaps)):
136 self.p_idxs.extend([p_idx]*N_mons)
137 self.m_idxs.append(slice(m_idx, m_idx+N_mons))
138 if getattr(self.hmaps[p_idx], "from_meq", False):
139 self.meq_idxs.all.add(m_idx)
140 if len(self.meq_idxs.all) > 2*len(self.meq_idxs.first_half):
141 self.meq_idxs.first_half.add(m_idx)
142 self.exps.extend(hmap)
143 self.cs.extend(hmap.values())
144 for exp in hmap:
145 if not exp: # space out A matrix with constants for mosek
146 row.append(m_idx)
147 col.append(0)
148 data.append(0)
149 for var in exp:
150 self.varlocs[var].append(m_idx)
151 m_idx += 1
152 self.p_idxs = np.array(self.p_idxs, "int32") # for later use as array
153 self.varidxs = {vk: i for i, vk in enumerate(self.varlocs)}
154 for j, (var, locs) in enumerate(self.varlocs.items()):
155 row.extend(locs)
156 col.extend([j]*len(locs))
157 data.extend(self.exps[i][var] for i in locs)
158 # A [mons, vks]: sparse array of each monomials' variables' exponents
159 self.A = CootMatrix(row, col, data)
161 # pylint: disable=too-many-statements, too-many-locals
162 def solve(self, solver=None, *, verbosity=1, gen_result=True, **kwargs):
163 """Solves a GeometricProgram and returns the solution.
165 Arguments
166 ---------
167 solver : str or function (optional)
168 By default uses a solver found during installation.
169 If "mosek_conif", "mosek_cli", or "cvxopt", uses that solver.
170 If a function, passes that function cs, A, p_idxs, and k.
171 verbosity : int (default 1)
172 If greater than 0, prints solver name and solve time.
173 **kwargs :
174 Passed to solver constructor and solver function.
177 Returns
178 -------
179 result : SolutionArray
180 """
181 solvername, solverfn = _get_solver(solver, kwargs)
182 solverargs = DEFAULT_SOLVER_KWARGS.get(solvername, {})
183 solverargs.update(kwargs)
184 solver_out = {}
186 if verbosity > 0:
187 print("Using solver '%s'" % solvername)
188 print(" for %i free variables" % len(self.varlocs))
189 print(" in %i posynomial inequalities." % len(self.k))
190 starttime = time()
191 infeasibility, original_stdout = None, sys.stdout
192 try:
193 sys.stdout = SolverLog(original_stdout, verbosity=verbosity-2)
194 solver_out = solverfn(c=self.cs, A=self.A, meq_idxs=self.meq_idxs,
195 k=self.k, p_idxs=self.p_idxs, **solverargs)
196 except Infeasible as e:
197 infeasibility = e
198 except Exception as e:
199 if isinstance(e, InvalidLicense):
200 raise InvalidLicense("license for solver \"%s\" is invalid."
201 % solvername) from e
202 raise UnknownInfeasible("Something unexpected went wrong.") from e
203 finally:
204 self.solve_log = "\n".join(sys.stdout)
205 sys.stdout = original_stdout
206 self.solver_out = solver_out
207 solver_out["solver"] = solvername
208 solver_out["soltime"] = time() - starttime
209 if verbosity > 0:
210 print("Solving took %.3g seconds." % solver_out["soltime"])
212 if infeasibility:
213 if isinstance(infeasibility, PrimalInfeasible):
214 msg = ("The model had no feasible points; "
215 "you may wish to relax some constraints or constants.")
216 elif isinstance(infeasibility, DualInfeasible):
217 msg = ("The model ran to an infinitely low cost;"
218 " bounding the right variables would prevent this.")
219 elif isinstance(infeasibility, UnknownInfeasible):
220 msg = "The solver failed for an unknown reason."
221 if verbosity > 0 and solver_out["soltime"] < 1:
222 print(msg + "\nSince this model solved in less than a second,"
223 " let's run `.debug()` automatically to check.\n`")
224 return self.model.debug(solver=solver)
225 msg += (" Running `.debug()` may pinpoint the trouble. You can"
226 " also try another solver, or increase the verbosity.")
227 raise infeasibility.__class__(msg) from infeasibility
229 if gen_result: # NOTE: SIDE EFFECTS
230 self._result = self.generate_result(solver_out,
231 verbosity=verbosity-2)
232 return self.result
233 # TODO: remove this "generate_result" closure
234 solver_out["generate_result"] = \
235 lambda: self.generate_result(solver_out, dual_check=False)
237 return solver_out
239 @property
240 def result(self):
241 "Creates and caches a result from the raw solver_out"
242 if not self._result:
243 self._result = self.generate_result(self.solver_out)
244 return self._result
246 def generate_result(self, solver_out, *, verbosity=0, dual_check=True):
247 "Generates a full SolutionArray and checks it."
248 if verbosity > 0:
249 soltime = solver_out["soltime"]
250 tic = time()
251 # result packing #
252 result = self._compile_result(solver_out) # NOTE: SIDE EFFECTS
253 if verbosity > 0:
254 print("Result packing took %.2g%% of solve time." %
255 ((time() - tic) / soltime * 100))
256 tic = time()
257 # solution checking #
258 try:
259 tol = SOLUTION_TOL.get(solver_out["solver"], 1e-5)
260 self.check_solution(result["cost"], solver_out['primal'],
261 solver_out["nu"], solver_out["la"], tol)
262 except Infeasible as chkerror:
263 chkwarn = str(chkerror)
264 if dual_check or ("Dual" not in chkwarn and "nu" not in chkwarn):
265 print("Solution check warning: %s" % chkwarn)
266 if verbosity > 0:
267 print("Solution checking took %.2g%% of solve time." %
268 ((time() - tic) / soltime * 100))
269 tic = time()
270 return result
272 def _generate_nula(self, solver_out):
273 if "nu" in solver_out:
274 # solver gave us monomial sensitivities, generate posynomial ones
275 solver_out["nu"] = nu = np.ravel(solver_out["nu"])
276 nu_by_posy = [nu[mi] for mi in self.m_idxs]
277 solver_out["la"] = la = np.array([sum(nup) for nup in nu_by_posy])
278 elif "la" in solver_out:
279 la = np.ravel(solver_out["la"])
280 if len(la) == len(self.hmaps) - 1:
281 # assume solver dropped the cost's sensitivity (always 1.0)
282 la = np.hstack(([1.0], la))
283 # solver gave us posynomial sensitivities, generate monomial ones
284 solver_out["la"] = la
285 z = np.log(self.cs) + self.A.dot(solver_out["primal"])
286 m_iss = [self.p_idxs == i for i in range(len(la))]
287 nu_by_posy = [la[p_i]*np.exp(z[m_is])/sum(np.exp(z[m_is]))
288 for p_i, m_is in enumerate(m_iss)]
289 solver_out["nu"] = np.hstack(nu_by_posy)
290 else:
291 raise RuntimeWarning("The dual solution was not returned.")
292 return la, nu_by_posy
294 def _compile_result(self, solver_out):
295 """Creates a result dict (as returned by solve() from solver output
297 This internal method is called from within the solve() method, unless
298 solver_out["status"] is not "optimal", in which case a RuntimeWarning
299 is raised prior to this method being called. In that case, users
300 may use this method to attempt to create a results dict from the
301 output of the failed solve.
303 Arguments
304 ---------
305 solver_out: dict
306 dict in format returned by solverfn within GeometricProgram.solve
308 Returns
309 -------
310 result: dict
311 dict in format returned by GeometricProgram.solve()
312 """
313 result = {"cost": float(solver_out["objective"])}
314 primal = solver_out["primal"]
315 if len(self.varlocs) != len(primal):
316 raise RuntimeWarning("The primal solution was not returned.")
317 result["freevariables"] = KeyDict(zip(self.varlocs, np.exp(primal)))
318 result["constants"] = KeyDict(self.substitutions)
319 result["variables"] = KeyDict(result["freevariables"])
320 result["variables"].update(result["constants"])
321 result["sensitivities"] = {"constraints": {}}
322 la, self.nu_by_posy = self._generate_nula(solver_out)
323 cost_senss = sum(nu_i*exp for (nu_i, exp) in zip(self.nu_by_posy[0],
324 self.cost.hmap))
325 self.v_ss = cost_senss.copy()
326 for las, nus, c in zip(la[1:], self.nu_by_posy[1:], self.hmaps[1:]):
327 while getattr(c, "parent", None) is not None:
328 c = c.parent
329 v_ss, c_senss = c.sens_from_dual(las, nus, result)
330 for vk, x in v_ss.items():
331 self.v_ss[vk] = x + self.v_ss.get(vk, 0)
332 while getattr(c, "generated_by", None) is not None:
333 c = c.generated_by
334 result["sensitivities"]["constraints"][c] = c_senss
335 # carry linked sensitivities over to their constants
336 for v in list(v for v in self.v_ss if v.gradients):
337 dlogcost_dlogv = self.v_ss.pop(v)
338 val = np.array(result["constants"][v])
339 for c, dv_dc in v.gradients.items():
340 with warnings.catch_warnings(): # skip pesky divide-by-zeros
341 warnings.simplefilter("ignore")
342 dlogv_dlogc = dv_dc * result["constants"][c]/val
343 before = self.v_ss.get(c, 0)
344 self.v_ss[c] = before + dlogcost_dlogv*dlogv_dlogc
345 if v in cost_senss:
346 if c in self.cost.varkeys:
347 dlogcost_dlogv = cost_senss.pop(v)
348 before = cost_senss.get(c, 0)
349 cost_senss[c] = before + dlogcost_dlogv*dlogv_dlogc
350 result["sensitivities"]["cost"] = cost_senss
351 result["sensitivities"]["variables"] = KeyDict(self.v_ss)
352 result["sensitivities"]["constants"] = \
353 result["sensitivities"]["variables"] # NOTE: backwards compat.
354 result["soltime"] = solver_out["soltime"]
355 return SolutionArray(result)
357 def check_solution(self, cost, primal, nu, la, tol, abstol=1e-20):
358 """Run checks to mathematically confirm solution solves this GP
360 Arguments
361 ---------
362 cost: float
363 cost returned by solver
364 primal: list
365 primal solution returned by solver
366 nu: numpy.ndarray
367 monomial lagrange multiplier
368 la: numpy.ndarray
369 posynomial lagrange multiplier
371 Raises
372 ------
373 Infeasible if any problems are found
374 """
375 A = self.A.tocsr()
376 def almost_equal(num1, num2):
377 "local almost equal test"
378 return (num1 == num2 or abs((num1 - num2) / (num1 + num2)) < tol
379 or abs(num1 - num2) < abstol)
380 # check primal sol #
381 primal_exp_vals = self.cs * np.exp(A.dot(primal)) # c*e^Ax
382 if not almost_equal(primal_exp_vals[self.m_idxs[0]].sum(), cost):
383 raise Infeasible("Primal solution computed cost did not match"
384 " solver-returned cost: %s vs %s." %
385 (primal_exp_vals[self.m_idxs[0]].sum(), cost))
386 for mi in self.m_idxs[1:]:
387 if primal_exp_vals[mi].sum() > 1 + tol:
388 raise Infeasible("Primal solution violates constraint: %s is "
389 "greater than 1" % primal_exp_vals[mi].sum())
390 # check dual sol #
391 # note: follows dual formulation in section 3.1 of
392 # http://web.mit.edu/~whoburg/www/papers/hoburg_phd_thesis.pdf
393 if not almost_equal(self.nu_by_posy[0].sum(), 1):
394 raise Infeasible("Dual variables associated with objective sum"
395 " to %s, not 1" % self.nu_by_posy[0].sum())
396 if any(nu < 0):
397 minnu = min(nu)
398 if minnu < -tol/1000:
399 raise Infeasible("Dual solution has negative entries as"
400 " large as %s." % minnu)
401 if any(np.abs(A.T.dot(nu)) > tol):
402 raise Infeasible("Sum of nu^T * A did not vanish.")
403 b = np.log(self.cs)
404 dual_cost = sum(
405 self.nu_by_posy[i].dot(
406 b[mi] - np.log(self.nu_by_posy[i]/la[i]))
407 for i, mi in enumerate(self.m_idxs) if la[i])
408 if not almost_equal(np.exp(dual_cost), cost):
409 raise Infeasible("Dual cost %s does not match primal cost %s"
410 % (np.exp(dual_cost), cost))
413def gen_meq_bounds(missingbounds, exps, meq_idxs): # pylint: disable=too-many-locals,too-many-branches
414 "Generate conditional monomial equality bounds"
415 meq_bounds = defaultdict(set)
416 for i in meq_idxs.first_half:
417 p_upper, p_lower, n_upper, n_lower = set(), set(), set(), set()
418 for key, x in exps[i].items():
419 if (key, "upper") in missingbounds:
420 if x > 0:
421 p_upper.add((key, "upper"))
422 else:
423 n_upper.add((key, "upper"))
424 if (key, "lower") in missingbounds:
425 if x > 0:
426 p_lower.add((key, "lower"))
427 else:
428 n_lower.add((key, "lower"))
429 # (consider x*y/z == 1)
430 # for a var (e.g. x) to be upper bounded by this monomial equality,
431 # - vars of the same sign/side (y) must be lower bounded
432 # - AND vars of the opposite sign/side (z) must be upper bounded
433 p_ub = n_lb = frozenset(n_upper).union(p_lower)
434 n_ub = p_lb = frozenset(p_upper).union(n_lower)
435 for keys, ub in ((p_upper, p_ub), (n_upper, n_ub)):
436 for key, _ in keys:
437 needed = ub.difference([(key, "lower")])
438 if needed:
439 meq_bounds[(key, "upper")].add(needed)
440 else:
441 del missingbounds[(key, "upper")]
442 for keys, lb in ((p_lower, p_lb), (n_lower, n_lb)):
443 for key, _ in keys:
444 needed = lb.difference([(key, "upper")])
445 if needed:
446 meq_bounds[(key, "lower")].add(needed)
447 else:
448 del missingbounds[(key, "lower")]
449 return meq_bounds
452def fulfill_meq_bounds(missingbounds, meq_bounds):
453 "Bounds variables with monomial equalities"
454 still_alive = True
455 while still_alive:
456 still_alive = False # if no changes are made, the loop exits
457 for bound in set(meq_bounds):
458 if bound not in missingbounds:
459 del meq_bounds[bound]
460 continue
461 for condition in meq_bounds[bound]:
462 if not any(bound in missingbounds for bound in condition):
463 del meq_bounds[bound]
464 del missingbounds[bound]
465 still_alive = True
466 break
467 for (var, bound) in meq_bounds:
468 boundstr = (", but would gain it from any of these sets: ")
469 for condition in list(meq_bounds[(var, bound)]):
470 meq_bounds[(var, bound)].remove(condition)
471 newcond = condition.intersection(missingbounds)
472 if newcond and not any(c.issubset(newcond)
473 for c in meq_bounds[(var, bound)]):
474 meq_bounds[(var, bound)].add(newcond)
475 boundstr += " or ".join(str(list(condition))
476 for condition in meq_bounds[(var, bound)])
477 missingbounds[(var, bound)] = boundstr