Coverage for /pythoncovmergedfiles/medio/medio/usr/local/lib/python3.8/site-packages/scipy/optimize/_shgo

1import numpy as np

2import copy

5class Complex:

6 def __init__(self, dim, func, func_args=(), symmetry=False, bounds=None,

7 g_cons=None, g_args=()):

8 self.dim = dim

9 self.bounds = bounds

10 self.symmetry = symmetry # TODO: Define the functions to be used

11 # here in init to avoid if checks

12 self.gen = 0

13 self.perm_cycle = 0

15 # Every cell is stored in a list of its generation,

16 # e.g., the initial cell is stored in self.H[0]

17 # 1st get new cells are stored in self.H[1] etc.

18 # When a cell is subgenerated it is removed from this list

20 self.H = [] # Storage structure of cells

21 # Cache of all vertices

22 self.V = VertexCache(func, func_args, bounds, g_cons, g_args)

24 # Generate n-cube here:

25 self.n_cube(dim, symmetry=symmetry)

27 # TODO: Assign functions to a the complex instead

28 if symmetry:

29 self.generation_cycle = 1

30 # self.centroid = self.C0()[-1].x

31 # self.C0.centroid = self.centroid

32 else:

33 self.add_centroid()

35 self.H.append([])

36 self.H[0].append(self.C0)

37 self.hgr = self.C0.homology_group_rank()

38 self.hgrd = 0 # Complex group rank differential

39 # self.hgr = self.C0.hg_n

41 # Build initial graph

42 self.graph_map()

44 self.performance = []

45 self.performance.append(0)

46 self.performance.append(0)

48 def __call__(self):

49 return self.H

51 def n_cube(self, dim, symmetry=False, printout=False):

52 """

53 Generate the simplicial triangulation of the N-D hypercube

54 containing 2**n vertices

55 """

56 origin = list(np.zeros(dim, dtype=int))

57 self.origin = origin

58 supremum = list(np.ones(dim, dtype=int))

59 self.supremum = supremum

61 # tuple versions for indexing

62 origintuple = tuple(origin)

63 supremumtuple = tuple(supremum)

65 x_parents = [origintuple]

67 if symmetry:

68 self.C0 = Simplex(0, 0, 0, self.dim) # Initial cell object

69 self.C0.add_vertex(self.V[origintuple])

71 i_s = 0

72 self.perm_symmetry(i_s, x_parents, origin)

73 self.C0.add_vertex(self.V[supremumtuple])

74 else:

75 self.C0 = Cell(0, 0, origin, supremum) # Initial cell object

76 self.C0.add_vertex(self.V[origintuple])

77 self.C0.add_vertex(self.V[supremumtuple])

79 i_parents = []

80 self.perm(i_parents, x_parents, origin)

82 if printout:

83 print("Initial hyper cube:")

84 for v in self.C0():

85 v.print_out()

87 def perm(self, i_parents, x_parents, xi):

88 # TODO: Cut out of for if outside linear constraint cutting planes

89 xi_t = tuple(xi)

91 # Construct required iterator

92 iter_range = [x for x in range(self.dim) if x not in i_parents]

94 for i in iter_range:

95 i2_parents = copy.copy(i_parents)

96 i2_parents.append(i)

97 xi2 = copy.copy(xi)

98 xi2[i] = 1

99 # Make new vertex list a hashable tuple

100 xi2_t = tuple(xi2)

101 # Append to cell

102 self.C0.add_vertex(self.V[xi2_t])

103 # Connect neighbors and vice versa

104 # Parent point

105 self.V[xi2_t].connect(self.V[xi_t])

106

107 # Connect all family of simplices in parent containers

108 for x_ip in x_parents:

109 self.V[xi2_t].connect(self.V[x_ip])

110

111 x_parents2 = copy.copy(x_parents)

112 x_parents2.append(xi_t)

113

114 # Permutate

115 self.perm(i2_parents, x_parents2, xi2)

116

117 def perm_symmetry(self, i_s, x_parents, xi):

118 # TODO: Cut out of for if outside linear constraint cutting planes

119 xi_t = tuple(xi)

120 xi2 = copy.copy(xi)

121 xi2[i_s] = 1

122 # Make new vertex list a hashable tuple

123 xi2_t = tuple(xi2)

124 # Append to cell

125 self.C0.add_vertex(self.V[xi2_t])

126 # Connect neighbors and vice versa

127 # Parent point

128 self.V[xi2_t].connect(self.V[xi_t])

129

130 # Connect all family of simplices in parent containers

131 for x_ip in x_parents:

132 self.V[xi2_t].connect(self.V[x_ip])

133

134 x_parents2 = copy.copy(x_parents)

135 x_parents2.append(xi_t)

136

137 i_s += 1

138 if i_s == self.dim:

139 return

140 # Permutate

141 self.perm_symmetry(i_s, x_parents2, xi2)

142

143 def add_centroid(self):

144 """Split the central edge between the origin and supremum of

145 a cell and add the new vertex to the complex"""

146 self.centroid = list(

147 (np.array(self.origin) + np.array(self.supremum)) / 2.0)

148 self.C0.add_vertex(self.V[tuple(self.centroid)])

149 self.C0.centroid = self.centroid

150

151 # Disconnect origin and supremum

152 self.V[tuple(self.origin)].disconnect(self.V[tuple(self.supremum)])

153

154 # Connect centroid to all other vertices

155 for v in self.C0():

156 self.V[tuple(self.centroid)].connect(self.V[tuple(v.x)])

157

158 self.centroid_added = True

159 return

160

161 # Construct incidence array:

162 def incidence(self):

163 if self.centroid_added:

164 self.structure = np.zeros([2 ** self.dim + 1, 2 ** self.dim + 1],

165 dtype=int)

166 else:

167 self.structure = np.zeros([2 ** self.dim, 2 ** self.dim],

168 dtype=int)

169

170 for v in self.HC.C0():

171 for v2 in v.nn:

172 self.structure[v.index, v2.index] = 1

173

174 return

175

176 # A more sparse incidence generator:

177 def graph_map(self):

178 """ Make a list of size 2**n + 1 where an entry is a vertex

179 incidence, each list element contains a list of indexes

180 corresponding to that entries neighbors"""

181

182 self.graph = [[v2.index for v2 in v.nn] for v in self.C0()]

183

184 # Graph structure method:

185 # 0. Capture the indices of the initial cell.

186 # 1. Generate new origin and supremum scalars based on current generation

187 # 2. Generate a new set of vertices corresponding to a new

188 # "origin" and "supremum"

189 # 3. Connected based on the indices of the previous graph structure

190 # 4. Disconnect the edges in the original cell

191

192 def sub_generate_cell(self, C_i, gen):

193 """Subgenerate a cell `C_i` of generation `gen` and

194 homology group rank `hgr`."""

195 origin_new = tuple(C_i.centroid)

196 centroid_index = len(C_i()) - 1

197

198 # If not gen append

199 try:

200 self.H[gen]

201 except IndexError:

202 self.H.append([])

203

204 # Generate subcubes using every extreme vertex in C_i as a supremum

205 # and the centroid of C_i as the origin

206 H_new = [] # list storing all the new cubes split from C_i

207 for i, v in enumerate(C_i()[:-1]):

208 supremum = tuple(v.x)

209 H_new.append(

210 self.construct_hypercube(origin_new, supremum, gen, C_i.hg_n))

211

212 for i, connections in enumerate(self.graph):

213 # Present vertex V_new[i]; connect to all connections:

214 if i == centroid_index: # Break out of centroid

215 break

216

217 for j in connections:

218 C_i()[i].disconnect(C_i()[j])

219

220 # Destroy the old cell

221 if C_i is not self.C0: # Garbage collector does this anyway; not needed

222 del C_i

223

224 # TODO: Recalculate all the homology group ranks of each cell

225 return H_new

226

227 def split_generation(self):

228 """

229 Run sub_generate_cell for every cell in the current complex self.gen

230 """

231 no_splits = False # USED IN SHGO

232 try:

233 for c in self.H[self.gen]:

234 if self.symmetry:

235 # self.sub_generate_cell_symmetry(c, self.gen + 1)

236 self.split_simplex_symmetry(c, self.gen + 1)

237 else:

238 self.sub_generate_cell(c, self.gen + 1)

239 except IndexError:

240 no_splits = True # USED IN SHGO

241

242 self.gen += 1

243 return no_splits # USED IN SHGO

244

245 def construct_hypercube(self, origin, supremum, gen, hgr,

246 printout=False):

247 """

248 Build a hypercube with triangulations symmetric to C0.

249

250 Parameters

251 ----------

252 origin : vec

253 supremum : vec (tuple)

254 gen : generation

255 hgr : parent homology group rank

256 """

257 # Initiate new cell

258 v_o = np.array(origin)

259 v_s = np.array(supremum)

260

261 C_new = Cell(gen, hgr, origin, supremum)

262 C_new.centroid = tuple((v_o + v_s) * .5)

263

264 # Build new indexed vertex list

265 V_new = []

266

267 for i, v in enumerate(self.C0()[:-1]):

268 v_x = np.array(v.x)

269 sub_cell_t1 = v_o - v_o * v_x

270 sub_cell_t2 = v_s * v_x

271

272 vec = sub_cell_t1 + sub_cell_t2

273

274 vec = tuple(vec)

275 C_new.add_vertex(self.V[vec])

276 V_new.append(vec)

277

278 # Add new centroid

279 C_new.add_vertex(self.V[C_new.centroid])

280 V_new.append(C_new.centroid)

281

282 # Connect new vertices #TODO: Thread into other loop; no need for V_new

283 for i, connections in enumerate(self.graph):

284 # Present vertex V_new[i]; connect to all connections:

285 for j in connections:

286 self.V[V_new[i]].connect(self.V[V_new[j]])

287

288 if printout:

289 print("A sub hyper cube with:")

290 print("origin: {}".format(origin))

291 print("supremum: {}".format(supremum))

292 for v in C_new():

293 v.print_out()

294

295 # Append the new cell to the to complex

296 self.H[gen].append(C_new)

297

298 return C_new

299

300 def split_simplex_symmetry(self, S, gen):

301 """

302 Split a hypersimplex S into two sub simplices by building a hyperplane

303 which connects to a new vertex on an edge (the longest edge in

304 dim = {2, 3}) and every other vertex in the simplex that is not

305 connected to the edge being split.

306

307 This function utilizes the knowledge that the problem is specified

308 with symmetric constraints

309

310 The longest edge is tracked by an ordering of the

311 vertices in every simplices, the edge between first and second

312 vertex is the longest edge to be split in the next iteration.

313 """

314 # If not gen append

315 try:

316 self.H[gen]

317 except IndexError:

318 self.H.append([])

319

320 # Find new vertex.

321 # V_new_x = tuple((np.array(C()[0].x) + np.array(C()[1].x)) / 2.0)

322 s = S()

323 firstx = s[0].x

324 lastx = s[-1].x

325 V_new = self.V[tuple((np.array(firstx) + np.array(lastx)) / 2.0)]

326

327 # Disconnect old longest edge

328 self.V[firstx].disconnect(self.V[lastx])

329

330 # Connect new vertices to all other vertices

331 for v in s[:]:

332 v.connect(self.V[V_new.x])

333

334 # New "lower" simplex

335 S_new_l = Simplex(gen, S.hg_n, self.generation_cycle,

336 self.dim)

337 S_new_l.add_vertex(s[0])

338 S_new_l.add_vertex(V_new) # Add new vertex

339 for v in s[1:-1]: # Add all other vertices

340 S_new_l.add_vertex(v)

341

342 # New "upper" simplex

343 S_new_u = Simplex(gen, S.hg_n, S.generation_cycle, self.dim)

344

345 # First vertex on new long edge

346 S_new_u.add_vertex(s[S_new_u.generation_cycle + 1])

347

348 for v in s[1:-1]: # Remaining vertices

349 S_new_u.add_vertex(v)

350

351 for k, v in enumerate(s[1:-1]): # iterate through inner vertices

352 if k == S.generation_cycle:

353 S_new_u.add_vertex(V_new)

354 else:

355 S_new_u.add_vertex(v)

356

357 S_new_u.add_vertex(s[-1]) # Second vertex on new long edge

358

359 self.H[gen].append(S_new_l)

360 self.H[gen].append(S_new_u)

361

362 return

363

364 # Plots

365 def plot_complex(self):

366 """

367 Here, C is the LIST of simplexes S in the

368 2- or 3-D complex

369

370 To plot a single simplex S in a set C, use e.g., [C[0]]

371 """

372 from matplotlib import pyplot

373 if self.dim == 2:

374 pyplot.figure()

375 for C in self.H:

376 for c in C:

377 for v in c():

378 if self.bounds is None:

379 x_a = np.array(v.x, dtype=float)

380 else:

381 x_a = np.array(v.x, dtype=float)

382 for i in range(len(self.bounds)):

383 x_a[i] = (x_a[i] * (self.bounds[i][1]

384 - self.bounds[i][0])

385 + self.bounds[i][0])

386

387 # logging.info('v.x_a = {}'.format(x_a))

388

389 pyplot.plot([x_a[0]], [x_a[1]], 'o')

390

391 xlines = []

392 ylines = []

393 for vn in v.nn:

394 if self.bounds is None:

395 xn_a = np.array(vn.x, dtype=float)

396 else:

397 xn_a = np.array(vn.x, dtype=float)

398 for i in range(len(self.bounds)):

399 xn_a[i] = (xn_a[i] * (self.bounds[i][1]

400 - self.bounds[i][0])

401 + self.bounds[i][0])

402

403 # logging.info('vn.x = {}'.format(vn.x))

404

405 xlines.append(xn_a[0])

406 ylines.append(xn_a[1])

407 xlines.append(x_a[0])

408 ylines.append(x_a[1])

409

410 pyplot.plot(xlines, ylines)

411

412 if self.bounds is None:

413 pyplot.ylim([-1e-2, 1 + 1e-2])

414 pyplot.xlim([-1e-2, 1 + 1e-2])

415 else:

416 pyplot.ylim(

417 [self.bounds[1][0] - 1e-2, self.bounds[1][1] + 1e-2])

418 pyplot.xlim(

419 [self.bounds[0][0] - 1e-2, self.bounds[0][1] + 1e-2])

420

421 pyplot.show()

422

423 elif self.dim == 3:

424 fig = pyplot.figure()

425 ax = fig.add_subplot(111, projection='3d')

426

427 for C in self.H:

428 for c in C:

429 for v in c():

430 x = []

431 y = []

432 z = []

433 # logging.info('v.x = {}'.format(v.x))

434 x.append(v.x[0])

435 y.append(v.x[1])

436 z.append(v.x[2])

437 for vn in v.nn:

438 x.append(vn.x[0])

439 y.append(vn.x[1])

440 z.append(vn.x[2])

441 x.append(v.x[0])

442 y.append(v.x[1])

443 z.append(v.x[2])

444 # logging.info('vn.x = {}'.format(vn.x))

445

446 ax.plot(x, y, z, label='simplex')

447

448 pyplot.show()

449 else:

450 print("dimension higher than 3 or wrong complex format")

451 return

452

453

454class VertexGroup:

455 def __init__(self, p_gen, p_hgr):

456 self.p_gen = p_gen # parent generation

457 self.p_hgr = p_hgr # parent homology group rank

458 self.hg_n = None

459 self.hg_d = None

460

461 # Maybe add parent homology group rank total history

462 # This is the sum off all previously split cells

463 # cumulatively throughout its entire history

464 self.C = []

465

466 def __call__(self):

467 return self.C

468

469 def add_vertex(self, V):

470 if V not in self.C:

471 self.C.append(V)

472

473 def homology_group_rank(self):

474 """

475 Returns the homology group order of the current cell

476 """

477 if self.hg_n is None:

478 self.hg_n = sum(1 for v in self.C if v.minimiser())

479

480 return self.hg_n

481

482 def homology_group_differential(self):

483 """

484 Returns the difference between the current homology group of the

485 cell and its parent group

486 """

487 if self.hg_d is None:

488 self.hgd = self.hg_n - self.p_hgr

489

490 return self.hgd

491

492 def polytopial_sperner_lemma(self):

493 """

494 Returns the number of stationary points theoretically contained in the

495 cell based information currently known about the cell

496 """

497 pass

498

499 def print_out(self):

500 """

501 Print the current cell to console

502 """

503 for v in self():

504 v.print_out()

505

506

507class Cell(VertexGroup):

508 """

509 Contains a cell that is symmetric to the initial hypercube triangulation

510 """

511

512 def __init__(self, p_gen, p_hgr, origin, supremum):

513 super().__init__(p_gen, p_hgr)

514

515 self.origin = origin

516 self.supremum = supremum

517 self.centroid = None # (Not always used)

518 # TODO: self.bounds

519

520

521class Simplex(VertexGroup):

522 """

523 Contains a simplex that is symmetric to the initial symmetry constrained

524 hypersimplex triangulation

525 """

526

527 def __init__(self, p_gen, p_hgr, generation_cycle, dim):

528 super().__init__(p_gen, p_hgr)

529

530 self.generation_cycle = (generation_cycle + 1) % (dim - 1)

531

532

533class Vertex:

534 def __init__(self, x, bounds=None, func=None, func_args=(), g_cons=None,

535 g_cons_args=(), nn=None, index=None):

536 self.x = x

537 self.order = sum(x)

538 x_a = np.array(x, dtype=float)

539 if bounds is not None:

540 for i, (lb, ub) in enumerate(bounds):

541 x_a[i] = x_a[i] * (ub - lb) + lb

542

543 # TODO: Make saving the array structure optional

544 self.x_a = x_a

545

546 # Note Vertex is only initiated once for all x so only

547 # evaluated once

548 if func is not None:

549 self.feasible = True

550 if g_cons is not None:

551 for g, args in zip(g_cons, g_cons_args):

552 if g(self.x_a, *args) < 0.0:

553 self.f = np.inf

554 self.feasible = False

555 break

556 if self.feasible:

557 self.f = func(x_a, *func_args)

558

559 if nn is not None:

560 self.nn = nn

561 else:

562 self.nn = set()

563

564 self.fval = None

565 self.check_min = True

566

567 # Index:

568 if index is not None:

569 self.index = index

570

571 def __hash__(self):

572 return hash(self.x)

573

574 def connect(self, v):

575 if v is not self and v not in self.nn:

576 self.nn.add(v)

577 v.nn.add(self)

578

579 if self.minimiser():

580 v._min = False

581 v.check_min = False

582

583 # TEMPORARY

584 self.check_min = True

585 v.check_min = True

586

587 def disconnect(self, v):

588 if v in self.nn:

589 self.nn.remove(v)

590 v.nn.remove(self)

591 self.check_min = True

592 v.check_min = True

593

594 def minimiser(self):

595 """Check whether this vertex is strictly less than all its neighbors"""

596 if self.check_min:

597 self._min = all(self.f < v.f for v in self.nn)

598 self.check_min = False

599

600 return self._min

601

602 def print_out(self):

603 print("Vertex: {}".format(self.x))

604 constr = 'Connections: '

605 for vc in self.nn:

606 constr += '{} '.format(vc.x)

607

608 print(constr)

609 print('Order = {}'.format(self.order))

610

611

612class VertexCache:

613 def __init__(self, func, func_args=(), bounds=None, g_cons=None,

614 g_cons_args=(), indexed=True):

615

616 self.cache = {}

617 self.func = func

618 self.g_cons = g_cons

619 self.g_cons_args = g_cons_args

620 self.func_args = func_args

621 self.bounds = bounds

622 self.nfev = 0

623 self.size = 0

624

625 if indexed:

626 self.index = -1

627

628 def __getitem__(self, x, indexed=True):

629 try:

630 return self.cache[x]

631 except KeyError:

632 if indexed:

633 self.index += 1

634 xval = Vertex(x, bounds=self.bounds,

635 func=self.func, func_args=self.func_args,

636 g_cons=self.g_cons,

637 g_cons_args=self.g_cons_args,

638 index=self.index)

639 else:

640 xval = Vertex(x, bounds=self.bounds,

641 func=self.func, func_args=self.func_args,

642 g_cons=self.g_cons,

643 g_cons_args=self.g_cons_args)

644

645 # logging.info("New generated vertex at x = {}".format(x))

646 # NOTE: Surprisingly high performance increase if logging is commented out

647 self.cache[x] = xval

648

649 # TODO: Check

650 if self.func is not None:

651 if self.g_cons is not None:

652 if xval.feasible:

653 self.nfev += 1

654 self.size += 1

655 else:

656 self.size += 1

657 else:

658 self.nfev += 1

659 self.size += 1

660

661 return self.cache[x]

Coverage for /pythoncovmergedfiles/medio/medio/usr/local/lib/python3.8/site-packages/scipy/optimize/_shgo_lib/triangulation.py: 11%

359 statements