1"Implements Model" 

2import numpy as np 

3from .costed import CostedConstraintSet 

4from ..nomials import Monomial 

5from .prog_factories import progify, solvify 

6from .gp import GeometricProgram 

7from .sgp import SequentialGeometricProgram 

8from ..small_scripts import mag 

9from ..tools.autosweep import autosweep_1d 

10from ..exceptions import InvalidGPConstraint 

11from .. import NamedVariables 

12from ..tools.docstring import expected_unbounded 

13from .set import add_meq_bounds 

14from ..exceptions import Infeasible 

15 

16 

17class Model(CostedConstraintSet): 

18 """Symbolic representation of an optimization problem. 

19 

20 The Model class is used both directly to create models with constants and 

21 sweeps, and indirectly inherited to create custom model classes. 

22 

23 Arguments 

24 --------- 

25 cost : Posynomial (optional) 

26 Defaults to `Monomial(1)`. 

27 

28 constraints : ConstraintSet or list of constraints (optional) 

29 Defaults to an empty list. 

30 

31 substitutions : dict (optional) 

32 This dictionary will be substituted into the problem before solving, 

33 and also allows the declaration of sweeps and linked sweeps. 

34 

35 Attributes with side effects 

36 ---------------------------- 

37 `program` is set during a solve 

38 `solution` is set at the end of a solve 

39 """ 

40 program = None 

41 solution = None 

42 

43 def __init__(self, cost=None, constraints=None, *args, **kwargs): 

44 setup_vars = None 

45 substitutions = kwargs.pop("substitutions", None) # reserved keyword 

46 if hasattr(self, "setup"): 

47 self.cost = None 

48 # lineage holds the (name, num) environment a model was created in, 

49 # including its own (name, num), and those of models above it 

50 with NamedVariables(self.__class__.__name__) as (self.lineage, 

51 setup_vars): 

52 args = tuple(arg for arg in [cost, constraints] 

53 if arg is not None) + args 

54 cs = self.setup(*args, **kwargs) # pylint: disable=no-member 

55 if (isinstance(cs, tuple) and len(cs) == 2 

56 and isinstance(cs[1], dict)): 

57 constraints, substitutions = cs 

58 else: 

59 constraints = cs 

60 cost = self.cost 

61 elif args and not substitutions: 

62 # backwards compatibility: substitutions as third argument 

63 substitutions, = args 

64 

65 cost = cost or Monomial(1) 

66 constraints = constraints or [] 

67 if setup_vars: 

68 # add all the vars created in .setup to the Model's varkeys 

69 # even if they aren't used in any constraints 

70 self.unique_varkeys = frozenset(v.key for v in setup_vars) 

71 CostedConstraintSet.__init__(self, cost, constraints, substitutions) 

72 docstr = self.__class__.__doc__ 

73 if self.lineage and docstr and "SKIP VERIFICATION" not in docstr: 

74 if "Unbounded" in docstr or "Bounded by" in docstr: 

75 self.verify_docstring() 

76 

77 gp = progify(GeometricProgram) 

78 solve = solvify(progify(GeometricProgram, "solve")) 

79 

80 sp = progify(SequentialGeometricProgram) 

81 localsolve = solvify(progify(SequentialGeometricProgram, "localsolve")) 

82 

83 def verify_docstring(self): # pylint:disable=too-many-locals,too-many-branches,too-many-statements 

84 "Verifies docstring bounds are sufficient but not excessive." 

85 err = "while verifying %s:\n" % self.__class__.__name__ 

86 bounded, meq_bounded = self.bounded.copy(), self.meq_bounded.copy() 

87 doc = self.__class__.__doc__ 

88 exp_unbounds = expected_unbounded(self, doc) 

89 unexp_bounds = bounded.intersection(exp_unbounds) 

90 if unexp_bounds: # anything bounded that shouldn't be? err! 

91 for direction in ["lower", "upper"]: 

92 badvks = [v for v, d in unexp_bounds if d == direction] 

93 if not badvks: 

94 continue 

95 badvks = ", ".join(str(v) for v in badvks) 

96 badvks += (" were" if len(badvks) > 1 else " was") 

97 err += (" %s %s-bounded; expected %s-unbounded" 

98 "\n" % (badvks, direction, direction)) 

99 raise ValueError(err) 

100 bounded.update(exp_unbounds) # if not, treat expected as bounded 

101 add_meq_bounds(bounded, meq_bounded) # and add more meqs 

102 self.missingbounds = {} # now let's figure out what's missing 

103 for bound in meq_bounded: # first add the un-dealt-with meq bounds 

104 for condition in list(meq_bounded[bound]): 

105 meq_bounded[bound].remove(condition) 

106 newcond = condition - bounded 

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

108 for c in meq_bounded[bound]): 

109 meq_bounded[bound].add(newcond) 

110 bsets = " or ".join(str(list(c)) for c in meq_bounded[bound]) 

111 self.missingbounds[bound] = (", but would gain it from any of" 

112 " these sets of bounds: " + bsets) 

113 # then add everything that's not in bounded 

114 if len(bounded)+len(self.missingbounds) != 2*len(self.varkeys): 

115 for key in self.varkeys: 

116 for bound in ("upper", "lower"): 

117 if (key, bound) not in bounded: 

118 if (key, bound) not in self.missingbounds: 

119 self.missingbounds[(key, bound)] = "" 

120 if self.missingbounds: # anything unbounded? err! 

121 boundstrs = "\n".join(" %s has no %s bound%s" % (v, b, x) 

122 for (v, b), x 

123 in self.missingbounds.items()) 

124 docstring = ("To fix this add the following to %s's" 

125 " docstring (you may not need it all):" 

126 " \n" % self.__class__.__name__) 

127 for direction in ["upper", "lower"]: 

128 mb = [k for (k, b) in self.missingbounds if b == direction] 

129 if mb: 

130 docstring += """ 

131%s Unbounded 

132--------------- 

133%s 

134""" % (direction.title(), ", ".join(set(k.name for k in mb))) 

135 raise ValueError(err + boundstrs + "\n\n" + docstring) 

136 

137 def sweep(self, sweeps, **solveargs): 

138 "Sweeps {var: values} pairs in sweeps. Returns swept solutions." 

139 sols = [] 

140 for sweepvar, sweepvals in sweeps.items(): 

141 original_val = self.substitutions.get(sweepvar, None) 

142 self.substitutions.update({sweepvar: ('sweep', sweepvals)}) 

143 try: 

144 sols.append(self.solve(**solveargs)) 

145 except InvalidGPConstraint: 

146 sols.append(self.localsolve(**solveargs)) 

147 if original_val: 

148 self.substitutions[sweepvar] = original_val 

149 else: 

150 del self.substitutions[sweepvar] 

151 return sols if len(sols) > 1 else sols[0] 

152 

153 def autosweep(self, sweeps, tol=0.01, samplepoints=100, **solveargs): 

154 """Autosweeps {var: (start, end)} pairs in sweeps to tol. 

155 

156 Returns swept and sampled solutions. 

157 The original simplex tree can be accessed at sol.bst 

158 """ 

159 sols = [] 

160 for sweepvar, sweepvals in sweeps.items(): 

161 sweepvar = self[sweepvar].key 

162 start, end = sweepvals 

163 bst = autosweep_1d(self, tol, sweepvar, [start, end], **solveargs) 

164 sols.append(bst.sample_at(np.linspace(start, end, samplepoints))) 

165 return sols if len(sols) > 1 else sols[0] 

166 

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

168 def debug(self, solver=None, verbosity=1, **solveargs): 

169 "Attempts to diagnose infeasible models." 

170 from .relax import ConstantsRelaxed, ConstraintsRelaxed 

171 from .bounded import Bounded 

172 

173 sol = None 

174 solveargs["solver"] = solver 

175 solveargs["verbosity"] = verbosity - 1 

176 solveargs["process_result"] = False 

177 

178 if verbosity: 

179 print("< DEBUGGING >") 

180 print("> Trying with bounded variables and relaxed constants:") 

181 

182 bounded = Bounded(self) 

183 if self.substitutions: 

184 tants = ConstantsRelaxed(bounded) 

185 feas = Model(tants.relaxvars.prod()**30 * self.cost, tants) 

186 else: 

187 feas = Model(self.cost, bounded) 

188 

189 try: 

190 try: 

191 sol = feas.solve(**solveargs) 

192 except InvalidGPConstraint: 

193 sol = feas.sp(use_pccp=False).localsolve(**solveargs) 

194 sol["boundedness"] = bounded.check_boundaries(sol, 

195 verbosity=verbosity) 

196 if self.substitutions: 

197 relaxed = get_relaxed([sol(r) for r in tants.relaxvars], 

198 tants.freedvars, 

199 min_return=0 if sol["boundedness"] else 1) 

200 if verbosity and relaxed: 

201 if sol["boundedness"]: 

202 print("and these constants relaxed:") 

203 else: 

204 print("\nSolves with these constants relaxed:") 

205 for (_, freed) in relaxed: 

206 print(" %s: relaxed from %-.4g to %-.4g" 

207 % (freed, mag(tants.constants[freed.key]), 

208 mag(sol(freed)))) 

209 print("") 

210 if verbosity > 0: 

211 print(">> Success!") 

212 except Infeasible: 

213 if verbosity > 0: 

214 print(">> Failure.") 

215 print("> Trying with relaxed constraints:") 

216 

217 try: 

218 traints = ConstraintsRelaxed(self) 

219 feas = Model(traints.relaxvars.prod()**30 * self.cost, traints) 

220 try: 

221 sol = feas.solve(**solveargs) 

222 except InvalidGPConstraint: 

223 sol = feas.sp(use_pccp=False).localsolve(**solveargs) 

224 relaxed_constraints = feas[0]["relaxed constraints"] 

225 relaxed = get_relaxed(sol(traints.relaxvars), 

226 range(len(relaxed_constraints))) 

227 if verbosity > 0 and relaxed: 

228 print("\nSolves with these constraints relaxed:") 

229 for relaxval, i in relaxed: 

230 relax_percent = "%i%%" % (0.5+(relaxval-1)*100) 

231 oldconstraint = traints.original_constraints[i] 

232 newconstraint = relaxed_constraints[i][0] 

233 subs = {traints.relaxvars[i]: relaxval} 

234 relaxdleft = newconstraint.left.sub(subs) 

235 relaxdright = newconstraint.right.sub(subs) 

236 print(" %3i: %5s relaxed, from %s %s %s \n" 

237 " to %s %s %s " 

238 % (i, relax_percent, oldconstraint.left, 

239 oldconstraint.oper, oldconstraint.right, 

240 relaxdleft, newconstraint.oper, relaxdright)) 

241 if verbosity > 0: 

242 print("\n>> Success!\n") 

243 except (ValueError, RuntimeWarning): 

244 if verbosity > 0: 

245 print(">> Failure\n") 

246 return sol 

247 

248 

249def get_relaxed(relaxvals, mapped_list, min_return=1): 

250 "Determines which relaxvars are considered 'relaxed'" 

251 sortrelaxed = sorted(zip(relaxvals, mapped_list), key=lambda x: -x[0]) 

252 # arbitrarily, 1.01 is the threshold below which we don't show slack 

253 mostrelaxed = max(sortrelaxed[0][0], 1.01) 

254 for i, (val, _) in enumerate(sortrelaxed): 

255 if i >= min_return and val <= 1.01 and (val-1) <= (mostrelaxed-1)/10: 

256 return sortrelaxed[:i] 

257 return sortrelaxed