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