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 #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 \sumi | \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]

@kinnala L above is some loss function. Can think of squared error for gaussian on something more complicated for something like quantile regression or another noise model.

I'm not look for the solver, just the accounting/set up of the basis functions. It is quite a lot of things to keep track of to match the boundary conditions, map to local coordinates, do the inference at specific points etc.

The nice thing about these kind of basis function approches is, for a fixed mesh and for a fixed set of data points (points to value the surface) you don't need to do anything with the fem after that. You just minimize a cost function w.r.t the coefficients.

[gdmcbain]

This reminds me of #593 from the author of smoothfit. Is that close? Following the discussion there, smoothfit does have code using skfem.

[gdmcbain]

There are sone preliminary experiments on minimization in #642, which is a way of having a certain class of nonlinear problems automatically solved in scikit-fem, viz. those that can be posed as constrained minimizations, with an objective, jacobian, hessian, and constraints. That issue is currently shelved as the scipy.optimize solver is weak. I was able to formulate and solve small problems elegantly, e.g. minimal surface, but thus far it's not competitive with hand-written Newton iteration.

[cottrell]

There are sone preliminary experiments on minimization in #642, which is a way of having a certain class of nonlinear problems automatically solved in scikit-fem, viz. those that can be posed as constrained minimizations, with an objective, jacobian, hessian, and constraints. That issue is currently shelved as the scipy.optimize solver is weak. I was able to formulate and solve small problems elegantly, e.g. minimal surface, but thus far it's not competitive with hand-written Newton iteration.

So I would completely off-load the optimization as another task (to be implemented however one wants). For example, one could combine multiple representations.

And for complicated non-linear stuff jax (or similar) can help a lot with the optimization.

Am looking through the examples you mentioned ... had not come across smoothfit yet either.

[cottrell]

Yep, it looks like smoothfit is doing this same thing with skfem:

https://github.com/nschloe/smoothfit/blob/main/src/smoothfit/main.py#L39

[gdmcbain]

Yeah, O. K.; so I guess it suffers the same blow-up in probes at L110?

[cottrell]

FWIW, this is some initial hack I've done to assemble what should be the basis function evaluations over an arbitrary set of "observation" points.

from pylab import *
import skfem as fem
from scipy.spatial import Delaunay
import numpy as np

def create_random_mesh(M_int=10, M_boundary=10, d=2, seed=1):
    np.random.seed(seed)

    # generate mesh points
    X_int = np.random.rand(M_int, d)

    def boundary_map(x):
        # around counter clockwise
        x0 = np.array([0, 0])[None, :]
        x1 = np.array([1, 0])[None, :]
        x2 = np.array([1, 1])[None, :]
        x3 = np.array([0, 1])[None, :]
        x = 4 * x[:, None]
        condlist = [x < 1, (1 <= x) & (x < 2), (2 <= x) & (x < 3), (3 <= x) & (x < 4)]
        funclist = [x * (x1 - x0), x1 + (x2 - x1) * (x - 1), x2 + (x3 - x2) * (x - 2), x3 + (x0 - x3) * (x - 3)]
        out = np.zeros((x.shape[0], 2))
        for c, f in zip(condlist, funclist):
            out = out + c * f
        return out

    X_boundary = boundary_map(np.linspace(0, 1, M_boundary))
    X = np.vstack([X_int, X_boundary])
    tri = Delaunay(X)
    mesh = fem.MeshTri(np.ascontiguousarray(tri.points.T), np.ascontiguousarray(tri.vertices.T))
    return mesh

mesh = create_random_mesh()
display(mesh)

# grid of observations for testing
nx = 200
ny = 200
def get_obs_grid(nx=nx, ny=ny):
    X, Y = np.meshgrid(np.linspace(0, 1, nx), np.linspace(0, 1, ny))
    return np.vstack([X.ravel(), Y.ravel()])

X = get_obs_grid()

# assemble a basis
element_degree = 1
n_elements_per_mesh_element = int((element_degree + 1) * (element_degree + 2) / 2)

if element_degree == 2:
    element = fem.ElementTriP2()
elif element_degree == 1:
    element = fem.ElementTriP1()
else:
    raise
basis = fem.CellBasis(mesh, element)
facet_basis = fem.FacetBasis(basis.mesh, basis.elem)

# assemble the basis evaluations
def my_mapping_affine_invF(mapping, x):
    """mapping is skfem.mapping.mapping_affine.MappingAffine instance"""
    return np.einsum('ijk,jkl->ikl', mapping.invA, x[:, None, :] - mapping.b[:, :, None])
    # return np.einsum('ijk,jl->ikl', invA, x) - np.einsum('ijk,jk->ik', invA, b)[:, :, None]  # alternative, could be better in some cases

def my_inside_mask(X):
    # (2, n_elem, n_obs) -> (n_elem, n_obs)
    return (X[0] >= 0) * (X[1] >= 0) * (1 - X[0] - X[1] >= 0)

X_local = my_mapping_affine_invF(basis.mapping, X)  # (2, 297, 30_000)
mask = my_inside_mask(X_local)  # you need this for masking the values since we are applying the inverse mapping everywhere?
# assert np.sum(mask, axis=0).max() == 1, 'each col should only be one element' # TODO WHY DOES THIS FAIL SOMETIMES?

# finite elements for each mesh element
phi = np.zeros((X_local.shape[1], X_local.shape[2], n_elements_per_mesh_element))
dphi = np.zeros((2, X_local.shape[1], X_local.shape[2], n_elements_per_mesh_element))
for ind in range(n_elements_per_mesh_element):
    phi_, dphi_ = basis.elem.lbasis(X_local, ind)
    phi[:, :, ind] = phi_
    dphi[:, :, :, ind] = dphi_

display(mesh.p.shape, mesh.t.shape, X.shape, X_local.shape, mask.shape, phi.shape, dphi.shape)

i_elem = 23
fig = figure(figsize=(20, 4))
ax = fig.subplots(1, n_elements_per_mesh_element)
for i in range(n_elements_per_mesh_element):
    z = phi[i_elem, :, i]
    x = X[0].reshape(nx, ny)
    y = X[1].reshape(nx, ny)
    z = z.reshape(nx, ny)
    mask_ = mask[i_elem].reshape(nx, ny)
    ax[i].imshow(z * mask_)

i_elem = 10
fig = figure(figsize=(20, 4))
ax = fig.subplots(1, n_elements_per_mesh_element)
for i in range(n_elements_per_mesh_element):
    z = phi[i_elem, :, i]
    x = X[0].reshape(nx, ny)
    y = X[1].reshape(nx, ny)
    z = z.reshape(nx, ny)
    mask_ = mask[i_elem].reshape(nx, ny)
    ax[i].imshow(z * mask_)

There are some things in the code logic starting at https://github.com/kinnala/scikit-fem/blob/master/skfem/assembly/basis/cell_basis.py#L146 that I don't really understand ... it looks like things are going to local than global coords and back again or something.

[kinnala]

In FEM, basis functions can only be evaluated locally. Thus, to evaluate the basis at a global point, we must first find the corresponding element (element_finder) and then the corresponding point at the reference triangle (invF). Using the point at the reference triangle you can finally evaluate the FEM basis. The gradient must finally be transformed back to the global coordinate system as using the tensor change of variables rule.

[kinnala]

Or maybe you are meaning what is going on inside element_finder? As stated before, the implementation of that method is suboptimal because it will perform invF in order to find out if a point is inside the element or not. It will further reduce the search space by fetching 5 candidate elements from scipy KDTree based on midpoint distance which will make the performance reasonable for a small number of points (even for large meshes), but it is still suboptimal for a large number of points.

An optimal implementation would be easy in a compiled language but here will require some numpy wizardry.

[kinnala]

Reading at your implementation, it will probably not work for a large mesh because it does not prune elements that are far away from the point and, thus, perform invF always for all elements.

[kinnala]

What I'm trying to say that this line

X_local = my_mapping_affine_invF(basis.mapping, X)  # (2, 297, 30_000)

will again blow your memory if you have 30 000 elements and 30 000 points.

So instead of optimizing for the number of points we went for the number of elements in order for it to work also in 3D cases where you usually have a significant number of elements.

[cottrell]

Reading at your implementation, it will probably not work for a large mesh because it does not prune elements that are far away from the point and, thus, perform invF always for all elements.

Yes, it will definitely not work for large mesh as I have not using the NN kdtree stuff yet. This is explicitly for the N_data >> N_mesh problem which only comes up in smoothing AFAIK.

I'm sure it can be improved using the kdtree NN tricks in the skfem code here. I haven't quite grokked whether or not the bit about the N_data that I did can be just used in the skfem code to avoid the problems as N_data gets large.

[kinnala]

You are right, that could be possible. I'll look into it.

[kinnala]

See #667, I'll finish it later.

[cottrell]

See #667, I'll finish it later.

Just checked out the branch and I am wondering what the output of probes is supposed to look like (in general not just on the branch). It seems like the dimensions are different from what I was expecting:

(from the above example)

[kinnala]

I'm not sure what is phi in the snippet, but I think 19 is the number of degrees-of-freedom since you are supposed to use it like E @ y where y is a DOF array.

[cottrell]

Yeah I think 19 is suppose to be the degrees of freedom but it is not what I expect or I don't know where to pull that from. phi is supposed to be the same thing. I think every element has 3 nodes or something.

I still don't really understand the logic of 5 nearest-neighbours trick ... I might have at one point but I've forgotten again. I suspect that when duplicate ix entries appear it is possibly not necessary.

[kinnala]

Each element has 3 nodes, yes, but globally some of those are shared with the neighbouring elements. Depends on what kind of finite element and mesh you are using, but for ElementTriP1 you have one DOF per vertex. You might be confusing local and global bases.

Your case is different, typically you have maybe 10 points and 1 million elements. Then it's less work to check 5 x 10 elements than all 1 million elements.

[kinnala]

Let me add that if you really want 3 DOFs per element then the correct object is ElementTriDG(ElementTriP1()) but that's unlikely unless you want to use a globally discontinuous basis.

[cottrell]

@kinnala Don't add the DOFs yet as I have no idea what is going on yet :) ... just trying to understand if probes is actually used by anyone. I am not sure it is correct at the moment.

Adding some notes for when I come back to this (I'll lose them otherwise):

I've added some debug (ongoing) notes in the branch below with a _v2 and _v3 version of the finder routine. v2 agrees with the original. I do not undestand how the intention with the argmax. It seems like there is a dimension missing as I would imagine the intention is to argmax over the 5 neighbours not the entire multiply replicated set of nodes (mesh.t).

https://github.com/cottrell/scikit-fem/tree/debug_finder

[kinnala]

You can grep over the examples and see that it's used by several of them. I don't use it personally but I think I understand the implementation well enough to have reasonable faith that it works as intended in most cases. I could probably also think of degenerate cases where it doesn't work (in particular, Mesh.element_finder) but I'm not too worried about that right now because people often tend to use reasonable meshes.

There are some tests for it under tests/test_basis.py and tests/test_mesh.py. If you have additional tests to propose, feel free to open a PR or another issue.

[kinnala]

Now that I look at the above Delaunay mesh, it doesn't actually look like what I usually think of a "reasonable FEM mesh", so it might well be a case which breaks the heuristic of using 5 nearest midpoints: for crappy meshes you probably either want to increase 5 to something larger or actually implement Mesh.element_finder properly using graph search.

[kinnala]

The reason why this (https://github.com/cottrell/scikit-fem/commit/9c3333c3db23d3c51ff88e58e24c520ff54a6fe3#diff-a18285e3df75b5c062df233ef0b63edd67fd9571acab48d72b37608b3cefe187R356) appears to replicate rows is the test case. Use instead say 5 points and 1e5 elements which is the original use case.

In my PR I dropped the KDtree optimization if there are more than 1e3 points which is also a rather heuristic limit but I guess the allocations become quite large after that. Dropping the optimization altogether would make the original use case impossible so unless someone implements this properly I cannot remove it even if it doesn't work with some nasty element shapes.

[kinnala]

One improvement I can think of right now is a check that if none of the 5 candidate elements contain the query point then falling back to the exhaustive search would make sense.

[gdmcbain]

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.

[kinnala]

One improvement I can think of right now is a check that if none of the 5 candidate elements contain the query point then falling back to the exhaustive search would make sense.

I made an attempt to implement such fallback behavior in #667 . I'm testing it with some randomized Delaunay meshes where I noticed that in about 0.1% of the cases the correct element is not among the candidates when comparing to scipy.spatial.Delaunay.find_simplex. Hoping for feedback.

[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.

[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.

The example is not meant to be a "good mesh" in any sense. Just testing robustness.

Data Assimilation is basically just inference/statistics/ML whatever you want to call it. The words "data assimilation" were more common in atmospheric sciences at one point. Typically it is something where one is combing model and data. So you have some physics (usually pde) and some noise/observation model for data. And you minimize some likelihood. For very simple cases where you have linear phyics and quadratic noise modesl (gaussians) the critical point equations.

This paper gives some great connections between splines, gaussian process regression and poisson equation https://arxiv.org/abs/2102.03306 that points to the simplest "smoothing" problem.

[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

https://github.com/kinnala/scikit-fem/blob/62882d61951e5283196cf22b9b9924b86f5b243b/docs/examples/ex15.py#L7

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

https://github.com/kinnala/scikit-fem/blob/d08af2a3f4607e28023656b1b6ad025b1a530aa1/docs/examples/ex22.py#L1

[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

https://github.com/kinnala/scikit-fem/blob/62882d61951e5283196cf22b9b9924b86f5b243b/docs/examples/ex15.py#L7

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

https://github.com/kinnala/scikit-fem/blob/d08af2a3f4607e28023656b1b6ad025b1a530aa1/docs/examples/ex22.py#L1

@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. https://github.com/cottrell/scikit-fem/commit/f32ca21198f30536fa3c70a7d4a447746d66ef98

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.