technologicat/python-wlsqm

Weighted least squares meshless interpolator and differentiator.

wlsqm

Weighted least squares meshless interpolator

Introduction

WLSQM (Weighted Least SQuares Meshless) is a fast and accurate meshless least-squares interpolator for Python, for scalar-valued data defined as point values on 1D, 2D and 3D point clouds.

Use cases include response surface modeling, and computing space derivatives of data known only as values at discrete points in space (this has applications in explicit algorithms for solving IBVPs). No grid or mesh is needed. No restriction is imposed on geometry other than "not degenerate", e.g. points in 2D should not all fall onto the same 1D line.

This is an independent implementation of the weighted least squares meshless algorithm described (in the 2nd order 2D case) in section 2.2.1 of Hong Wang (2012), Evolutionary Design Optimization with Nash Games and Hybridized Mesh/Meshless Methods in Computational Fluid Dynamics, Jyväskylä Studies in Computing 162, University of Jyväskylä. ISBN 978-951-39-5007-1 (PDF)

This implementation is targeted for high performance in a single-node environment, such as a laptop. Cython is used to accelerate the low-level routines. The main target is the x86_64 architecture, but any 64-bit architecture should be fine with the appropriate compiler option changes to setup.py.

Currently automated unit tests are missing; this is an area that is likely to be improved. Otherwise the code is already rather stable; any major new features are unlikely to be added, and the API is considered stable.

Features

• Given scalar data values on a set of points in 1D, 2D or 3D, construct a piecewise polynomial global surrogate model (a.k.a. response surface), using up to 4th order polynomials.

• Sliced arrays are supported for input, both for the geometry (points) and data (function values).

• Obtain any derivative of the model, up to the order of the polynomial. Derivatives at each "local model reference point" xi are directly available as DOFs of the solution. Derivatives at any other point can be automatically interpolated from the model. Differentiation of polynomials has been hardcoded to obtain high performance.

• Knowns. At the model reference point xi, the function value and/or any of the derivatives can be specified as knowns. The knowns are internally automatically eliminated (making the equation system smaller) and only the unknowns are fitted. The function value itself may also be unknown, which is useful for implementing Neumann BCs in a PDE (IBVP) solving context.

• Selectable weighting method for the fitting error, to support different use cases:

  • uniform (wlsqm.fitter.defs.WEIGHT_UNIFORM), for best overall fit for function values
  • emphasize points closer to xi (wlsqm.fitter.defs.WEIGHT_CENTER), to improve derivatives at the reference point xi by reducing the influence of points far away from the reference point.

• Sensitivity data of solution DOFs (on the data values at points other than the reference in the local neighborhood) can be optionally computed.

• Expert mode with separate prepare and solve stages, for faster fitting of many data sets using the same geometry. Also performs global model patching, using the set of local models fitted.

  CAVEAT: wlsqm.fitter.expert.ExpertSolver instances are not currently pickleable or copyable. This is a known limitation that may (or may not) change in the future.

  It is nevertheless recommended to use ExpertSolver, since this allows for easy simultaneous solving of many local models (in parallel), automatic global model patching, and reuse of problem matrices when the geometry of the point cloud does not change.

• Speed:

  • Performance-critical parts are implemented in Cython, and the GIL is released during computation.

  • LAPACK is used directly via SciPy's Cython-level bindings (see the ntasks parameter in various API functions in wlsqm). This is especially useful when many (1e4 or more) local models are being fitted, as the solver loop does not require holding the GIL.

  • OpenMP is used for parallelization over the independent local problems (also in the linear solver step).

  • The polynomial evaluation code has been manually optimized to reduce the number of FLOPS required.

    In 1D, the Horner form is used. The 2D and 3D cases use a symmetric form that extends the 1D Horner form into multiple dimensions (see wlsqm/fitter/polyeval.pyx for details). The native FMA (fused multiply-add) instruction of the CPU is used in the evaluation to further reduce FLOPS required, and to improve accuracy (utilizing the fact it rounds only once).

• Accuracy:

  • Problem matrices are preconditioned by a symmetry-preserving scaling algorithm (D. Ruiz 2001; exact reference given in wlsqm/utils/lapackdrivers.pyx) to obtain best possible accuracy from the direct linear solver. This is critical especially for high-order fits.
  • The fitting procedure optionally accomodates an internal iterative refinement loop to mitigate the effect of roundoff.
  • FMA, as mentioned above.

Installation

From PyPI (wlsqm v0.1.1+)

Install as user:

pip install wlsqm --user

Install as admin:

sudo pip install wlsqm

From GitHub

As user:

git clone https://github.com/Technologicat/python-wlsqm.git
cd python-wlsqm
python setup.py install --user

As admin, change the last command to

sudo python setup.py install

Documentation

For usage examples, see examples/wlsqm_example.py for a tour, and examples/expertsolver_example.py for a minimal example concentrating specifically on ExpertSolver.

For the technical details, see the docstrings and comments in the code itself.

Mathematics documented at:

• https://yousource.it.jyu.fi/jjrandom2/freya/trees/master/docs [dead link, relevant files mirrored below]

where the relevant files are [mirrored locally on GitHub]:

• wlsqm.pdf (old documentation for the old pure-Python version of WLSQM included in FREYA, plus the sensitivity calculation)
• eulerflow.pdf (clearer presentation of the original version, but without the sensitivity calculation)
• wlsqm_gen.pdf (theory diff on how to make a version that handles also missing function values; also why WLSQM works and some analysis of its accuracy)

The documentation is slightly out of date; see TODO for details on what needs updating and how.

Experiencing crashes?

Check that you are loading the same BLAS your LAPACK and SciPy link against:

shopt -s globstar
ldd /usr/lib/**/*lapack*.so | grep blas
ldd $(dirname $(python -c "import scipy; print(scipy.__file__)"))/linalg/cython_lapack.so | grep blas

In Debian-based Linux, you can change the active BLAS implementation by:

sudo update-alternatives --config libblas.so
sudo update-alternatives --config libblas.so.3

This may (or may not) be different from what NumPy links against:

ldd $(dirname $(python -c "import numpy; print(numpy.__file__)"))/core/multiarray.so | grep blas

WLSQM itself does not link against LAPACK or BLAS; it utilizes the cython_lapack module of SciPy.

Dependencies

• NumPy
• SciPy
• Cython
• Matplotlib (for usage examples)

License

BSD. Copyright 2016-2017 Juha Jeronen and University of Jyväskylä.

Acknowledgement

This work was financially supported by the Jenny and Antti Wihuri Foundation.

python-wlsqm's Issues

build error

brilliant work！but I have some problems
python3.6
pip install wlsqm

wlsqm/utils/lapackdrivers.pyx:1740:38: Converting to Python object not allowed without gil

Error compiling Cython file:
------------------------------------------------------------
...
    dgesvd( &jobu, &jobvt, &m, &n, A, &m, S, <double*>0, &ldu, <double*>0, &ldvt, &worksize, &lwork, &info )
    lwork  = <int>worksize
    cdef double* p_work = <double*>malloc( lwork*sizeof(double) )

    # LDA(A) = n
    dgesvd( &jobu, &jobvt, &m, &n, A, &m, S, <double*>0, &ldu, <double*>0, &ldvt, p_work, &lwork, &info )
                                         ^
------------------------------------------------------------

wlsqm/utils/lapackdrivers.pyx:1740:42: Converting to Python object not allowed without gil

Error compiling Cython file:
------------------------------------------------------------
...
    dgesvd( &jobu, &jobvt, &m, &n, A, &m, S, <double*>0, &ldu, <double*>0, &ldvt, &worksize, &lwork, &info )
    lwork  = <int>worksize
    cdef double* p_work = <double*>malloc( lwork*sizeof(double) )

    # LDA(A) = n
    dgesvd( &jobu, &jobvt, &m, &n, A, &m, S, <double*>0, &ldu, <double*>0, &ldvt, p_work, &lwork, &info )
                                            ^
------------------------------------------------------------

wlsqm/utils/lapackdrivers.pyx:1740:45: Converting to Python object not allowed without gil

Error compiling Cython file:
------------------------------------------------------------
...
    dgesvd( &jobu, &jobvt, &m, &n, A, &m, S, <double*>0, &ldu, <double*>0, &ldvt, &worksize, &lwork, &info )
    lwork  = <int>worksize
    cdef double* p_work = <double*>malloc( lwork*sizeof(double) )

    # LDA(A) = n
    dgesvd( &jobu, &jobvt, &m, &n, A, &m, S, <double*>0, &ldu, <double*>0, &ldvt, p_work, &lwork, &info )
                                                        ^
------------------------------------------------------------

wlsqm/utils/lapackdrivers.pyx:1740:57: Converting to Python object not allowed without gil

Error compiling Cython file:
------------------------------------------------------------
...
    dgesvd( &jobu, &jobvt, &m, &n, A, &m, S, <double*>0, &ldu, <double*>0, &ldvt, &worksize, &lwork, &info )
    lwork  = <int>worksize
    cdef double* p_work = <double*>malloc( lwork*sizeof(double) )

    # LDA(A) = n
    dgesvd( &jobu, &jobvt, &m, &n, A, &m, S, <double*>0, &ldu, <double*>0, &ldvt, p_work, &lwork, &info )
                                                              ^
------------------------------------------------------------

wlsqm/utils/lapackdrivers.pyx:1740:63: Converting to Python object not allowed without gil

Error compiling Cython file:
------------------------------------------------------------
...
    dgesvd( &jobu, &jobvt, &m, &n, A, &m, S, <double*>0, &ldu, <double*>0, &ldvt, &worksize, &lwork, &info )
    lwork  = <int>worksize
    cdef double* p_work = <double*>malloc( lwork*sizeof(double) )

    # LDA(A) = n
    dgesvd( &jobu, &jobvt, &m, &n, A, &m, S, <double*>0, &ldu, <double*>0, &ldvt, p_work, &lwork, &info )
                                                                          ^
------------------------------------------------------------

wlsqm/utils/lapackdrivers.pyx:1740:75: Converting to Python object not allowed without gil

Error compiling Cython file:
------------------------------------------------------------
...
    dgesvd( &jobu, &jobvt, &m, &n, A, &m, S, <double*>0, &ldu, <double*>0, &ldvt, &worksize, &lwork, &info )
    lwork  = <int>worksize
    cdef double* p_work = <double*>malloc( lwork*sizeof(double) )

    # LDA(A) = n
    dgesvd( &jobu, &jobvt, &m, &n, A, &m, S, <double*>0, &ldu, <double*>0, &ldvt, p_work, &lwork, &info )
                                                                                 ^
------------------------------------------------------------

wlsqm/utils/lapackdrivers.pyx:1740:82: Converting to Python object not allowed without gil

Error compiling Cython file:
------------------------------------------------------------
...
    dgesvd( &jobu, &jobvt, &m, &n, A, &m, S, <double*>0, &ldu, <double*>0, &ldvt, &worksize, &lwork, &info )
    lwork  = <int>worksize
    cdef double* p_work = <double*>malloc( lwork*sizeof(double) )

    # LDA(A) = n
    dgesvd( &jobu, &jobvt, &m, &n, A, &m, S, <double*>0, &ldu, <double*>0, &ldvt, p_work, &lwork, &info )
                                                                                         ^
------------------------------------------------------------

wlsqm/utils/lapackdrivers.pyx:1740:90: Converting to Python object not allowed without gil

Error compiling Cython file:
------------------------------------------------------------
...
    dgesvd( &jobu, &jobvt, &m, &n, A, &m, S, <double*>0, &ldu, <double*>0, &ldvt, &worksize, &lwork, &info )
    lwork  = <int>worksize
    cdef double* p_work = <double*>malloc( lwork*sizeof(double) )

    # LDA(A) = n
    dgesvd( &jobu, &jobvt, &m, &n, A, &m, S, <double*>0, &ldu, <double*>0, &ldvt, p_work, &lwork, &info )
                                                                                                 ^
------------------------------------------------------------

wlsqm/utils/lapackdrivers.pyx:1740:98: Converting to Python object not allowed without gil
Traceback (most recent call last):
  File "setup.py", line 220, in <module>
    gdb_debug = debug ),
  File "/home/duan/miniconda3/envs/cv/lib/python3.6/site-packages/Cython/Build/Dependencies.py", line 1097, in cythonize
    cythonize_one(*args)
  File "/home/duan/miniconda3/envs/cv/lib/python3.6/site-packages/Cython/Build/Dependencies.py", line 1220, in cythonize_one
    raise CompileError(None, pyx_file)
Cython.Compiler.Errors.CompileError: wlsqm/utils/lapackdrivers.pyx

use python3.4
successful installed, but

>>> import wlsqm
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/home/duan/PycharmProjects/vessel/python-wlsqm/wlsqm/__init__.py", line 26, in <module>
    from .fitter.defs   import *  # definitions (constants) (common)
ImportError: No module named 'wlsqm.fitter.defs'
>>> 

Python 3.6 build error

Fails to compile under Python 3.6. See #2.

High priority now that Python 3.4 official support has ended (in March 2019), and many new-ish installed pythons are at least 3.6.

Update contact email

The contact email in setup.py is outdated. PyPI uses it to provide the link to email the maintainer. Update it for the next release.

ValueError: Buffer and memoryview are not contiguous in the same dimension.

I try to use this code to form the derivates of my meshless scattered data cloud of CFD.
x is the ndarray[37817,2].
F is the value of x ndarray[37817,].

when input the data in the expertsolver_example.py.
the error occurs:

/home/wjq/MeshlessDerivatives/main.py
Traceback (most recent call last):
File "/home/wjq/MeshlessDerivatives/main.py", line 183, in
main()
File "/home/wjq/MeshlessDerivatives/main.py", line 177, in main
X,Y,Z = project_onto_regular_grid_2D(x, F, fit_order=2, nk=50)
File "/home/wjq/MeshlessDerivatives/main.py", line 92, in project_onto_regular_grid_2D
solver.prepare( xi=x, xk=x[hoods] ) # generate problem matrices from the geometry of the point cloud
File "wlsqm/fitter/expert.pyx", line 399, in wlsqm.fitter.expert.ExpertSolver.prepare
ValueError: Buffer and memoryview are not contiguous in the same dimension.

Process finished with exit code 1

Thanks!