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 

13 

14 

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

16 

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

18class SequentialGeometricProgram: 

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

20 

21 Arguments 

22 --------- 

23 cost : Posynomial 

24 Objective to minimize when solving 

25 constraints : list of Constraint or SignomialConstraint objects 

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

27 verbosity : int (optional) 

28 Currently has no effect: SequentialGeometricPrograms don't know 

29 anything new after being created, unlike GeometricPrograms. 

30 

31 Attributes with side effects 

32 ---------------------------- 

33 `gps` is set during a solve 

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

35 

36 Examples 

37 -------- 

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

39 # minimize 

40 x, 

41 [ # subject to 

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

43 y/10 # <= 1 

44 ]) 

45 >>> gp.solve() 

46 """ 

47 gps = solver_outs = _results = result = model = None 

48 with NamedVariables("RelaxPCCP"): 

49 slack = Variable("C") 

50 

51 def __init__(self, cost, model, substitutions, *, 

52 use_pccp=True, pccp_penalty=2e2, checkbounds=True): 

53 self.pccp_penalty = pccp_penalty 

54 if cost.any_nonpositive_cs: 

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

56 

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

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

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

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

61 if not use_pccp: 

62 self.slack = 1 

63 else: 

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

65 cost *= self.slack**pccp_penalty 

66 self.approxconstraints = [] 

67 self.sgpvks = set() 

68 x0 = KeyDict(substitutions) 

69 x0.varkeys = model.varkeys # for string access and so forth 

70 for cs in model.flat(): 

71 try: 

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

73 raise InvalidGPConstraint(cs) 

74 if not isinstance(cs, PosynomialInequality): 

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

76 self.gpconstraints.append(cs) 

77 except InvalidGPConstraint: 

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

79 raise InvalidSGPConstraint(cs) 

80 self.sgpconstraints.append(cs) 

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

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

83 constraint.generated_by = cs 

84 self.approxconstraints.append(constraint) 

85 self.sgpvks.update(constraint.vks) 

86 if not self.sgpconstraints: 

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

88 

89SequentialGeometricPrograms should only be created with Models containing 

90Signomial Constraints, since Models without Signomials have global 

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

92 self._gp = GeometricProgram( 

93 cost, self.approxconstraints + self.gpconstraints, 

94 substitutions, checkbounds=checkbounds) 

95 self._gp.x0 = x0 

96 self.a_idxs = defaultdict(list) 

97 cost_mons = self._gp.k[0] 

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

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

100 if cost_mons <= m_idx <= sp_mons: 

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

102 

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

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

105 iteration_limit=50, **solveargs): 

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

107 

108 Arguments 

109 --------- 

110 solver : str or function (optional) 

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

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

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

114 verbosity : int (optional) 

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

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

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

118 x0 : dict (optional) 

119 Initial location to approximate signomials about. 

120 reltol : float 

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

122 consecutive solve's objective values. 

123 iteration_limit : int 

124 Maximum GP iterations allowed. 

125 mutategp: boolean 

126 Prescribes whether to mutate the previously generated GP 

127 or to create a new GP with every solve. 

128 **solveargs : 

129 Passed to solver function. 

130 

131 Returns 

132 ------- 

133 result : dict 

134 A dictionary containing the translated solver result. 

135 """ 

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

137 starttime = time() 

138 if verbosity > 0: 

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

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

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

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

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

144 prevcost, cost, rel_improvement = None, None, None 

145 while rel_improvement is None or rel_improvement > reltol: 

146 prevcost = cost 

147 if len(self.gps) > iteration_limit: 

148 raise Infeasible( 

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

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

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

152 gp = self.gp(x0, cleanx0=True) 

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

154 if verbosity > 1: 

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

156 if verbosity > 2: 

157 print("===============") 

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

159 gen_result=False, **solveargs) 

160 self.solver_outs.append(solver_out) 

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

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

163 if verbosity > 2: 

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

165 self._results.append(result) 

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

167 elif verbosity > 1: 

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

169 if prevcost is None: 

170 continue 

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

172 if cost*(1 - EPS) > prevcost + EPS and verbosity > -1: 

173 pywarnings.warn( 

174 "SGP not convergent: Cost rose by %.2g%% on GP solve %i." 

175 " Details can be found in `m.program.results` or by" 

176 " solving at a higher verbosity. Note that convergence" 

177 " is not guaranteed for models with SignomialEqualities." 

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

179 rel_improvement = cost = None 

180 # solved successfully! 

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

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

183 if verbosity > 1: 

184 print() 

185 if verbosity > 0: 

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

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

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

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

190 if excess_slack <= EPS: 

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

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

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

194 slackconstraint = self.gpconstraints[0] 

195 del self.result["sensitivities"]["constraints"][slackconstraint] 

196 elif verbosity > -1: 

197 pywarnings.warn( 

198 "Final solution let signomial constraints slacken by" 

199 " %.2g%%. Calling .localsolve with a higher" 

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

201 " final slack if the model is solvable with less. If" 

202 " you think it might not be, check by solving with " 

203 "`use_pccp=False, x0=(this model's final solution)`.\n" 

204 % (100*excess_slack, self.pccp_penalty)) 

205 return self.result 

206 

207 @property 

208 def results(self): 

209 "Creates and caches results from the raw solver_outs" 

210 if not self._results: 

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

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

213 return self._results 

214 

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

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

217 if not x0: 

218 return self._gp # return last generated 

219 if not cleanx0: 

220 x0 = KeyDict(x0) 

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

222 p_idx = 0 

223 for sgpc in self.sgpconstraints: 

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

225 approxc = self.approxconstraints[p_idx] 

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

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

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

229 self._gp.hmaps[p_idx] = hmap 

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

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

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

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

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

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

236 try: # modify a particular A entry 

237 row_idx = a_idxs.pop() 

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

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

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

241 except IndexError: # numbers of exps increased 

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

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

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

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

246 for row_idx in a_idxs: # number of exps decreased 

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

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

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

250 return self._gp