Brief tutorial for the PyAPD library¶

Installation via pip¶

This should be as easy as running the following command.

Note that PyAPD relies on

• pykeops -- a kernel operations library which ensures we are able to always use near 100% of GPU power at our disposal.
• torch -- the machine learning library that is GPU friendly, we use it for defining the model, storing data and linking between torch tensors and PyKeOps Lazy Tensors
• pytorch-minimize -- for optimisation routines, which are designed to mimic scipy while remaining on the GPU.
• matplotlib -- for plotting
• numpy -- for some of the data manipulation, (using torch is preferred though and going forward we may get rid of numpy altogether).
• scipy -- pytorch-minimize calls scipy and we manually override one of the pytorch-minimize modules to ensure that error tolerance for gradients can be vectorial, that is different for each component of the gradient.

In the near future we will most likely add pyvista to the list of dependencies to have internal 3D plotting capabilities.

These dependencies are installed automatically when installing PyAPD.

Example usage with brief description¶

We begin by importing the package and defining an instance of an apd_system, which is the central object of the library.

By default, we assemble a $D=2$ dimensional system with $N=10$ grains occupying a square domain $[0,1]^D$.

From the GPU perspective, other important parameters are

• apd1.dt defining the data type, either torch.float32 (single precision) or torch.float64 (double precision). Single precision will be generally much faster as most GPUs are geared towards single precision arithmetic. Unfortunately for our application, when $N$ is large, this can lead to precision loss, we advise caution.
• apd1.device which by default gets set to "cuda" if torch.cuda.is_available() else "cpu", meaning that we work with GPUs when they are available, otherwise we rever to a cpu.

While these parameters can be adjusted for an already defined apd_system, there may be some bugs in the implementation, so a switch to a new instance of apd_system is advised if one wants to switch to a different device or a different data type.

We also have several parameters that are important for fitting optimal APDs and generating synthetic microstructure.

• apd1.error_tolerance is the relative error tolerance we are happy with with respect to volumes of grains, by default set to $\varepsilon = 0.01$. If we want $i$th grain to have volume $v_i$, this means that we are happy to accept that a grain, $L_i$, has volume satisfying either (option error_wrt_each_grain = True in solvers (to be discused below)):

or (option error_wrt_each_grain = False): $$ \frac{\big||L_i|- v_i\big|}{\min_j v_j} \leq \varepsilon. $$

In terms of seed points, they are stored at apd1.X and their LazyTensor counterpart is stored at apd1.x. If we do not pass them explicitly, they are chosen on random with an internal call to apd1.set_X which also relies on the parameter:

• apd1.radius_of_exclusion is used when generating random seed points and defines a relative radius around an already sampled seed ensuring no other seed point is generated inside it. The radius shrinks with the number of points and is given by

where $R_{{\rm roe}}$ is apd1.radius_of_exclusion. This feature was introduced in relation to precision loss that we observed for very large, $N = 10,000$, especially when $D=2$. This feature can be disabled by setting apd1.radius_of_exclusion = 0.

In terms of anisotropy matrices, they are stored at apd1.As and their LazyTensor counterpart at apd1.a. If we do not pass them explicitly, they are chosen on random with an internal call to apd1.set_As which relies on two parameters:

• apd1.ani_thres determines the allowed level of anisotropy each grain can exhibit and is used when generating synthetic microstructure. If ``apd1.ani_thres = 0, the anisotropy is switched off and have ${\mathbf A}_i = {\rm Id}$ and thus we generate a power diagram. If on the other handapd1.ani_thres = 1`, then any level of anistropy is allowed, i.e. in principle one of the axes can be arbitrarily small with respect to others.
• apd1.det_constraint determines whether we constrain the anisotropy matrices to satisfy ${\rm det} \mathbf{A}_i = 1$. It is set to True be default and we recommend keeping it this way.

We further note that

apd1.seed determines whether we work with a specific random seed for sampling seed points and anisotropy matrices, which can be useful for reproducibility purposes. By default we have it set to -1 which is meant to mean that we do not fix a random seed.

In terms of weights, they are stored at apd1.W and their LazyTensor counterpart at apd1.w. If apd1.heuristic_W = False then they are initial set to $W = 0$, whereas if apd1.heuristic_W = True, they are set to the heuristic guess discussed in

• Teferra, Kirubel, and David J. Rowenhorst. "Direct parameter estimation for generalised balanced power diagrams." Philosophical Magazine Letters 98.2 (2018): 79-87.

Target masses are stored at apd1.target_masses and if they are not passed explicitly, they are by default set to uniform masses, which is $1/N$. Note that if we pass explicitly target masses, they should always add up to 1 and will be internally weighted by apd1.box_size which corresponds to the size of the domain.

The final parameter to discuss is apd1.pixel_size_prefactor, which regulates how many pixels we should create. We work with rectilinear pixels/voxels with the length of each pixel/voxel in the $i$th dimension given by $\Omega_i / M_i$, where $\Omega_i$ is the length of the domain in the $i$th direction.

We use a heuristic formula which, first of all, assumes that $M = M_i$ for all $i$. Secondly, we want each pixel to less than $\varepsilon \min_j v_j$, which can be written as $$ |\Omega|\left(\frac{1}{M^D}\right) \leq \varepsilon \min_j v_j. $$ Rearranging, we get $$ M \geq \left(\frac{|\Omega|}{\varepsilon \min_j v_j}\right) $$ which we translate to $$ M = C\left(\frac{|\Omega|}{\varepsilon \min_j v_j}\right) $$ and $C$ is precisely apd1.pixel_size_prefactor. We recommend keeping at at least $C=1$ and increase if we struggle to hit desired tolerance. Also, setting $C$ at high value is often recommended for plotting purposes,

By default, we set apd1.optimality = False to indicate that we have not optimised the APD yet

A specific graphical insight into the model comes from trying to plot seed points $X_i$ with ellipses associated with $\mathbf A_i$, which shows the preferred anisotropy

Pixel method -- computing an APD, finding optimal APDs, Lloyds algorithm¶

We are now in a position to run through the capabilities of the library.

With apd1.pixel_params specified, we can proceed to assembling the pixels, which will be stored at apd1.Y and their LazyTensor counterparts at apd1.y.

In preparation for future updates with regards to adaptable pixel sizes, we also store size of each pixel separate at apd1.PS

As we can see, we have $M^D$ pixels

For a fixed triplet of seed points $\{X\}_{i=1}^N$, anisotropy matrices $\{\mathbf A_i\}_{i=1}^N$ and weights $\{w_i\}_{i=1}^N$, we can assemble the APD as follows.

This is an assignment vector assign each pixel to a grain:

With APD computed, we can plot it as follows.

We can further check whether the current APD is optimal with respect to target masses:

and find the set of optimal weights by calling:

Let us plot it.

Since this is a random run and there are not that many grains, it is likely that some grains will be disjoint. To correct for it, we can run a generalised Lloyds algorithm as follows.

If the resulting APD still has disjoint grains, just run a few more iterations of Lloyds. In any case, we see that this is an interesting approximate space filling algorithm, since now the ellipses are nicely spread out.

One further comment is that for plotting we may wish to increase the pixel parameter $M$, which can be done quickly as follows.

What's worth noting is that this only needs to happen at the plotting stage -- for finding optimal set of weights, we are typically alright with much lower values of this parameter.

A big example to finish off¶

The method is rather fast, even on the T4 GPU from Google Colab, which we can demonstrate quickly as follows.

It is fast in 3D too.