Coverage for gpkit/constraints/gp.py: 90%

336 statements  

« prev     ^ index     » next       coverage.py v7.4.0, created at 2024-01-07 22:15 -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) 

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 # pylint: disable=too-few-public-methods 

25 "Class to hold MonoEqualityIndexes" 

26 

27 def __init__(self): 

28 self.all = set() 

29 self.first_half = set() 

30 

31 

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 

54 

55 

56class GeometricProgram: 

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

58 """Standard mathematical representation of a GP. 

59 

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 

64 

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 

77 

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) 

98 

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 

125 

126 def gen(self): 

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

128 

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 

137 

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) 

171 

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. 

175 

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. 

188 

189 

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

199 

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 

222 

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

227 

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 

248 

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 

254 

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 

261 

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 

289 

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 

311 

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 

342 

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) 

391 

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 

395 

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 

406 

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

450 

451 

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 

489 

490 

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