Coverage for /usr/lib/python3/dist-packages/scipy/optimize/_shgo_lib/_complex.py: 3%
666 statements
« prev ^ index » next coverage.py v7.9.1, created at 2025-06-14 15:55 +0200
1"""Base classes for low memory simplicial complex structures."""
2import copy
3import logging
4import itertools
5import decimal
6from functools import cache
8import numpy
10from ._vertex import (VertexCacheField, VertexCacheIndex)
13class Complex:
14 """
15 Base class for a simplicial complex described as a cache of vertices
16 together with their connections.
18 Important methods:
19 Domain triangulation:
20 Complex.triangulate, Complex.split_generation
21 Triangulating arbitrary points (must be traingulable,
22 may exist outside domain):
23 Complex.triangulate(sample_set)
24 Converting another simplicial complex structure data type to the
25 structure used in Complex (ex. OBJ wavefront)
26 Complex.convert(datatype, data)
28 Important objects:
29 HC.V: The cache of vertices and their connection
30 HC.H: Storage structure of all vertex groups
32 Parameters
33 ----------
34 dim : int
35 Spatial dimensionality of the complex R^dim
36 domain : list of tuples, optional
37 The bounds [x_l, x_u]^dim of the hyperrectangle space
38 ex. The default domain is the hyperrectangle [0, 1]^dim
39 Note: The domain must be convex, non-convex spaces can be cut
40 away from this domain using the non-linear
41 g_cons functions to define any arbitrary domain
42 (these domains may also be disconnected from each other)
43 sfield :
44 A scalar function defined in the associated domain f: R^dim --> R
45 sfield_args : tuple
46 Additional arguments to be passed to `sfield`
47 vfield :
48 A scalar function defined in the associated domain
49 f: R^dim --> R^m
50 (for example a gradient function of the scalar field)
51 vfield_args : tuple
52 Additional arguments to be passed to vfield
53 symmetry : None or list
54 Specify if the objective function contains symmetric variables.
55 The search space (and therefore performance) is decreased by up to
56 O(n!) times in the fully symmetric case.
58 E.g. f(x) = (x_1 + x_2 + x_3) + (x_4)**2 + (x_5)**2 + (x_6)**2
60 In this equation x_2 and x_3 are symmetric to x_1, while x_5 and
61 x_6 are symmetric to x_4, this can be specified to the solver as:
63 symmetry = [0, # Variable 1
64 0, # symmetric to variable 1
65 0, # symmetric to variable 1
66 3, # Variable 4
67 3, # symmetric to variable 4
68 3, # symmetric to variable 4
69 ]
71 constraints : dict or sequence of dict, optional
72 Constraints definition.
73 Function(s) ``R**n`` in the form::
75 g(x) <= 0 applied as g : R^n -> R^m
76 h(x) == 0 applied as h : R^n -> R^p
78 Each constraint is defined in a dictionary with fields:
80 type : str
81 Constraint type: 'eq' for equality, 'ineq' for inequality.
82 fun : callable
83 The function defining the constraint.
84 jac : callable, optional
85 The Jacobian of `fun` (only for SLSQP).
86 args : sequence, optional
87 Extra arguments to be passed to the function and Jacobian.
89 Equality constraint means that the constraint function result is to
90 be zero whereas inequality means that it is to be
91 non-negative.constraints : dict or sequence of dict, optional
92 Constraints definition.
93 Function(s) ``R**n`` in the form::
95 g(x) <= 0 applied as g : R^n -> R^m
96 h(x) == 0 applied as h : R^n -> R^p
98 Each constraint is defined in a dictionary with fields:
100 type : str
101 Constraint type: 'eq' for equality, 'ineq' for inequality.
102 fun : callable
103 The function defining the constraint.
104 jac : callable, optional
105 The Jacobian of `fun` (unused).
106 args : sequence, optional
107 Extra arguments to be passed to the function and Jacobian.
109 Equality constraint means that the constraint function result is to
110 be zero whereas inequality means that it is to be non-negative.
112 workers : int optional
113 Uses `multiprocessing.Pool <multiprocessing>`) to compute the field
114 functions in parrallel.
115 """
116 def __init__(self, dim, domain=None, sfield=None, sfield_args=(),
117 symmetry=None, constraints=None, workers=1):
118 self.dim = dim
120 # Domains
121 self.domain = domain
122 if domain is None:
123 self.bounds = [(float(0), float(1.0)), ] * dim
124 else:
125 self.bounds = domain
126 self.symmetry = symmetry
127 # here in init to avoid if checks
129 # Field functions
130 self.sfield = sfield
131 self.sfield_args = sfield_args
133 # Process constraints
134 # Constraints
135 # Process constraint dict sequence:
136 if constraints is not None:
137 self.min_cons = constraints
138 self.g_cons = []
139 self.g_args = []
140 if (type(constraints) is not tuple) and (type(constraints)
141 is not list):
142 constraints = (constraints,)
144 for cons in constraints:
145 if cons['type'] in ('ineq'):
146 self.g_cons.append(cons['fun'])
147 try:
148 self.g_args.append(cons['args'])
149 except KeyError:
150 self.g_args.append(())
151 self.g_cons = tuple(self.g_cons)
152 self.g_args = tuple(self.g_args)
153 else:
154 self.g_cons = None
155 self.g_args = None
157 # Homology properties
158 self.gen = 0
159 self.perm_cycle = 0
161 # Every cell is stored in a list of its generation,
162 # ex. the initial cell is stored in self.H[0]
163 # 1st get new cells are stored in self.H[1] etc.
164 # When a cell is sub-generated it is removed from this list
166 self.H = [] # Storage structure of vertex groups
168 # Cache of all vertices
169 if (sfield is not None) or (self.g_cons is not None):
170 # Initiate a vertex cache and an associated field cache, note that
171 # the field case is always initiated inside the vertex cache if an
172 # associated field scalar field is defined:
173 if sfield is not None:
174 self.V = VertexCacheField(field=sfield, field_args=sfield_args,
175 g_cons=self.g_cons,
176 g_cons_args=self.g_args,
177 workers=workers)
178 elif self.g_cons is not None:
179 self.V = VertexCacheField(field=sfield, field_args=sfield_args,
180 g_cons=self.g_cons,
181 g_cons_args=self.g_args,
182 workers=workers)
183 else:
184 self.V = VertexCacheIndex()
186 self.V_non_symm = [] # List of non-symmetric vertices
188 def __call__(self):
189 return self.H
191 # %% Triangulation methods
192 def cyclic_product(self, bounds, origin, supremum, centroid=True):
193 """Generate initial triangulation using cyclic product"""
194 # Define current hyperrectangle
195 vot = tuple(origin)
196 vut = tuple(supremum) # Hyperrectangle supremum
197 self.V[vot]
198 vo = self.V[vot]
199 yield vo.x
200 self.V[vut].connect(self.V[vot])
201 yield vut
202 # Cyclic group approach with second x_l --- x_u operation.
204 # These containers store the "lower" and "upper" vertices
205 # corresponding to the origin or supremum of every C2 group.
206 # It has the structure of `dim` times embedded lists each containing
207 # these vertices as the entire complex grows. Bounds[0] has to be done
208 # outside the loops before we have symmetric containers.
209 # NOTE: This means that bounds[0][1] must always exist
210 C0x = [[self.V[vot]]]
211 a_vo = copy.copy(list(origin))
212 a_vo[0] = vut[0] # Update aN Origin
213 a_vo = self.V[tuple(a_vo)]
214 # self.V[vot].connect(self.V[tuple(a_vo)])
215 self.V[vot].connect(a_vo)
216 yield a_vo.x
217 C1x = [[a_vo]]
218 # C1x = [[self.V[tuple(a_vo)]]]
219 ab_C = [] # Container for a + b operations
221 # Loop over remaining bounds
222 for i, x in enumerate(bounds[1:]):
223 # Update lower and upper containers
224 C0x.append([])
225 C1x.append([])
226 # try to access a second bound (if not, C1 is symmetric)
227 try:
228 # Early try so that we don't have to copy the cache before
229 # moving on to next C1/C2: Try to add the operation of a new
230 # C2 product by accessing the upper bound
231 x[1]
232 # Copy lists for iteration
233 cC0x = [x[:] for x in C0x[:i + 1]]
234 cC1x = [x[:] for x in C1x[:i + 1]]
235 for j, (VL, VU) in enumerate(zip(cC0x, cC1x)):
236 for k, (vl, vu) in enumerate(zip(VL, VU)):
237 # Build aN vertices for each lower-upper pair in N:
238 a_vl = list(vl.x)
239 a_vu = list(vu.x)
240 a_vl[i + 1] = vut[i + 1]
241 a_vu[i + 1] = vut[i + 1]
242 a_vl = self.V[tuple(a_vl)]
244 # Connect vertices in N to corresponding vertices
245 # in aN:
246 vl.connect(a_vl)
248 yield a_vl.x
250 a_vu = self.V[tuple(a_vu)]
251 # Connect vertices in N to corresponding vertices
252 # in aN:
253 vu.connect(a_vu)
255 # Connect new vertex pair in aN:
256 a_vl.connect(a_vu)
258 # Connect lower pair to upper (triangulation
259 # operation of a + b (two arbitrary operations):
260 vl.connect(a_vu)
261 ab_C.append((vl, a_vu))
263 # Update the containers
264 C0x[i + 1].append(vl)
265 C0x[i + 1].append(vu)
266 C1x[i + 1].append(a_vl)
267 C1x[i + 1].append(a_vu)
269 # Update old containers
270 C0x[j].append(a_vl)
271 C1x[j].append(a_vu)
273 # Yield new points
274 yield a_vu.x
276 # Try to connect aN lower source of previous a + b
277 # operation with a aN vertex
278 ab_Cc = copy.copy(ab_C)
280 for vp in ab_Cc:
281 b_v = list(vp[0].x)
282 ab_v = list(vp[1].x)
283 b_v[i + 1] = vut[i + 1]
284 ab_v[i + 1] = vut[i + 1]
285 b_v = self.V[tuple(b_v)] # b + vl
286 ab_v = self.V[tuple(ab_v)] # b + a_vl
287 # Note o---o is already connected
288 vp[0].connect(ab_v) # o-s
289 b_v.connect(ab_v) # s-s
291 # Add new list of cross pairs
292 ab_C.append((vp[0], ab_v))
293 ab_C.append((b_v, ab_v))
295 except IndexError:
296 cC0x = C0x[i]
297 cC1x = C1x[i]
298 VL, VU = cC0x, cC1x
299 for k, (vl, vu) in enumerate(zip(VL, VU)):
300 # Build aN vertices for each lower-upper pair in N:
301 a_vu = list(vu.x)
302 a_vu[i + 1] = vut[i + 1]
303 # Connect vertices in N to corresponding vertices
304 # in aN:
305 a_vu = self.V[tuple(a_vu)]
306 # Connect vertices in N to corresponding vertices
307 # in aN:
308 vu.connect(a_vu)
309 # Connect new vertex pair in aN:
310 # a_vl.connect(a_vu)
311 # Connect lower pair to upper (triangulation
312 # operation of a + b (two arbitrary operations):
313 vl.connect(a_vu)
314 ab_C.append((vl, a_vu))
315 C0x[i + 1].append(vu)
316 C1x[i + 1].append(a_vu)
317 # Yield new points
318 a_vu.connect(self.V[vut])
319 yield a_vu.x
320 ab_Cc = copy.copy(ab_C)
321 for vp in ab_Cc:
322 if vp[1].x[i] == vut[i]:
323 ab_v = list(vp[1].x)
324 ab_v[i + 1] = vut[i + 1]
325 ab_v = self.V[tuple(ab_v)] # b + a_vl
326 # Note o---o is already connected
327 vp[0].connect(ab_v) # o-s
329 # Add new list of cross pairs
330 ab_C.append((vp[0], ab_v))
332 # Clean class trash
333 try:
334 del C0x
335 del cC0x
336 del C1x
337 del cC1x
338 del ab_C
339 del ab_Cc
340 except UnboundLocalError:
341 pass
343 # Extra yield to ensure that the triangulation is completed
344 if centroid:
345 vo = self.V[vot]
346 vs = self.V[vut]
347 # Disconnect the origin and supremum
348 vo.disconnect(vs)
349 # Build centroid
350 vc = self.split_edge(vot, vut)
351 for v in vo.nn:
352 v.connect(vc)
353 yield vc.x
354 return vc.x
355 else:
356 yield vut
357 return vut
359 def triangulate(self, n=None, symmetry=None, centroid=True,
360 printout=False):
361 """
362 Triangulate the initial domain, if n is not None then a limited number
363 of points will be generated
365 Parameters
366 ----------
367 n : int, Number of points to be sampled.
368 symmetry :
370 Ex. Dictionary/hashtable
371 f(x) = (x_1 + x_2 + x_3) + (x_4)**2 + (x_5)**2 + (x_6)**2
373 symmetry = symmetry[0]: 0, # Variable 1
374 symmetry[1]: 0, # symmetric to variable 1
375 symmetry[2]: 0, # symmetric to variable 1
376 symmetry[3]: 3, # Variable 4
377 symmetry[4]: 3, # symmetric to variable 4
378 symmetry[5]: 3, # symmetric to variable 4
379 }
380 centroid : bool, if True add a central point to the hypercube
381 printout : bool, if True print out results
383 NOTES:
384 ------
385 Rather than using the combinatorial algorithm to connect vertices we
386 make the following observation:
388 The bound pairs are similar a C2 cyclic group and the structure is
389 formed using the cartesian product:
391 H = C2 x C2 x C2 ... x C2 (dim times)
393 So construct any normal subgroup N and consider H/N first, we connect
394 all vertices within N (ex. N is C2 (the first dimension), then we move
395 to a left coset aN (an operation moving around the defined H/N group by
396 for example moving from the lower bound in C2 (dimension 2) to the
397 higher bound in C2. During this operation connection all the vertices.
398 Now repeat the N connections. Note that these elements can be connected
399 in parrallel.
400 """
401 # Inherit class arguments
402 if symmetry is None:
403 symmetry = self.symmetry
404 # Build origin and supremum vectors
405 origin = [i[0] for i in self.bounds]
406 self.origin = origin
407 supremum = [i[1] for i in self.bounds]
409 self.supremum = supremum
411 if symmetry is None:
412 cbounds = self.bounds
413 else:
414 cbounds = copy.copy(self.bounds)
415 for i, j in enumerate(symmetry):
416 if i is not j:
417 # pop second entry on second symmetry vars
418 cbounds[i] = [self.bounds[symmetry[i]][0]]
419 # Sole (first) entry is the sup value and there is no
420 # origin:
421 cbounds[i] = [self.bounds[symmetry[i]][1]]
422 if (self.bounds[symmetry[i]] is not
423 self.bounds[symmetry[j]]):
424 logging.warning(f"Variable {i} was specified as "
425 f"symmetetric to variable {j}, however"
426 f", the bounds {i} ="
427 f" {self.bounds[symmetry[i]]} and {j}"
428 f" ="
429 f" {self.bounds[symmetry[j]]} do not "
430 f"match, the mismatch was ignored in "
431 f"the initial triangulation.")
432 cbounds[i] = self.bounds[symmetry[j]]
434 if n is None:
435 # Build generator
436 self.cp = self.cyclic_product(cbounds, origin, supremum, centroid)
437 for i in self.cp:
438 i
440 try:
441 self.triangulated_vectors.append((tuple(self.origin),
442 tuple(self.supremum)))
443 except (AttributeError, KeyError):
444 self.triangulated_vectors = [(tuple(self.origin),
445 tuple(self.supremum))]
447 else:
448 # Check if generator already exists
449 try:
450 self.cp
451 except (AttributeError, KeyError):
452 self.cp = self.cyclic_product(cbounds, origin, supremum,
453 centroid)
455 try:
456 while len(self.V.cache) < n:
457 next(self.cp)
458 except StopIteration:
459 try:
460 self.triangulated_vectors.append((tuple(self.origin),
461 tuple(self.supremum)))
462 except (AttributeError, KeyError):
463 self.triangulated_vectors = [(tuple(self.origin),
464 tuple(self.supremum))]
466 if printout:
467 # for v in self.C0():
468 # v.print_out()
469 for v in self.V.cache:
470 self.V[v].print_out()
472 return
474 def refine(self, n=1):
475 if n is None:
476 try:
477 self.triangulated_vectors
478 self.refine_all()
479 return
480 except AttributeError as ae:
481 if str(ae) == "'Complex' object has no attribute " \
482 "'triangulated_vectors'":
483 self.triangulate(symmetry=self.symmetry)
484 return
485 else:
486 raise
488 nt = len(self.V.cache) + n # Target number of total vertices
489 # In the outer while loop we iterate until we have added an extra `n`
490 # vertices to the complex:
491 while len(self.V.cache) < nt: # while loop 1
492 try: # try 1
493 # Try to access triangulated_vectors, this should only be
494 # defined if an initial triangulation has already been
495 # performed:
496 self.triangulated_vectors
497 # Try a usual iteration of the current generator, if it
498 # does not exist or is exhausted then produce a new generator
499 try: # try 2
500 next(self.rls)
501 except (AttributeError, StopIteration, KeyError):
502 vp = self.triangulated_vectors[0]
503 self.rls = self.refine_local_space(*vp, bounds=self.bounds)
504 next(self.rls)
506 except (AttributeError, KeyError):
507 # If an initial triangulation has not been completed, then
508 # we start/continue the initial triangulation targeting `nt`
509 # vertices, if nt is greater than the initial number of
510 # vertices then the `refine` routine will move back to try 1.
511 self.triangulate(nt, self.symmetry)
512 return
514 def refine_all(self, centroids=True):
515 """Refine the entire domain of the current complex."""
516 try:
517 self.triangulated_vectors
518 tvs = copy.copy(self.triangulated_vectors)
519 for i, vp in enumerate(tvs):
520 self.rls = self.refine_local_space(*vp, bounds=self.bounds)
521 for i in self.rls:
522 i
523 except AttributeError as ae:
524 if str(ae) == "'Complex' object has no attribute " \
525 "'triangulated_vectors'":
526 self.triangulate(symmetry=self.symmetry, centroid=centroids)
527 else:
528 raise
530 # This adds a centroid to every new sub-domain generated and defined
531 # by self.triangulated_vectors, in addition the vertices ! to complete
532 # the triangulation
533 return
535 def refine_local_space(self, origin, supremum, bounds, centroid=1):
536 # Copy for later removal
537 origin_c = copy.copy(origin)
538 supremum_c = copy.copy(supremum)
540 # Initiate local variables redefined in later inner `for` loop:
541 vl, vu, a_vu = None, None, None
543 # Change the vector orientation so that it is only increasing
544 s_ov = list(origin)
545 s_origin = list(origin)
546 s_sv = list(supremum)
547 s_supremum = list(supremum)
548 for i, vi in enumerate(s_origin):
549 if s_ov[i] > s_sv[i]:
550 s_origin[i] = s_sv[i]
551 s_supremum[i] = s_ov[i]
553 vot = tuple(s_origin)
554 vut = tuple(s_supremum) # Hyperrectangle supremum
556 vo = self.V[vot] # initiate if doesn't exist yet
557 vs = self.V[vut]
558 # Start by finding the old centroid of the new space:
559 vco = self.split_edge(vo.x, vs.x) # Split in case not centroid arg
561 # Find set of extreme vertices in current local space
562 sup_set = copy.copy(vco.nn)
563 # Cyclic group approach with second x_l --- x_u operation.
565 # These containers store the "lower" and "upper" vertices
566 # corresponding to the origin or supremum of every C2 group.
567 # It has the structure of `dim` times embedded lists each containing
568 # these vertices as the entire complex grows. Bounds[0] has to be done
569 # outside the loops before we have symmetric containers.
570 # NOTE: This means that bounds[0][1] must always exist
572 a_vl = copy.copy(list(vot))
573 a_vl[0] = vut[0] # Update aN Origin
574 if tuple(a_vl) not in self.V.cache:
575 vo = self.V[vot] # initiate if doesn't exist yet
576 vs = self.V[vut]
577 # Start by finding the old centroid of the new space:
578 vco = self.split_edge(vo.x, vs.x) # Split in case not centroid arg
580 # Find set of extreme vertices in current local space
581 sup_set = copy.copy(vco.nn)
582 a_vl = copy.copy(list(vot))
583 a_vl[0] = vut[0] # Update aN Origin
584 a_vl = self.V[tuple(a_vl)]
585 else:
586 a_vl = self.V[tuple(a_vl)]
588 c_v = self.split_edge(vo.x, a_vl.x)
589 c_v.connect(vco)
590 yield c_v.x
591 Cox = [[vo]]
592 Ccx = [[c_v]]
593 Cux = [[a_vl]]
594 ab_C = [] # Container for a + b operations
595 s_ab_C = [] # Container for symmetric a + b operations
597 # Loop over remaining bounds
598 for i, x in enumerate(bounds[1:]):
599 # Update lower and upper containers
600 Cox.append([])
601 Ccx.append([])
602 Cux.append([])
603 # try to access a second bound (if not, C1 is symmetric)
604 try:
605 t_a_vl = list(vot)
606 t_a_vl[i + 1] = vut[i + 1]
608 # New: lists are used anyway, so copy all
609 # %%
610 # Copy lists for iteration
611 cCox = [x[:] for x in Cox[:i + 1]]
612 cCcx = [x[:] for x in Ccx[:i + 1]]
613 cCux = [x[:] for x in Cux[:i + 1]]
614 # Try to connect aN lower source of previous a + b
615 # operation with a aN vertex
616 ab_Cc = copy.copy(ab_C) # NOTE: We append ab_C in the
617 # (VL, VC, VU) for-loop, but we use the copy of the list in the
618 # ab_Cc for-loop.
619 s_ab_Cc = copy.copy(s_ab_C)
621 # Early try so that we don't have to copy the cache before
622 # moving on to next C1/C2: Try to add the operation of a new
623 # C2 product by accessing the upper bound
624 if tuple(t_a_vl) not in self.V.cache:
625 # Raise error to continue symmetric refine
626 raise IndexError
627 t_a_vu = list(vut)
628 t_a_vu[i + 1] = vut[i + 1]
629 if tuple(t_a_vu) not in self.V.cache:
630 # Raise error to continue symmetric refine:
631 raise IndexError
633 for vectors in s_ab_Cc:
634 # s_ab_C.append([c_vc, vl, vu, a_vu])
635 bc_vc = list(vectors[0].x)
636 b_vl = list(vectors[1].x)
637 b_vu = list(vectors[2].x)
638 ba_vu = list(vectors[3].x)
640 bc_vc[i + 1] = vut[i + 1]
641 b_vl[i + 1] = vut[i + 1]
642 b_vu[i + 1] = vut[i + 1]
643 ba_vu[i + 1] = vut[i + 1]
645 bc_vc = self.V[tuple(bc_vc)]
646 bc_vc.connect(vco) # NOTE: Unneeded?
647 yield bc_vc
649 # Split to centre, call this centre group "d = 0.5*a"
650 d_bc_vc = self.split_edge(vectors[0].x, bc_vc.x)
651 d_bc_vc.connect(bc_vc)
652 d_bc_vc.connect(vectors[1]) # Connect all to centroid
653 d_bc_vc.connect(vectors[2]) # Connect all to centroid
654 d_bc_vc.connect(vectors[3]) # Connect all to centroid
655 yield d_bc_vc.x
656 b_vl = self.V[tuple(b_vl)]
657 bc_vc.connect(b_vl) # Connect aN cross pairs
658 d_bc_vc.connect(b_vl) # Connect all to centroid
660 yield b_vl
661 b_vu = self.V[tuple(b_vu)]
662 bc_vc.connect(b_vu) # Connect aN cross pairs
663 d_bc_vc.connect(b_vu) # Connect all to centroid
665 b_vl_c = self.split_edge(b_vu.x, b_vl.x)
666 bc_vc.connect(b_vl_c)
668 yield b_vu
669 ba_vu = self.V[tuple(ba_vu)]
670 bc_vc.connect(ba_vu) # Connect aN cross pairs
671 d_bc_vc.connect(ba_vu) # Connect all to centroid
673 # Split the a + b edge of the initial triangulation:
674 os_v = self.split_edge(vectors[1].x, ba_vu.x) # o-s
675 ss_v = self.split_edge(b_vl.x, ba_vu.x) # s-s
676 b_vu_c = self.split_edge(b_vu.x, ba_vu.x)
677 bc_vc.connect(b_vu_c)
678 yield os_v.x # often equal to vco, but not always
679 yield ss_v.x # often equal to bc_vu, but not always
680 yield ba_vu
681 # Split remaining to centre, call this centre group
682 # "d = 0.5*a"
683 d_bc_vc = self.split_edge(vectors[0].x, bc_vc.x)
684 d_bc_vc.connect(vco) # NOTE: Unneeded?
685 yield d_bc_vc.x
686 d_b_vl = self.split_edge(vectors[1].x, b_vl.x)
687 d_bc_vc.connect(vco) # NOTE: Unneeded?
688 d_bc_vc.connect(d_b_vl) # Connect dN cross pairs
689 yield d_b_vl.x
690 d_b_vu = self.split_edge(vectors[2].x, b_vu.x)
691 d_bc_vc.connect(vco) # NOTE: Unneeded?
692 d_bc_vc.connect(d_b_vu) # Connect dN cross pairs
693 yield d_b_vu.x
694 d_ba_vu = self.split_edge(vectors[3].x, ba_vu.x)
695 d_bc_vc.connect(vco) # NOTE: Unneeded?
696 d_bc_vc.connect(d_ba_vu) # Connect dN cross pairs
697 yield d_ba_vu
699 # comb = [c_vc, vl, vu, a_vl, a_vu,
700 # bc_vc, b_vl, b_vu, ba_vl, ba_vu]
701 comb = [vl, vu, a_vu,
702 b_vl, b_vu, ba_vu]
703 comb_iter = itertools.combinations(comb, 2)
704 for vecs in comb_iter:
705 self.split_edge(vecs[0].x, vecs[1].x)
706 # Add new list of cross pairs
707 ab_C.append((d_bc_vc, vectors[1], b_vl, a_vu, ba_vu))
708 ab_C.append((d_bc_vc, vl, b_vl, a_vu, ba_vu)) # = prev
710 for vectors in ab_Cc:
711 bc_vc = list(vectors[0].x)
712 b_vl = list(vectors[1].x)
713 b_vu = list(vectors[2].x)
714 ba_vl = list(vectors[3].x)
715 ba_vu = list(vectors[4].x)
716 bc_vc[i + 1] = vut[i + 1]
717 b_vl[i + 1] = vut[i + 1]
718 b_vu[i + 1] = vut[i + 1]
719 ba_vl[i + 1] = vut[i + 1]
720 ba_vu[i + 1] = vut[i + 1]
721 bc_vc = self.V[tuple(bc_vc)]
722 bc_vc.connect(vco) # NOTE: Unneeded?
723 yield bc_vc
725 # Split to centre, call this centre group "d = 0.5*a"
726 d_bc_vc = self.split_edge(vectors[0].x, bc_vc.x)
727 d_bc_vc.connect(bc_vc)
728 d_bc_vc.connect(vectors[1]) # Connect all to centroid
729 d_bc_vc.connect(vectors[2]) # Connect all to centroid
730 d_bc_vc.connect(vectors[3]) # Connect all to centroid
731 d_bc_vc.connect(vectors[4]) # Connect all to centroid
732 yield d_bc_vc.x
733 b_vl = self.V[tuple(b_vl)]
734 bc_vc.connect(b_vl) # Connect aN cross pairs
735 d_bc_vc.connect(b_vl) # Connect all to centroid
736 yield b_vl
737 b_vu = self.V[tuple(b_vu)]
738 bc_vc.connect(b_vu) # Connect aN cross pairs
739 d_bc_vc.connect(b_vu) # Connect all to centroid
740 yield b_vu
741 ba_vl = self.V[tuple(ba_vl)]
742 bc_vc.connect(ba_vl) # Connect aN cross pairs
743 d_bc_vc.connect(ba_vl) # Connect all to centroid
744 self.split_edge(b_vu.x, ba_vl.x)
745 yield ba_vl
746 ba_vu = self.V[tuple(ba_vu)]
747 bc_vc.connect(ba_vu) # Connect aN cross pairs
748 d_bc_vc.connect(ba_vu) # Connect all to centroid
749 # Split the a + b edge of the initial triangulation:
750 os_v = self.split_edge(vectors[1].x, ba_vu.x) # o-s
751 ss_v = self.split_edge(b_vl.x, ba_vu.x) # s-s
752 yield os_v.x # often equal to vco, but not always
753 yield ss_v.x # often equal to bc_vu, but not always
754 yield ba_vu
755 # Split remaining to centre, call this centre group
756 # "d = 0.5*a"
757 d_bc_vc = self.split_edge(vectors[0].x, bc_vc.x)
758 d_bc_vc.connect(vco) # NOTE: Unneeded?
759 yield d_bc_vc.x
760 d_b_vl = self.split_edge(vectors[1].x, b_vl.x)
761 d_bc_vc.connect(vco) # NOTE: Unneeded?
762 d_bc_vc.connect(d_b_vl) # Connect dN cross pairs
763 yield d_b_vl.x
764 d_b_vu = self.split_edge(vectors[2].x, b_vu.x)
765 d_bc_vc.connect(vco) # NOTE: Unneeded?
766 d_bc_vc.connect(d_b_vu) # Connect dN cross pairs
767 yield d_b_vu.x
768 d_ba_vl = self.split_edge(vectors[3].x, ba_vl.x)
769 d_bc_vc.connect(vco) # NOTE: Unneeded?
770 d_bc_vc.connect(d_ba_vl) # Connect dN cross pairs
771 yield d_ba_vl
772 d_ba_vu = self.split_edge(vectors[4].x, ba_vu.x)
773 d_bc_vc.connect(vco) # NOTE: Unneeded?
774 d_bc_vc.connect(d_ba_vu) # Connect dN cross pairs
775 yield d_ba_vu
776 c_vc, vl, vu, a_vl, a_vu = vectors
778 comb = [vl, vu, a_vl, a_vu,
779 b_vl, b_vu, ba_vl, ba_vu]
780 comb_iter = itertools.combinations(comb, 2)
781 for vecs in comb_iter:
782 self.split_edge(vecs[0].x, vecs[1].x)
784 # Add new list of cross pairs
785 ab_C.append((bc_vc, b_vl, b_vu, ba_vl, ba_vu))
786 ab_C.append((d_bc_vc, d_b_vl, d_b_vu, d_ba_vl, d_ba_vu))
787 ab_C.append((d_bc_vc, vectors[1], b_vl, a_vu, ba_vu))
788 ab_C.append((d_bc_vc, vu, b_vu, a_vl, ba_vl))
790 for j, (VL, VC, VU) in enumerate(zip(cCox, cCcx, cCux)):
791 for k, (vl, vc, vu) in enumerate(zip(VL, VC, VU)):
792 # Build aN vertices for each lower-upper C3 group in N:
793 a_vl = list(vl.x)
794 a_vu = list(vu.x)
795 a_vl[i + 1] = vut[i + 1]
796 a_vu[i + 1] = vut[i + 1]
797 a_vl = self.V[tuple(a_vl)]
798 a_vu = self.V[tuple(a_vu)]
799 # Note, build (a + vc) later for consistent yields
800 # Split the a + b edge of the initial triangulation:
801 c_vc = self.split_edge(vl.x, a_vu.x)
802 self.split_edge(vl.x, vu.x) # Equal to vc
803 # Build cN vertices for each lower-upper C3 group in N:
804 c_vc.connect(vco)
805 c_vc.connect(vc)
806 c_vc.connect(vl) # Connect c + ac operations
807 c_vc.connect(vu) # Connect c + ac operations
808 c_vc.connect(a_vl) # Connect c + ac operations
809 c_vc.connect(a_vu) # Connect c + ac operations
810 yield c_vc.x
811 c_vl = self.split_edge(vl.x, a_vl.x)
812 c_vl.connect(vco)
813 c_vc.connect(c_vl) # Connect cN group vertices
814 yield c_vl.x
815 # yield at end of loop:
816 c_vu = self.split_edge(vu.x, a_vu.x)
817 c_vu.connect(vco)
818 # Connect remaining cN group vertices
819 c_vc.connect(c_vu) # Connect cN group vertices
820 yield c_vu.x
822 a_vc = self.split_edge(a_vl.x, a_vu.x) # is (a + vc) ?
823 a_vc.connect(vco)
824 a_vc.connect(c_vc)
826 # Storage for connecting c + ac operations:
827 ab_C.append((c_vc, vl, vu, a_vl, a_vu))
829 # Update the containers
830 Cox[i + 1].append(vl)
831 Cox[i + 1].append(vc)
832 Cox[i + 1].append(vu)
833 Ccx[i + 1].append(c_vl)
834 Ccx[i + 1].append(c_vc)
835 Ccx[i + 1].append(c_vu)
836 Cux[i + 1].append(a_vl)
837 Cux[i + 1].append(a_vc)
838 Cux[i + 1].append(a_vu)
840 # Update old containers
841 Cox[j].append(c_vl) # !
842 Cox[j].append(a_vl)
843 Ccx[j].append(c_vc) # !
844 Ccx[j].append(a_vc) # !
845 Cux[j].append(c_vu) # !
846 Cux[j].append(a_vu)
848 # Yield new points
849 yield a_vc.x
851 except IndexError:
852 for vectors in ab_Cc:
853 ba_vl = list(vectors[3].x)
854 ba_vu = list(vectors[4].x)
855 ba_vl[i + 1] = vut[i + 1]
856 ba_vu[i + 1] = vut[i + 1]
857 ba_vu = self.V[tuple(ba_vu)]
858 yield ba_vu
859 d_bc_vc = self.split_edge(vectors[1].x, ba_vu.x) # o-s
860 yield ba_vu
861 d_bc_vc.connect(vectors[1]) # Connect all to centroid
862 d_bc_vc.connect(vectors[2]) # Connect all to centroid
863 d_bc_vc.connect(vectors[3]) # Connect all to centroid
864 d_bc_vc.connect(vectors[4]) # Connect all to centroid
865 yield d_bc_vc.x
866 ba_vl = self.V[tuple(ba_vl)]
867 yield ba_vl
868 d_ba_vl = self.split_edge(vectors[3].x, ba_vl.x)
869 d_ba_vu = self.split_edge(vectors[4].x, ba_vu.x)
870 d_ba_vc = self.split_edge(d_ba_vl.x, d_ba_vu.x)
871 yield d_ba_vl
872 yield d_ba_vu
873 yield d_ba_vc
874 c_vc, vl, vu, a_vl, a_vu = vectors
875 comb = [vl, vu, a_vl, a_vu,
876 ba_vl,
877 ba_vu]
878 comb_iter = itertools.combinations(comb, 2)
879 for vecs in comb_iter:
880 self.split_edge(vecs[0].x, vecs[1].x)
882 # Copy lists for iteration
883 cCox = Cox[i]
884 cCcx = Ccx[i]
885 cCux = Cux[i]
886 VL, VC, VU = cCox, cCcx, cCux
887 for k, (vl, vc, vu) in enumerate(zip(VL, VC, VU)):
888 # Build aN vertices for each lower-upper pair in N:
889 a_vu = list(vu.x)
890 a_vu[i + 1] = vut[i + 1]
892 # Connect vertices in N to corresponding vertices
893 # in aN:
894 a_vu = self.V[tuple(a_vu)]
895 yield a_vl.x
896 # Split the a + b edge of the initial triangulation:
897 c_vc = self.split_edge(vl.x, a_vu.x)
898 self.split_edge(vl.x, vu.x) # Equal to vc
899 c_vc.connect(vco)
900 c_vc.connect(vc)
901 c_vc.connect(vl) # Connect c + ac operations
902 c_vc.connect(vu) # Connect c + ac operations
903 c_vc.connect(a_vu) # Connect c + ac operations
904 yield (c_vc.x)
905 c_vu = self.split_edge(vu.x,
906 a_vu.x) # yield at end of loop
907 c_vu.connect(vco)
908 # Connect remaining cN group vertices
909 c_vc.connect(c_vu) # Connect cN group vertices
910 yield (c_vu.x)
912 # Update the containers
913 Cox[i + 1].append(vu)
914 Ccx[i + 1].append(c_vu)
915 Cux[i + 1].append(a_vu)
917 # Update old containers
918 s_ab_C.append([c_vc, vl, vu, a_vu])
920 yield a_vu.x
922 # Clean class trash
923 try:
924 del Cox
925 del Ccx
926 del Cux
927 del ab_C
928 del ab_Cc
929 except UnboundLocalError:
930 pass
932 try:
933 self.triangulated_vectors.remove((tuple(origin_c),
934 tuple(supremum_c)))
935 except ValueError:
936 # Turn this into a logging warning?
937 pass
938 # Add newly triangulated vectors:
939 for vs in sup_set:
940 self.triangulated_vectors.append((tuple(vco.x), tuple(vs.x)))
942 # Extra yield to ensure that the triangulation is completed
943 if centroid:
944 vcn_set = set()
945 c_nn_lists = []
946 for vs in sup_set:
947 # Build centroid
948 c_nn = self.vpool(vco.x, vs.x)
949 try:
950 c_nn.remove(vcn_set)
951 except KeyError:
952 pass
953 c_nn_lists.append(c_nn)
955 for c_nn in c_nn_lists:
956 try:
957 c_nn.remove(vcn_set)
958 except KeyError:
959 pass
961 for vs, c_nn in zip(sup_set, c_nn_lists):
962 # Build centroid
963 vcn = self.split_edge(vco.x, vs.x)
964 vcn_set.add(vcn)
965 try: # Shouldn't be needed?
966 c_nn.remove(vcn_set)
967 except KeyError:
968 pass
969 for vnn in c_nn:
970 vcn.connect(vnn)
971 yield vcn.x
972 else:
973 pass
975 yield vut
976 return
978 def refine_star(self, v):
979 """Refine the star domain of a vertex `v`."""
980 # Copy lists before iteration
981 vnn = copy.copy(v.nn)
982 v1nn = []
983 d_v0v1_set = set()
984 for v1 in vnn:
985 v1nn.append(copy.copy(v1.nn))
987 for v1, v1nn in zip(vnn, v1nn):
988 vnnu = v1nn.intersection(vnn)
990 d_v0v1 = self.split_edge(v.x, v1.x)
991 for o_d_v0v1 in d_v0v1_set:
992 d_v0v1.connect(o_d_v0v1)
993 d_v0v1_set.add(d_v0v1)
994 for v2 in vnnu:
995 d_v1v2 = self.split_edge(v1.x, v2.x)
996 d_v0v1.connect(d_v1v2)
997 return
999 @cache
1000 def split_edge(self, v1, v2):
1001 v1 = self.V[v1]
1002 v2 = self.V[v2]
1003 # Destroy original edge, if it exists:
1004 v1.disconnect(v2)
1005 # Compute vertex on centre of edge:
1006 try:
1007 vct = (v2.x_a - v1.x_a) / 2.0 + v1.x_a
1008 except TypeError: # Allow for decimal operations
1009 vct = (v2.x_a - v1.x_a) / decimal.Decimal(2.0) + v1.x_a
1011 vc = self.V[tuple(vct)]
1012 # Connect to original 2 vertices to the new centre vertex
1013 vc.connect(v1)
1014 vc.connect(v2)
1015 return vc
1017 def vpool(self, origin, supremum):
1018 vot = tuple(origin)
1019 vst = tuple(supremum)
1020 # Initiate vertices in case they don't exist
1021 vo = self.V[vot]
1022 vs = self.V[vst]
1024 # Remove origin - supremum disconnect
1026 # Find the lower/upper bounds of the refinement hyperrectangle
1027 bl = list(vot)
1028 bu = list(vst)
1029 for i, (voi, vsi) in enumerate(zip(vot, vst)):
1030 if bl[i] > vsi:
1031 bl[i] = vsi
1032 if bu[i] < voi:
1033 bu[i] = voi
1035 # NOTE: This is mostly done with sets/lists because we aren't sure
1036 # how well the numpy arrays will scale to thousands of
1037 # dimensions.
1038 vn_pool = set()
1039 vn_pool.update(vo.nn)
1040 vn_pool.update(vs.nn)
1041 cvn_pool = copy.copy(vn_pool)
1042 for vn in cvn_pool:
1043 for i, xi in enumerate(vn.x):
1044 if bl[i] <= xi <= bu[i]:
1045 pass
1046 else:
1047 try:
1048 vn_pool.remove(vn)
1049 except KeyError:
1050 pass # NOTE: Not all neigbouds are in initial pool
1051 return vn_pool
1053 def vf_to_vv(self, vertices, simplices):
1054 """
1055 Convert a vertex-face mesh to a vertex-vertex mesh used by this class
1057 Parameters
1058 ----------
1059 vertices : list
1060 Vertices
1061 simplices : list
1062 Simplices
1063 """
1064 if self.dim > 1:
1065 for s in simplices:
1066 edges = itertools.combinations(s, self.dim)
1067 for e in edges:
1068 self.V[tuple(vertices[e[0]])].connect(
1069 self.V[tuple(vertices[e[1]])])
1070 else:
1071 for e in simplices:
1072 self.V[tuple(vertices[e[0]])].connect(
1073 self.V[tuple(vertices[e[1]])])
1074 return
1076 def connect_vertex_non_symm(self, v_x, near=None):
1077 """
1078 Adds a vertex at coords v_x to the complex that is not symmetric to the
1079 initial triangulation and sub-triangulation.
1081 If near is specified (for example; a star domain or collections of
1082 cells known to contain v) then only those simplices containd in near
1083 will be searched, this greatly speeds up the process.
1085 If near is not specified this method will search the entire simplicial
1086 complex structure.
1088 Parameters
1089 ----------
1090 v_x : tuple
1091 Coordinates of non-symmetric vertex
1092 near : set or list
1093 List of vertices, these are points near v to check for
1094 """
1095 if near is None:
1096 star = self.V
1097 else:
1098 star = near
1099 # Create the vertex origin
1100 if tuple(v_x) in self.V.cache:
1101 if self.V[v_x] in self.V_non_symm:
1102 pass
1103 else:
1104 return
1106 self.V[v_x]
1107 found_nn = False
1108 S_rows = []
1109 for v in star:
1110 S_rows.append(v.x)
1112 S_rows = numpy.array(S_rows)
1113 A = numpy.array(S_rows) - numpy.array(v_x)
1114 # Iterate through all the possible simplices of S_rows
1115 for s_i in itertools.combinations(range(S_rows.shape[0]),
1116 r=self.dim + 1):
1117 # Check if connected, else s_i is not a simplex
1118 valid_simplex = True
1119 for i in itertools.combinations(s_i, r=2):
1120 # Every combination of vertices must be connected, we check of
1121 # the current iteration of all combinations of s_i are
1122 # connected we break the loop if it is not.
1123 if ((self.V[tuple(S_rows[i[1]])] not in
1124 self.V[tuple(S_rows[i[0]])].nn)
1125 and (self.V[tuple(S_rows[i[0]])] not in
1126 self.V[tuple(S_rows[i[1]])].nn)):
1127 valid_simplex = False
1128 break
1130 S = S_rows[tuple([s_i])]
1131 if valid_simplex:
1132 if self.deg_simplex(S, proj=None):
1133 valid_simplex = False
1135 # If s_i is a valid simplex we can test if v_x is inside si
1136 if valid_simplex:
1137 # Find the A_j0 value from the precalculated values
1138 A_j0 = A[tuple([s_i])]
1139 if self.in_simplex(S, v_x, A_j0):
1140 found_nn = True
1141 # breaks the main for loop, s_i is the target simplex:
1142 break
1144 # Connect the simplex to point
1145 if found_nn:
1146 for i in s_i:
1147 self.V[v_x].connect(self.V[tuple(S_rows[i])])
1148 # Attached the simplex to storage for all non-symmetric vertices
1149 self.V_non_symm.append(self.V[v_x])
1150 # this bool value indicates a successful connection if True:
1151 return found_nn
1153 def in_simplex(self, S, v_x, A_j0=None):
1154 """Check if a vector v_x is in simplex `S`.
1156 Parameters
1157 ----------
1158 S : array_like
1159 Array containing simplex entries of vertices as rows
1160 v_x :
1161 A candidate vertex
1162 A_j0 : array, optional,
1163 Allows for A_j0 to be pre-calculated
1165 Returns
1166 -------
1167 res : boolean
1168 True if `v_x` is in `S`
1169 """
1170 A_11 = numpy.delete(S, 0, 0) - S[0]
1172 sign_det_A_11 = numpy.sign(numpy.linalg.det(A_11))
1173 if sign_det_A_11 == 0:
1174 # NOTE: We keep the variable A_11, but we loop through A_jj
1175 # ind=
1176 # while sign_det_A_11 == 0:
1177 # A_11 = numpy.delete(S, ind, 0) - S[ind]
1178 # sign_det_A_11 = numpy.sign(numpy.linalg.det(A_11))
1180 sign_det_A_11 = -1 # TODO: Choose another det of j instead?
1181 # TODO: Unlikely to work in many cases
1183 if A_j0 is None:
1184 A_j0 = S - v_x
1186 for d in range(self.dim + 1):
1187 det_A_jj = (-1)**d * sign_det_A_11
1188 # TODO: Note that scipy might be faster to add as an optional
1189 # dependency
1190 sign_det_A_j0 = numpy.sign(numpy.linalg.det(numpy.delete(A_j0, d,
1191 0)))
1192 # TODO: Note if sign_det_A_j0 == then the point is coplanar to the
1193 # current simplex facet, so perhaps return True and attach?
1194 if det_A_jj == sign_det_A_j0:
1195 continue
1196 else:
1197 return False
1199 return True
1201 def deg_simplex(self, S, proj=None):
1202 """Test a simplex S for degeneracy (linear dependence in R^dim).
1204 Parameters
1205 ----------
1206 S : np.array
1207 Simplex with rows as vertex vectors
1208 proj : array, optional,
1209 If the projection S[1:] - S[0] is already
1210 computed it can be added as an optional argument.
1211 """
1212 # Strategy: we test all combination of faces, if any of the
1213 # determinants are zero then the vectors lie on the same face and is
1214 # therefore linearly dependent in the space of R^dim
1215 if proj is None:
1216 proj = S[1:] - S[0]
1218 # TODO: Is checking the projection of one vertex against faces of other
1219 # vertices sufficient? Or do we need to check more vertices in
1220 # dimensions higher than 2?
1221 # TODO: Literature seems to suggest using proj.T, but why is this
1222 # needed?
1223 if numpy.linalg.det(proj) == 0.0: # TODO: Repalace with tolerance?
1224 return True # Simplex is degenerate
1225 else:
1226 return False # Simplex is not degenerate