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) 

14 

15 

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

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

18 

19 

20class MonoEqualityIndexes: 

21 "Class to hold MonoEqualityIndexes" 

22 

23 def __init__(self): 

24 self.all = set() 

25 self.first_half = set() 

26 

27 

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 

49 

50 

51class GeometricProgram: 

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

53 """Standard mathematical representation of a GP. 

54 

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 

59 

60 Examples 

61 -------- 

62 >>> gp = gpkit.geometric_program.GeometricProgram( 

63 # minimize 

64 x, 

65 [ # subject to 

66 1/x # <= 1, implicitly 

67 ], {}) 

68 >>> gp.solve() 

69 """ 

70 _result = solve_log = solver_out = model = v_ss = nu_by_posy = None 

71 

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 ValueError("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) 

91 

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 

118 

119 def gen(self): 

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

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

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

123 # m_idxs [mons]: monomial indices of each posynomial 

124 self.m_idxs = [] 

125 # p_idxs [mons]: posynomial index of each monomial 

126 self.p_idxs = [] 

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

128 self.cs, self.exps = [], [] 

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

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

131 # meq_idxs: {all indices of equality mons} and {just the first halves} 

132 self.meq_idxs = MonoEqualityIndexes() 

133 m_idx = 0 

134 row, col, data = [], [], [] 

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

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

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

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

139 self.meq_idxs.all.add(m_idx) 

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

141 self.meq_idxs.first_half.add(m_idx) 

142 self.exps.extend(hmap) 

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

144 for exp in hmap: 

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

146 row.append(m_idx) 

147 col.append(0) 

148 data.append(0) 

149 for var in exp: 

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

151 m_idx += 1 

152 self.p_idxs = np.array(self.p_idxs, "int32") # for later use as array 

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

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

155 row.extend(locs) 

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

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

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

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

160 

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

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

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

164 

165 Arguments 

166 --------- 

167 solver : str or function (optional) 

168 By default uses a solver found during installation. 

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

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

171 verbosity : int (default 1) 

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

173 **kwargs : 

174 Passed to solver constructor and solver function. 

175 

176 

177 Returns 

178 ------- 

179 result : SolutionArray 

180 """ 

181 solvername, solverfn = _get_solver(solver, kwargs) 

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

183 solverargs.update(kwargs) 

184 solver_out = {} 

185 

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

190 starttime = time() 

191 infeasibility, original_stdout = None, sys.stdout 

192 try: 

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

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

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

196 except Infeasible as e: 

197 infeasibility = e 

198 except Exception as e: 

199 if isinstance(e, InvalidLicense): 

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

201 % solvername) from e 

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

203 finally: 

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

205 sys.stdout = original_stdout 

206 self.solver_out = solver_out 

207 solver_out["solver"] = solvername 

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

209 if verbosity > 0: 

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

211 

212 if infeasibility: 

213 if isinstance(infeasibility, PrimalInfeasible): 

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

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

216 elif isinstance(infeasibility, DualInfeasible): 

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

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

219 elif isinstance(infeasibility, UnknownInfeasible): 

220 msg = "The solver failed for an unknown reason." 

221 if verbosity > 0 and solver_out["soltime"] < 1: 

222 print(msg + "\nSince this model solved in less than a second," 

223 " let's run `.debug()` automatically to check.\n`") 

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

225 msg += (" Running `.debug()` may pinpoint the trouble. You can" 

226 " also try another solver, or increase the verbosity.") 

227 raise infeasibility.__class__(msg) from infeasibility 

228 

229 if gen_result: # NOTE: SIDE EFFECTS 

230 self._result = self.generate_result(solver_out, 

231 verbosity=verbosity-2) 

232 return self.result 

233 # TODO: remove this "generate_result" closure 

234 solver_out["generate_result"] = \ 

235 lambda: self.generate_result(solver_out, dual_check=False) 

236 

237 return solver_out 

238 

239 @property 

240 def result(self): 

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

242 if not self._result: 

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

244 return self._result 

245 

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

247 "Generates a full SolutionArray and checks it." 

248 if verbosity > 0: 

249 soltime = solver_out["soltime"] 

250 tic = time() 

251 # result packing # 

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

253 if verbosity > 0: 

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

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

256 tic = time() 

257 # solution checking # 

258 try: 

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

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

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

262 except Infeasible as chkerror: 

263 chkwarn = str(chkerror) 

264 if dual_check or ("Dual" not in chkwarn and "nu" not in chkwarn): 

265 print("Solution check warning: %s" % chkwarn) 

266 if verbosity > 0: 

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

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

269 tic = time() 

270 return result 

271 

272 def _generate_nula(self, solver_out): 

273 if "nu" in solver_out: 

274 # solver gave us monomial sensitivities, generate posynomial ones 

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

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

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

278 elif "la" in solver_out: 

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

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

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

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

283 # solver gave us posynomial sensitivities, generate monomial ones 

284 solver_out["la"] = la 

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

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

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

288 for p_i, m_is in enumerate(m_iss)] 

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

290 else: 

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

292 return la, nu_by_posy 

293 

294 def _compile_result(self, solver_out): 

295 """Creates a result dict (as returned by solve() from solver output 

296 

297 This internal method is called from within the solve() method, unless 

298 solver_out["status"] is not "optimal", in which case a RuntimeWarning 

299 is raised prior to this method being called. In that case, users 

300 may use this method to attempt to create a results dict from the 

301 output of the failed solve. 

302 

303 Arguments 

304 --------- 

305 solver_out: dict 

306 dict in format returned by solverfn within GeometricProgram.solve 

307 

308 Returns 

309 ------- 

310 result: dict 

311 dict in format returned by GeometricProgram.solve() 

312 """ 

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

314 primal = solver_out["primal"] 

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

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

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

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

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

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

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

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

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

324 self.cost.hmap)) 

325 self.v_ss = cost_senss.copy() 

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

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

328 c = c.parent 

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

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

331 self.v_ss[vk] = x + self.v_ss.get(vk, 0) 

332 while getattr(c, "generated_by", None) is not None: 

333 c = c.generated_by 

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

335 # carry linked sensitivities over to their constants 

336 for v in list(v for v in self.v_ss if v.gradients): 

337 dlogcost_dlogv = self.v_ss.pop(v) 

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

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

340 with warnings.catch_warnings(): # skip pesky divide-by-zeros 

341 warnings.simplefilter("ignore") 

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

343 before = self.v_ss.get(c, 0) 

344 self.v_ss[c] = before + dlogcost_dlogv*dlogv_dlogc 

345 if v in cost_senss: 

346 if c in self.cost.varkeys: 

347 dlogcost_dlogv = cost_senss.pop(v) 

348 before = cost_senss.get(c, 0) 

349 cost_senss[c] = before + dlogcost_dlogv*dlogv_dlogc 

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

351 result["sensitivities"]["variables"] = KeyDict(self.v_ss) 

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

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

354 result["soltime"] = solver_out["soltime"] 

355 return SolutionArray(result) 

356 

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

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

359 

360 Arguments 

361 --------- 

362 cost: float 

363 cost returned by solver 

364 primal: list 

365 primal solution returned by solver 

366 nu: numpy.ndarray 

367 monomial lagrange multiplier 

368 la: numpy.ndarray 

369 posynomial lagrange multiplier 

370 

371 Raises 

372 ------ 

373 Infeasible if any problems are found 

374 """ 

375 A = self.A.tocsr() 

376 def almost_equal(num1, num2): 

377 "local almost equal test" 

378 return (num1 == num2 or abs((num1 - num2) / (num1 + num2)) < tol 

379 or abs(num1 - num2) < abstol) 

380 # check primal sol # 

381 primal_exp_vals = self.cs * np.exp(A.dot(primal)) # c*e^Ax 

382 if not almost_equal(primal_exp_vals[self.m_idxs[0]].sum(), cost): 

383 raise Infeasible("Primal solution computed cost did not match" 

384 " solver-returned cost: %s vs %s." % 

385 (primal_exp_vals[self.m_idxs[0]].sum(), cost)) 

386 for mi in self.m_idxs[1:]: 

387 if primal_exp_vals[mi].sum() > 1 + tol: 

388 raise Infeasible("Primal solution violates constraint: %s is " 

389 "greater than 1" % primal_exp_vals[mi].sum()) 

390 # check dual sol # 

391 # note: follows dual formulation in section 3.1 of 

392 # http://web.mit.edu/~whoburg/www/papers/hoburg_phd_thesis.pdf 

393 if not almost_equal(self.nu_by_posy[0].sum(), 1): 

394 raise Infeasible("Dual variables associated with objective sum" 

395 " to %s, not 1" % self.nu_by_posy[0].sum()) 

396 if any(nu < 0): 

397 minnu = min(nu) 

398 if minnu < -tol/1000: 

399 raise Infeasible("Dual solution has negative entries as" 

400 " large as %s." % minnu) 

401 if any(np.abs(A.T.dot(nu)) > tol): 

402 raise Infeasible("Sum of nu^T * A did not vanish.") 

403 b = np.log(self.cs) 

404 dual_cost = sum( 

405 self.nu_by_posy[i].dot( 

406 b[mi] - np.log(self.nu_by_posy[i]/la[i])) 

407 for i, mi in enumerate(self.m_idxs) if la[i]) 

408 if not almost_equal(np.exp(dual_cost), cost): 

409 raise Infeasible("Dual cost %s does not match primal cost %s" 

410 % (np.exp(dual_cost), cost)) 

411 

412 

413def gen_meq_bounds(missingbounds, exps, meq_idxs): # pylint: disable=too-many-locals,too-many-branches 

414 "Generate conditional monomial equality bounds" 

415 meq_bounds = defaultdict(set) 

416 for i in meq_idxs.first_half: 

417 p_upper, p_lower, n_upper, n_lower = set(), set(), set(), set() 

418 for key, x in exps[i].items(): 

419 if (key, "upper") in missingbounds: 

420 if x > 0: 

421 p_upper.add((key, "upper")) 

422 else: 

423 n_upper.add((key, "upper")) 

424 if (key, "lower") in missingbounds: 

425 if x > 0: 

426 p_lower.add((key, "lower")) 

427 else: 

428 n_lower.add((key, "lower")) 

429 # (consider x*y/z == 1) 

430 # for a var (e.g. x) to be upper bounded by this monomial equality, 

431 # - vars of the same sign/side (y) must be lower bounded 

432 # - AND vars of the opposite sign/side (z) must be upper bounded 

433 p_ub = n_lb = frozenset(n_upper).union(p_lower) 

434 n_ub = p_lb = frozenset(p_upper).union(n_lower) 

435 for keys, ub in ((p_upper, p_ub), (n_upper, n_ub)): 

436 for key, _ in keys: 

437 needed = ub.difference([(key, "lower")]) 

438 if needed: 

439 meq_bounds[(key, "upper")].add(needed) 

440 else: 

441 del missingbounds[(key, "upper")] 

442 for keys, lb in ((p_lower, p_lb), (n_lower, n_lb)): 

443 for key, _ in keys: 

444 needed = lb.difference([(key, "upper")]) 

445 if needed: 

446 meq_bounds[(key, "lower")].add(needed) 

447 else: 

448 del missingbounds[(key, "lower")] 

449 return meq_bounds 

450 

451 

452def fulfill_meq_bounds(missingbounds, meq_bounds): 

453 "Bounds variables with monomial equalities" 

454 still_alive = True 

455 while still_alive: 

456 still_alive = False # if no changes are made, the loop exits 

457 for bound in set(meq_bounds): 

458 if bound not in missingbounds: 

459 del meq_bounds[bound] 

460 continue 

461 for condition in meq_bounds[bound]: 

462 if not any(bound in missingbounds for bound in condition): 

463 del meq_bounds[bound] 

464 del missingbounds[bound] 

465 still_alive = True 

466 break 

467 for (var, bound) in meq_bounds: 

468 boundstr = (", but would gain it from any of these sets: ") 

469 for condition in list(meq_bounds[(var, bound)]): 

470 meq_bounds[(var, bound)].remove(condition) 

471 newcond = condition.intersection(missingbounds) 

472 if newcond and not any(c.issubset(newcond) 

473 for c in meq_bounds[(var, bound)]): 

474 meq_bounds[(var, bound)].add(newcond) 

475 boundstr += " or ".join(str(list(condition)) 

476 for condition in meq_bounds[(var, bound)]) 

477 missingbounds[(var, bound)] = boundstr