1#TODO: add depth-culling to simplify 

2# clean up factor generation 

3 

4# pylint: skip-file 

5import string 

6from collections import defaultdict, namedtuple, Counter 

7from gpkit.nomials import Monomial, Posynomial, Variable 

8from gpkit.nomials.map import NomialMap 

9from gpkit.small_scripts import mag 

10from gpkit.small_classes import FixedScalar, HashVector 

11from gpkit.exceptions import DimensionalityError 

12from gpkit.repr_conventions import unitstr as get_unitstr 

13from gpkit.varkey import VarKey 

14import numpy as np 

15 

16 

17Tree = namedtuple("Tree", ["key", "value", "branches"]) 

18Transform = namedtuple("Transform", ["factor", "power", "origkey"]) 

19def is_factor(key): 

20 return (isinstance(key, Transform) and key.power == 1) 

21def is_power(key): 

22 return (isinstance(key, Transform) and key.power != 1) 

23 

24def get_free_vks(posy, solution): 

25 "Returns all free vks of a given posynomial for a given solution" 

26 return set(vk for vk in posy.vks if vk not in solution["constants"]) 

27 

28def get_model_breakdown(solution): 

29 breakdowns = {"|sensitivity|": 0} 

30 for modelname, senss in solution["sensitivities"]["models"].items(): 

31 senss = abs(senss) # for those monomial equalities 

32 *namespace, name = modelname.split(".") 

33 subbd = breakdowns 

34 subbd["|sensitivity|"] += senss 

35 for parent in namespace: 

36 if parent not in subbd: 

37 subbd[parent] = {parent: {}} 

38 subbd = subbd[parent] 

39 if "|sensitivity|" not in subbd: 

40 subbd["|sensitivity|"] = 0 

41 subbd["|sensitivity|"] += senss 

42 subbd[name] = {"|sensitivity|": senss} 

43 # print(breakdowns["HyperloopSystem"]["|sensitivity|"]) 

44 breakdowns = {"|sensitivity|": 0} 

45 for constraint, senss in solution["sensitivities"]["constraints"].items(): 

46 senss = abs(senss) # for those monomial 

47 if senss <= 1e-5: 

48 continue 

49 subbd = breakdowns 

50 subbd["|sensitivity|"] += senss 

51 for parent in constraint.lineagestr().split("."): 

52 if parent == "": 

53 continue 

54 if parent not in subbd: 

55 subbd[parent] = {} 

56 subbd = subbd[parent] 

57 if "|sensitivity|" not in subbd: 

58 subbd["|sensitivity|"] = 0 

59 subbd["|sensitivity|"] += senss 

60 # treat vectors as namespace 

61 constraint = constraint.str_without({"unnecessary lineage", "units", ":MAGIC:"+constraint.lineagestr()}) 

62 subbd[constraint] = {"|sensitivity|": senss} 

63 for vk in solution.vks: 

64 if vk not in solution["sensitivities"]["variables"]: 

65 continue 

66 senss = abs(solution["sensitivities"]["variables"][vk]) 

67 if hasattr(senss, "shape"): 

68 senss = np.nansum(senss) 

69 if senss <= 1e-5: 

70 continue 

71 subbd = breakdowns 

72 subbd["|sensitivity|"] += senss 

73 for parent in vk.lineagestr().split("."): 

74 if parent == "": 

75 continue 

76 if parent not in subbd: 

77 subbd[parent] = {} 

78 subbd = subbd[parent] 

79 if "|sensitivity|" not in subbd: 

80 subbd["|sensitivity|"] = 0 

81 subbd["|sensitivity|"] += senss 

82 # treat vectors as namespace (indexing vectors above) 

83 vk = vk.str_without({"lineage"}) + get_valstr(vk, solution, " = %s").replace(", fixed", "") 

84 subbd[vk] = {"|sensitivity|": senss} 

85 # TODO: track down in a live-solve environment why this isn't the same 

86 # print(breakdowns["HyperloopSystem"]["|sensitivity|"]) 

87 return breakdowns 

88 

89def crawl_modelbd(bd, name="Model"): 

90 tree = Tree(name, bd.pop("|sensitivity|"), []) 

91 for subname, subtree in sorted(bd.items(), 

92 key=lambda kv: -kv[1]["|sensitivity|"]): 

93 tree.branches.append(crawl_modelbd(subtree, subname)) 

94 return tree 

95 

96def divide_out_vk(vk, pow, lt, gt): 

97 hmap = NomialMap({HashVector({vk: 1}): 1.0}) 

98 hmap.units = vk.units 

99 var = Monomial(hmap)**pow 

100 lt, gt = lt/var, gt/var 

101 lt.ast = gt.ast = None 

102 return lt, gt 

103 

104# @profile 

105def get_breakdowns(solution): 

106 """Returns {key: (lt, gt, constraint)} for breakdown constrain in solution. 

107 

108 A breakdown constraint is any whose "gt" contains a single free variable. 

109 

110 (At present, monomial constraints check both sides as "gt") 

111 """ 

112 breakdowns = defaultdict(list) 

113 beatout = defaultdict(set) 

114 for constraint, senss in sorted(solution["sensitivities"]["constraints"].items(), key=lambda kv: (-abs(kv[1]), str(kv[0]))): 

115 if abs(senss) <= 1e-5: # only tight-ish ones 

116 continue 

117 if constraint.oper == ">=": 

118 gt, lt = (constraint.left, constraint.right) 

119 elif constraint.oper == "<=": 

120 lt, gt = (constraint.left, constraint.right) 

121 elif constraint.oper == "=": 

122 if senss > 0: # l_over_r is more sensitive - see nomials/math.py 

123 lt, gt = (constraint.left, constraint.right) 

124 else: # r_over_l is more sensitive - see nomials/math.py 

125 gt, lt = (constraint.left, constraint.right) 

126 if lt.any_nonpositive_cs or len(gt.hmap) > 1: 

127 continue # no signomials # TODO: approximate signomials at sol 

128 pos_gtvks = {vk for vk, pow in gt.exp.items() if pow > 0} 

129 if len(pos_gtvks) > 1: 

130 pos_gtvks &= get_free_vks(gt, solution) # remove constants 

131 for vk in gt.vks: 

132 if vk.name == "podturnaroundtime": 

133 print(constraint, senss, len(pos_gtvks)) 

134 if len(pos_gtvks) == 1: 

135 chosenvk, = pos_gtvks 

136 breakdowns[chosenvk].append((lt, gt, constraint)) 

137 for constraint, senss in sorted(solution["sensitivities"]["constraints"].items(), key=lambda kv: (-abs(kv[1]), str(kv[0]))): 

138 if abs(senss) <= 1e-5: # only tight-ish ones 

139 continue 

140 if constraint.oper == ">=": 

141 gt, lt = (constraint.left, constraint.right) 

142 elif constraint.oper == "<=": 

143 lt, gt = (constraint.left, constraint.right) 

144 elif constraint.oper == "=": 

145 if senss > 0: # l_over_r is more sensitive - see nomials/math.py 

146 lt, gt = (constraint.left, constraint.right) 

147 else: # r_over_l is more sensitive - see nomials/math.py 

148 gt, lt = (constraint.left, constraint.right) 

149 if lt.any_nonpositive_cs or len(gt.hmap) > 1: 

150 continue # no signomials # TODO: approximate signomials at sol 

151 pos_gtvks = {vk for vk, pow in gt.exp.items() if pow > 0} 

152 if len(pos_gtvks) > 1: 

153 pos_gtvks &= get_free_vks(gt, solution) # remove constants 

154 if len(pos_gtvks) != 1: # we'll choose our favorite vk 

155 for vk, pow in gt.exp.items(): 

156 if pow < 0: # remove all non-positive 

157 lt, gt = divide_out_vk(vk, pow, lt, gt) 

158 # bring over common factors from lt 

159 lt_pows = defaultdict(set) 

160 for exp in lt.hmap: 

161 for vk, pow in exp.items(): 

162 lt_pows[vk].add(pow) 

163 for vk, pows in lt_pows.items(): 

164 if len(pows) == 1: 

165 pow, = pows 

166 if pow < 0: # ...but only if they're positive 

167 lt, gt = divide_out_vk(vk, pow, lt, gt) 

168 # don't choose something that's already been broken down 

169 candidatevks = {vk for vk in gt.vks if vk not in breakdowns} 

170 if candidatevks: 

171 chosenvk, *_ = sorted(candidatevks, key=lambda vk: (-gt.exp[vk]*solution["sensitivities"]["variablerisk"].get(vk, 0), str(vk))) 

172 for vk in gt.vks: 

173 if vk is not chosenvk: 

174 lt, gt = divide_out_vk(vk, pow, lt, gt) 

175 breakdowns[chosenvk].append((lt, gt, constraint)) 

176 breakdowns = dict(breakdowns) # remove the defaultdict-ness 

177 

178 prevlen = None 

179 while len(BASICALLY_FIXED_VARIABLES) != prevlen: 

180 prevlen = len(BASICALLY_FIXED_VARIABLES) 

181 for key in breakdowns: 

182 if key not in BASICALLY_FIXED_VARIABLES: 

183 get_fixity(key, breakdowns, solution, BASICALLY_FIXED_VARIABLES) 

184 return breakdowns 

185 

186BASICALLY_FIXED_VARIABLES = set() 

187 

188 

189def get_fixity(key, bd, solution, basically_fixed=set(), visited=set()): 

190 lt, gt, _ = bd[key][0] 

191 free_vks = get_free_vks(lt, solution).union(get_free_vks(gt, solution)) 

192 for vk in free_vks: 

193 if vk is key or vk in BASICALLY_FIXED_VARIABLES: 

194 continue # currently checking or already checked 

195 if vk not in bd: 

196 return # a very free variable, can't even be broken down 

197 if vk in visited: 

198 return # tried it before, implicitly it didn't work out 

199 # maybe it's basically fixed? 

200 visited.add(key) 

201 get_fixity(vk, bd, solution, basically_fixed, visited) 

202 if vk not in BASICALLY_FIXED_VARIABLES: 

203 return # ...well, we tried 

204 basically_fixed.add(key) 

205 

206SOLCACHE = {} 

207def solcache(solution, key): # replaces solution(key) 

208 if key not in SOLCACHE: 

209 SOLCACHE[key] = solution(key) 

210 return SOLCACHE[key] 

211 

212# @profile # ~84% of total last check # TODO: remove 

213def crawl(key, bd, solution, basescale=1, permissivity=2, verbosity=0, 

214 visited_bdkeys=None, gone_negative=False): 

215 "Returns the tree of breakdowns of key in bd, sorting by solution's values" 

216 if key in bd: 

217 # TODO: do multiple if sensitivities are quite close? 

218 composition, keymon, constraint = bd[key][0] 

219 elif isinstance(key, Posynomial): 

220 composition = key 

221 keymon = None 

222 else: 

223 raise TypeError("the `key` argument must be a VarKey or Posynomial.") 

224 

225 if visited_bdkeys is None: 

226 visited_bdkeys = set() 

227 if verbosity == 1: 

228 solution.set_necessarylineage() 

229 if verbosity: 

230 indent = verbosity-1 # HACK: a bit of overloading, here 

231 keyvalstr = "%s (%s)" % (key.str_without(["unnecessary lineage", "units"]), 

232 get_valstr(key, solution)) 

233 print(" "*indent + keyvalstr + ", which breaks down further") 

234 indent += 1 

235 orig_subtree = subtree = [] 

236 tree = Tree(key, basescale, subtree) 

237 visited_bdkeys.add(key) 

238 if keymon is None: 

239 scale = solution(key)/basescale 

240 else: 

241 interesting_vks = {key} 

242 subkey, = interesting_vks 

243 power = keymon.exp[subkey] 

244 boring_vks = set(keymon.vks) - interesting_vks 

245 scale = solution(key)**power/basescale 

246 # TODO: make method that can handle both kinds of transforms 

247 if power != 1 or boring_vks or mag(keymon.c) != 1 or keymon.units != key.units: 

248 units = 1 

249 exp = HashVector() 

250 for vk in interesting_vks: 

251 exp[vk] = keymon.exp[vk] 

252 if vk.units: 

253 units *= vk.units**keymon.exp[vk] 

254 subhmap = NomialMap({exp: 1}) 

255 try: 

256 subhmap.units = None if units == 1 else units 

257 except DimensionalityError: 

258 # pints was unable to divide a unit by itself bc 

259 # it has terrible floating-point errors. 

260 # so let's assume it isn't dimensionless 

261 # even though it probably is 

262 subhmap.units = units 

263 freemon = Monomial(subhmap) 

264 factor = Monomial(keymon/freemon) 

265 scale = scale * solution(factor) 

266 if factor != 1: 

267 factor = factor**(-1/power) # invert the transform 

268 factor.ast = None 

269 if verbosity: 

270 keyvalstr = "%s (%s)" % (factor.str_without(["unnecessary lineage", "units"]), 

271 get_valstr(factor, solution)) 

272 print(" "*indent + "(with a factor of " + keyvalstr + " )") 

273 subsubtree = [] 

274 transform = Transform(factor, 1, keymon) 

275 orig_subtree.append(Tree(transform, basescale, subsubtree)) 

276 orig_subtree = subsubtree 

277 if power != 1: 

278 if verbosity: 

279 print(" "*indent + "(with a power of %.2g )" % power) 

280 subsubtree = [] 

281 transform = Transform(1, 1/power, keymon) # inverted bc it's on the gt side 

282 orig_subtree.append(Tree(transform, basescale, subsubtree)) 

283 orig_subtree = subsubtree 

284 if verbosity: 

285 if keymon is not None: 

286 print(" "*indent + "in: " + constraint.str_without(["units", "lineage"])) 

287 print(" "*indent + "by:") 

288 indent += 1 

289 

290 # TODO: use ast_parsing instead of chop? 

291 monsols = [solcache(solution, mon) for mon in composition.chop()] # ~20% of total last check # TODO: remove 

292 parsed_monsols = [getattr(mon, "value", mon) for mon in monsols] 

293 monvals = [float(mon/scale) for mon in parsed_monsols] # ~10% of total last check # TODO: remove 

294 # sort by value, preserving order in case of value tie 

295 sortedmonvals = sorted(zip(monvals, range(len(monvals)), 

296 composition.chop()), reverse=True) 

297 for scaledmonval, _, mon in sortedmonvals: 

298 if not scaledmonval: 

299 continue 

300 subtree = orig_subtree # return to the original subtree 

301 

302 # time for some filtering 

303 interesting_vks = mon.vks 

304 potential_filters = [ 

305 {vk for vk in interesting_vks if vk not in bd}, 

306 mon.vks - get_free_vks(mon, solution), 

307 {vk for vk in interesting_vks if vk in BASICALLY_FIXED_VARIABLES} 

308 ] 

309 if scaledmonval < 1 - permissivity: # skip breakdown filter 

310 potential_filters = potential_filters[1:] 

311 for filter in potential_filters: 

312 if interesting_vks - filter: # don't remove the last one 

313 interesting_vks = interesting_vks - filter 

314 # if filters weren't enough and permissivity is high enough, sort! 

315 if len(interesting_vks) > 1 and permissivity > 1: 

316 best_vks = sorted((vk for vk in interesting_vks if vk in bd), 

317 key=lambda vk: 

318 # TODO: without exp: "most strongly broken-down component" 

319 # but it could use nus (or v_ss) to say 

320 # "breakdown which the solution is most sensitive to" 

321 # ...right now it's in-between 

322 (-mon.exp[vk]*abs(solution["sensitivities"]["constraints"][bd[vk][0][2]]), 

323 str(bd[vk][0][0]))) # ~5% of total last check # TODO: remove 

324 # changing to str(vk) above does some odd stuff 

325 if best_vks: 

326 interesting_vks = set([best_vks[0]]) 

327 boring_vks = mon.vks - interesting_vks 

328 

329 subkey = None 

330 if len(interesting_vks) == 1: 

331 subkey, = interesting_vks 

332 if subkey in visited_bdkeys and len(sortedmonvals) == 1: 

333 continue # don't even go there 

334 if subkey not in bd: 

335 power = 1 # no need for a transform 

336 else: 

337 power = mon.exp[subkey] 

338 if power < 0 and gone_negative: 

339 subkey = None # don't breakdown another negative 

340 

341 if subkey is None: 

342 power = 1 

343 if scaledmonval > 1 - permissivity and not boring_vks: 

344 boring_vks = interesting_vks 

345 interesting_vks = set() 

346 if not interesting_vks: 

347 # prioritize showing some boring_vks as if they were "free" 

348 if len(boring_vks) == 1: 

349 interesting_vks = boring_vks 

350 boring_vks = set() 

351 else: 

352 for vk in list(boring_vks): 

353 if vk.units and not vk.units.dimensionless: 

354 interesting_vks.add(vk) 

355 boring_vks.remove(vk) 

356 

357 if interesting_vks and (boring_vks or mag(mon.c) != 1): 

358 units = 1 

359 exp = HashVector() 

360 for vk in interesting_vks: 

361 exp[vk] = mon.exp[vk] 

362 if vk.units: 

363 units *= vk.units**mon.exp[vk] 

364 subhmap = NomialMap({exp: 1}) 

365 subhmap.units = None if units is 1 else units 

366 freemon = Monomial(subhmap) 

367 factor = mon/freemon # autoconvert... 

368 if (factor.units is None and isinstance(factor, FixedScalar) 

369 and abs(factor.value - 1) <= 1e-4): 

370 factor = 1 # minor fudge to clear numerical inaccuracies 

371 if factor != 1 : 

372 factor.ast = None 

373 if verbosity: 

374 keyvalstr = "%s (%s)" % (factor.str_without(["unnecessary lineage", "units"]), 

375 get_valstr(factor, solution)) 

376 print(" "*indent + "(with a factor of %s )" % keyvalstr) 

377 subsubtree = [] 

378 transform = Transform(factor, 1, mon) 

379 subtree.append(Tree(transform, scaledmonval, subsubtree)) 

380 subtree = subsubtree 

381 mon = freemon # simplifies units 

382 if power != 1: 

383 if verbosity: 

384 print(" "*indent + "(with a power of %.2g )" % power) 

385 subsubtree = [] 

386 transform = Transform(1, power, mon) 

387 subtree.append(Tree(transform, scaledmonval, subsubtree)) 

388 subtree = subsubtree 

389 mon = mon**(1/power) 

390 mon.ast = None 

391 # TODO: make minscale an argument - currently an arbitrary 0.01 

392 if (subkey is not None and subkey not in visited_bdkeys 

393 and subkey in bd and scaledmonval > 0.05): 

394 if verbosity: 

395 verbosity = indent + 1 # slight hack 

396 subsubtree = crawl(subkey, bd, solution, scaledmonval, 

397 permissivity, verbosity, set(visited_bdkeys), 

398 gone_negative) 

399 subtree.append(subsubtree) 

400 else: 

401 if verbosity: 

402 keyvalstr = "%s (%s)" % (mon.str_without(["unnecessary lineage", "units"]), 

403 get_valstr(mon, solution)) 

404 print(" "*indent + keyvalstr) 

405 subtree.append(Tree(mon, scaledmonval, [])) 

406 if verbosity == 1: 

407 solution.set_necessarylineage(clear=True) 

408 return tree 

409 

410SYMBOLS = string.ascii_uppercase + string.ascii_lowercase 

411for ambiguous_symbol in "lILT": 

412 SYMBOLS = SYMBOLS.replace(ambiguous_symbol, "") 

413 

414def get_spanstr(legend, length, label, leftwards, solution): 

415 "Returns span visualization, collapsing labels to symbols" 

416 if label is None: 

417 return " "*length 

418 spacer, lend, rend = "│", "┯", "┷" 

419 if isinstance(label, Transform): 

420 spacer, lend, rend = "╎", "╤", "╧" 

421 if label.power != 1: 

422 spacer = " " 

423 lend = rend = "^" if label.power > 0 else "/" 

424 # remove origkeys so they collide in the legends dictionary 

425 label = Transform(label.factor, label.power, None) 

426 if label.power == 1 and len(str(label.factor)) == 1: 

427 legend[label] = str(label.factor) 

428 

429 if label not in legend: 

430 legend[label] = SYMBOLS[len(legend)] 

431 shortname = legend[label] 

432 

433 if length <= 1: 

434 return shortname 

435 shortside = int(max(0, length - 2)/2) 

436 longside = int(max(0, length - 3)/2) 

437 if leftwards: 

438 if length == 2: 

439 return lend + shortname 

440 return lend + spacer*shortside + shortname + spacer*longside + rend 

441 else: 

442 if length == 2: 

443 return shortname + rend 

444 # HACK: no corners on long rightwards - only used for depth 0 

445 return "┃"*(longside+1) + shortname + "┃"*(shortside+1) 

446 

447def simplify(tree, extent, solution, collapse, depth=0, justsplit=False): 

448 # TODO: add vertical simplification? 

449 key, val, branches = tree 

450 if collapse: # collapse Transforms with power 1 

451 while any(isinstance(branch.key, Transform) and branch.key.power > 0 for branch in branches): 

452 newbranches = [] 

453 for branch in branches: 

454 # isinstance(branch.key, Transform) and branch.key.power > 0 

455 if isinstance(branch.key, Transform) and branch.key.power > 0: 

456 newbranches.extend(branch.branches) 

457 else: 

458 newbranches.append(branch) 

459 branches = newbranches 

460 

461 scale = extent/val 

462 values = [b.value for b in branches] 

463 bkey_indexs = {} 

464 for i, b in enumerate(branches): 

465 k = get_keystr(b.key, solution) 

466 if isinstance(b.key, Transform): 

467 if len(b.branches) == 1: 

468 k = get_keystr(b.branches[0].key, solution) 

469 if k in bkey_indexs: 

470 values[bkey_indexs[k]] += values[i] 

471 values[i] = None 

472 else: 

473 bkey_indexs[k] = i 

474 if any(v is None for v in values): 

475 branches, values = zip(*((b, v) for b, v in zip(branches, values) if v is not None)) 

476 branches = list(branches) 

477 values = list(values) 

478 extents = [int(round(scale*v)) for v in values] 

479 surplus = extent - sum(extents) 

480 for i, b in enumerate(branches): 

481 if isinstance(b.key, Transform): 

482 subscale = extents[i]/b.value 

483 if not any(round(subscale*subv) for _, subv, _ in b.branches): 

484 extents[i] = 0 # transform with no worthy heirs: misc it 

485 if not any(extents): 

486 return Tree(key, extent, []) 

487 if not all(extents): # create a catch-all 

488 branches = branches.copy() 

489 miscvkeys, miscval = [], 0 

490 for subextent in reversed(extents): 

491 if not subextent or (branches[-1].value < miscval and surplus < 0): 

492 extents.pop() 

493 k, v, _ = branches.pop() 

494 if isinstance(k, Transform): 

495 k = k.origkey # TODO: this is the only use of origkey - remove it 

496 if isinstance(k, tuple): 

497 vkeys = [(-kv[1], str(kv[0]), kv[0]) for kv in k] 

498 if not isinstance(k, tuple): 

499 vkeys = [(-v, str(k), k)] 

500 miscvkeys += vkeys 

501 surplus -= (round(scale*(miscval + v)) 

502 - round(scale*miscval) - subextent) 

503 miscval += v 

504 misckeys = tuple(k for _, _, k in sorted(miscvkeys)) 

505 branches.append(Tree(misckeys, miscval, [])) 

506 extents.append(int(round(scale*miscval))) 

507 if surplus: 

508 sign = int(np.sign(surplus)) 

509 bump_priority = sorted((ext, sign*b.value, i) for i, (b, ext) 

510 in enumerate(zip(branches, extents))) 

511 while surplus: 

512 extents[bump_priority.pop()[-1]] += sign 

513 surplus -= sign 

514 

515 # simplify based on how we're branching 

516 branchfactor = len([ext for ext in extents if ext]) - 1 

517 if depth and not isinstance(key, Transform): 

518 if extent == 1 or branchfactor >= max(extent-2, 1): 

519 # if we'd branch too much, stop 

520 return Tree(key, extent, []) 

521 if collapse and not branchfactor and not justsplit: 

522 # if we didn't just split and aren't about to, collapse 

523 return simplify(branches[0], extent, solution, collapse, 

524 depth=depth+1, justsplit=False) 

525 if branchfactor: 

526 justsplit = True 

527 elif not isinstance(key, Transform): # justsplit passes through transforms 

528 justsplit = False 

529 

530 tree = Tree(key, extent, []) 

531 for branch, subextent in zip(branches, extents): 

532 if subextent: 

533 branch = simplify(branch, subextent, solution, collapse, 

534 depth=depth+1, justsplit=justsplit) 

535 if (collapse and is_power(branch.key) 

536 and all(is_power(b.key) for b in branch.branches)): 

537 # collapse stacked powers 

538 power = branch.key.power 

539 for b in branch.branches: 

540 key = Transform(1, power*b.key.power, None) 

541 if key.power == 1: # powers canceled, collapse both 

542 tree.branches.extend(b.branches) 

543 else: # collapse this level 

544 tree.branches.append(Tree(key, b.value, b.branches)) 

545 else: 

546 tree.branches.append(branch) 

547 return tree 

548 

549def layer(map, tree, maxdepth, depth=0): 

550 "Turns the tree into a 2D-array" 

551 if depth <= maxdepth: 

552 if len(map) <= depth: 

553 map.append([]) 

554 key, extent, branches = tree 

555 map[depth].append((key, extent)) 

556 if not branches: 

557 branches = [Tree(None, extent, [])] # pad it out 

558 for branch in branches: 

559 layer(map, branch, maxdepth, depth+1) 

560 return map 

561 

562def plumb(tree, depth=0): 

563 "Finds maximum depth of a tree" 

564 maxdepth = depth 

565 for branch in tree.branches: 

566 maxdepth = max(maxdepth, plumb(branch, depth+1)) 

567 return maxdepth 

568 

569# @profile # ~16% of total last check # TODO: remove 

570def graph(tree, solution, extent=None, maxdepth=None, showlegend=False, maxwidth=110): 

571 "Prints breakdown" 

572 if maxdepth is None: 

573 maxdepth = plumb(tree) 

574 if extent is None: # auto-zoom-in from 20 

575 prev_extent = None 

576 extent = 20 

577 while prev_extent != extent: 

578 subtree = simplify(tree, extent, solution, collapse=not showlegend) 

579 mt = layer([], subtree, maxdepth) 

580 prev_extent = extent 

581 extent = min(extent, max(*(4*len(at_depth) for at_depth in mt))) 

582 else: 

583 subtree = simplify(tree, extent, solution, collapse=not showlegend) 

584 mt = layer([], subtree, maxdepth) 

585 legend = {} 

586 chararray = np.full((len(mt), extent), "", "object") 

587 for depth, elements_at_depth in enumerate(mt): 

588 row = "" 

589 for i, (element, length) in enumerate(elements_at_depth): 

590 leftwards = depth > 0 and length > 2 

591 row += get_spanstr(legend, length, element, leftwards, solution) 

592 chararray[depth, :] = list(row) 

593 

594 solution.set_necessarylineage() 

595 # Format depth=0 

596 A_key, = [key for key, value in legend.items() if value == "A"] 

597 A_str = get_keystr(A_key, solution) 

598 prefix = "" 

599 if isinstance(A_key, VarKey) and A_key.necessarylineage: 

600 prefix = A_key.lineagestr() 

601 A_valstr = get_valstr(A_key, solution, into="(%s)") 

602 fmt = "{0:>%s}" % (max(len(A_str), len(A_valstr)) + 3) 

603 for j, entry in enumerate(chararray[0,:]): 

604 if entry == "A": 

605 chararray[0,j] = fmt.format(A_str + "╺┫") 

606 chararray[0,j+1] = fmt.format(A_valstr + " ┃") 

607 else: 

608 chararray[0,j] = fmt.format(entry) 

609 # Format depths 1+ 

610 labeled = set() 

611 new_legend = {} 

612 for pos in range(extent): 

613 for depth in reversed(range(1,len(mt))): 

614 value = chararray[depth, pos] 

615 if value not in SYMBOLS: 

616 continue 

617 key, = [k for k, val in legend.items() if val == value] 

618 if getattr(key, "vks", None) and len(key.vks) == 1 and all(vk in new_legend for vk in key.vks): 

619 key, = key.vks 

620 if key not in new_legend and (isinstance(key, tuple) or (depth != len(mt) - 1 and chararray[depth+1, pos] != " ")): 

621 new_legend[key] = SYMBOLS[len(new_legend)] 

622 if key in new_legend: 

623 chararray[depth, pos] = new_legend[key] 

624 if isinstance(key, tuple) and not isinstance(key, Transform): 

625 chararray[depth, pos] = "*" + chararray[depth, pos] 

626 if showlegend: 

627 continue 

628 tryup, trydn = True, True 

629 span = 0 

630 if not showlegend and is_power(key): 

631 chararray[depth, pos] = "^" if key.power > 0 else "/" 

632 continue 

633 

634 keystr = get_keystr(key, solution, prefix) 

635 valuestr = get_valstr(key, solution, into=" (%s)") 

636 if keystr in labeled: 

637 valuestr = "" 

638 if not showlegend: 

639 fmt = "{0:<%s}" % max(len(keystr) + 3, len(valuestr) + 2) 

640 else: 

641 fmt = "{0:<1}" 

642 while tryup or trydn: 

643 span += 1 

644 if tryup: 

645 if pos - span < 0: 

646 tryup = False 

647 else: 

648 upchar = chararray[depth, pos-span] 

649 if upchar == "│": 

650 chararray[depth, pos-span] = fmt.format("┃") 

651 elif upchar == "┯": 

652 chararray[depth, pos-span] = fmt.format("┓") 

653 else: 

654 tryup = False 

655 if trydn: 

656 if pos + span >= extent: 

657 trydn = False 

658 else: 

659 dnchar = chararray[depth, pos+span] 

660 if dnchar == "│": 

661 chararray[depth, pos+span] = fmt.format("┃") 

662 elif dnchar == "┷": 

663 chararray[depth, pos+span] = fmt.format("┛") 

664 else: 

665 trydn = False 

666 #TODO: make submodels show up with this; bd should be an argument 

667 if showlegend and (key in bd or (hasattr(key, "vks") and key.vks and any(vk in bd for vk in key.vks))): 

668 linkstr = "┣┉" 

669 else: 

670 linkstr = "┣╸" 

671 if not isinstance(key, FixedScalar): 

672 labeled.add(keystr) 

673 if span > 1 and (not showlegend or pos + 2 >= extent or chararray[depth, pos+1] == "┃"): 

674 chararray[depth, pos+1] = fmt.format(chararray[depth, pos+1].rstrip() + valuestr) 

675 elif showlegend: 

676 keystr += valuestr 

677 chararray[depth, pos] = fmt.format(linkstr + keystr) 

678 # Rotate and print 

679 toowiderows = [] 

680 rows = chararray.T.tolist() 

681 if not showlegend: # remove according to maxwidth 

682 for i, orig_row in enumerate(rows): 

683 depth_occluded = -1 

684 width = None 

685 row = orig_row.copy() 

686 row.append("") 

687 while width is None or width > maxwidth: 

688 row = row[:-1] 

689 rowstr = " " + "".join(row).rstrip() 

690 width = len(rowstr) 

691 depth_occluded += 1 

692 if depth_occluded: 

693 previous_is_pow = orig_row[-depth_occluded-1] in "^/" 

694 if abs(depth_occluded) + 1 + previous_is_pow < len(orig_row): 

695 strdepth = len(" " + "".join(orig_row[:-depth_occluded])) 

696 toowiderows.append((strdepth, i)) 

697 rowstrs = [" " + "".join(row).rstrip() for row in rows] 

698 for depth_occluded, i in sorted(toowiderows, reverse=True): 

699 if len(rowstrs[i]) <= depth_occluded: 

700 continue # already occluded 

701 if "┣" == rowstrs[i][depth_occluded]: 

702 pow = 0 

703 while rowstrs[i][depth_occluded-pow-1] in "^/": 

704 pow += 1 

705 rowstrs[i] = rowstrs[i][:depth_occluded-pow] 

706 connected = "^┃┓┛┣╸" 

707 for dir in [-1, 1]: 

708 idx = i + dir 

709 while (0 <= idx < len(rowstrs) 

710 and len(rowstrs[idx]) > depth_occluded 

711 and rowstrs[idx][depth_occluded] 

712 and rowstrs[idx][depth_occluded] in connected): 

713 while rowstrs[idx][depth_occluded-pow-1] in "^/": 

714 pow += 1 

715 rowstrs[idx] = rowstrs[idx][:depth_occluded-pow] 

716 idx += dir 

717 vertstr = "\n".join(rowstr.rstrip() for rowstr in rowstrs) 

718 print() 

719 print(vertstr) 

720 print() 

721 legend = new_legend 

722 

723 if showlegend: # create and print legend 

724 legend_lines = [] 

725 for key, shortname in sorted(legend.items(), key=lambda kv: kv[1]): 

726 legend_lines.append(legend_entry(key, shortname, solution)) 

727 maxlens = [max(len(el) for el in col) for col in zip(*legend_lines)] 

728 fmts = ["{0:<%s}" % L for L in maxlens] 

729 for line in legend_lines: 

730 line = "".join(fmt.format(cell) 

731 for fmt, cell in zip(fmts, line) if cell).rstrip() 

732 print(" " + line) 

733 

734 solution.set_necessarylineage(clear=True) 

735 

736def legend_entry(key, shortname, solution): 

737 "Returns list of legend elements" 

738 operator = note = "" 

739 keystr = valuestr = " " 

740 operator = "= " if shortname else " + " 

741 if is_factor(key): 

742 operator = " ×" 

743 key = key.factor 

744 free, quasifixed = False, False 

745 if any(vk not in BASICALLY_FIXED_VARIABLES 

746 for vk in get_free_vks(key, solution)): 

747 note = " [free factor]" 

748 if is_power(key): 

749 valuestr = " ^%.3g" % key.power 

750 else: 

751 valuestr = get_valstr(key, solution, into=" "+operator+"%s") 

752 if not isinstance(key, FixedScalar): 

753 keystr = get_keystr(key, solution) 

754 return ["%-4s" % shortname, keystr, valuestr, note] 

755 

756def get_keystr(key, solution, prefix=""): 

757 if key is solution.costposy: 

758 out = "Cost" 

759 elif hasattr(key, "str_without"): 

760 out = key.str_without({"unnecessary lineage", "units", ":MAGIC:"+prefix}) 

761 elif isinstance(key, tuple): 

762 out = "[%i terms]" % len(key) 

763 else: 

764 out = str(key) 

765 return out if len(out) <= 67 else out[:66]+"…" 

766 

767def get_valstr(key, solution, into="%s"): 

768 "Returns formatted string of the value of key in solution." 

769 # get valuestr 

770 try: 

771 value = solution(key) 

772 except (ValueError, TypeError): 

773 try: 

774 value = sum(solution(subkey) for subkey in key) 

775 except (ValueError, TypeError): 

776 return " " 

777 if isinstance(value, FixedScalar): 

778 value = value.value 

779 if 1e3 <= mag(value) < 1e6: 

780 valuestr = "{:,.0f}".format(mag(value)) 

781 else: 

782 valuestr = "%-.3g" % mag(value) 

783 # get unitstr 

784 if hasattr(key, "unitstr"): 

785 unitstr = key.unitstr() 

786 else: 

787 try: 

788 if hasattr(value, "units"): 

789 value.ito_reduced_units() 

790 except DimensionalityError: 

791 pass 

792 unitstr = get_unitstr(value) 

793 if unitstr[:2] == "1/": 

794 unitstr = "/" + unitstr[2:] 

795 if key in solution["constants"] or (hasattr(key, "vks") and key.vks and all(vk in solution["constants"] for vk in key.vks)): 

796 unitstr += ", fixed" 

797 return into % (valuestr + unitstr) 

798 

799 

800import pickle 

801from gpkit import ureg 

802ureg.define("pax = 1") 

803ureg.define("paxkm = km") 

804ureg.define("trip = 1") 

805 

806print("STARTING...") 

807from gpkit.tests.helpers import StdoutCaptured 

808 

809import difflib 

810 

811permissivity = 2 

812 

813# sol = pickle.load(open("solar.p", "rb")) 

814# bd = get_breakdowns(sol) 

815# mbd = get_model_breakdown(sol) 

816# mtree = crawl_modelbd(mbd) 

817# graph(mtree, sol, showlegend=True, extent=20) 

818# graph(mtree.branches[0].branches[0].branches[0], sol, showlegend=False, extent=20) 

819# 

820# tree = crawl(sol.costposy, bd, sol, permissivity=2, verbosity=0) 

821# graph(tree, sol) 

822# 

823# 

824# # key, = [vk for vk in bd if "Wing.Planform.b" in str(vk)] 

825# # tree = crawl(key, bd, sol, permissivity=2, verbosity=1) 

826# # graph(tree, sol) 

827# # key, = [vk for vk in bd if "Aircraft.Fuselage.R[0,0]" in str(vk)] 

828# # tree = crawl(key, bd, sol, permissivity=2, verbosity=1) 

829# # graph(tree, sol) 

830# # key, = [vk for vk in bd if "Mission.Climb.AircraftDrag.CD[0]" in str(vk)] 

831# # tree = crawl(key, bd, sol, permissivity=2, verbosity=1) 

832# # graph(tree, sol) 

833# 

834# 

835# keys = sorted(bd.keys(), key=str) 

836# 

837# # with StdoutCaptured("solarbreakdowns.log"): 

838# # graph(mtree, sol, showlegend=False) 

839# # graph(mtree, sol, showlegend=True) 

840# # tree = crawl(sol.costposy, bd, sol, permissivity=permissivity) 

841# # graph(tree, sol) 

842# # for key in keys: 

843# # tree = crawl(key, bd, sol, permissivity=permissivity) 

844# # graph(tree, sol) 

845# 

846# with StdoutCaptured("solarbreakdowns.log.new"): 

847# graph(mtree, sol, showlegend=False) 

848# graph(mtree, sol, showlegend=True) 

849# tree = crawl(sol.costposy, bd, sol, permissivity=permissivity) 

850# graph(tree, sol) 

851# for key in keys: 

852# tree = crawl(key, bd, sol, permissivity=permissivity) 

853# try: 

854# graph(tree, sol) 

855# except: 

856# raise ValueError(key) 

857# 

858# with open("solarbreakdowns.log", "r") as original: 

859# with open("solarbreakdowns.log.new", "r") as new: 

860# diff = difflib.unified_diff( 

861# original.readlines(), 

862# new.readlines(), 

863# fromfile="original", 

864# tofile="new", 

865# ) 

866# for line in diff: 

867# print(line[:-1]) 

868# 

869# print("SOLAR DONE") 

870 

871 

872sol = pickle.load(open("bd.p", "rb")) 

873bd = get_breakdowns(sol) 

874mbd = get_model_breakdown(sol) 

875mtree = crawl_modelbd(mbd) 

876 

877# graph(mtree, sol, showlegend=False) 

878 

879key, = [vk for vk in bd if "podtripenergy[1,0]" in str(vk)] 

880tree = crawl(key, bd, sol, permissivity=2, verbosity=1) 

881graph(tree, sol, extent=40) 

882tree = crawl(sol.costposy, bd, sol, permissivity=2, verbosity=1) 

883graph(tree, sol) 

884# key, = [vk for vk in bd if "statorlaminationmassperpolepair" in str(vk)] 

885# tree = crawl(key, bd, sol, permissivity=2, verbosity=1) 

886# graph(tree, sol) 

887# key, = [vk for vk in bd if "numberofcellsperstring" in str(vk)] 

888# tree = crawl(key, bd, sol, permissivity=2, verbosity=1) 

889# graph(tree, sol) 

890 

891keys = sorted((key for key in bd.keys() if not key.idx or len(key.shape) == 1), 

892 key=lambda k: k.str_without(excluded={})) 

893 

894# with StdoutCaptured("breakdowns.log"): 

895# graph(mtree, sol, showlegend=False) 

896# graph(mtree.branches[0].branches[1], sol, showlegend=False) 

897# graph(mtree, sol, showlegend=True) 

898# tree = crawl(sol.costposy, bd, sol, permissivity=permissivity) 

899# graph(tree, sol) 

900# for key in keys: 

901# tree = crawl(key, bd, sol, permissivity=permissivity) 

902# graph(tree, sol) 

903 

904with StdoutCaptured("breakdowns.log.new"): 

905 graph(mtree, sol, showlegend=False) 

906 graph(mtree.branches[0].branches[1], sol, showlegend=False) 

907 graph(mtree, sol, showlegend=True) 

908 tree = crawl(sol.costposy, bd, sol, permissivity=permissivity) 

909 graph(tree, sol) 

910 for key in keys: 

911 tree = crawl(key, bd, sol, permissivity=permissivity) 

912 try: 

913 graph(tree, sol) 

914 except: 

915 raise ValueError(key) 

916 

917with open("breakdowns.log", "r") as original: 

918 with open("breakdowns.log.new", "r") as new: 

919 diff = difflib.unified_diff( 

920 original.readlines(), 

921 new.readlines(), 

922 fromfile="original", 

923 tofile="new", 

924 ) 

925 for line in diff: 

926 print(line[:-1]) 

927 

928print("DONE")