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

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)

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

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

23class MonoEqualityIndexes:

24 "Class to hold MonoEqualityIndexes"

26 def __init__(self):

27 self.all = set()

28 self.first_half = set()

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

54class GeometricProgram:

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

56 """Standard mathematical representation of a GP.

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

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

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)

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