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.constraints.gp.GeometricProgram(
63 # minimize
64 x,
65 [ # subject to
66 x >= 1,
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 TypeError("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.
122 k [posys]: number of monomials (rows of A) present in each constraint
123 m_idxs [mons]: monomial indices of each posynomial
124 p_idxs [mons]: posynomial index of each monomial
125 cs, exps [mons]: coefficient and exponents of each monomial
126 varlocs: {vk: monomial indices of each variables' location}
127 meq_idxs: {all indices of equality mons} and {the first index of each}
128 varidxs: {vk: which column corresponds to it in A}
129 A [mons, vks]: sparse array of each monomials' variables' exponents
131 """
132 self.k = [len(hmap) for hmap in self.hmaps]
133 self.m_idxs, self.p_idxs, self.cs, self.exps = [], [], [], []
134 self.varkeys = self.varlocs = defaultdict(list)
135 self.meq_idxs = MonoEqualityIndexes()
136 m_idx = 0
137 row, col, data = [], [], []
138 for p_idx, (N_mons, hmap) in enumerate(zip(self.k, self.hmaps)):
139 self.p_idxs.extend([p_idx]*N_mons)
140 self.m_idxs.append(slice(m_idx, m_idx+N_mons))
141 if getattr(self.hmaps[p_idx], "from_meq", False):
142 self.meq_idxs.all.add(m_idx)
143 if len(self.meq_idxs.all) > 2*len(self.meq_idxs.first_half):
144 self.meq_idxs.first_half.add(m_idx)
145 self.exps.extend(hmap)
146 self.cs.extend(hmap.values())
147 for exp in hmap:
148 if not exp: # space out A matrix with constants for mosek
149 row.append(m_idx)
150 col.append(0)
151 data.append(0)
152 for var in exp:
153 self.varlocs[var].append(m_idx)
154 m_idx += 1
155 self.p_idxs = np.array(self.p_idxs, "int32") # to use array equalities
156 self.varidxs = {vk: i for i, vk in enumerate(self.varlocs)}
157 for j, (var, locs) in enumerate(self.varlocs.items()):
158 row.extend(locs)
159 col.extend([j]*len(locs))
160 data.extend(self.exps[i][var] for i in locs)
161 self.A = CootMatrix(row, col, data)
163 # pylint: disable=too-many-statements, too-many-locals
164 def solve(self, solver=None, *, verbosity=1, gen_result=True, **kwargs):
165 """Solves a GeometricProgram and returns the solution.
167 Arguments
168 ---------
169 solver : str or function (optional)
170 By default uses a solver found during installation.
171 If "mosek_conif", "mosek_cli", or "cvxopt", uses that solver.
172 If a function, passes that function cs, A, p_idxs, and k.
173 verbosity : int (default 1)
174 If greater than 0, prints solver name and solve time.
175 gen_result : bool (default True)
176 If True, makes a human-readable SolutionArray from solver output.
177 **kwargs :
178 Passed to solver constructor and solver function.
181 Returns
182 -------
183 SolutionArray (or dict if gen_result is False)
184 """
185 solvername, solverfn = _get_solver(solver, kwargs)
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))
191 solverargs = DEFAULT_SOLVER_KWARGS.get(solvername, {})
192 solverargs.update(kwargs)
193 starttime = time()
194 solver_out, infeasibility, original_stdout = {}, None, sys.stdout
195 try:
196 sys.stdout = SolverLog(original_stdout, verbosity=verbosity-2)
197 solver_out = solverfn(c=self.cs, A=self.A, meq_idxs=self.meq_idxs,
198 k=self.k, p_idxs=self.p_idxs, **solverargs)
199 except Infeasible as e:
200 infeasibility = e
201 except InvalidLicense as e:
202 raise InvalidLicense("license for solver \"%s\" is invalid."
203 % solvername) from e
204 except Exception as e:
205 raise UnknownInfeasible("Something unexpected went wrong.") from e
206 finally:
207 self.solve_log = "\n".join(sys.stdout)
208 sys.stdout = original_stdout
209 self.solver_out = solver_out
211 solver_out["solver"] = solvername
212 solver_out["soltime"] = time() - starttime
213 if verbosity > 0:
214 print("Solving took %.3g seconds." % solver_out["soltime"])
216 if infeasibility:
217 if isinstance(infeasibility, PrimalInfeasible):
218 msg = ("The model had no feasible points; "
219 "you may wish to relax some constraints or constants.")
220 elif isinstance(infeasibility, DualInfeasible):
221 msg = ("The model ran to an infinitely low cost;"
222 " bounding the right variables would prevent this.")
223 elif isinstance(infeasibility, UnknownInfeasible):
224 msg = "The solver failed for an unknown reason."
225 if (verbosity > 0 and solver_out["soltime"] < 1
226 and hasattr(self, "model")): # fast, top-level model
227 print(msg + "\nSince the model solved in less than a second,"
228 " let's run `.debug()` to analyze what happened.\n`")
229 return self.model.debug(solver=solver)
230 # else, raise a clarifying error
231 msg += (" Running `.debug()` may pinpoint the trouble. You can"
232 " also try another solver, or increase the verbosity.")
233 raise infeasibility.__class__(msg) from infeasibility
235 if not gen_result:
236 return solver_out
237 # else, generate a human-readable SolutionArray
238 self._result = self.generate_result(solver_out, verbosity=verbosity-2)
239 return self.result
241 @property
242 def result(self):
243 "Creates and caches a result from the raw solver_out"
244 if not self._result:
245 self._result = self.generate_result(self.solver_out)
246 return self._result
248 def generate_result(self, solver_out, *, verbosity=0, dual_check=True):
249 "Generates a full SolutionArray and checks it."
250 if verbosity > 0:
251 soltime = solver_out["soltime"]
252 tic = time()
253 # result packing #
254 result = self._compile_result(solver_out) # NOTE: SIDE EFFECTS
255 if verbosity > 0:
256 print("Result packing took %.2g%% of solve time." %
257 ((time() - tic) / soltime * 100))
258 tic = time()
259 # solution checking #
260 try:
261 tol = SOLUTION_TOL.get(solver_out["solver"], 1e-5)
262 self.check_solution(result["cost"], solver_out['primal'],
263 solver_out["nu"], solver_out["la"], tol)
264 except Infeasible as chkerror:
265 chkwarn = str(chkerror)
266 if not ("Dual" in chkwarn and not dual_check):
267 print("Solution check warning: %s" % chkwarn)
268 if verbosity > 0:
269 print("Solution checking took %.2g%% of solve time." %
270 ((time() - tic) / soltime * 100))
271 return result
273 def _generate_nula(self, solver_out):
274 if "nu" in solver_out:
275 # solver gave us monomial sensitivities, generate posynomial ones
276 solver_out["nu"] = nu = np.ravel(solver_out["nu"])
277 nu_by_posy = [nu[mi] for mi in self.m_idxs]
278 solver_out["la"] = la = np.array([sum(nup) for nup in nu_by_posy])
279 elif "la" in solver_out:
280 la = np.ravel(solver_out["la"])
281 if len(la) == len(self.hmaps) - 1:
282 # assume solver dropped the cost's sensitivity (always 1.0)
283 la = np.hstack(([1.0], la))
284 # solver gave us posynomial sensitivities, generate monomial ones
285 solver_out["la"] = la
286 z = np.log(self.cs) + self.A.dot(solver_out["primal"])
287 m_iss = [self.p_idxs == i for i in range(len(la))]
288 nu_by_posy = [la[p_i]*np.exp(z[m_is])/sum(np.exp(z[m_is]))
289 for p_i, m_is in enumerate(m_iss)]
290 solver_out["nu"] = np.hstack(nu_by_posy)
291 else:
292 raise RuntimeWarning("The dual solution was not returned.")
293 return la, nu_by_posy
295 def _compile_result(self, solver_out):
296 result = {"cost": float(solver_out["objective"])}
297 primal = solver_out["primal"]
298 if len(self.varlocs) != len(primal):
299 raise RuntimeWarning("The primal solution was not returned.")
300 result["freevariables"] = KeyDict(zip(self.varlocs, np.exp(primal)))
301 result["constants"] = KeyDict(self.substitutions)
302 result["variables"] = KeyDict(result["freevariables"])
303 result["variables"].update(result["constants"])
304 result["sensitivities"] = {"constraints": {}}
305 la, self.nu_by_posy = self._generate_nula(solver_out)
306 cost_senss = sum(nu_i*exp for (nu_i, exp) in zip(self.nu_by_posy[0],
307 self.cost.hmap))
308 self.v_ss = cost_senss.copy()
309 for las, nus, c in zip(la[1:], self.nu_by_posy[1:], self.hmaps[1:]):
310 while getattr(c, "parent", None) is not None:
311 c = c.parent
312 v_ss, c_senss = c.sens_from_dual(las, nus, result)
313 for vk, x in v_ss.items():
314 self.v_ss[vk] = x + self.v_ss.get(vk, 0)
315 while getattr(c, "generated_by", None) is not None:
316 c = c.generated_by
317 result["sensitivities"]["constraints"][c] = c_senss
318 # carry linked sensitivities over to their constants
319 for v in list(v for v in self.v_ss if v.gradients):
320 dlogcost_dlogv = self.v_ss.pop(v)
321 val = np.array(result["constants"][v])
322 for c, dv_dc in v.gradients.items():
323 with warnings.catch_warnings(): # skip pesky divide-by-zeros
324 warnings.simplefilter("ignore")
325 dlogv_dlogc = dv_dc * result["constants"][c]/val
326 before = self.v_ss.get(c, 0)
327 self.v_ss[c] = before + dlogcost_dlogv*dlogv_dlogc
328 if v in cost_senss:
329 if c in self.cost.varkeys:
330 dlogcost_dlogv = cost_senss.pop(v)
331 before = cost_senss.get(c, 0)
332 cost_senss[c] = before + dlogcost_dlogv*dlogv_dlogc
333 result["sensitivities"]["cost"] = cost_senss
334 result["sensitivities"]["variables"] = KeyDict(self.v_ss)
335 result["sensitivities"]["constants"] = \
336 result["sensitivities"]["variables"] # NOTE: backwards compat.
337 result["soltime"] = solver_out["soltime"]
338 return SolutionArray(result)
340 def check_solution(self, cost, primal, nu, la, tol, abstol=1e-20):
341 """Run checks to mathematically confirm solution solves this GP
343 Arguments
344 ---------
345 cost: float
346 cost returned by solver
347 primal: list
348 primal solution returned by solver
349 nu: numpy.ndarray
350 monomial lagrange multiplier
351 la: numpy.ndarray
352 posynomial lagrange multiplier
354 Raises
355 ------
356 Infeasible if any problems are found
357 """
358 A = self.A.tocsr()
359 def almost_equal(num1, num2):
360 "local almost equal test"
361 return (num1 == num2 or abs((num1 - num2) / (num1 + num2)) < tol
362 or abs(num1 - num2) < abstol)
363 # check primal sol #
364 primal_exp_vals = self.cs * np.exp(A.dot(primal)) # c*e^Ax
365 if not almost_equal(primal_exp_vals[self.m_idxs[0]].sum(), cost):
366 raise Infeasible("Primal solution computed cost did not match"
367 " solver-returned cost: %s vs %s." %
368 (primal_exp_vals[self.m_idxs[0]].sum(), cost))
369 for mi in self.m_idxs[1:]:
370 if primal_exp_vals[mi].sum() > 1 + tol:
371 raise Infeasible("Primal solution violates constraint: %s is "
372 "greater than 1" % primal_exp_vals[mi].sum())
373 # check dual sol #
374 # note: follows dual formulation in section 3.1 of
375 # http://web.mit.edu/~whoburg/www/papers/hoburg_phd_thesis.pdf
376 if not almost_equal(self.nu_by_posy[0].sum(), 1):
377 raise Infeasible("Dual variables associated with objective sum"
378 " to %s, not 1" % self.nu_by_posy[0].sum())
379 if any(nu < 0):
380 minnu = min(nu)
381 if minnu < -tol/1000:
382 raise Infeasible("Dual solution has negative entries as"
383 " large as %s." % minnu)
384 if any(np.abs(A.T.dot(nu)) > tol):
385 raise Infeasible("Dual: sum of nu^T * A did not vanish.")
386 b = np.log(self.cs)
387 dual_cost = sum(
388 self.nu_by_posy[i].dot(
389 b[mi] - np.log(self.nu_by_posy[i]/la[i]))
390 for i, mi in enumerate(self.m_idxs) if la[i])
391 if not almost_equal(np.exp(dual_cost), cost):
392 raise Infeasible("Dual cost %s does not match primal cost %s"
393 % (np.exp(dual_cost), cost))
396def gen_meq_bounds(missingbounds, exps, meq_idxs): # pylint: disable=too-many-locals,too-many-branches
397 "Generate conditional monomial equality bounds"
398 meq_bounds = defaultdict(set)
399 for i in meq_idxs.first_half:
400 p_upper, p_lower, n_upper, n_lower = set(), set(), set(), set()
401 for key, x in exps[i].items():
402 if (key, "upper") in missingbounds:
403 if x > 0:
404 p_upper.add((key, "upper"))
405 else:
406 n_upper.add((key, "upper"))
407 if (key, "lower") in missingbounds:
408 if x > 0:
409 p_lower.add((key, "lower"))
410 else:
411 n_lower.add((key, "lower"))
412 # (consider x*y/z == 1)
413 # for a var (e.g. x) to be upper bounded by this monomial equality,
414 # - vars of the same sign (y) must be lower bounded
415 # - AND vars of the opposite sign (z) must be upper bounded
416 p_ub = n_lb = frozenset(n_upper).union(p_lower)
417 n_ub = p_lb = frozenset(p_upper).union(n_lower)
418 for keys, ub in ((p_upper, p_ub), (n_upper, n_ub)):
419 for key, _ in keys:
420 needed = ub.difference([(key, "lower")])
421 if needed:
422 meq_bounds[(key, "upper")].add(needed)
423 else:
424 del missingbounds[(key, "upper")]
425 for keys, lb in ((p_lower, p_lb), (n_lower, n_lb)):
426 for key, _ in keys:
427 needed = lb.difference([(key, "upper")])
428 if needed:
429 meq_bounds[(key, "lower")].add(needed)
430 else:
431 del missingbounds[(key, "lower")]
432 return meq_bounds
435def fulfill_meq_bounds(missingbounds, meq_bounds):
436 "Bounds variables with monomial equalities"
437 still_alive = True
438 while still_alive:
439 still_alive = False # if no changes are made, the loop exits
440 for bound in set(meq_bounds):
441 if bound not in missingbounds:
442 del meq_bounds[bound]
443 continue
444 for condition in meq_bounds[bound]:
445 if not any(bound in missingbounds for bound in condition):
446 del meq_bounds[bound]
447 del missingbounds[bound]
448 still_alive = True
449 break
450 for (var, bound) in meq_bounds:
451 boundstr = (", but would gain it from any of these sets: ")
452 for condition in list(meq_bounds[(var, bound)]):
453 meq_bounds[(var, bound)].remove(condition)
454 newcond = condition.intersection(missingbounds)
455 if newcond and not any(c.issubset(newcond)
456 for c in meq_bounds[(var, bound)]):
457 meq_bounds[(var, bound)].add(newcond)
458 boundstr += " or ".join(str(list(condition))
459 for condition in meq_bounds[(var, bound)])
460 missingbounds[(var, bound)] = boundstr