1"""Test substitution capability across gpkit""" 

2import pickle 

3import unittest 

4import numpy as np 

5import numpy.testing as npt 

6from ad import adnumber, ADV 

7import gpkit 

8from gpkit import SignomialsEnabled, NamedVariables 

9from gpkit import Variable, VectorVariable, Model, Signomial 

10from gpkit.small_scripts import mag 

11from gpkit.tests.helpers import run_tests 

12from gpkit.exceptions import UnboundedGP, DimensionalityError 

13 

14# pylint: disable=invalid-name,attribute-defined-outside-init,unused-variable 

15 

16 

17class TestNomialSubs(unittest.TestCase): 

18 """Test substitution for nomial-family objects""" 

19 

20 def test_vectorized_linked(self): 

21 class VectorLinked(Model): 

22 "simple vectorized link" 

23 def setup(self): 

24 self.y = y = Variable("y", 1) 

25 

26 def vectorlink(c): 

27 "linked vector function" 

28 if isinstance(c[y], ADV): 

29 return np.array(c[y])+adnumber([1, 2, 3]) 

30 return c[y]+np.array([1, 2, 3]) 

31 self.x = x = VectorVariable(3, "x", vectorlink) 

32 m = VectorLinked() 

33 self.assertEqual(m.substitutions[m.x[0].key](m.substitutions), 2) 

34 self.assertEqual(m.gp().substitutions[m.x[0].key], 2) 

35 self.assertEqual(m.gp().substitutions[m.x[1].key], 3) 

36 self.assertEqual(m.gp().substitutions[m.x[2].key], 4) 

37 

38 def test_numeric(self): 

39 """Basic substitution of numeric value""" 

40 x = Variable("x") 

41 p = x**2 

42 self.assertEqual(p.sub({x: 3}), 9) 

43 self.assertEqual(p.sub({x.key: 3}), 9) 

44 self.assertEqual(p.sub({"x": 3}), 9) 

45 

46 def test_dimensionless_units(self): 

47 x = Variable('x', 3, 'ft') 

48 y = Variable('y', 1, 'm') 

49 if x.units is not None: 

50 # units are enabled 

51 self.assertAlmostEqual((x/y).value, 0.9144) 

52 

53 def test_vector(self): 

54 x = Variable("x") 

55 y = Variable("y") 

56 z = VectorVariable(2, "z") 

57 p = x*y*z 

58 self.assertTrue(all(p.sub({x: 1, "y": 2}) == 2*z)) 

59 self.assertTrue(all(p.sub({x: 1, y: 2, "z": [1, 2]}) == 

60 z.sub({z: [2, 4]}))) 

61 self.assertRaises(ValueError, z.sub, {z: [1, 2, 3]}) 

62 

63 xvec = VectorVariable(3, "x", "m") 

64 xs = xvec[:2].sum() 

65 for x_ in ["x", xvec]: 

66 self.assertAlmostEqual(mag(xs.sub({x_: [1, 2, 3]}).c), 3.0) 

67 

68 def test_variable(self): 

69 """Test special single-argument substitution for Variable""" 

70 x = Variable('x') 

71 y = Variable('y') 

72 m = x*y**2 

73 self.assertEqual(x.sub(3), 3) 

74 # make sure x was not mutated 

75 self.assertEqual(x, Variable('x')) 

76 self.assertNotEqual(x.sub(3), Variable('x')) 

77 # also make sure the old way works 

78 self.assertEqual(x.sub({x: 3}), 3) 

79 # and for vectors 

80 xvec = VectorVariable(3, 'x') 

81 self.assertEqual(xvec[1].sub(3), 3) 

82 

83 def test_signomial(self): 

84 """Test Signomial substitution""" 

85 D = Variable('D', units="N") 

86 x = Variable('x', units="N") 

87 y = Variable('y', units="N") 

88 a = Variable('a') 

89 with SignomialsEnabled(): 

90 sc = (a*x + (1 - a)*y - D) 

91 subbed = sc.sub({a: 0.1}) 

92 self.assertEqual(subbed, 0.1*x + 0.9*y - D) 

93 self.assertTrue(isinstance(subbed, Signomial)) 

94 subbed = sc.sub({a: 2.0}) 

95 self.assertTrue(isinstance(subbed, Signomial)) 

96 self.assertEqual(subbed, 2*x - y - D) 

97 _ = a.sub({a: -1}).value # fix monomial assumptions 

98 

99 

100class TestModelSubs(unittest.TestCase): 

101 """Test substitution for Model objects""" 

102 

103 def test_bad_gp_sub(self): 

104 x = Variable("x") 

105 y = Variable("y") 

106 m = Model(x, [y >= 1], {y: x}) 

107 with self.assertRaises(TypeError): 

108 m.solve() 

109 

110 def test_quantity_sub(self): 

111 if gpkit.units: 

112 x = Variable("x", 1, "cm") 

113 y = Variable("y", 1) 

114 self.assertEqual(x.sub({x: 1*gpkit.units.m}).c.magnitude, 100) 

115 # NOTE: uncomment the below if requiring Quantity substitutions 

116 # self.assertRaises(ValueError, x.sub, x, 1) 

117 self.assertRaises(DimensionalityError, x.sub, {x: 1*gpkit.ureg.N}) 

118 self.assertRaises(DimensionalityError, y.sub, {y: 1*gpkit.ureg.N}) 

119 v = gpkit.VectorVariable(3, "v", "cm") 

120 subbed = v.sub({v: [1, 2, 3]*gpkit.ureg.m}) 

121 self.assertEqual([z.c.magnitude for z in subbed], [100, 200, 300]) 

122 v = VectorVariable(1, "v", "km") 

123 v_min = VectorVariable(1, "v_min", "km") 

124 m = Model(v.prod(), [v >= v_min], 

125 {v_min: [2*gpkit.units("nmi")]}) 

126 cost = m.solve(verbosity=0)["cost"] 

127 self.assertAlmostEqual(cost/3.704, 1.0) 

128 m = Model(v.prod(), [v >= v_min], 

129 {v_min: np.array([2])*gpkit.units("nmi")}) 

130 cost = m.solve(verbosity=0)["cost"] 

131 self.assertAlmostEqual(cost/3.704, 1.0) 

132 

133 def test_phantoms(self): 

134 x = Variable("x") 

135 x_ = Variable("x", 1, lineage=[("test", 0)]) 

136 xv = VectorVariable(2, "x", [1, 1], lineage=[("vec", 0)]) 

137 m = Model(x, [x >= x_, x_ == xv.prod()]) 

138 m.solve(verbosity=0) 

139 with self.assertRaises(ValueError): 

140 _ = m.substitutions["x"] 

141 with self.assertRaises(KeyError): 

142 _ = m.substitutions["y"] 

143 with self.assertRaises(ValueError): 

144 _ = m["x"] 

145 self.assertIn(x, m.variables_byname("x")) 

146 self.assertIn(x_, m.variables_byname("x")) 

147 

148 def test_persistence(self): 

149 x = gpkit.Variable("x") 

150 y = gpkit.Variable("y") 

151 ymax = gpkit.Variable("y_{max}", 0.1) 

152 

153 with gpkit.SignomialsEnabled(): 

154 m = gpkit.Model(x, [x >= 1-y, y <= ymax]) 

155 m.substitutions[ymax] = 0.2 

156 self.assertAlmostEqual(m.localsolve(verbosity=0)["cost"], 0.8, 3) 

157 m = gpkit.Model(x, [x >= 1-y, y <= ymax]) 

158 with self.assertRaises(UnboundedGP): # from unbounded ymax 

159 m.localsolve(verbosity=0) 

160 m = gpkit.Model(x, [x >= 1-y, y <= ymax]) 

161 m.substitutions[ymax] = 0.1 

162 self.assertAlmostEqual(m.localsolve(verbosity=0)["cost"], 0.9, 3) 

163 

164 def test_united_sub_sweep(self): 

165 A = Variable("A", "USD") 

166 h = Variable("h", "USD/count") 

167 Q = Variable("Q", "count") 

168 Y = Variable("Y", "USD") 

169 m = Model(Y, [Y >= h*Q + A/Q]) 

170 m.substitutions.update({A: 500*gpkit.units("USD"), 

171 h: 35*gpkit.units("USD"), 

172 Q: ("sweep", [50, 100, 500])}) 

173 firstcost = m.solve(verbosity=0)["cost"][0] 

174 self.assertAlmostEqual(1760/firstcost, 1, 5) 

175 

176 def test_skipfailures(self): 

177 x = Variable("x") 

178 x_min = Variable("x_{min}", [1, 2]) 

179 

180 m = Model(x, [x <= 1, x >= x_min]) 

181 sol = m.solve(verbosity=0, skipsweepfailures=True) 

182 sol.table() 

183 self.assertEqual(len(sol), 1) 

184 

185 with self.assertRaises(RuntimeWarning): 

186 sol = m.solve(verbosity=0, skipsweepfailures=False) 

187 

188 m.substitutions[x_min][1][0] = 5 # so no sweeps solve 

189 with self.assertRaises(RuntimeWarning): 

190 sol = m.solve(verbosity=0, skipsweepfailures=True) 

191 

192 def test_vector_sweep(self): 

193 """Test sweep involving VectorVariables""" 

194 x = Variable("x") 

195 x_min = Variable("x_min", 1) 

196 y = VectorVariable(2, "y") 

197 m = Model(x, [x >= y.prod()]) 

198 m.substitutions.update({y: ('sweep', [[2, 3], [5, 7], [9, 11]])}) 

199 a = m.solve(verbosity=0)["cost"] 

200 b = [6, 15, 27, 14, 35, 63, 22, 55, 99] 

201 self.assertTrue(all(abs(a-b)/(a+b) < 1e-7)) 

202 x_min = Variable("x_min", 1) # constant to check array indexing 

203 m = Model(x, [x >= y.prod(), x >= x_min]) 

204 m.substitutions.update({y: ('sweep', [[2, 3], [5, 7, 11]])}) 

205 sol = m.solve(verbosity=0) 

206 a = sol["cost"] 

207 b = [10, 15, 14, 21, 22, 33] 

208 self.assertTrue(all(abs(a-b)/(a+b) < 1e-7)) 

209 self.assertEqual(sol["constants"][x_min], 1) 

210 for i, bi in enumerate(b): 

211 self.assertEqual(sol.atindex(i)["constants"][x_min], 1) 

212 ai = m.solution.atindex(i)["cost"] 

213 self.assertTrue(abs(ai-bi)/(ai+bi) < 1e-7) 

214 m = Model(x, [x >= y.prod()]) 

215 m.substitutions.update({y: ('sweep', [[2, 3, 9], [5, 7, 11]])}) 

216 self.assertRaises(ValueError, m.solve, verbosity=0) 

217 m.substitutions.update({y: [2, ("sweep", [3, 5])]}) 

218 a = m.solve(verbosity=0)["cost"] 

219 b = [6, 10] 

220 self.assertTrue(all(abs(a-b)/(a+b) < 1e-7)) 

221 # create a numpy float array, then insert a sweep element 

222 m.substitutions.update({y: [2, 3]}) 

223 m.substitutions.update({y[1]: ("sweep", [3, 5])}) 

224 a = m.solve(verbosity=0)["cost"] 

225 self.assertTrue(all(abs(a-b)/(a+b) < 1e-7)) 

226 

227 def test_calcconst(self): 

228 x = Variable("x", "hours") 

229 t_day = Variable("t_{day}", 12, "hours") 

230 t_night = Variable("t_{night}", lambda c: 24 - c[t_day], "hours") 

231 _ = pickle.dumps(t_night) 

232 m = Model(x, [x >= t_day, x >= t_night]) 

233 sol = m.solve(verbosity=0) 

234 self.assertAlmostEqual(sol(t_night)/gpkit.ureg.hours, 12) 

235 m.substitutions.update({t_day: ("sweep", [6, 8, 9, 13])}) 

236 sol = m.solve(verbosity=0) 

237 npt.assert_allclose(sol["sensitivities"]["variables"][t_day], 

238 [-1/3, -0.5, -0.6, +1], 1e-5) 

239 self.assertEqual(len(sol["cost"]), 4) 

240 npt.assert_allclose([float(l) for l in 

241 (sol(t_day) + sol(t_night))/gpkit.ureg.hours], 24) 

242 

243 def test_vector_init(self): 

244 N = 6 

245 Weight = 50000 

246 xi_dist = 6*Weight/float(N)*( 

247 (np.array(range(1, N+1)) - .5/float(N))/float(N) - 

248 (np.array(range(1, N+1)) - .5/float(N))**2/float(N)**2) 

249 

250 xi = VectorVariable(N, "xi", xi_dist, "N", "Constant Thrust per Bin") 

251 P = Variable("P", "N", "Total Power") 

252 phys_constraints = [P >= xi.sum()] 

253 objective = P 

254 eqns = phys_constraints 

255 m = Model(objective, eqns) 

256 sol = m.solve(verbosity=0) 

257 a, b = sol("xi"), xi_dist*gpkit.ureg.N 

258 self.assertTrue(all(abs(a-b)/(a+b) < 1e-7)) 

259 

260 # pylint: disable=too-many-locals 

261 def test_model_composition_units(self): 

262 class Above(Model): 

263 """A simple upper bound on x 

264 

265 Lower Unbounded 

266 --------------- 

267 x 

268 """ 

269 def setup(self): 

270 x = self.x = Variable("x", "ft") 

271 x_max = Variable("x_{max}", 1, "yard") 

272 self.cost = 1/x 

273 return [x <= x_max] 

274 

275 class Below(Model): 

276 """A simple lower bound on x 

277 

278 Upper Unbounded 

279 --------------- 

280 x 

281 """ 

282 def setup(self): 

283 x = self.x = Variable("x", "m") 

284 x_min = Variable("x_{min}", 1, "cm") 

285 self.cost = x 

286 return [x >= x_min] 

287 

288 a, b = Above(), Below() 

289 concatm = Model(a.cost*b.cost, [a, b]) 

290 concat_cost = concatm.solve(verbosity=0)["cost"] 

291 almostequal = self.assertAlmostEqual 

292 yard, cm = gpkit.ureg("yard"), gpkit.ureg("cm") 

293 ft, meter = gpkit.ureg("ft"), gpkit.ureg("m") 

294 if not isinstance(a["x"].key.units, str): 

295 almostequal(a.solve(verbosity=0)["cost"], ft/yard, 5) 

296 almostequal(b.solve(verbosity=0)["cost"], cm/meter, 5) 

297 almostequal(cm/yard, concat_cost, 5) 

298 NamedVariables.reset_modelnumbers() 

299 a1, b1 = Above(), Below() 

300 self.assertEqual(a1["x"].key.lineage, (("Above", 0),)) 

301 m = Model(a1["x"], [a1, b1, b1["x"] == a1["x"]]) 

302 sol = m.solve(verbosity=0) 

303 if not isinstance(a1["x"].key.units, str): 

304 almostequal(sol["cost"], cm/ft, 5) 

305 a1, b1 = Above(), Below() 

306 self.assertEqual(a1["x"].key.lineage, (("Above", 1),)) 

307 m = Model(b1["x"], [a1, b1, b1["x"] == a1["x"]]) 

308 sol = m.solve(verbosity=0) 

309 if not isinstance(b1["x"].key.units, str): 

310 almostequal(sol["cost"], cm/meter, 5) 

311 self.assertIn(a1["x"], sol["variables"]) 

312 self.assertIn(b1["x"], sol["variables"]) 

313 self.assertNotIn(a["x"], sol["variables"]) 

314 self.assertNotIn(b["x"], sol["variables"]) 

315 

316 def test_getkey(self): 

317 class Top(Model): 

318 """Some high level model 

319 

320 Upper Unbounded 

321 --------------- 

322 y 

323 """ 

324 def setup(self): 

325 y = self.y = Variable('y') 

326 s = Sub() 

327 sy = s["y"] 

328 self.cost = y 

329 return [s, y >= sy, sy >= 1] 

330 

331 class Sub(Model): 

332 """A simple sub model 

333 

334 Upper Unbounded 

335 --------------- 

336 y 

337 """ 

338 def setup(self): 

339 y = self.y = Variable('y') 

340 self.cost = y 

341 return [y >= 2] 

342 

343 sol = Top().solve(verbosity=0) 

344 self.assertAlmostEqual(sol['cost'], 2) 

345 

346 def test_model_recursion(self): 

347 class Top(Model): 

348 """Some high level model 

349 

350 Upper Unbounded 

351 --------------- 

352 x 

353 

354 """ 

355 def setup(self): 

356 sub = Sub() 

357 x = self.x = Variable("x") 

358 self.cost = x 

359 return sub, [x >= sub["y"], sub["y"] >= 1] 

360 

361 class Sub(Model): 

362 """A simple sub model 

363 

364 Upper Unbounded 

365 --------------- 

366 y 

367 

368 """ 

369 def setup(self): 

370 y = self.y = Variable('y') 

371 self.cost = y 

372 return [y >= 2] 

373 

374 sol = Top().solve(verbosity=0) 

375 self.assertAlmostEqual(sol['cost'], 2) 

376 

377 def test_vector_sub(self): 

378 x = VectorVariable(3, "x") 

379 y = VectorVariable(3, "y") 

380 ymax = VectorVariable(3, "ymax") 

381 

382 with SignomialsEnabled(): 

383 # issue1077 links to a case that failed for SPs only 

384 m = Model(x.prod(), [x + y >= 1, y <= ymax]) 

385 

386 m.substitutions["ymax"] = [0.3, 0.5, 0.8] 

387 m.localsolve(verbosity=0) 

388 

389 def test_spsubs(self): 

390 x = Variable("x", 5) 

391 y = Variable("y", lambda c: 2*c[x]) 

392 z = Variable("z") 

393 w = Variable("w") 

394 

395 with SignomialsEnabled(): 

396 cnstr = [z + w >= y*x, w <= y] 

397 

398 m = Model(z, cnstr) 

399 m.localsolve(verbosity=0) 

400 self.assertTrue(m.substitutions["y"], "__call__") 

401 

402class TestNomialMapSubs(unittest.TestCase): 

403 "Tests substitutions of nomialmaps" 

404 def test_monomial_sub(self): 

405 z = Variable("z") 

406 w = Variable("w") 

407 

408 with self.assertRaises(ValueError): 

409 z.hmap.sub({z.key: w.key}, varkeys=z.vks) 

410 

411 def test_subinplace_zero(self): 

412 z = Variable("z") 

413 w = Variable("w") 

414 

415 p = 2*w + z*w + 2 

416 

417 self.assertEqual(p.sub({z: -2}), 2) 

418 

419TESTS = [TestNomialSubs, TestModelSubs, TestNomialMapSubs] 

420 

421if __name__ == "__main__": # pragma: no cover 

422 run_tests(TESTS)