Add minimal Python OOP translation

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

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
+}

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
+}

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"

sies/__init__.py

@@ -0,0 +1,3 @@
+from .shape import C2Boundary
+from .acq import MConfig
+from .pde import SmallInclusions, ElectricFish

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]

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
+

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

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]))