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

333 statements  

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) 

17 

18 

19DEFAULT_SOLVER_KWARGS = {"cvxopt": {"kktsolver": "ldl"}} 

20SOLUTION_TOL = {"cvxopt": 1e-3, "mosek_cli": 1e-4, "mosek_conif": 1e-3} 

21 

22 

23class MonoEqualityIndexes: 

24 "Class to hold MonoEqualityIndexes" 

25 

26 def __init__(self): 

27 self.all = set() 

28 self.first_half = set() 

29 

30 

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 

52 

53 

54class GeometricProgram: 

55 # pylint: disable=too-many-instance-attributes 

56 """Standard mathematical representation of a GP. 

57 

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 

62 

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 

75 

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) 

96 

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 

123 

124 def gen(self): 

125 """Generates nomial and solve data (A, p_idxs) from posynomials. 

126 

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 

135 

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) 

169 

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. 

173 

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. 

186 

187 

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)) 

197 

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 

220 

221 solver_out["solver"] = solvername 

222 solver_out["soltime"] = time() - starttime 

223 if verbosity > 0: 

224 print("Solving took %.3g seconds." % solver_out["soltime"]) 

225 

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 

246 

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 

252 

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 

259 

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 

287 

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 

309 

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 

340 

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) 

389 

390 def check_solution(self, cost, primal, nu, la, tol, abstol=1e-20): 

391 """Run checks to mathematically confirm solution solves this GP 

392 

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 

403 

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)) 

446 

447 

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 

485 

486 

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