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

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 TypeError("substitution {%s: %s} has invalid value type" 

82 " %s." % (key, sub, type(sub))) 

83 cost_hmap = cost.hmap.sub(self.substitutions, cost.vks) 

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 

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 

130 

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) 

162 

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. 

166 

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. 

179 

180 

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

190 

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 

210 

211 solver_out["solver"] = solvername 

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

213 if verbosity > 0: 

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

215 

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 

234 

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 

240 

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 

247 

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 

272 

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 

294 

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 gpv_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 gpv_ss[vk] = x + gpv_ss.get(vk, 0) 

315 while getattr(c, "generated_by", 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 gpv_ss if v.gradients): 

320 dlogcost_dlogv = gpv_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 = gpv_ss.get(c, 0) 

327 gpv_ss[c] = before + dlogcost_dlogv*dlogv_dlogc 

328 if v in cost_senss: 

329 if c in self.cost.vks: 

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

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

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

337 result["soltime"] = solver_out["soltime"] 

338 return SolutionArray(result) 

339 

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

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

342 

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 

353 

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

394 

395 

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 

433 

434 

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