1"""Implement the SequentialGeometricProgram class""" 

2import warnings as pywarnings 

3from time import time 

4from collections import defaultdict 

5import numpy as np 

6from ..exceptions import (InvalidGPConstraint, Infeasible, UnnecessarySGP, 

7 InvalidPosynomial, InvalidSGPConstraint) 

8from ..keydict import KeyDict 

9from ..nomials import Variable 

10from .gp import GeometricProgram 

11from ..nomials import PosynomialInequality, Posynomial 

12from .. import NamedVariables 

13from ..small_scripts import appendsolwarning, initsolwarning 

14 

15 

16EPS = 1e-6 # 1 +/- this is used in a few relative differences 

17 

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

19class SequentialGeometricProgram: 

20 """Prepares a collection of signomials for a SP solve. 

21 

22 Arguments 

23 --------- 

24 cost : Posynomial 

25 Objective to minimize when solving 

26 constraints : list of Constraint or SignomialConstraint objects 

27 Constraints to maintain when solving (implicitly Signomials <= 1) 

28 verbosity : int (optional) 

29 Currently has no effect: SequentialGeometricPrograms don't know 

30 anything new after being created, unlike GeometricPrograms. 

31 

32 Attributes with side effects 

33 ---------------------------- 

34 `gps` is set during a solve 

35 `result` is set at the end of a solve 

36 

37 Examples 

38 -------- 

39 >>> gp = gpkit.geometric_program.SequentialGeometricProgram( 

40 # minimize 

41 x, 

42 [ # subject to 

43 1/x - y/x, # <= 1, implicitly 

44 y/10 # <= 1 

45 ]) 

46 >>> gp.solve() 

47 """ 

48 gps = solver_outs = _results = result = model = None 

49 with NamedVariables("RelaxPCCP"): 

50 slack = Variable("C") 

51 

52 def __init__(self, cost, model, substitutions, 

53 *, use_pccp=True, pccp_penalty=2e2, checkbounds=True, **_): 

54 self.pccp_penalty = pccp_penalty 

55 if cost.any_nonpositive_cs: 

56 raise InvalidPosynomial("""an SGP's cost must be Posynomial 

57 

58 The equivalent of a Signomial objective can be constructed by constraining 

59 a dummy variable `z` to be greater than the desired Signomial objective `s` 

60 (z >= s) and then minimizing that dummy variable.""") 

61 self.gpconstraints, self.sgpconstraints = [], [] 

62 if not use_pccp: 

63 self.slack = 1 

64 else: 

65 self.gpconstraints.append(self.slack >= 1) 

66 cost *= self.slack**pccp_penalty 

67 self.approxconstraints = [] 

68 self.sgpvks = set() 

69 x0 = KeyDict(substitutions) 

70 x0.vks = model.vks # for string access and so forth 

71 for cs in model.flat(): 

72 try: 

73 if not hasattr(cs, "as_hmapslt1"): 

74 raise InvalidGPConstraint(cs) 

75 if not isinstance(cs, PosynomialInequality): 

76 cs.as_hmapslt1(substitutions) # gp-compatible? 

77 self.gpconstraints.append(cs) 

78 except InvalidGPConstraint: 

79 if not hasattr(cs, "as_gpconstr"): 

80 raise InvalidSGPConstraint(cs) 

81 self.sgpconstraints.append(cs) 

82 for hmaplt1 in cs.as_gpconstr(x0).as_hmapslt1({}): 

83 constraint = (Posynomial(hmaplt1) <= self.slack) 

84 constraint.generated_by = cs 

85 self.approxconstraints.append(constraint) 

86 self.sgpvks.update(constraint.vks) 

87 if not self.sgpconstraints: 

88 raise UnnecessarySGP("""Model valid as a Geometric Program. 

89 

90SequentialGeometricPrograms should only be created with Models containing 

91Signomial Constraints, since Models without Signomials have global 

92solutions and can be solved with 'Model.solve()'.""") 

93 self._gp = GeometricProgram( 

94 cost, self.approxconstraints + self.gpconstraints, 

95 substitutions, checkbounds=checkbounds) 

96 self._gp.x0 = x0 

97 self.a_idxs = defaultdict(list) 

98 cost_mons = self._gp.k[0] 

99 sp_mons = sum(self._gp.k[:1+len(self.approxconstraints)]) 

100 for row_idx, m_idx in enumerate(self._gp.A.row): 

101 if cost_mons <= m_idx <= sp_mons: 

102 self.a_idxs[self._gp.p_idxs[m_idx]].append(row_idx) 

103 

104 # pylint: disable=too-many-locals,too-many-branches,too-many-statements 

105 def localsolve(self, solver=None, *, verbosity=1, x0=None, reltol=1e-4, 

106 iteration_limit=50, **solveargs): 

107 """Locally solves a SequentialGeometricProgram and returns the solution. 

108 

109 Arguments 

110 --------- 

111 solver : str or function (optional) 

112 By default uses one of the solvers found during installation. 

113 If set to "mosek", "mosek_cli", or "cvxopt", uses that solver. 

114 If set to a function, passes that function cs, A, p_idxs, and k. 

115 verbosity : int (optional) 

116 If greater than 0, prints solve time and number of iterations. 

117 Each GP is created and solved with verbosity one less than this, so 

118 if greater than 1, prints solver name and time for each GP. 

119 x0 : dict (optional) 

120 Initial location to approximate signomials about. 

121 reltol : float 

122 Iteration ends when this is greater than the distance between two 

123 consecutive solve's objective values. 

124 iteration_limit : int 

125 Maximum GP iterations allowed. 

126 mutategp: boolean 

127 Prescribes whether to mutate the previously generated GP 

128 or to create a new GP with every solve. 

129 **solveargs : 

130 Passed to solver function. 

131 

132 Returns 

133 ------- 

134 result : dict 

135 A dictionary containing the translated solver result. 

136 """ 

137 self.gps, self.solver_outs, self._results = [], [], [] 

138 starttime = time() 

139 if verbosity > 0: 

140 print("Starting a sequence of GP solves") 

141 print(" for %i free variables" % len(self.sgpvks)) 

142 print(" in %i locally-GP constraints" % len(self.sgpconstraints)) 

143 print(" and for %i free variables" % len(self._gp.varlocs)) 

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

145 prevcost, cost, rel_improvement = None, None, None 

146 while rel_improvement is None or rel_improvement > reltol: 

147 prevcost = cost 

148 if len(self.gps) > iteration_limit: 

149 raise Infeasible( 

150 "Unsolved after %s iterations. Check `m.program.results`;" 

151 " if they're converging, try `.localsolve(...," 

152 " iteration_limit=NEWLIMIT)`." % len(self.gps)) 

153 gp = self.gp(x0, cleanx0=(len(self.gps) >= 1)) # clean the first x0 

154 self.gps.append(gp) # NOTE: SIDE EFFECTS 

155 if verbosity > 1: 

156 print("\nGP Solve %i" % len(self.gps)) 

157 if verbosity > 2: 

158 print("===============") 

159 solver_out = gp.solve(solver, verbosity=verbosity-1, 

160 gen_result=False, **solveargs) 

161 self.solver_outs.append(solver_out) 

162 cost = float(solver_out["objective"]) 

163 x0 = dict(zip(gp.varlocs, np.exp(solver_out["primal"]))) 

164 if verbosity > 2: 

165 result = gp.generate_result(solver_out, verbosity=verbosity-3) 

166 self._results.append(result) 

167 print(result.table(self.sgpvks)) 

168 elif verbosity > 1: 

169 print("Solved cost was %.4g." % cost) 

170 if prevcost is None: 

171 continue 

172 rel_improvement = (prevcost - cost)/(prevcost + cost) 

173 if cost/prevcost >= 1 + 10*EPS: 

174 pywarnings.warn( 

175 "SGP not convergent: Cost rose by %.2g%% (%.6g to %.6g) on" 

176 " GP solve %i. Details can be found in `m.program.results`" 

177 " or by solving at a higher verbosity. Note convergence" 

178 " is not guaranteed for models with SignomialEqualities." 

179 % (100*(cost - prevcost)/prevcost, prevcost, cost, len(self.gps))) 

180 rel_improvement = cost = None 

181 # solved successfully! 

182 self.result = gp.generate_result(solver_out, verbosity=verbosity-3) 

183 self.result["soltime"] = time() - starttime 

184 if verbosity > 1: 

185 print() 

186 if verbosity > 0: 

187 print("Solving took %.3g seconds and %i GP solves." 

188 % (self.result["soltime"], len(self.gps))) 

189 if hasattr(self.slack, "key"): 

190 initsolwarning(self.result, "Slack Non-GP Constraints") 

191 excess_slack = self.result["variables"][self.slack.key] - 1 # pylint: disable=no-member 

192 if excess_slack > EPS: 

193 msg = ("Final PCCP solution let non-GP constraints slacken by" 

194 " %.2g%%." % (100*excess_slack)) 

195 appendsolwarning(msg, (1 + excess_slack), self.result, 

196 "Slack Non-GP Constraints") 

197 if verbosity > -1: 

198 print(msg + 

199 " Calling .localsolve(pccp_penalty=...) with a" 

200 " higher `pccp_penalty` (it was %.3g this time) will" 

201 " reduce slack if the model is solvable with less. To" 

202 " verify that the slack is needed, generate an SGP with" 

203 " `use_pccp=False` and start it from this model's" 

204 " solution: e.g. `m.localsolve(use_pccp=False, x0=" 

205 "m.solution[\"variables\"])`." % self.pccp_penalty) 

206 del self.result["freevariables"][self.slack.key] # pylint: disable=no-member 

207 del self.result["variables"][self.slack.key] # pylint: disable=no-member 

208 del self.result["sensitivities"]["variables"][self.slack.key] # pylint: disable=no-member 

209 slcon = self.gpconstraints[0] 

210 slconsenss = self.result["sensitivities"]["constraints"][slcon] 

211 del self.result["sensitivities"]["constraints"][slcon] 

212 # TODO: create constraint in RelaxPCCP namespace 

213 self.result["sensitivities"]["models"][""] -= slconsenss 

214 if not self.result["sensitivities"]["models"][""]: 

215 del self.result["sensitivities"]["models"][""] 

216 return self.result 

217 

218 @property 

219 def results(self): 

220 "Creates and caches results from the raw solver_outs" 

221 if not self._results: 

222 self._results = [gp.generate_result(s_o, dual_check=False) 

223 for gp, s_o in zip(self.gps, self.solver_outs)] 

224 return self._results 

225 

226 def gp(self, x0=None, *, cleanx0=False): 

227 "Update self._gp for x0 and return it." 

228 if not x0: 

229 return self._gp # return last generated 

230 if not cleanx0: 

231 cleanedx0 = KeyDict() 

232 cleanedx0.vks = self._gp.x0.vks 

233 cleanedx0.update(x0) 

234 x0 = cleanedx0 

235 self._gp.x0.update({vk: x0[vk] for vk in self.sgpvks if vk in x0}) 

236 p_idx = 0 

237 for sgpc in self.sgpconstraints: 

238 for hmaplt1 in sgpc.as_gpconstr(self._gp.x0).as_hmapslt1({}): 

239 approxc = self.approxconstraints[p_idx] 

240 approxc.unsubbed = [Posynomial(hmaplt1)/self.slack] 

241 p_idx += 1 # p_idx=0 is the cost; sp constraints are after it 

242 hmap, = approxc.as_hmapslt1(self._gp.substitutions) 

243 self._gp.hmaps[p_idx] = hmap 

244 m_idx = self._gp.m_idxs[p_idx].start 

245 a_idxs = list(self.a_idxs[p_idx]) # A's entries we can modify 

246 for i, (exp, c) in enumerate(hmap.items()): 

247 self._gp.exps[m_idx + i] = exp 

248 self._gp.cs[m_idx + i] = c 

249 for var, x in exp.items(): 

250 try: # modify a particular A entry 

251 row_idx = a_idxs.pop() 

252 self._gp.A.row[row_idx] = m_idx + i 

253 self._gp.A.col[row_idx] = self._gp.varidxs[var] 

254 self._gp.A.data[row_idx] = x 

255 except IndexError: # numbers of exps increased 

256 self.a_idxs[p_idx].append(len(self._gp.A.row)) 

257 self._gp.A.row.append(m_idx + i) 

258 self._gp.A.col.append(self._gp.varidxs[var]) 

259 self._gp.A.data.append(x) 

260 for row_idx in a_idxs: # number of exps decreased 

261 self._gp.A.row[row_idx] = 0 # zero out this entry 

262 self._gp.A.col[row_idx] = 0 

263 self._gp.A.data[row_idx] = 0 

264 return self._gp