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