Periodicity

Given $G$ a group of isometries, the periodicity parameter takes the cost $c(x,y)=\underset{g\in G}{min}|x-g(y){|}^2$: the transport takes place on the quotient space ${\mathbb{R}}^d/G$. It is accepted by plan, w2, ot_solve, OtProblem and PowerDiagram.

For instance, to compute a plan inside a Klein bottle:

flip = np.diag( [ 1, -1 ] )
p = sdot.plan( sdot.SumOfDiracs( positions ), sdot.Box(),
    periodicity = [
        ( flip, [ 1, 1 ] ),
        [ 0, 1 ]
    ]
)

The group

periodicity lists the generators of the group: their inverses and compositions are added. Each one is an isometry $x\mapsto Mx+V$, $M$ orthogonal:

a generator	the map
a vector V	the translation $x\mapsto x+V$
a pair ( M, V )	$x\mapsto Mx+V$: a rotation, a symmetry, a glide, a screw...
a $(d+1{)}^2$ homogeneous matrix	the same, as one matrix

f, g = sdot.SumOfDiracs( positions ), sdot.Box()
quarter = np.array( [ [ 0, -1 ], [ 1, 0 ] ] )

sdot.plan( f, g, periodicity = [ [ 1, 0 ], [ 0, 1 ] ] )                 # a torus
sdot.plan( f, g, periodicity = [ ( quarter, [ 0, 0 ] ) ] )              # a cone: a quarter turn about a corner
sdot.plan( f, g, periodicity = [ True, False ] )                        # a cylinder along x
sdot.plan( f, g, periodicity = [ [ 1, 1 ] ] )                           # a cylinder along any direction

Translations that are not independent are reduced to a basis of their lattice: [ [ 1, 0 ], [ 2, 0 ] ] is the period $1$, and [ [ 1, 0 ], [ 1.1, 0 ] ] the period $0.1$. A group that is not discrete is refused: [ [ 1, 0 ], [ √2, 0 ] ], a turn by one radian, a fivefold turn with a translation. So is a generator that is the identity.

The target is the content of a fundamental domain: its images by the group must not overlap. This is checked for every group, the first time a diagram is built with it: an image of the support that holds a ball inside the support raises an error naming the element of the group. Images that only touch the support are fine. The sides of the support that the group does not glue stay walls.

Two shortcuts read the support of the target:

• True, or one boolean per axis ( [ True, False ]: a cylinder along $x$ ): one translation per side of an axis-aligned box;
• target.domain: each facet of the support is glued onto the opposite parallel facet, by the translation that sends one onto the other. On a parallelogram, that translation runs along a side, not across it. The two facets must be translates of each other, and translations that are not independent are reduced to a basis of their lattice: the three pairs of sides of a hexagon give a lattice of rank 2.

The fundamental domain of a target

The group alone does not say where the fundamental domain is. It is the support of the target, as the distribution declares it. The values of the density are never read to find it. A density that is not constant is repeated as it is, and its images are summed:

target	its fundamental domain
Box, Polytope, ConvexHull	the polytope itself
Image, SplineGrid	the block of its pixels, placed by frame and origin ( a parallelogram for an oblique frame: glue it with target.domain )
Mesh	the box of its nodes, without the margin of its support
SumOfGaussians	support_box = ( lo, hi ), which is then required: the box of the tails depends on the gaussians and cannot be a period

With True, the booleans or target.domain, the periods are read from this domain. With explicit generators, the domain is where the diracs are folded and where the walls are found, and it must lie in a fundamental domain of the group ( checked as above ).

Visualization

A cell is computed whole, around its dirac as given. To check a quotient by eye, draw the cells, then their copies by the group around the identity, fainter: the copies must glue without a seam. The domain is dashed.

p = sdot.plan( sdot.SumOfDiracs( positions ), sdot.Box(), periodicity = True )

viz = sdot.Visualizer()
p.add_to_viz( viz )    # the cells, then their copies by the first ring of the group, fainter

# more rings, more visible, and counted in the zoom of the page
p.add_to_viz( viz, nb_periodicity_rings = 2, opacity_periodicity_rings = 0.5, fit_periodicity_rings = True )
viz.write_html( "torus.html" )

The same by hand, with the pieces it is made of:

pd = p.power_diagram()
for k, a in enumerate( pd.periodicity_rings( 1 ) ):   # the elements of the first ring, as affine maps [ M | V ], the identity first
    pd.cells.moved( a ).add_to_viz( viz, opacity = 1.0 if k == 0 else 0.25, frames = k == 0 )   # ( the copies out of the zoom )

periodicity_rings( k ) gives the elements of the group up to the ring k: ring 1 holds the elements whose image of the domain touches the domain ( the 8 neighbours of a square on a torus ), ring 2 those touching ring 1, and so on. Cell.moved moves a cell by one of them, or each cell by its own.

Cells are where their diracs are

The diracs can lie anywhere. Below, three of them are given outside the domain, one period away or more, and one has a large mass. With ring = 0 ( nb_periodicity_rings = 0 ) the cells are shown alone: those of the three diracs lie outside the domain, away from the others, and the big one crosses a side. With ring = 1 their copies fill the gaps, and the quotient appears.

Open in a separate page ↗

The barycenters and the moments are those of these cells, next to the diracs as they were given. A Lloyd step $x_i\leftarrow b_i$ is then a short move, never a jump across the domain, and the gradient $2m_i(x_i-b_i)$ of the cost ( cost_and_position_grad ) stays continuous when a dirac crosses a side, whatever the group does to it there.

One dirac moving to the right on a torus ( the slider: its abscissa x ). Its cell follows it whole, out of the domain, while its copy by the translation $(-1,0)$ comes in on the left side. Nothing jumps when $x$ passes $1$:

Open in a separate page ↗

Everything is differentiable as without periodicity: with respect to the positions, the masses of the diracs and the parameters of the target.

A gallery of quotients

The rest of the page runs Lloyd's iterations on several quotients. Each dirac moves to the barycenter of its cell until the cells are as regular as the space allows, so the result shows the geometry of the quotient itself. The slider at the top of each viewer goes through the iterations.

pb = sdot.OtProblem( sdot.SumOfDiracs( x ), target, periodicity = generators )
for it in range( 50 ):
    pb.source = sdot.SumOfDiracs( x )
    x = np.asarray( pb.solve().barycenters )      # next to the given diracs, whatever the side they crossed

A torus

periodicity = True: one translation per side of the box.

Open in a separate page ↗

A cone

periodicity = [ ( quarter, [ 0, 0 ] ) ], with quarter the rotation by 90° about the origin. The group is finite ( four elements ) and has no translation at all. The quotient is a cone whose apex is the corner of the square. The target is an Image whose density peaks at the apex, and the cells there wrap around it, each one spread over the four quarters.

Open in a separate page ↗

A wallpaper group

periodicity = [ ( quarter, [ 0, 0 ] ), [ 1, 0 ] ]: quarter turns and translations generate the wallpaper group p4, a pattern that repeats with a fourfold rotational symmetry. Its fundamental domain is a square of side $1/2$, and its copies around it are turned by multiples of 90°.

Open in a separate page ↗

A hexagonal lattice

periodicity = [ [ 1, 0 ], [ 1/2, √3/2 ] ]: translations again, but the second one is not orthogonal to the first. A box of height $\sqrt{3}/2$ is glued to itself with a shift of half its width, like bricks. This is the lattice of the hexagonal packing, and 16 diracs settle into regular hexagons, which a square torus could not hold.

Open in a separate page ↗

A hexagon glued by target.domain

The same lattice, from another fundamental domain: a regular hexagon whose opposite sides are glued, read from the support by periodicity = target.domain. The cells end as regular hexagons themselves.

t = np.arange( 6 ) * np.pi / 3
target = sdot.ConvexHull( np.stack( [ np.cos( t ), np.sin( t ) ], -1 ) )
pb = sdot.OtProblem( sdot.SumOfDiracs( x ), target, periodicity = target.domain )

Open in a separate page ↗

A cylinder along any direction

periodicity = [ True, False ] gives a cylinder along $x$, but the booleans only read the sides of an axis-aligned box. Any other cylinder is a translation: here the support is a ConvexHull, a $2\times 1$ rectangle turned by 30°, and periodicity = [ 2 u ] glues it along its long side, $u=(\cos 30°,\sin 30°)$. Its two other sides stay walls.

Open in a separate page ↗

In 3D

Nothing above is specific to the plane. A screw, for example, moves up one step while turning by a quarter about the vertical axis:

quarter_z = np.array( [ [ 0, -1, 0 ], [ 1, 0, 0 ], [ 0, 0, 1 ] ] )
p = sdot.plan( sdot.SumOfDiracs( x ), sdot.Box(), periodicity = [ ( quarter_z, [ 1, 0, 1 ] ), [ 1, 0, 0 ], [ 0, 1, 0 ] ] )

Under the hood

There is one cell per dirac. It is built around its dirac, in the plane, and is not cut at the sides of the support. A cell cannot leave the Voronoi cell of its dirac among its own images, since they all have the same weight.

Internally, the diracs are first folded onto the support by the group ( sdot.wrap does the same thing by hand ). This changes nothing in what is computed, since a dirac and its images are the same point of the quotient, and the cells come back next to the diracs as given ( an isometry ). What it buys is compactness: the differences between the diracs stay within a period, so the images that can cut a cell are a short, fixed list for all the cells ( 9 translations on a 2D torus, however far apart the given diracs are ), and each cell starts from one box, the support widened by half a period.

The planes come from the diracs and from their images by the elements of the group that come close to the support. The tree of the diracs is traversed once per element, with its boxes moved. The density is integrated as the sum of its images, each cutting the cell into its own pieces. In 2D, a torus traverses the tree 9 times instead of once, but the pruning drops most of the traversals right away.

On a GPU

The diagrams ( PowerDiagram: measures, moments, their derivatives ) work on a CUDA device for any group. The solve on a CUDA device handles a group of translations ( a torus, a cylinder, a sheared lattice ) with a constant density ( a Box, a Polytope, an Image of equal values ), gaussians or a Mesh, in 2D and 3D. A group with rotations or symmetries, or an Image that is not a constant, is solved on the CPU ( LOOM_DEVICE=cpu ).

Gaussians and meshes

A SumOfGaussians is wrapped: each gaussian is summed over its images, its tails included, and its support_box is the fundamental domain ( see Gaussians ). A Mesh repeats the box of its nodes, without the margin of its support ( see Mesh ). Both copy themselves by the elements of the group that come near the cells, instead of being moved image by image in the kernels as an Image or a Polytope are.

What is not supported yet

• 1D ( the circle is another algorithm ), the "limits" Newton step ( "trials" is used ), the width continuation of an image ( the solve starts from the sharp image ).

See the roadmap. The figures are written by docs/guide/ground-metrics/periodicity_figures.py.