diff --git a/README.md b/README.md index cb2fee6..22d47d1 100644 --- a/README.md +++ b/README.md @@ -29,3 +29,13 @@ The library is written in Matlab object-oriented programming language and is org For usage of the package, go to the path examples/ and run the demo scripts. Your suggestions and remarks are welcome. Author: Han Wang, Email: han.wang@ens.fr +## Python package + +Une traduction minimale en Python se trouve dans le dossier `sies`. Vous pouvez installer ce petit paquet en mode développement avec `uv`. + +```bash +uv pip install -e . +``` + +Puis lancez les notebooks dans `examples/` pour voir quelques démos. + diff --git a/examples/eit_circle_demo.ipynb b/examples/eit_circle_demo.ipynb new file mode 100644 index 0000000..42f5deb --- /dev/null +++ b/examples/eit_circle_demo.ipynb @@ -0,0 +1,12 @@ +{ + "cells": [ + {"cell_type": "markdown", "metadata": {}, "source": ["# EIT sur un cercle"]}, + {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["import sys, os\n", "sys.path.append(os.path.abspath('..'))\n"]}, + {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["import numpy as np\n", "import matplotlib.pyplot as plt\n", "from sies import C2Boundary, MConfig, ElectricFish\n", "theta = np.linspace(0, 2*np.pi, 16, endpoint=False)\n", "src_rcv = np.vstack([np.cos(theta), np.sin(theta)]) * 2.0\n", "obj_theta = np.linspace(0, 2*np.pi, 50, endpoint=False)\n", "obj = C2Boundary(np.vstack([0.5*np.cos(obj_theta), 0.5*np.sin(obj_theta)]))\n", "cfg = MConfig(sources=src_rcv, receivers=src_rcv)\n", "cfg.Omega0 = C2Boundary(src_rcv)\n", "fish = ElectricFish(obj, [1.0], [1.0], cfg)\n", "plt.figure();\n", "plt.plot(obj.points[0], obj.points[1], label='objet');\n", "plt.scatter(src_rcv[0], src_rcv[1], c='r', label='capteurs');\n", "plt.axis('equal');\n", "plt.legend();\n"]}, + {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["def potential(x, source):\n", " return -np.log(np.linalg.norm(x - source[:, None], axis=0)) / (2*np.pi)\n", "grid_x, grid_y = np.mgrid[-1:1:100j, -1:1:100j]\n", "phi = potential(np.vstack([grid_x.ravel(), grid_y.ravel()]), src_rcv[:,0]).reshape(grid_x.shape)\n", "plt.figure();\n", "plt.pcolormesh(grid_x, grid_y, phi, shading='auto');\n", "plt.plot(obj.points[0], obj.points[1], 'k');\n", "plt.axis('equal');\n"]}, + {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["G = fish.grammatrix\n", "G"]} + ], + "metadata": {"kernelspec": {"display_name": "Python 3", "language": "python", "name": "python3"}, "language_info": {"name": "python", "version": "3.10"}}, + "nbformat": 4, + "nbformat_minor": 4 +} diff --git a/examples/electric_fish_demo.ipynb b/examples/electric_fish_demo.ipynb new file mode 100644 index 0000000..b8e9fe6 --- /dev/null +++ b/examples/electric_fish_demo.ipynb @@ -0,0 +1,51 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# ElectricFish demo" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import sys, os\n", + "sys.path.append(os.path.abspath(\"..\"))\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "from sies import C2Boundary, MConfig, ElectricFish\n", + "theta = np.linspace(0, 2 * np.pi, 10, endpoint=False)\n", + "points = np.vstack([np.cos(theta), np.sin(theta)])\n", + "boundary = C2Boundary(points)\n", + "cfg = MConfig(sources=np.array([[0.0], [0.0]]), receivers=points)\n", + "cfg.Omega0 = boundary\n", + "fish = ElectricFish(boundary, [1.0], [1.0], cfg)\n", + "fish.grammatrix" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "name": "python", + "version": "3.10" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} diff --git a/pyproject.toml b/pyproject.toml new file mode 100644 index 0000000..d35ab3b --- /dev/null +++ b/pyproject.toml @@ -0,0 +1,12 @@ +[project] +name = "sies" +version = "0.0.1" +description = "Python utilities for shape identification in electro-sensing" +authors = [{name = "OpenAI", email = "codex@openai.com"}] +readme = "README.md" +requires-python = ">=3.8" +dependencies = ["numpy"] + +[build-system] +requires = ["setuptools>=61"] +build-backend = "setuptools.build_meta" diff --git a/sies/__init__.py b/sies/__init__.py new file mode 100644 index 0000000..d4f5759 --- /dev/null +++ b/sies/__init__.py @@ -0,0 +1,3 @@ +from .shape import C2Boundary +from .acq import MConfig +from .pde import SmallInclusions, ElectricFish diff --git a/sies/acq.py b/sies/acq.py new file mode 100644 index 0000000..007b5df --- /dev/null +++ b/sies/acq.py @@ -0,0 +1,52 @@ +import numpy as np + +class MConfig: + """Acquisition configuration with sources and receivers. + + Parameters + ---------- + sources : array_like + Source coordinates of shape ``(2, Ns)``. + receivers : array_like + Receiver coordinates of shape ``(2, Nr)``. + center : array_like, optional + Reference center of the system. + """ + + def __init__(self, sources, receivers, center=None): + self.src_prv = [np.asarray(sources, dtype=float)] + self.rcv_prv = [np.asarray(receivers, dtype=float)] + self.Ng = 1 + if center is None: + self.center = np.mean(np.hstack([self.src_prv[0], self.rcv_prv[0]]), axis=1) + else: + self.center = np.asarray(center, dtype=float) + + @property + def Ns(self): + return self.src_prv[0].shape[1] + + @property + def Nr(self): + return self.rcv_prv[0].shape[1] + + @property + def Ns_total(self): + return self.Ns * self.Ng + + @property + def Nr_total(self): + return self.Nr * self.Ng + + @property + def data_dim(self): + return self.Ns * self.Ng * self.Nr + + def group(self, g): + return self.src_prv[g], self.rcv_prv[g] + + def src(self, sidx): + return self.src_prv[0][:, sidx - 1] + + def rcv(self, s): + return self.rcv_prv[0] diff --git a/sies/pde.py b/sies/pde.py new file mode 100644 index 0000000..64bca58 --- /dev/null +++ b/sies/pde.py @@ -0,0 +1,44 @@ +import numpy as np +from .shape import C2Boundary +from .acq import MConfig + +class SmallInclusions: + """Base class for PDEs with small inclusions.""" + + def __init__(self, inclusions, cfg): + if not isinstance(cfg, MConfig): + raise TypeError("cfg must be an instance of MConfig") + self.cfg = cfg + self.D = [] + self.nb_incls = 0 + if isinstance(inclusions, C2Boundary): + self.add_inclusion(inclusions) + else: + for inc in inclusions: + self.add_inclusion(inc) + + def add_inclusion(self, inc): + if not isinstance(inc, C2Boundary): + raise TypeError("inclusion must be C2Boundary") + self.D.append(inc) + self.nb_incls += 1 + + +class ElectricFish(SmallInclusions): + """Simplified electric fish model.""" + + def __init__(self, D, cnd, pmtt, cfg, step_bem=1): + super().__init__(D, cfg) + self.Omega = cfg.Omega0 if hasattr(cfg, "Omega0") else cfg.rcv_prv[0] + self.impd = getattr(cfg, "impd", 0.0) + self.cnd = np.asarray(cnd, dtype=float) + self.pmtt = np.asarray(pmtt, dtype=float) + self.step_bem1 = step_bem + self.step_bem2 = 1 + self.Psi = np.eye(self.Omega.nb_points)[::step_bem] + + @property + def grammatrix(self): + sigma = np.ones(self.Psi.shape[0]) + return self.Psi.T @ np.diag(sigma) @ self.Psi + diff --git a/sies/shape.py b/sies/shape.py new file mode 100644 index 0000000..a004b2e --- /dev/null +++ b/sies/shape.py @@ -0,0 +1,47 @@ +import numpy as np + +class C2Boundary: + """Simple C2 boundary representation. + + Parameters + ---------- + points : array_like + Array of shape ``(2, N)`` representing boundary points. + name : str, optional + Name of the boundary. + """ + + def __init__(self, points, name=""): + self.points = np.asarray(points, dtype=float) + if self.points.shape[0] != 2: + raise ValueError("Points must have shape (2, N)") + self.name = name + self.center_of_mass = self.points.mean(axis=1) + self.nb_points = self.points.shape[1] + self._compute_geometry() + + def _compute_geometry(self): + diff = np.roll(self.points, -1, axis=1) - self.points + norm = np.linalg.norm(diff, axis=0) + norm[norm == 0] = 1.0 + self.tvec = diff + self.normal = np.vstack((-diff[1], diff[0])) / norm + self.tvec_norm = norm + self.sigma = 2 * np.pi / self.nb_points * norm + + @property + def box(self): + d = self.points - self.center_of_mass[:, None] + w = d[0].max() - d[0].min() + h = d[1].max() - d[1].min() + return np.array([w, h]) + + @property + def diameter(self): + d = self.points - self.center_of_mass[:, None] + return 2 * np.sqrt((d ** 2).sum(axis=0)).max() + + def isinside(self, x): + x = np.asarray(x).reshape(2) + r = np.linalg.norm(x - self.center_of_mass) + return r < self.diameter / 2 diff --git a/tests/test_electric_fish.py b/tests/test_electric_fish.py new file mode 100644 index 0000000..210a457 --- /dev/null +++ b/tests/test_electric_fish.py @@ -0,0 +1,24 @@ +import os, sys +sys.path.insert(0, os.path.abspath(os.path.join(os.path.dirname(__file__), ".."))) +import numpy as np +from sies import C2Boundary, MConfig, ElectricFish + + +def test_boundary_properties(): + theta = np.linspace(0, 2 * np.pi, 20, endpoint=False) + points = np.vstack([np.cos(theta), np.sin(theta)]) + boundary = C2Boundary(points) + assert boundary.nb_points == 20 + assert np.isclose(boundary.diameter, 2, atol=1e-2) + + +def test_electric_fish_grammatrix(): + theta = np.linspace(0, 2 * np.pi, 10, endpoint=False) + points = np.vstack([np.cos(theta), np.sin(theta)]) + boundary = C2Boundary(points) + cfg = MConfig(sources=np.array([[0.0], [0.0]]), receivers=points) + cfg.Omega0 = boundary + fish = ElectricFish(boundary, [1.0], [1.0], cfg) + G = fish.grammatrix + assert G.shape[0] == fish.Psi.shape[1] + assert np.allclose(G, np.eye(G.shape[0]))