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) 

16 

17 

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

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

20 

21 

22class MonoEqualityIndexes: 

23 "Class to hold MonoEqualityIndexes" 

24 

25 def __init__(self): 

26 self.all = set() 

27 self.first_half = set() 

28 

29 

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 

51 

52 

53class GeometricProgram: 

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

55 """Standard mathematical representation of a GP. 

56 

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 

61 

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 

74 

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] + list(hmapgen) 

92 self.gen() # Generate various maps into the posy- and monomials 

93 if checkbounds: 

94 self.check_bounds(err_on_missing_bounds=True) 

95 

96 def check_bounds(self, *, err_on_missing_bounds=False): 

97 "Checks if any variables are unbounded, through equality constraints." 

98 missingbounds = {} 

99 for var, locs in self.varlocs.items(): 

100 upperbound, lowerbound = False, False 

101 for i in locs: 

102 if i not in self.meq_idxs.all: 

103 if self.exps[i][var] > 0: # pylint:disable=simplifiable-if-statement 

104 upperbound = True 

105 else: 

106 lowerbound = True 

107 if upperbound and lowerbound: 

108 break 

109 if not upperbound: 

110 missingbounds[(var, "upper")] = "." 

111 if not lowerbound: 

112 missingbounds[(var, "lower")] = "." 

113 if not missingbounds: 

114 return {} # all bounds found in inequalities 

115 meq_bounds = gen_meq_bounds(missingbounds, self.exps, self.meq_idxs) 

116 fulfill_meq_bounds(missingbounds, meq_bounds) 

117 if missingbounds and err_on_missing_bounds: 

118 raise UnboundedGP( 

119 "\n\n".join("%s has no %s bound%s" % (v, b, x) 

120 for (v, b), x in missingbounds.items())) 

121 return missingbounds 

122 

123 def gen(self): 

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

125 

126 k [posys]: number of monomials (rows of A) present in each constraint 

127 m_idxs [mons]: monomial indices of each posynomial 

128 p_idxs [mons]: posynomial index of each monomial 

129 cs, exps [mons]: coefficient and exponents of each monomial 

130 varlocs: {vk: monomial indices of each variables' location} 

131 meq_idxs: {all indices of equality mons} and {the first index of each} 

132 varidxs: {vk: which column corresponds to it in A} 

133 A [mons, vks]: sparse array of each monomials' variables' exponents 

134 

135 """ 

136 self.k = [len(hmap) for hmap in self.hmaps] 

137 self.m_idxs, self.p_idxs, self.cs, self.exps = [], [], [], [] 

138 self.varkeys = self.varlocs = defaultdict(list) 

139 self.meq_idxs = MonoEqualityIndexes() 

140 m_idx = 0 

141 row, col, data = [], [], [] 

142 for p_idx, (N_mons, hmap) in enumerate(zip(self.k, self.hmaps)): 

143 self.p_idxs.extend([p_idx]*N_mons) 

144 self.m_idxs.append(slice(m_idx, m_idx+N_mons)) 

145 if getattr(self.hmaps[p_idx], "from_meq", False): 

146 self.meq_idxs.all.add(m_idx) 

147 if len(self.meq_idxs.all) > 2*len(self.meq_idxs.first_half): 

148 self.meq_idxs.first_half.add(m_idx) 

149 self.exps.extend(hmap) 

150 self.cs.extend(hmap.values()) 

151 for exp in hmap: 

152 if not exp: # space out A matrix with constants for mosek 

153 row.append(m_idx) 

154 col.append(0) 

155 data.append(0) 

156 for var in exp: 

157 self.varlocs[var].append(m_idx) 

158 m_idx += 1 

159 self.p_idxs = np.array(self.p_idxs, "int32") # to use array equalities 

160 self.varidxs = {vk: i for i, vk in enumerate(self.varlocs)} 

161 self.choicevaridxs = {vk: i for i, vk in enumerate(self.varlocs) 

162 if vk.choices} 

163 for j, (var, locs) in enumerate(self.varlocs.items()): 

164 row.extend(locs) 

165 col.extend([j]*len(locs)) 

166 data.extend(self.exps[i][var] for i in locs) 

167 self.A = CootMatrix(row, col, data) 

168 

169 # pylint: disable=too-many-statements, too-many-locals 

170 def solve(self, solver=None, *, verbosity=1, gen_result=True, **kwargs): 

171 """Solves a GeometricProgram and returns the solution. 

172 

173 Arguments 

174 --------- 

175 solver : str or function (optional) 

176 By default uses a solver found during installation. 

177 If "mosek_conif", "mosek_cli", or "cvxopt", uses that solver. 

178 If a function, passes that function cs, A, p_idxs, and k. 

179 verbosity : int (default 1) 

180 If greater than 0, prints solver name and solve time. 

181 gen_result : bool (default True) 

182 If True, makes a human-readable SolutionArray from solver output. 

183 **kwargs : 

184 Passed to solver constructor and solver function. 

185 

186 

187 Returns 

188 ------- 

189 SolutionArray (or dict if gen_result is False) 

190 """ 

191 solvername, solverfn = _get_solver(solver, kwargs) 

192 if verbosity > 0: 

193 print("Using solver '%s'" % solvername) 

194 print(" for %i free variables" % len(self.varlocs)) 

195 print(" in %i posynomial inequalities." % len(self.k)) 

196 

197 solverargs = DEFAULT_SOLVER_KWARGS.get(solvername, {}) 

198 solverargs.update(kwargs) 

199 if self.choicevaridxs and solvername == "mosek_conif": 

200 solverargs["choicevaridxs"] = self.choicevaridxs 

201 self.integersolve = True 

202 starttime = time() 

203 solver_out, infeasibility, original_stdout = {}, None, sys.stdout 

204 try: 

205 sys.stdout = SolverLog(original_stdout, verbosity=verbosity-2) 

206 solver_out = solverfn(c=self.cs, A=self.A, meq_idxs=self.meq_idxs, 

207 k=self.k, p_idxs=self.p_idxs, **solverargs) 

208 except Infeasible as e: 

209 infeasibility = e 

210 except InvalidLicense as e: 

211 raise InvalidLicense("license for solver \"%s\" is invalid." 

212 % solvername) from e 

213 except Exception as e: 

214 raise UnknownInfeasible("Something unexpected went wrong.") from e 

215 finally: 

216 self.solve_log = "\n".join(sys.stdout) 

217 sys.stdout = original_stdout 

218 self.solver_out = solver_out 

219 

220 solver_out["solver"] = solvername 

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

222 if verbosity > 0: 

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

224 

225 if infeasibility: 

226 if isinstance(infeasibility, PrimalInfeasible): 

227 msg = ("The model had no feasible points; " 

228 "you may wish to relax some constraints or constants.") 

229 elif isinstance(infeasibility, DualInfeasible): 

230 msg = ("The model ran to an infinitely low cost;" 

231 " bounding the right variables would prevent this.") 

232 elif isinstance(infeasibility, UnknownInfeasible): 

233 msg = ("Solver failed for an unknown reason. Relaxing" 

234 " constraints/constants, bounding variables, or" 

235 " using a different solver might fix it.") 

236 if (verbosity > 0 and solver_out["soltime"] < 1 

237 and hasattr(self, "model")): # fast, top-level model 

238 print(msg + "\nSince the model solved in less than a second," 

239 " let's run `.debug()` to analyze what happened.\n`") 

240 return self.model.debug(solver=solver) 

241 # else, raise a clarifying error 

242 msg += (" Running `.debug()` or increasing verbosity may pinpoint" 

243 " the trouble.") 

244 raise infeasibility.__class__(msg) from infeasibility 

245 

246 if not gen_result: 

247 return solver_out 

248 # else, generate a human-readable SolutionArray 

249 self._result = self.generate_result(solver_out, verbosity=verbosity-2) 

250 return self.result 

251 

252 @property 

253 def result(self): 

254 "Creates and caches a result from the raw solver_out" 

255 if not self._result: 

256 self._result = self.generate_result(self.solver_out) 

257 return self._result 

258 

259 def generate_result(self, solver_out, *, verbosity=0, dual_check=True): 

260 "Generates a full SolutionArray and checks it." 

261 if verbosity > 0: 

262 soltime = solver_out["soltime"] 

263 tic = time() 

264 # result packing # 

265 result = self._compile_result(solver_out) # NOTE: SIDE EFFECTS 

266 if verbosity > 0: 

267 print("Result packing took %.2g%% of solve time." % 

268 ((time() - tic) / soltime * 100)) 

269 tic = time() 

270 # solution checking # 

271 initsolwarning(result, "Solution Inconsistency") 

272 try: 

273 tol = SOLUTION_TOL.get(solver_out["solver"], 1e-5) 

274 self.check_solution(result["cost"], solver_out['primal'], 

275 solver_out["nu"], solver_out["la"], tol) 

276 except Infeasible as chkerror: 

277 msg = str(chkerror) 

278 if not ("Dual" in msg and not dual_check): 

279 appendsolwarning(msg, None, result, "Solution Inconsistency") 

280 if verbosity > -4: 

281 print("Solution check warning: %s" % msg) 

282 if verbosity > 0: 

283 print("Solution checking took %.2g%% of solve time." % 

284 ((time() - tic) / soltime * 100)) 

285 return result 

286 

287 def _generate_nula(self, solver_out): 

288 if "nu" in solver_out: 

289 # solver gave us monomial sensitivities, generate posynomial ones 

290 solver_out["nu"] = nu = np.ravel(solver_out["nu"]) 

291 nu_by_posy = [nu[mi] for mi in self.m_idxs] 

292 solver_out["la"] = la = np.array([sum(nup) for nup in nu_by_posy]) 

293 elif "la" in solver_out: 

294 la = np.ravel(solver_out["la"]) 

295 if len(la) == len(self.hmaps) - 1: 

296 # assume solver dropped the cost's sensitivity (always 1.0) 

297 la = np.hstack(([1.0], la)) 

298 # solver gave us posynomial sensitivities, generate monomial ones 

299 solver_out["la"] = la 

300 z = np.log(self.cs) + self.A.dot(solver_out["primal"]) 

301 m_iss = [self.p_idxs == i for i in range(len(la))] 

302 nu_by_posy = [la[p_i]*np.exp(z[m_is])/sum(np.exp(z[m_is])) 

303 for p_i, m_is in enumerate(m_iss)] 

304 solver_out["nu"] = np.hstack(nu_by_posy) 

305 else: 

306 raise RuntimeWarning("The dual solution was not returned.") 

307 return la, nu_by_posy 

308 

309 def _compile_result(self, solver_out): 

310 result = {"cost": float(solver_out["objective"])} 

311 primal = solver_out["primal"] 

312 if len(self.varlocs) != len(primal): 

313 raise RuntimeWarning("The primal solution was not returned.") 

314 result["freevariables"] = KeyDict(zip(self.varlocs, np.exp(primal))) 

315 result["constants"] = KeyDict(self.substitutions) 

316 result["variables"] = KeyDict(result["freevariables"]) 

317 result["variables"].update(result["constants"]) 

318 result["soltime"] = solver_out["soltime"] 

319 if self.integersolve: 

320 result["choicevariables"] = KeyDict( \ 

321 {k: v for k, v in result["freevariables"].items() 

322 if k in self.choicevaridxs}) 

323 result["warnings"] = {"No Dual Solution": [(\ 

324 "This model has the discretized choice variables" 

325 " %s and hence no dual solution. You can fix those variables" 

326 " to their optimal value and get sensitivities to the resulting" 

327 " continuous problem by updating your model's substitions with" 

328 " `sol[\"choicevariables\"]`." 

329 % sorted(self.choicevaridxs.keys()), self.choicevaridxs)]} 

330 return SolutionArray(result) 

331 if self.choicevaridxs: 

332 result["warnings"] = {"Freed Choice Variables": [(\ 

333 "This model has the discretized choice variables" 

334 " %s, but since the '%s' solver doesn't support discretization" 

335 " they were treated as continuous variables." 

336 % (sorted(self.choicevaridxs.keys()), solver_out["solver"]), 

337 self.choicevaridxs)]} # TODO: choicevaridxs seems unnecessary 

338 

339 result["sensitivities"] = {"constraints": {}} 

340 la, self.nu_by_posy = self._generate_nula(solver_out) 

341 cost_senss = sum(nu_i*exp for (nu_i, exp) in zip(self.nu_by_posy[0], 

342 self.cost.hmap)) 

343 gpv_ss = cost_senss.copy() 

344 m_senss = defaultdict(float) 

345 absv_ss = {vk: abs(x) for vk, x in cost_senss.items()} 

346 for las, nus, c in zip(la[1:], self.nu_by_posy[1:], self.hmaps[1:]): 

347 while getattr(c, "parent", None) is not None: 

348 c = c.parent 

349 v_ss, c_senss = c.sens_from_dual(las, nus, result) 

350 for vk, x in v_ss.items(): 

351 gpv_ss[vk] = x + gpv_ss.get(vk, 0) 

352 absv_ss[vk] = abs(x) + absv_ss.get(vk, 0) 

353 while getattr(c, "generated_by", None): 

354 c = c.generated_by 

355 result["sensitivities"]["constraints"][c] = c_senss 

356 m_senss[lineagestr(c)] += abs(c_senss) 

357 result["sensitivities"]["models"] = dict(m_senss) 

358 # carry linked sensitivities over to their constants 

359 for v in list(v for v in gpv_ss if v.gradients): 

360 dlogcost_dlogv = gpv_ss.pop(v) 

361 dlogcost_dlogabsv = absv_ss.pop(v) 

362 val = np.array(result["constants"][v]) 

363 for c, dv_dc in v.gradients.items(): 

364 with pywarnings.catch_warnings(): # skip pesky divide-by-zeros 

365 pywarnings.simplefilter("ignore") 

366 dlogv_dlogc = dv_dc * result["constants"][c]/val 

367 gpv_ss[c] = gpv_ss.get(c, 0) + dlogcost_dlogv*dlogv_dlogc 

368 absv_ss[c] = (absv_ss.get(c, 0) 

369 + abs(dlogcost_dlogabsv*dlogv_dlogc)) 

370 if v in cost_senss: 

371 if c in self.cost.vks: # TODO: seems unnecessary 

372 dlogcost_dlogv = cost_senss.pop(v) 

373 before = cost_senss.get(c, 0) 

374 cost_senss[c] = before + dlogcost_dlogv*dlogv_dlogc 

375 # add fixed variables sensitivities to models 

376 for vk, senss in gpv_ss.items(): 

377 m_senss[lineagestr(vk)] += abs(senss) 

378 result["sensitivities"]["cost"] = cost_senss 

379 result["sensitivities"]["variables"] = KeyDict(gpv_ss) 

380 result["sensitivities"]["variablerisk"] = KeyDict(absv_ss) 

381 result["sensitivities"]["constants"] = \ 

382 result["sensitivities"]["variables"] # NOTE: backwards compat. 

383 return SolutionArray(result) 

384 

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

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

387 

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 

398 

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

441 

442 

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 

480 

481 

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