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