Coverage for gpkit/constraints/model.py: 77%

142 statements  

« prev     ^ index     » next       coverage.py v7.4.0, created at 2024-01-07 22:56 -0500

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 ..tools.autosweep import autosweep_1d 

9from ..exceptions import InvalidGPConstraint 

10from .. import NamedVariables 

11from ..tools.docstring import expected_unbounded 

12from .relax import ConstantsRelaxed, ConstraintsRelaxed 

13from .bounded import Bounded 

14from .set import add_meq_bounds 

15from ..exceptions import Infeasible 

16 

17 

18class Model(CostedConstraintSet): 

19 """Symbolic representation of an optimization problem. 

20 

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

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

23 

24 Arguments 

25 --------- 

26 cost : Posynomial (optional) 

27 Defaults to `Monomial(1)`. 

28 

29 constraints : ConstraintSet or list of constraints (optional) 

30 Defaults to an empty list. 

31 

32 substitutions : dict (optional) 

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

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

35 

36 Attributes with side effects 

37 ---------------------------- 

38 `program` is set during a solve 

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

40 """ 

41 program = None 

42 solution = None 

43 

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

45 #pylint: disable=keyword-arg-before-vararg 

46 setup_vars = None 

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

48 if hasattr(self, "setup"): 

49 self.cost = None 

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

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

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

53 setup_vars): 

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

55 if arg is not None) + args 

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

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

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

59 constraints, substitutions = cs 

60 else: 

61 constraints = cs 

62 cost = self.cost 

63 elif args and not substitutions: 

64 # backwards compatibility: substitutions as third argument 

65 substitutions, = args 

66 

67 cost = cost or Monomial(1) 

68 constraints = constraints or [] 

69 if setup_vars: 

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

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

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

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

74 docstr = self.__class__.__doc__ 

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

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

77 self.verify_docstring() 

78 

79 gp = progify(GeometricProgram) 

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

81 

82 sp = progify(SequentialGeometricProgram) 

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

84 

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

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

87 err = f"while verifying {self.__class__.__name__}:\n" 

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

89 doc = self.__class__.__doc__ 

90 exp_unbounds = expected_unbounded(self, doc) 

91 unexp_bounds = bounded.intersection(exp_unbounds) 

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

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

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

95 if not badvks: 

96 continue 

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

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

99 err += (f" {badvks} {direction}-bounded; expected " 

100 f"{direction}-unbounded\n") 

101 raise ValueError(err) 

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

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

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

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

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

107 meq_bounded[bound].remove(condition) 

108 newcond = condition - bounded 

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

110 for c in meq_bounded[bound]): 

111 meq_bounded[bound].add(newcond) 

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

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

114 " these sets of bounds: " + bsets) 

115 # then add everything that's not in bounded 

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

117 for key in self.varkeys: 

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

119 if (key, bound) not in bounded: 

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

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

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

123 boundstrs = "\n".join(f" {v} has no {b} bound{x}" 

124 for (v, b), x 

125 in self.missingbounds.items()) 

126 docstring = ("To fix this add the following to " 

127 f"{self.__class__.__name__}'s docstring " 

128 "(you may not need it all):\n") 

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

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

131 if mb: 

132 #pylint: disable=consider-using-f-string 

133 docstring += """ 

134%s Unbounded 

135--------------- 

136%s 

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

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

139 

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

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

142 sols = [] 

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

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

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

146 try: 

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

148 except InvalidGPConstraint: 

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

150 if original_val: 

151 self.substitutions[sweepvar] = original_val 

152 else: 

153 del self.substitutions[sweepvar] 

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

155 

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

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

158 

159 Returns swept and sampled solutions. 

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

161 """ 

162 sols = [] 

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

164 sweepvar = self[sweepvar].key 

165 start, end = sweepvals 

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

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

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

169 

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

171 "Attempts to diagnose infeasible models." 

172 

173 sol = None 

174 solveargs["solver"] = solver 

175 solveargs["verbosity"] = verbosity - 1 

176 solveargs["process_result"] = False 

177 

178 bounded = Bounded(self) 

179 tants = ConstantsRelaxed(bounded) 

180 if tants.relaxvars.size: 

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

182 else: 

183 feas = Model(self.cost, bounded) 

184 

185 try: 

186 try: 

187 sol = feas.solve(**solveargs) 

188 except InvalidGPConstraint: 

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

190 # limited results processing 

191 bounded.check_boundaries(sol) 

192 tants.check_relaxed(sol) 

193 except Infeasible: 

194 if verbosity: 

195 print("<DEBUG> Model is not feasible with relaxed constants" 

196 " and bounded variables.") 

197 traints = ConstraintsRelaxed(self) 

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

199 try: 

200 try: 

201 sol = feas.solve(**solveargs) 

202 except InvalidGPConstraint: 

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

204 # limited results processing 

205 traints.check_relaxed(sol) 

206 except Infeasible: 

207 print("<DEBUG> Model is not feasible with bounded constraints.") 

208 if sol and verbosity: 

209 warnings = sol.table(tables=["warnings"]).split("\n")[3:-2] 

210 if warnings: 

211 print("<DEBUG> Model is feasible with these modifications:") 

212 print("\n" + "\n".join(warnings) + "\n") 

213 else: 

214 print("<DEBUG> Model seems feasible without modification," 

215 " or only needs relaxations of less than 1%." 

216 " Check the returned solution for details.") 

217 return sol