Skip to content

Transport_program

Jump to bottom
Bruno Levy edited this page Jul 25, 2022 · 24 revisions

How to program your own semi-discrete optimal transport at home

Prerequisites

Before starting, if you do not have it already, you will need to install the Graphite software (and its WarpDrive plugin, that is included in the distribution):

  • if you are under Windows, install pre-compiled package (instructions here)
  • if you are under Linux or MacOS, install Graphite from sources ([instructions here)[https://github.com/BrunoLevy/GraphiteThree/wiki#installing])

Programming with Graphite

In this tutorial, we will use the LUA programming language. LUA is a scripting language similar to Python, if you know Python already you will (nearly) find yourself at home. More information about LUA here.

Why LUA rather than Python ?

  • LUA is smaller/faster/more elegantly designed (in my own subjective opinion). A more objective measure: it has a super small footprint (in terms of source size, package size). It is directly included in Graphite (nothing to install).

  • If you really prefer Python (that is much more popular than LUA and that has many existing packages), there will be a Python version of this tutorial soon (Graphite also supports Python through the "gompy" plugin).

Now you can take a quick look at the Scripting Graphite tutorial to see how the editor works (see in particular automatic completion and help bubbles that can save you some time).

Step 1: Laguerre diagram

We shall now see how to create a Laguerre diagram in Graphite and how to visualize it.

Create a square

The first thing we need to do is defining the domain. We create a unit square, using the create_square() method of the Shapes interface. Then we need to triangulate the domain (this will simply split the square into two triangles in our case).

scene_graph.clear()
Omega = scene_graph.create_object('OGF::MeshGrob')
Omega.rename('Omega')
Omega.I.Shapes.create_square()
Omega.I.Surface.triangulate()
Omega.visible = false

Sample it with random points

Now we create a point set. In this example, we will use a set of N=100 points picked randomly in our unit square:

N = 100 -- Number of points
Omega.I.Points.sample_surface(
   {nb_points=N,Lloyd_iter=0,Newton_iter=0}
)
points = scene_graph.objects.points

Compute Laguerre diagram

And finally, we can compute the Laguerre diagram (restricted to the domain). We first create a new mesh to store it. Laguerre diagrams are defined from a pointset and a vector of "weights". Here we use (for now) a vector of weights with all weights set to 0 (default value set by NL.create_vector()).

RVD = scene_graph.create_object('OGF::MeshGrob','RVD')
weight   = NL.create_vector(N)

Then we can compute the Laguerre diagram. The WarpDrive plugin adds a special Transport interface solely designed for scripting (it does not appear in the menus). It has a function to compute the Laguerre diagram:

points.I.Transport.compute_Laguerre_diagram(
   Omega, weight, RVD, 'EULER_2D'
)

Display the cells in different colors

And finally, we can set the graphic attributes of the computed Laguerre diagram to display the cells in different colors:

RVD.shader.painting='ATTRIBUTE'
RVD.shader.attribute='facets.chart'
RVD.shader.colormap = 'plasma;false;732;false;false;;'
RVD.shader.autorange()

The complete program is available here. You can also load it from the Examples... menu of the text editor (Examples.../WarpDrive/advanced/Transport_2d_01.lua).

Change the weight of a cell in a Laguerre diagram

Let us now see what happens when we change the weight of one cell in the Laguerre diagram. For that, we write a function compute(), at the beginning of our file as follows:

function compute()
   weight[0] = weight[0] + 0.01
   OT.compute_Laguerre_diagram(Omega, weight, RVD, 'EULER_2D')   
   RVD.shader.autorange()
   RVD.update()
end

Now we want to create a button that will call our function compute() each time we push it. At the end of our file, we add:

OT_dialog = {} 
OT_dialog.visible = true
OT_dialog.name = 'Transport' 
OT_dialog.x = 100
OT_dialog.y = 400
OT_dialog.w = 150
OT_dialog.h = 200
OT_dialog.width = 400

function OT_dialog.draw_window()
   if imgui.Button('Compute',-1,-1) then
      compute()
   end
end

graphite_main_window.add_module(OT_dialog)

This declares a new Module to the graphic user interface of Graphite. See the Scripting Graphite tutorial for more details. Each time the Compute button is pressed, our function compute() is called. It increases a bit the weight of one of the points and recomputes the Laguerre diagram, so that you can see the impact of modifying the weight of a single point.

The complete program is available here. You can also load it from the Examples... menu of the text editor (Examples.../WarpDrive/advanced/Transport_2d_02.lua).

Run the program (<F5>), then push the compute button multiple times, you will see one of the cells growing bigger and bigger. It will eventually nibble the entire domain (all the other cells become eventually empty, yes a Laguerre diagram can have empty cells).

Translating a Laguerre diagram by just changing the weight

Now you may think about a Laguerre diagram as a Voronoi diagram plus tuning buttons (the weights). By changing the tuning buttons, one may increase or decrease the size of the associated Laguerre cells. In fact there is more: did you know that you can translate a Laguerre diagram by an arbitrary vector just by changing the weights ?

Let us take a close look at the common boundary between two cells i and j, associated with points xi and xj. A point x that is on this common boundary satisfies the following equation:

$$d2(x,xi) - wi = d2(x,xj) - wj$$

where wi and wj are the weights associated with points xi and xj and where d2(.,.) denotes the squared distance between two points.

Side note: One can notice that this equation only involves linear terms in the coordinates of x (the squared lengh of x appears on both sides of the equation and cancel out), hence this common boundary is a straight line.

Suppose that wi = wj = 0, then what you have is the Voronoi diagram of xi and xj. Imagine now that you want to translate the whole diagram by an arbitrary 2D vector T. Is is possible to do so just by changing wi and wj ?

Try to derive it on your own before reading further

In other words, how can we set wi and wj in such a way that the set of points x that satisfy:

$$d2(x,xi+T) - wi = d2(x,xj+T) - wj$$

corresponds to the same straight line as the edge between the Voroonoi cells of xi and xj ?

It can be easily checked that setting wi = -2(T.xi) and wj = -2(T.xj) works (where . denotes the dot product).

OK, let us now program it !

First thing we need to do is to group a little bit our points in the lower-left corner of the square, so that we will still see them when we will translate them. Right after we create the point set, we declare a mesh editor (E = points.I.Editor), find the attribute that stores the point coordinates (coords = E.find_attribute('vertices.point')) and modify the coordinates:

  E = points.I.Editor
  coords = E.find_attribute('vertices.point')
  for i=0,N-1 do
     coords[3*i]   = coords[3*i]/2.0
     coords[3*i+1] = coords[3*i+1]/2.0
  end

Then, we will keep the same graphic interface (the big Compute button), and change the compute() function as follows:

  Tx = 0.0
  Ty = 0.0
  function compute()
     Tx = Tx + 0.05
     Ty = Ty + 0.05
     for i=0,N-1 do
        weight[i] = - 2.0 * Tx * coords[3*i]
                    - 2.0 * Ty * coords[3*i+1]   
     end
     OT.compute_Laguerre_diagram(Omega, weight, RVD, 'EULER_2D')   
     RVD.shader.autorange()
     RVD.update()
   end

The complete program is available here. You can also load it from the Examples... menu of the text editor (Examples.../WarpDrive/advanced/Transport_2d_03.lua).

Run the program (<F5>), then push the compute button multiple times, you will see the whole diagram moving towards the upper-right corner.

What we have learnt so far:

  • how to compute a Laguerre diagram in Graphite
  • by changing the weights of a Laguerre diagram, one can change the area of the cells
  • by changing the weights of a Laguerre diagram, one can translate it at a arbitrary position

Semi-discrete optimal transport without step-length control

Semi-discrete Optimal Transport gives a way of computing the weights in such a way that all the cells have prescribed areas. As before, we will start from a random distribution of points, and from their Voronoi diagram (all weights set to zero). Then we will iteratively update the weights until all the Laguerre cells have the prescribed areas (it will be the same, nu_i=1/N', for all cell iin the present example). For now, we want to do one updating step each time theCompute` button is pressed.

We just need to change the compute() function of our previous program, to make it compute one Newton step. We know (detailed derivation here) that this means solving a discrete Poisson system. Warpdrive has a function to assemble the Poisson system for you (we will see lated how to do that on our own), as follows:

  function compute()
     local L = NL.create_matrix(N,N) -- P1 Laplacian of Laguerre cells
     local b = NL.create_vector(N)   -- right-hand side
     OT.compute_Laguerre_cells_P1_Laplacian(
        Omega, weight, L, b, 'EULER_2D'
     )
     ...
  end

This creates a sparse matrix L with the discrete Laplacian corresponding to the current Laguerre diagram, and initializes the vector b with the areas of the Laguerre cells. See also the scripting tutorial about linear algebra in Lua.

Then we need to initialize the right-hand side b of our Poisson system. Each component of b is supposed to be equal to the target area (nu_i=1/N in our case) minus the current area of the Laguerre cell. Hence we update b as follows (still in the same compute() function):

   ...
   for i=0,N-1 do
     local nu_i = 1.0/N  -- desired area for Laguerre cell i
     b[i] = nu_i - b[i]  -- rhs = desired area - actual area
   end
   ...

Then we solve the linear system:

   ...
   local p = NL.create_vector(N)   -- Newton step vector
   L.solve_symmetric(b,p)          -- solve for p in Lp=b
   ...

This gives us a Newton step vector p. We still do not know how much along p we should descend, for now we descend by a fixed parameter alpha=1/8 (more on this later):

   ...
   local alpha = 1.0/8.0           -- Steplength (constant for now)
   NL.blas.axpy(alpha, p, weight)  -- weight = weight + alpha * p
   ...

Putting everything together, our new compute() function looks like that:

function compute()
   -- compute L(Laplacian) and b(init. with Laguerre cells areas)
   local L = NL.create_matrix(N,N) -- P1 Laplacian of Laguerre cells
   local b = NL.create_vector(N)   -- right-hand side
   OT.compute_Laguerre_cells_P1_Laplacian(
      Omega, weight, L, b, 'EULER_2D'
   )
   
   for i=0,N-1 do
     local nu_i = 1.0/N  -- desired area for Laguerre cell i
     b[i] = nu_i - b[i]  -- rhs = desired area - actual area
   end

   local p = NL.create_vector(N)   -- Newton step vector
   L.solve_symmetric(b,p)          -- solve for p in Lp=b
   local alpha = 1.0/8.0           -- Steplength (constant for now)
   NL.blas.axpy(alpha, p, weight)  -- weight = weight + alpha * p

   OT.compute_Laguerre_diagram(Omega, weight, RVD, 'EULER_2D')   
   RVD.shader.autorange()
   RVD.update()
end

The complete program is available here. You can also load it from the Examples... menu of the text editor (Examples.../WarpDrive/advanced/Transport_2d_04.lua).

  • Run the program (<F5>), then push the compute button multiple times, you will see that it converges to a Laguerre diagram with the same area for all the cells.
  • Set shrink_points to true, run the program then push the compute button multiple times. This will show a more extreme change of cells area. As you can see, in a Laguerre diagram, a point does not necessarily belong to its cell (unlike Voronoi diagrams).

Clone this wiki locally