Reorganize README and clarify multivariate case.

README.md

@@ -2,31 +2,34 @@ EasterEig
 =========
 [![License: GPL v3](https://img.shields.io/badge/License-GPLv3-blue.svg)](https://www.gnu.org/licenses/gpl-3.0) [![CI-Ubuntu](https://github.com/nennigb/EasterEig/actions/workflows/ci-ubuntu.yml/badge.svg)](https://github.com/nennigb/EasterEig/actions/workflows/ci-ubuntu.yml) [![pypi release](https://img.shields.io/pypi/v/eastereig.svg)](https://pypi.org/project/eastereig/)
 
-
-
-Consider a parametric eigenvalue problem depending on a parameter $$\nu$$. This arises for instance in
+Consider a **parametric eigenvalue problem** depending on one scalar $\nu$ or given vector $\boldsymbol\nu =(\nu_1,\nu_2,\ldots,\nu_N) \in \mathbb{C}^N$ of paramaters. This arises for instance in
 
   - waveguides, where the _wavenumber_ (eigenvalue) depends on the frequency (parameter)
-  - waveguides with absorbing materials on the wall, where _modal attenuation_ (eigenvalue imaginary part) depends on the liner properties like impedance, density (parameter)
-  - structural dynamics with a randomly varying parameter, where the _resonances frequencies_ (eigenvalue) depend on the parameter
+  - waveguides with absorbing materials on the wall, where _modal attenuation_ (eigenvalue imaginary part) depends on the liner properties like impedance, admittance, density (parameter)
+  - structural dynamics with a randomly varying parameter, where the _resonances frequencies_ (eigenvalue) depend on for instance of material parameters like Young modulus or density
   - ...
 
-Exceptional points (EP) of non-Hermitian systems correspond to particular values of the parameter leading to defective eigenvalue.
-At EP, both eigenvalues and eigenvectors are merging. 
+The aim of this package is to **reconstruct** the eigenvalue loci and to **locate** exceptional points (EPs). The EPs in non-Hermitian systems correspond to particular values of the parameters leading to defective eigenvalue. At EPs, both eigenvalues and eigenvectors are merging.
 
-The aim of this package is to **locate** exceptional points and to **reconstruct** the eigenvalue loci. The theoretical part of this work is described in [1], as for the location of _exceptional points_ and illustrated in [2] for eigenvalues reconstruction in structural dynamics.
+![Riemann surfaces around an EP2.](figures/EP2_light.png)  ![Riemann surfaces around an EP3.](figures/EP3_light.png)
 
-The method requires the computation of successive derivatives of two selected eigenvalues with respect to the parameter so that, after recombination, regular functions can be constructed. This algebraic manipulation enables
- * exceptional points (EP) localization, using standard root-finding algorithms; 
- * computation of the associated Puiseux series up to an arbitrary order.
+The theoretical parts of this work are described in [1] for the location of _exceptional points_ and in [2] for eigenvalues reconstruction. The extension to several parameters is presented in [3].
+
+The method requires the computation of successive derivatives of some selected eigenvalues with respect to the parameter so that, after recombination, regular functions can be constructed. This algebraic manipulation enables
+ * Exceptional points localization, using standard root-finding algorithms; 
+ * Computation of the associated Puiseux series up to an arbitrary order.
 
-This representation, which is associated with the topological structure of Riemann surfaces, allows to efficiently approximate the selected pair in a certain neighborhood of the EP.
+This representation, which is associated with the topological structure of Riemann surfaces, allows to efficiently approximate the selected eigenpair in a certain neighborhood of the EP.
 
 To use this package :
 
-  1. an access to the **operator derivative** must be possible
-  2. the eigenvalue problem must be recast into
-  		\[ \mathbf{L} (\lambda(\nu), \nu) \mathbf{x} (\nu) =\mathbf{0}  \]
+  1. An access to the **operator derivative** with respect to $\boldsymbol\nu$ is required
+  2. The parametric eigenvalue problem must have the form
+$$\mathbf{L} (\lambda(\boldsymbol\nu), \boldsymbol\nu) \mathbf{x} (\boldsymbol\nu) =\mathbf{0},$$
+where, for a given vector $\boldsymbol\nu$ which contains $N$ independent complex-valued parameters, $\lambda(\boldsymbol\nu)$ is an eigenvalue and $\mathbf{x}(\boldsymbol\nu)\neq \mathbf{0}$ is the associated right eigenvector. Here the matrix $\mathbf{L}$ admits the decomposition
+$$\mathbf{L} (\lambda, \boldsymbol\nu)  =\sum_{i \geq 0} f_i(\lambda) \mathbf{K}_i(\boldsymbol\nu)$$
+where $f_i$ is a polynomial function and matrices $\mathbf{K}_i$
+are supposed to be an analytic function of the parameters vector $\boldsymbol\nu$.
 
 The matrices of discrete operators can be either of numpy type for _full_, scipy type for _sparse_ or petsc mpiaij type for _sparse parallel_ matrices.
 
@@ -38,117 +41,72 @@ References
    [1] B. Nennig and E. Perrey-Debain. A high order continuation method to locate exceptional points and to compute Puiseux series with applications to acoustic waveguides. J. Comp. Phys., 109425, (2020). [[doi](https://dx.doi.org/10.1016/j.jcp.2020.109425)]; [[open access](https://arxiv.org/abs/1909.11579)]
 
    [2] M. Ghienne and B. Nennig. Beyond the limitations of perturbation methods for real random eigenvalue problems using Exceptional Points and analytic continuation. Journal of Sound and vibration, (2020). [[doi](https://doi.org/10.1016/j.jsv.2020.115398)]; [[open access](https://hal.archives-ouvertes.fr/hal-02536849)]
+
+   [3] B. Nennig, E. Perrey-Debain, Martin Ghienne. Localization of high order exceptional points : applications to acoustic waveguides. 16ème Congrès Français d’Acoustique, Apr 2022, Marseille, France.  [[open access](https://hal.science/hal-03838635v1/document)]
 
+Install 
+--------
 
-Basic workflow and class hierarchy
-----------------------------------
-
-`eastereig` provides several top level classes:
-
-  1. **OP class**, defines operators of your problem
-  2. **Eig class**, handles eigenvalues, their derivatives and reconstruction
-  3. **CharPol class**, combines several eigenvalues and their derivatives to reconstruction a part of the characteristic polynomial
-  4. **EP class**, combines Eig object to locate EP and compute Puiseux series
-  5. **Loci class**, stores numerical value of eigenvalues loci and allows easy Riemann surface plotting
-
-Dependencies
--------------
-
-`eastereig` is based on numpy (full) and scipy (sparse) for most internal computation and can handle _large_ parallel sparse matrices thanks to **optional** import of [petsc4py](https://petsc4py.readthedocs.io/en/stable/install.html) (and mumps), 
+`eastereig` is based on numpy (full) and scipy (sparse) for most internal computation and can handle _large_ parallel sparse matrices thanks to **optional** import of [petsc4py (>=3.20)](https://petsc.org/release/petsc4py/install.html) (and mumps), 
 [slepc4py](https://slepc4py.readthedocs.io/en/stable/install.html) and
 and [mpi4py](https://mpi4py.readthedocs.io/en/stable/install.html). As non-hermitian problems involve complex-valued eigenvalues, computations are realized with complex arithmetic and the **complex petsc version** is expected.
 [sympy](https://docs.sympy.org) is used for formal manipulation of multivariate polynomials.
 Riemann surface can also be plotted using the `Loci` class either with `matplotlib` or with [`pyvista`](https://github.com/pyvista/pyvista) and and `pyvistaqt` (optional).
 
-> **Remarks :**
-> To run an example with petsc (parallel), you need to run python with `mpirun` (or `mpiexec`). For instance, to run a program with 2 proc
-> `mpirun -n 2 python myprog.py`
-
-Install 
---------
+Before installing `eastereig`, you'll need python (tested for v >= 3.8), and to install manually optional dependencies you need:
+* The python packages `pypolsys` (optional homotopy EP solver) available from https://github.com/nennigb/pypolsys
+* The python packages `pyvista` and `pyvistaqt` (optional Riemann surfaces plotting)
+* The python package `petsc4py` and `slepc4py` (optional sparse parallel matrices surpport)
+* A fortran compiler (optional, tested with gfortran)
+**other dependencies will be installed automatically**.
 
-Before installing `eastereig`, you'll need: 
-* python (tested for v >= 3.5);
-* python packages: numpy, setuptools, wheel
-* pip (optional).
-* fortran compiler (optional)
-Note that on ubuntu, you will need to use `pip3` instead of `pip` and `python3` instead of `python`. Please see the steps given in the continous integration script [workflows](.github/workflows/ci-ubuntu.yml).
-
-
-By default, the fortan evaluation of multivariate polynomial is desactivated. To enable it, set the environnement variable: `EASTEREIG_USE_FPOLY=True`. On ubuntu like system, run
+By default, the fortan evaluation of multivariate polynomial is deactivated. To enable it, you need a fortran compiler and to set the environment variable: `EASTEREIG_USE_FPOLY=True`. On ubuntu like system, run
 ```console
 export EASTEREIG_USE_FPOLY=True
 ```
-
-### Using pip (preferred)
-Consider using `pip` over custom script (rationale [here](https://pip.pypa.io/en/stable/reference/pip_install/)). 
-
-You can install `eastereig` either from pypi (main releases only):
+Then, you can install `eastereig` either from pypi (main releases only):
 ```console
 pip install eastereig [--user]
 ``` 
-or from github:
-```console
-pip install path/to/EeasterEig-version.tar.gz [--user]
-```
-or in _editable_ mode if you want to modify the sources
-```console
-pip install -e path/to/EeasterEig
-```
-
-### Using python setuptools
-Go to root folder.
-and run:
-```console
-python setup.py install [--user]
-```
-
-To get the lastest updates (dev releases), run: 
+or from github (more updates):
 ```console
-python setup.py develop [--user]
+pip install path/to/EasterEig-version.tar.gz [--user]
 ```
+Note that on ubuntu, you will need to use `pip3` instead of `pip` and `python3` instead of `python`. Please see the steps given in the continous integration script [workflows](.github/workflows/ci-ubuntu.yml).
+We recommend installing `eastereig` in virtual environemment unless you know what you are doing.
 
 Running tests
 -------------
-Tests are handled with doctest. 
+Tests are handled with `doctest` and with `unittest`
 
 To execute the full test suite, run :
 ```console
 python -m eastereig
 ```
 
-Documentation
---------------
 
-## Generate documentation 
-Run: 
-```console
-pdoc3 --html --force --config latex_math=True  eastereig [--skip-errors]
-```
-N.B: The doctring are compatible with several Auto-generate API documentation, like pdoc.
-This notably allows to see latex includes.
+Basic workflow and class hierarchy
+----------------------------------
+`eastereig` provides several top level classes:
 
-## Generate class diagram
-Run: 
-```console
-pyreverse -s0 eastereig -m yes -f ALL
-dot -Tsvg classes.dot -o classes.svg
-```
-N.B: Class diagram generation is done using `pyreverse` (installed with pylint and spyder).
+  1. **OP class**, defines operators of your problem
+  2. **Eig class**, handles eigenvalues, their derivatives and reconstruction
+  3. **CharPol class**, combines several eigenvalues and their derivatives to reconstruction a part of the characteristic polynomial
+  4. **EP class**, combines Eig object to locate EP and compute Puiseux series
+  5. **Loci class**, stores numerical value of eigenvalues loci and allows easy Riemann surface plotting
 
-## Generate documentation and class diagram
-Both aspects are included in the `makedoc.py' script. So, just run :
-```console
-python ./makedoc.py
-```
 
 Getting started
 ---------------
-
 Several working examples are available in `./examples/` folder
 
   1. Acoustic waveguide with an impedance boundary condition (with the different supported linear libraries)
   2. 3-dof toy model of a structure with one random parameter (with numpy)
+  3. 3-dof toy with two parameters leading to EP3
+
+> **Remarks :**
+> To run an example with petsc (parallel), you need to run python with `mpirun` (or `mpiexec`). For instance, to run a program with 2 proc
+> `mpirun -n 2 python myprog.py`
 
 To get started, the first step is to define your problem. Basically it means to link the discrete operators (matrices) and their derivatives to the `eastereig` OP class.
 The problem has to be recast in the following form:
@@ -163,11 +121,11 @@ K = [K0, K1, K2].
 ```python
 dK = [dK0, dK1, dK2].
 ```
-Finally, **the functions** that returns derivatives with respect to $$\lambda$$ are stored in 'flda'
+Finally, **the functions** that returns derivatives with respect to $\lambda$ are stored in 'flda'
 ```python
 flda = [None, ee.lda_func.Lda, ee.lda_func.Lda2].
 ```
-Basic linear and quadratic dependency are defined in the module `lda_func`. Others dependencies can be easily implemented; provided that the appropriate eigenvalue solver is also implemented). The `None` keyword is used when there is no dependency to the eigenvalue, e. g. $$\mathbf{K}_0$$.
+Basic linear and quadratic dependency are defined in the module `lda_func`. Others dependencies can be easily implemented; provided that the appropriate eigenvalue solver is also implemented). The `None` keyword is used when there is no dependency to the eigenvalue, e. g. $\mathbf{K}_0$.
 
 This formulation is used to automatically compute (i) the successive derivatives of the operator and (ii) the RHS associated to the bordered matrix.
 
@@ -177,7 +135,7 @@ These variables are defined by creating **a subclass** that inherits from the ea
 import eastereig as ee
 
 class MyOP(ee.OP):
-    """Create a subclass of the OP class to describe your problem."""
+    """Create a subclass of the OP class to describe your problem with scalar parameter."""
 
     def __init__(self):
         """Initialize the problem."""
@@ -220,6 +178,31 @@ class MyOP(ee.OP):
 		return dM1
 ```
 
+Documentation
+--------------
+
+## Generate documentation 
+Run: 
+```console
+pdoc3 --html --force --config latex_math=True  eastereig [--skip-errors]
+```
+N.B: The doctring are compatible with several Auto-generate API documentation, like pdoc.
+This notably allows to see latex includes.
+
+## Generate class diagram
+Run: 
+```console
+pyreverse -s0 eastereig -m yes -f ALL
+dot -Tsvg classes.dot -o classes.svg
+```
+N.B: Class diagram generation is done using `pyreverse` (installed with pylint and spyder).
+
+## Generate documentation and class diagram
+Both aspects are included in the `makedoc.py' script. So, just run :
+```console
+python ./makedoc.py
+```
+
 How to contribute ?
 -------------------