Coverage for gpkit/constraints/model.py: 77%
142 statements
« prev ^ index » next coverage.py v7.4.0, created at 2024-01-03 16:49 -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 .set import add_meq_bounds
13from ..exceptions import Infeasible
16class Model(CostedConstraintSet):
17 """Symbolic representation of an optimization problem.
19 The Model class is used both directly to create models with constants and
20 sweeps, and indirectly inherited to create custom model classes.
22 Arguments
23 ---------
24 cost : Posynomial (optional)
25 Defaults to `Monomial(1)`.
27 constraints : ConstraintSet or list of constraints (optional)
28 Defaults to an empty list.
30 substitutions : dict (optional)
31 This dictionary will be substituted into the problem before solving,
32 and also allows the declaration of sweeps and linked sweeps.
34 Attributes with side effects
35 ----------------------------
36 `program` is set during a solve
37 `solution` is set at the end of a solve
38 """
39 program = None
40 solution = None
42 def __init__(self, cost=None, constraints=None, *args, **kwargs):
43 setup_vars = None
44 substitutions = kwargs.pop("substitutions", None) # reserved keyword
45 if hasattr(self, "setup"):
46 self.cost = None
47 # lineage holds the (name, num) environment a model was created in,
48 # including its own (name, num), and those of models above it
49 with NamedVariables(self.__class__.__name__) as (self.lineage,
50 setup_vars):
51 args = tuple(arg for arg in [cost, constraints]
52 if arg is not None) + args
53 cs = self.setup(*args, **kwargs) # pylint: disable=no-member
54 if (isinstance(cs, tuple) and len(cs) == 2
55 and isinstance(cs[1], dict)):
56 constraints, substitutions = cs
57 else:
58 constraints = cs
59 cost = self.cost
60 elif args and not substitutions:
61 # backwards compatibility: substitutions as third argument
62 substitutions, = args
64 cost = cost or Monomial(1)
65 constraints = constraints or []
66 if setup_vars:
67 # add all the vars created in .setup to the Model's varkeys
68 # even if they aren't used in any constraints
69 self.unique_varkeys = frozenset(v.key for v in setup_vars)
70 CostedConstraintSet.__init__(self, cost, constraints, substitutions)
71 docstr = self.__class__.__doc__
72 if self.lineage and docstr and "SKIP VERIFICATION" not in docstr:
73 if "Unbounded" in docstr or "Bounded by" in docstr:
74 self.verify_docstring()
76 gp = progify(GeometricProgram)
77 solve = solvify(progify(GeometricProgram, "solve"))
79 sp = progify(SequentialGeometricProgram)
80 localsolve = solvify(progify(SequentialGeometricProgram, "localsolve"))
82 def verify_docstring(self): # pylint:disable=too-many-locals,too-many-branches,too-many-statements
83 "Verifies docstring bounds are sufficient but not excessive."
84 err = "while verifying %s:\n" % self.__class__.__name__
85 bounded, meq_bounded = self.bounded.copy(), self.meq_bounded.copy()
86 doc = self.__class__.__doc__
87 exp_unbounds = expected_unbounded(self, doc)
88 unexp_bounds = bounded.intersection(exp_unbounds)
89 if unexp_bounds: # anything bounded that shouldn't be? err!
90 for direction in ["lower", "upper"]:
91 badvks = [v for v, d in unexp_bounds if d == direction]
92 if not badvks:
93 continue
94 badvks = ", ".join(str(v) for v in badvks)
95 badvks += (" were" if len(badvks) > 1 else " was")
96 err += (" %s %s-bounded; expected %s-unbounded"
97 "\n" % (badvks, direction, direction))
98 raise ValueError(err)
99 bounded.update(exp_unbounds) # if not, treat expected as bounded
100 add_meq_bounds(bounded, meq_bounded) # and add more meqs
101 self.missingbounds = {} # now let's figure out what's missing
102 for bound in meq_bounded: # first add the un-dealt-with meq bounds
103 for condition in list(meq_bounded[bound]):
104 meq_bounded[bound].remove(condition)
105 newcond = condition - bounded
106 if newcond and not any(c.issubset(newcond)
107 for c in meq_bounded[bound]):
108 meq_bounded[bound].add(newcond)
109 bsets = " or ".join(str(list(c)) for c in meq_bounded[bound])
110 self.missingbounds[bound] = (", but would gain it from any of"
111 " these sets of bounds: " + bsets)
112 # then add everything that's not in bounded
113 if len(bounded)+len(self.missingbounds) != 2*len(self.varkeys):
114 for key in self.varkeys:
115 for bound in ("upper", "lower"):
116 if (key, bound) not in bounded:
117 if (key, bound) not in self.missingbounds:
118 self.missingbounds[(key, bound)] = ""
119 if self.missingbounds: # anything unbounded? err!
120 boundstrs = "\n".join(" %s has no %s bound%s" % (v, b, x)
121 for (v, b), x
122 in self.missingbounds.items())
123 docstring = ("To fix this add the following to %s's"
124 " docstring (you may not need it all):"
125 " \n" % self.__class__.__name__)
126 for direction in ["upper", "lower"]:
127 mb = [k for (k, b) in self.missingbounds if b == direction]
128 if mb:
129 docstring += """
130%s Unbounded
131---------------
132%s
133""" % (direction.title(), ", ".join(set(k.name for k in mb)))
134 raise ValueError(err + boundstrs + "\n\n" + docstring)
136 def sweep(self, sweeps, **solveargs):
137 "Sweeps {var: values} pairs in sweeps. Returns swept solutions."
138 sols = []
139 for sweepvar, sweepvals in sweeps.items():
140 original_val = self.substitutions.get(sweepvar, None)
141 self.substitutions.update({sweepvar: ('sweep', sweepvals)})
142 try:
143 sols.append(self.solve(**solveargs))
144 except InvalidGPConstraint:
145 sols.append(self.localsolve(**solveargs))
146 if original_val:
147 self.substitutions[sweepvar] = original_val
148 else:
149 del self.substitutions[sweepvar]
150 return sols if len(sols) > 1 else sols[0]
152 def autosweep(self, sweeps, tol=0.01, samplepoints=100, **solveargs):
153 """Autosweeps {var: (start, end)} pairs in sweeps to tol.
155 Returns swept and sampled solutions.
156 The original simplex tree can be accessed at sol.bst
157 """
158 sols = []
159 for sweepvar, sweepvals in sweeps.items():
160 sweepvar = self[sweepvar].key
161 start, end = sweepvals
162 bst = autosweep_1d(self, tol, sweepvar, [start, end], **solveargs)
163 sols.append(bst.sample_at(np.linspace(start, end, samplepoints)))
164 return sols if len(sols) > 1 else sols[0]
166 def debug(self, solver=None, verbosity=1, **solveargs):
167 "Attempts to diagnose infeasible models."
168 from .relax import ConstantsRelaxed, ConstraintsRelaxed
169 from .bounded import Bounded
171 sol = None
172 solveargs["solver"] = solver
173 solveargs["verbosity"] = verbosity - 1
174 solveargs["process_result"] = False
176 bounded = Bounded(self)
177 tants = ConstantsRelaxed(bounded)
178 if tants.relaxvars.size:
179 feas = Model(tants.relaxvars.prod()**30 * self.cost, tants)
180 else:
181 feas = Model(self.cost, bounded)
183 try:
184 try:
185 sol = feas.solve(**solveargs)
186 except InvalidGPConstraint:
187 sol = feas.sp(use_pccp=False).localsolve(**solveargs)
188 # limited results processing
189 bounded.check_boundaries(sol)
190 tants.check_relaxed(sol)
191 except Infeasible:
192 if verbosity:
193 print("<DEBUG> Model is not feasible with relaxed constants"
194 " and bounded variables.")
195 traints = ConstraintsRelaxed(self)
196 feas = Model(traints.relaxvars.prod()**30 * self.cost, traints)
197 try:
198 try:
199 sol = feas.solve(**solveargs)
200 except InvalidGPConstraint:
201 sol = feas.sp(use_pccp=False).localsolve(**solveargs)
202 # limited results processing
203 traints.check_relaxed(sol)
204 except Infeasible:
205 print("<DEBUG> Model is not feasible with bounded constraints.")
206 if sol and verbosity:
207 warnings = sol.table(tables=["warnings"]).split("\n")[3:-2]
208 if warnings:
209 print("<DEBUG> Model is feasible with these modifications:")
210 print("\n" + "\n".join(warnings) + "\n")
211 else:
212 print("<DEBUG> Model seems feasible without modification,"
213 " or only needs relaxations of less than 1%."
214 " Check the returned solution for details.")
215 return sol