Examples providing your own mesh and retriving the basis functions?

[cottrell]

Are there any examples for supplying your own mesh (for example from scipy Delaunay) and extracting basis functions of some degree?

For example you might evaluate the basis functions at a set of points for fitting against data as a kind of smoothing spline.

[kinnala]

There is an example on defining the mesh points and triangles in the README. It is not specific to scipy Delaunay though. I think there are several examples on interpolating the finite element solution at specific points. One of them is here: https://github.com/kinnala/scikit-fem/blob/master/docs/examples/ex12.py#L39

[kinnala]

I think scipy Delaunay is something like this:

points = np.array([[0, 0], [0, 1.1], [1, 0], [1, 1]])
from scipy.spatial import Delaunay
tri = Delaunay(points)
from skfem import MeshTri
mesh = MeshTri(points.T, tri.simplices.T)

[cottrell]

@kinnala Thanks. Yes that works and I am up to the point of trying to understand why this is blowing up memory under the hood.

In [17]: import numpy as np 
    ...: import skfem as fem 
    ...: from scipy.spatial import Delaunay  
    ...: X_knots = np.random.randn(200, 2) 
    ...: X = np.random.randn(32_000, 2) 
    ...:  
    ...: tri = Delaunay(X_knots) 
    ...: mesh = fem.MeshTri(tri.points.T, tri.simplices.T) 
    ...: Vh = fem.Basis(mesh, fem.ElementTriP2())                                                                                                                                                                  

In [18]: Vh.probes(X.T)                                                                                                                                                                                            
---------------------------------------------------------------------------
MemoryError                               Traceback (most recent call last)
<ipython-input-18-ced6a0ac1a22> in <module>
----> 1 Vh.probes(X.T)

~/anaconda3/envs/38/lib/python3.8/site-packages/skfem/assembly/basis/cell_basis.py in probes(self, x)
    144         """
    145 
--> 146         cells = self.mesh.element_finder(mapping=self.mapping)(*x)
    147         pts = self.mapping.invF(x[:, :, np.newaxis], tind=cells)
    148         phis = np.array(

~/anaconda3/envs/38/lib/python3.8/site-packages/skfem/mesh/mesh.py in finder(x, y)
   1088         def finder(x, y):
   1089             ix = tree.query(np.array([x, y]).T, 5)[1].flatten()
-> 1090             X = mapping.invF(np.array([x, y])[:, None], ix)
   1091             inside = (
   1092                 (X[0] >= 0) *

~/anaconda3/envs/38/lib/python3.8/site-packages/skfem/mapping/mapping_affine.py in invF(self, x, tind)
    138             invA, b = self.invA[:, :, tind], self.b[:, tind]
    139 
--> 140         y = (x.T - b.T).T
    141 
    142         return np.einsum('ijk,jkl->ikl', invA, y)

MemoryError: Unable to allocate 76.3 GiB for an array with shape (32000, 160000, 2) and data type float64

In [19]: %debug                                                                                                                                                                                                    
> /home/cottrell/anaconda3/envs/38/lib/python3.8/site-packages/skfem/mapping/mapping_affine.py(140)invF()
    138             invA, b = self.invA[:, :, tind], self.b[:, tind]
    139 
--> 140         y = (x.T - b.T).T
    141 
    142         return np.einsum('ijk,jkl->ikl', invA, y)

ipdb> x.shape                                                                                                                                                                                                      
(2, 1, 32000)
ipdb> bb = b                                                                                                                                                                                                       
ipdb> bb.shape                                                                                                                                                                                                     
(2, 160000)

I think it should just be evaluating the basis functions and returning something that is sparse 30k x 800ish which is fine. So I obviously must not yet underestand what I am using :)

[kinnala]

I don't use these features very often but I think CellBasis.probes creates some sort of interpolation operator instead of just interpolating at specific points. Maybe CellBasis.interpolator is more useful for your use case?

[kinnala]

Alright, seems that interpolator uses probes nowadays so no help there. So you want to interpolate at a large number of points at the same time? Could look into further optimizing that specific case unless you want to use a for-loop. This is not a very natural thing to do in the context of finite elements because you must find the correct element for each point separately which will be a nontrivial process. That part is not very well optimized at the moment, it uses scipy KDTree to find candidate elements and then tests those via inverse reference mapping.

[gdmcbain]

This is not a very natural thing to do in the context of finite elements because you must find the correct element for each point separately which will be a nontrivial process.

This is true, but is it not so much about finite elements as about unstructured meshes? Wouldn't it equally afflict finite volumes? The blow-up is in element_finder which only refers to the mesh, not the basis.

While, I agree, it's not a very natural thing to do in finite elements, it is a very natural thing to want to do in applications that otherwise have recourse to finite elements; three examples:

• interpolating between meshes
• particle-tracking
• method of characteristics method of characteristics #600

I believe FreeFEM has good implementations of the first and last of these, probably using a common facility for interpolation. It is C++ rather than Python and I've refrained from looking at the implementation anyway as it's GPL and so couldn't be copied here.

One thing that might be exploited in each of the three examples is that one has a good initial guess for the cell at second and subsequent steps; e.g. the particle-tracking and finding the foot of a characteristic requires time-stepping and a step us likely to be within a cell or to a neighbour through a particular face. Similarly for interpolating between meshes if the new is deformed from the old. Perhaps this is more natural in a finite volume mesh data structure in which the facets of cells are known, but we do already have that: t2f.

Another thing might be that the 32 000 probes are embarrassingly parallel. (I wouldn't think that that would ease the memory blow-up though, depending on your architecture.)

Is your application like one of the above three?

[gdmcbain]

Of course if you do have a scipy Delaunay object you can use find_simplex. It would be interesting to see how that compares for speed and memory.

[kinnala]

Interpolating between meshes could be further abstracted using something similar to MortarBasis and a library such as libsupermesh: https://zenodo.org/record/1316942#.YOk2onX7Td4

[cottrell]

@kinnala Ah ok, a for loop is fine for now to try something out ... If I understand correctly there is some interaction effect going on.

It's been over a decade I've done anything with FE solveres but on trying to remember this stuff it seems interesting that the computational flows are quite different for the data-assimilation problem vs galerkin DE solving. I might be way off but am thinking of the variational log-likelihood setting as "data-loss" + "smoothing-loss" + boundary conditions as similar to non-linear forced poisson eq. I suspect the weird flow comes in because the data loss requires evaluation at (potentially) off-mesh points but only once (I am considering a fixed mesh in this context).

Thanks for the help!

[kinnala]

Of course if you do have a scipy Delaunay object you can use find_simplex. It would be interesting to see how that compares for speed and memory.

One could try to hack this method so that the user can pass another query function, e.g., something based on find_simplex: https://github.com/kinnala/scikit-fem/blob/master/skfem/mesh/mesh_tri_1.py#L323

Edit: Maybe it's simpler to inherit MeshTri1 and then override element_finder.

Edit 2: Reading find_simplex source code, it utilizes some property of Delaunay triangulation and traverses the Voronoi graph from facet to facet in order to find the correct element.

[kinnala]

Actually, to me this part of scipy.spatial.Delaunay does not seem to use any Delaunay-specific property: https://github.com/scipy/scipy/blob/master/scipy/spatial/qhull.pyx#L1323

So I wonder if that could be used for any mesh if we just start from an arbitrary edge.

Edit: Shame that there is not any wrapper to call it directly. I think calling cdef requires creating a wrapper and recompiling scipy.

[cottrell]

As I'm reading through things, the other issue I'm seeing is that I'm not sure yet how to create the Laplacian operator directly as I think the problems being dealt with here are typically linear DE and cast to the weak form via integration by parts (this stuff is foggy memory to me so not sure I'm making sense).

The data assimilation problem is typically something minimization like this (like a smoothing spline)
$$
\text{argmin}{u} \left(\sum{i=1}^m L(u(x_i), y_i) + \lambda \int | \Delta u(x) |^2\text{d}x\right)
$$

Even though it is is non-linear, if the grid is fixed I think you can just evaluate the basis functions over the $x_i$ once, and also you can form the Laplacian operator once and they you have a fairly simple for for the gradient of the loss something like

$$
\text{argmin}{A} \left(\sum{i=1}^m L(B_i A, y_i) + \lambda \sum_i | \sum_{j} D_{ij} A_j |^2\right)
$$

Unless I'm missing some trick I think this is a different kind of operational flow using the same basis and operator discretization tools

[kinnala]

I have no idea what that is, but in finite elements a discrete Laplacian is typically the matrix corresponding to -D:

@BilinearForm
def laplacian(u, v, w):
    from skfem.helpers import grad, dot
    return dot(grad(u), grad(v))

D = laplacian.assemble(basis)  # create a Basis object and pass it here

[kinnala]

What is L above?

[kinnala]

And yes, scikit-fem does not automagically solve nonlinear equations (like, e.g., FEniCS does) and the linearization for resolving any nonlinearities must be implemented by the user.

[cottrell]

I'm not clear on the example. Is such poor quality representative of the untended application ? I didn't know what "data assimilation " was so looked at Wikipedia. Actually i read https://fr.wikipedia.org/wiki/Assimilation_de_donn%C3%A9es rather than the English page as it has an example and there the mesh is as regular as can be. But I assume i'm misunderstanding the example and the application. Could we get some clarification on that? No point chasing something that's not going to be used; e.g. efficiency on pathological meshes when all practical meshes are going to be fine.

On the other hand, the example is Delaunay. If this holds in the application, why not exploit the Delaunay simplex finder? But is it Delaunay in the intended application ? Clarification please.

And regarding Delaunay, I just used that as a quick example. I think for smoothing applications where one expects the scale of the problem to somehow align with intensity (number) of observations it could make sense to use some delauny mesh.

I think the most important thing to keep in mind for ML/smoothing/data assimilation is that the one is less constrained by numerical stability concerns (as with PDE) and more trying to live within a parametric budget. So if I have millions of points on a 2d grid, I might lay down some finite-budget mesh either based on some prior knowledge or based on data density. But my goal is a) control over-fitting and b) use a limited number of dofs.

[cottrell]

No point chasing something that's not going to be used; e.g. efficiency on pathological meshes when all practical meshes are going to be fine.

I agree but it seems to be a rather small change if I got it correctly: just verify that the correct element was found and bypass the KDTree optimization if not.

One thing I would check though is that points are not duplicated in that ix slicing. I think a np.unique there solves all of this no?

[cottrell]

Also just to confirm and clarify what I would be using this stuff for in case it helps. Suppose I assemble some basis and I have 100 data points corresponding to observations of values of the function. I would assemble this data-matrix "A" once and optimize that plus other terms somehow.

In [17]: A = basis.probes(X[:,:100]).shape; A.shape                                                                                                                  
Out[17]: (100, 19)
In [28]: theta.shape                                                                                                                                                                                               
Out[28]: (19,)
In [29]: y.shape                                                                                                                                                                                                   
Out[29]: (100,)
In [30]: np.sum((A @ theta - y) ** 2)  # this is a term in some objective function

And this example is just some gaussian noise case. But you can have more complicated loss functions. But the reason one turns to FEM is that the work of matching boundaries, handling deriviates etc is quite a bit of accounting to write oneself.

I think a lot of the ML community has been focussed on "learning the knots" (the mesh) which turns the problem into a fully non-linear problem (for example RBF networks) and you don't have any advantage in pre-computing these projection matrices since you're basis functions themselves are changing with each iteration. I am trying to stick with fixed meshes for now to aim for faster calibration times.

[gdmcbain]

Thanks for the Raket paper, I'll take a look; I'm quite at home with Green's functions, so that'll help.

(You'll see in #629 & #632 that the adjoint of the operator returned by the .probes method that we've been discussing above is the discrete approximation of Dirac's δ; hence the solution of a boundary value problem with it as right-hand side is the discrete Green's function.)

If horrid meshes aren't going to be the daily grind, I'd suggest solving an easier case first; it does make a big difference in the finite element method. If you do have the luxury of placing the points of the mesh wherever you like within the region, use something like Gmsh, dmsh, adaptmesh, or triangle, to place them reasonably.

[kinnala]

No point chasing something that's not going to be used; e.g. efficiency on pathological meshes when all practical meshes are going to be fine.

I agree but it seems to be a rather small change if I got it correctly: just verify that the correct element was found and bypass the KDTree optimization if not.

One thing I would check though is that points are not duplicated in that ix slicing. I think a np.unique there solves all of this no?

That might handle the part with lots of input points but not the part with pathological meshes. Sometimes the correct element simply is not part of ix if there are many slit elements.

[gdmcbain]

Sorry, I just read Raket (2021) and understood it (it's elementary) but didn't see the connexion to the present issue. The only domain considered there is 0 < x < 1, for which one could just use MeshLine as in

scikit-fem/docs/examples/ex15.py

Line 7 in 62882d6

m = MeshLine(np.linspace(0, 1, 10))

I assume the domains relevant to the application are harder than that. What are they like? How are they specified?

I think a lot of the ML community has been focussed on "learning the knots" (the mesh)

I wonder whether this is like automatic mesh adaptation, as in

scikit-fem/docs/examples/ex22.py

Line 1 in d08af2a

r"""Adaptive Poisson equation.

[cottrell]

Sorry, I just read Raket (2021) and understood it (it's elementary) but didn't see the connexion to the present issue. The only domain considered there is 0 < x < 1, for which one could just use MeshLine as in

scikit-fem/docs/examples/ex15.py

Line 7 in 62882d6

m = MeshLine(np.linspace(0, 1, 10))

I assume the domains relevant to the application are harder than that. What are they like? How are they specified?

I think a lot of the ML community has been focussed on "learning the knots" (the mesh)

I wonder whether this is like automatic mesh adaptation, as in

scikit-fem/docs/examples/ex22.py

Line 1 in d08af2a

r"""Adaptive Poisson equation.

@gdmcbain Yes, that paper is just a nice pedagogical linkage between various concepts. I think it was a course originally. The relevant application is super simple. You can of course do data assimilation with any PDE, SPDE etc.

There are packages like this https://github.com/alan-turing-institute/stat-fem and there are some papers in the README (https://arxiv.org/abs/1905.06391) for example (have not read it myself). Reason I am attacted to the set up with skfem is the pure python-ness of it so I can easily install or without any heavy extra stuff I do not need in images.

If the small branch pr I just pushed is correct it basically solves the large number of points problem for me. It seems to give same results and passes the tests locally. cottrell@f32ca21

I would also forget about the horrid mesh I provided as that is just some bad example I created quickly.

I think I am nearly there with my understanding of the layout of things. Will try to assemble a nice working example at the end of this.

[gdmcbain]

O. K., great. This is sounding good and very interesting. Girolami et al looks like the reference I needed. Thanks very much.

I see that statFEM depends on Firedrake. That's an impressive package but it is difficult to install. It was very similar considerations that led me from there to here.

On the branch pr, I hope I'm not holding things up by not having reviewed #667 #668; if so I'll get onto that tomorrow.

[kinnala]

No hurry, I'll look into combining #667 and #668 at some point. #668 is missing MeshTet1.element_finder and some additional tests and fixes that I added in #667. I'll try to also evaluate which approach gives better performance.

[cottrell]

Another thing that I'm wondering is if the intention is to include the data term as a bilinear form (with a sum of diracs)? Seems like it would match the flow more ... trying to see if there are examples like this. It might be that probes is still the main target for this kind of flow.

One thing I'm also looking for is if there is an established way of doing probes for a derivative or integral of the basis functions. For example if you have observational data about the derivatives instead of the the actual values of the solution.

[gdmcbain]

include the data term as a bilinear form (with a sum of diracs)? Seems like it would match the flow more ... trying to see if there are examples like this. It might be that probes is still the main target for this kind of flow.

Close. A Dirac point source or force is a linear form (RHS, rather than LHS for the bilinear forms). It's actually given by the transpose of CellBasis.probes which is here implemented in

https://scikit-fem.readthedocs.io/en/latest/listofexamples.html#example-38-point-source

scikit-fem/docs/examples/ex38.py

Line 46 in d08af2a

b = basis.point_source(source)

But is this exactly what's wanted in your problem? I'm not sure.

[gdmcbain]

One thing I'm also looking for is if there is an established way of doing probes for a derivative or integral of the basis functions. For example if you have observational data about the derivatives instead of the the actual values of the solution.

No, there isn't an established way but it's an interesting idea. If it's a functional it should be possible to express it as a matrix (ndarray) like that returned by CellBasis.probes. This matrix is really not a matrix like the discretization of a linear form, it's a stack of row vectors, one per probe-point.

For integrals there is

scikit-fem/skfem/assembly/form/functional.py

Line 10 in d08af2a

class Functional(Form):

Really a point-probe is mathematically a functional too but because its expression involves Dirac's δ it needs to be treated specially. Integrals over domains or boundaries or parts thereof are typically much better behaved and can just use Functional; I use them for things like lift coefficient, kinetic energy, &c.

[cottrell]

include the data term as a bilinear form (with a sum of diracs)? Seems like it would match the flow more ... trying to see if there are examples like this. It might be that probes is still the main target for this kind of flow.

Close. A Dirac point source or force is a linear form (RHS, rather than LHS for the bilinear forms). It's actually given by the transpose of CellBasis.probes which is here implemented in

https://scikit-fem.readthedocs.io/en/latest/listofexamples.html#example-38-point-source

scikit-fem/docs/examples/ex38.py

Line 46 in d08af2a

b = basis.point_source(source)

But is this exactly what's wanted in your problem? I'm not sure.

So I am just kind of remembering (and misremembering) how all this fits together now. I was writing down the attached variational formulation and it made me realize that yes, we could think of a "dirac" basis for the point evalutions as mentioned in the docstring, but it is probably not super useful to literally implement a Form in the code against a Dirac basis. I mean, that is exactly what probes is doing but I'm not sure expressing the probes evaluation as some Form against this particular basis is ever useful for the current solver flows.

... unless of course there are examples with representations in different bases (perhaps the boudary examples) that I am not aware of yet.

[gdmcbain]

What about splitting the right-hand side? The term (δ(x − s) u (x), h (x)) is a bilinear form and so could go on the left, the term (δ(x − s) y (x), h (x)) is a linear form and stays on the right.

[gdmcbain]

Since (δ(x − s) u (x), h (x)) = u (s) h (s), the ij-th element of that discretized bilinear form would be the product of the values of the i-th and j-th basis functions at the probe-point s. But I'll leave it at that till I've read Girolami et al; I'm sure they would've sorted all this out already.