Is there an inverse of basis.interpolate()?

[gatling-nrl]

Just that... Is there an inverse of basis.interpolate()?

[kinnala]

In theory it is L2 projection, i.e. skfem.utils.projection. However, I'm still not completely happy about skfem.utils.projection. I'm working on a new version which will hopefully also be useful for this.

This is what I have in mind, proj() is the new function, and has similar syntax to asm():

basis = CellBasis(mesh, ElementTriP1())

# w is like the last parameter to forms
y1 = proj(lambda w: w.x[0], basis)

# take derivative with gradient averaging
x = mesh.p[0]
y2 = proj(lambda w: w['u'].grad[0], basis, u=x)

# real derivative
y3 = proj(lambda w: w['u'].grad[0], basis.with_element(ElementTriP0()), u=basis.interpolate(x))

# inverse
y4 = proj(lambda w: w['u'], basis, u=basis.interpolate(x))
# y4 should be almost equal to x

[gatling-nrl]

but I don't think there is a need to change the names again

I strongly agree. I also understand better your thoughts: refine() to modify and existing object in place and refined() to return a new object and leave the original reference alone.

Frequently I am left wondering which of those two things is happening when I call numpy functions. (Sadly, even familiar ones...) So I often have to resort to google, once again. Nevertheless, I prefer verbs and clear docstrings... and very much prefer backward compatibility over consistent or pedantic naming.

[gatling-nrl]

Let me add that the name proj above is somewhat arbitrary, I just wanted to demonstrate the functionality, it might be possible to implement in projection in a backwards compatible manner.

Can your new ideas not be implemented in AbstractBasis.proj(ect)? It is just such an important (I think) feature of skfem that it seems unintuitve to have it relegated to something called utils in the first place. Isn't it a first class type operation in the semantics? Can you even do the projecting without the AbstractBasis object?

[kinnala]

Here is a simple draft implementation for proj above:

def proj(fun, basis, **kwargs):
    from skfem.models import mass
    M = mass.assemble(basis)
    f = LinearForm(lambda v, w: fun(w) * v).assemble(basis, **kwargs)
    return solve(M, f)

As you can see, it is much simpler than the implementation of projection and takes care of almost all special cases supported by projection. Adding this to skfem.utils is technically not a backwards incompatible change (in the spirit of semantic versioning) since we are adding another function. Only removing or changing functionality is backwards incompatible.

[kinnala]

Let me add that the name proj above is somewhat arbitrary, I just wanted to demonstrate the functionality, it might be possible to implement in projection in a backwards compatible manner.

Can your new ideas not be implemented in AbstractBasis.proj(ect)? It is just such an important (I think) feature of skfem that it seems unintuitve to have it relegated to something called utils in the first place. Isn't it a first class type operation in the semantics? Can you even do the projecting without the AbstractBasis object?

Implementation inside AbstractBasis would be something like:

class AbstractBasis:
    ...
    def project(self, fun, **kwargs):
        from skfem.models import mass
        from skfem.utils import solve
        from skfem.assembly import LinearForm
        M = mass.assemble(basis)
        f = LinearForm(lambda v, w: fun(w) * v).assemble(self, **kwargs)
        return solve(M, f)

[kinnala]

This is how the snippets would look like:

basis = CellBasis(mesh, ElementTriP1())

# w is like the last parameter to forms
y1 = basis.project(lambda w: w.x[0])

# take derivative with gradient averaging
x = mesh.p[0]
y2 = basis.project(lambda w: w['u'].grad[0], u=x)

# real derivative
basis0 = basis.with_element(ElementTriP0())
y3 = basis0.project(lambda w: w['u'].grad[0], u=basis.interpolate(x))

# inverse
y4 = basis.project(lambda w: w['u'], u=basis.interpolate(x))
# y4 should be almost equal to x

I think it kind of makes sense.

[kinnala]

Alright, to answer the original question, with the existing projection it would like this:

import numpy as np
from skfem import *
m = MeshTri()
basis = Basis(m, ElementTriP1())
x1 = m.p[0]  # the function f(x) = x
x1_interp = basis.interpolate(x1)
x2 = projection(LinearForm(lambda v, w: x1_interp * v), basis)
print(np.linalg.norm(x1 - x2))

which will print 2.6766507790745728e-16.

[gatling-nrl]

$ wmic cpu get caption, deviceid, name, numberofcores, maxclockspeed, status
Caption                                 DeviceID  MaxClockSpeed  Name                                      NumberOfCores  Status
Intel64 Family 6 Model 158 Stepping 12  CPU0      3600           Intel(R) Core(TM) i9-9900K CPU @ 3.60GHz  8              OK

[kinnala]

This shows I'm using OpenBLAS and not MKL:

In [2]: numpy.show_config()
blas_mkl_info:
  NOT AVAILABLE
blis_info:
  NOT AVAILABLE
openblas_info:
    libraries = ['lapack', 'lapacke', 'blas', 'cblas', 'lapack', 'lapacke', 'blas', 'cblas']
    library_dirs = ['/nix/store/5gq47dyyzjz4p2h603drb2srs66hr2fr-lapack-3/lib', '/nix/store/2cz8569qf7khb1cqdli88y1ywwyqw5b3-blas-3/lib']
    language = c
    define_macros = [('HAVE_CBLAS', None)]
    runtime_library_dirs = ['/nix/store/5gq47dyyzjz4p2h603drb2srs66hr2fr-lapack-3/lib', '/nix/store/2cz8569qf7khb1cqdli88y1ywwyqw5b3-blas-3/lib']
blas_opt_info:
    libraries = ['lapack', 'lapacke', 'blas', 'cblas', 'lapack', 'lapacke', 'blas', 'cblas']
    library_dirs = ['/nix/store/5gq47dyyzjz4p2h603drb2srs66hr2fr-lapack-3/lib', '/nix/store/2cz8569qf7khb1cqdli88y1ywwyqw5b3-blas-3/lib']
    language = c
    define_macros = [('HAVE_CBLAS', None)]
    runtime_library_dirs = ['/nix/store/5gq47dyyzjz4p2h603drb2srs66hr2fr-lapack-3/lib', '/nix/store/2cz8569qf7khb1cqdli88y1ywwyqw5b3-blas-3/lib']
lapack_mkl_info:
  NOT AVAILABLE
openblas_lapack_info:
    libraries = ['lapack', 'lapacke', 'blas', 'cblas', 'lapack', 'lapacke', 'blas', 'cblas']
    library_dirs = ['/nix/store/5gq47dyyzjz4p2h603drb2srs66hr2fr-lapack-3/lib', '/nix/store/2cz8569qf7khb1cqdli88y1ywwyqw5b3-blas-3/lib']
    language = c
    define_macros = [('HAVE_CBLAS', None)]
    runtime_library_dirs = ['/nix/store/5gq47dyyzjz4p2h603drb2srs66hr2fr-lapack-3/lib', '/nix/store/2cz8569qf7khb1cqdli88y1ywwyqw5b3-blas-3/lib']
lapack_opt_info:
    libraries = ['lapack', 'lapacke', 'blas', 'cblas', 'lapack', 'lapacke', 'blas', 'cblas']
    library_dirs = ['/nix/store/5gq47dyyzjz4p2h603drb2srs66hr2fr-lapack-3/lib', '/nix/store/2cz8569qf7khb1cqdli88y1ywwyqw5b3-blas-3/lib']
    language = c
    define_macros = [('HAVE_CBLAS', None)]
    runtime_library_dirs = ['/nix/store/5gq47dyyzjz4p2h603drb2srs66hr2fr-lapack-3/lib', '/nix/store/2cz8569qf7khb1cqdli88y1ywwyqw5b3-blas-3/lib']

Moreover,

In [3]: numpy.__version__
Out[3]: '1.20.2'

In [4]: import scipy

In [5]: scipy.__version__
Out[5]: '1.6.1'

[gatling-nrl]

np.show_config()
blas_info:
    libraries = ['cblas', 'blas', 'cblas', 'blas', 'cblas', 'blas']
    library_dirs = ['C:/Users/ggatling/Miniconda3\\Library\\lib']
    include_dirs = ['C:/Users/ggatling/Miniconda3\\Library\\include']
    language = f77
    define_macros = [('HAVE_CBLAS', None)]
blas_opt_info:
    define_macros = [('NO_ATLAS_INFO', 1), ('HAVE_CBLAS', None)]
    libraries = ['cblas', 'blas', 'cblas', 'blas', 'cblas', 'blas']
    library_dirs = ['C:/Users/ggatling/Miniconda3\\Library\\lib']
    include_dirs = ['C:/Users/ggatling/Miniconda3\\Library\\include']
    language = f77
lapack_info:
    libraries = ['lapack', 'blas', 'lapack', 'blas']
    library_dirs = ['C:/Users/ggatling/Miniconda3\\Library\\lib']
    language = f77
lapack_opt_info:
    libraries = ['lapack', 'blas', 'lapack', 'blas', 'cblas', 'blas', 'cblas', 'blas', 'cblas', 'blas']
    library_dirs = ['C:/Users/ggatling/Miniconda3\\Library\\lib']
    language = f77
    define_macros = [('NO_ATLAS_INFO', 1), ('HAVE_CBLAS', None)]
    include_dirs = ['C:/Users/ggatling/Miniconda3\\Library\\include']
Supported SIMD extensions in this NumPy install:
    baseline = SSE,SSE2,SSE3
    found = SSSE3,SSE41,POPCNT,SSE42,AVX,F16C,FMA3,AVX2
    not found = AVX512F,AVX512CD,AVX512_SKX,AVX512_CLX,AVX512_CNL

np.__version__
'1.21.2'
scipy.__version__
'1.7.1'

[kinnala]

It seems that SuperLU is included in scipy GitHub repository so basically the sparse solver version depends on the scipy version.

[kinnala]

Alright, I know that the version of numpy in conda is optimized for performance and will most likely use some specific compilation flags or different BLAS version. I'm using a build script from nixpkgs which is just a really basic build: https://github.com/NixOS/nixpkgs/blob/master/pkgs/development/python-modules/scipy/default.nix

Also the version number difference could explain it.