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