diff --git a/examples/conductivity_demo.py b/examples/conductivity_demo.py new file mode 100644 index 0000000..52e3426 --- /dev/null +++ b/examples/conductivity_demo.py @@ -0,0 +1,95 @@ +import marimo + +__generated_with = "0.1.0" +app = marimo.App() + + +@app.cell +def __(): + import marimo as mo + import numpy as np + import matplotlib.pyplot as plt + import sies.shape as shape + import sies.acq as acq + import sies.pde as pde + import sies.asymp as asymp + return acq, asymp, mo, np, pde, plt, shape + + +@app.cell +def __(np, shape): + # Definition of small inclusions + # B = shape.Triangle(0.5, np.pi/3, 2**10, 10) + B = shape.Ellipse(0.5, 0.25, 2**10) # Using ellipse for easier visualization first + + # Multiple inclusions + D = [B + [0.1, 0.1]] + cnd = [10.0] + pmtt = [1.0] + return B, D, cnd, pmtt + + +@app.cell +def __(D, acq, np, pde, plt): + # Set up an environment for experience + # Neutrality: surprisingly, this has a better conditioning + cfg = acq.Coincided([0, 0], 10, 50, viewmode=(1, np.pi/16, 2*np.pi), grouped=False, neutCoeff=[1, -1], neutRad=0.01) + + P = pde.Conductivity_R2(D, [10.0], [1.0], cfg) + + plt.figure(figsize=(6, 6)) + P.plot() + plt.axis('equal') + plt.title("Inclusion and Acquisition System") + plt.gca() + return P, cfg + + +@app.cell +def __(P, np): + # Simulation of the MSR data + freqlist = np.linspace(0, 100 * np.pi, 5) + data = P.data_simulation(freqlist) + return data, freqlist + + +@app.cell +def __(D, asymp, cnd, freqlist, pmtt): + # Compute first the theoretical value of CGPT + ord_val = 1 + M_theo = [] + for _f in freqlist: + _lamb = asymp.lambda_contrast(cnd, pmtt, _f) + M_theo.append(asymp.theoretical_CGPT(D, _lamb, ord_val)) + return M_theo, ord_val + + +@app.cell +def __(M_theo, P, data, freqlist, np, ord_val): + # Reconstruct CGPT and show error + nlvl = 0.01 # Add some noise + data_noisy = P.add_white_noise(data, nlvl) + + K = max(1, ord_val) + out = P.reconstruct_CGPT(data_noisy["MSR_noisy"], K, maxiter=100000, tol=1e-10, symmode=True, method='lsqr') + + print("Relative error between theoretical and reconstructed CGPT matrix at different frequencies:") + for _f_idx, _f in enumerate(freqlist): + _toto = out["CGPT"][_f_idx][:2*ord_val, :2*ord_val] + _theo = M_theo[_f_idx] + + _err = np.linalg.norm(_theo - _toto, 'fro') / np.linalg.norm(_theo, 'fro') + + # SVD error (as in Matlab demo) + _u0, _s0, _vh0 = np.linalg.svd(_theo) + _u1, _s1, _vh1 = np.linalg.svd(_toto) + _sv0 = _s0[0] / _s0[1] + _sv1 = _s1[0] / _s1[1] + _errsvd = np.abs(_sv0 - _sv1) / np.abs(_sv1) + + print(f"Frequency: {_f:.2f}, error: {_err:.4f}, error of sv: {_errsvd:.4f}") + return K, data_noisy, out + + +if __name__ == "__main__": + app.run() diff --git a/pyproject.toml b/pyproject.toml new file mode 100644 index 0000000..c1acfea --- /dev/null +++ b/pyproject.toml @@ -0,0 +1,20 @@ +[project] +name = "sies" +version = "0.1.0" +description = "Shape identification in electro-sensing (Python port)" +readme = "README.md" +requires-python = ">=3.10" +dependencies = [ + "numpy", + "scipy", + "matplotlib", + "marimo", +] + +[build-system] +requires = ["setuptools>=61"] +build-backend = "setuptools.build_meta" + +[tool.setuptools.packages.find] +where = ["."] +include = ["sies*"] diff --git a/sies/__init__.py b/sies/__init__.py new file mode 100644 index 0000000..9f88688 --- /dev/null +++ b/sies/__init__.py @@ -0,0 +1,4 @@ +from .shape import Shape, Triangle, Ellipse, Flower +from .acq import AcquisitionConfig, Concentric, Coincided +from .pde import SmallInclusions, Conductivity_R2 +from .asymp import lambda_contrast, theoretical_CGPT diff --git a/sies/acq/__init__.py b/sies/acq/__init__.py new file mode 100644 index 0000000..65c4d7a --- /dev/null +++ b/sies/acq/__init__.py @@ -0,0 +1,270 @@ +import numpy as np +import matplotlib.pyplot as plt + +class AcquisitionConfig: + """ + Abstract class for the configuration of acquisition system. + + An acquisition system consists of sources and receivers. This class + contains only the geometrical properties, such as the position of + sources and receivers, but no physical properties like frequency. + + Parameters + ---------- + src_prv : list of ndarray + Coordinates of sources by group. + rcv_prv : list of ndarray + Coordinates of receivers by group. + center : array_like, optional + Reference center of the measurement system. + """ + + def __init__(self, src_prv, rcv_prv, center=None): + self.src_prv = [np.asarray(s) for s in src_prv] + self.rcv_prv = [np.asarray(r) for r in rcv_prv] + self.Ng = len(self.src_prv) + self.center = np.asarray(center).flatten() if center is not None else np.array([0, 0]) + + self.Ns = self.src_prv[0].shape[1] if self.Ng > 0 else 0 + self.Nr = self.rcv_prv[0].shape[1] if self.Ng > 0 else 0 + + self.Ns_total = self.Ns * self.Ng + self.Nr_total = self.Nr * self.Ng + self.data_dim = self.Ns * self.Ng * self.Nr + + def group(self, g): + """ + Get coordinates of sources and receivers for a specific group. + + Parameters + ---------- + g : int + Group index (0-indexed). + + Returns + ------- + tuple + (src, rcv) coordinates for the group. + """ + return self.src_prv[g], self.rcv_prv[g] + + def src_query(self, s): + """ + For a global source index, get group index and local index. + + Parameters + ---------- + s : int + Global source index (0-indexed). + + Returns + ------- + tuple + (gid, sid) group index and local source index. + """ + if s >= self.Ns_total or s < 0: + raise IndexError("Source index out of range") + gid = s // self.Ns + sid = s % self.Ns + return gid, sid + + def src(self, sidx=None): + """ + Get coordinates of specific sources. + + Parameters + ---------- + sidx : int or list of int, optional + Index or indices of sources. If None, all sources are returned. + + Returns + ------- + ndarray + Coordinates of requested sources. + """ + if sidx is None: + sidx = np.arange(self.Ns_total) + + indices = np.atleast_1d(sidx) + val = np.zeros((2, len(indices))) + for i, s in enumerate(indices): + gid, sid = self.src_query(s) + val[:, i] = self.src_prv[gid][:, sid] + + return val if not isinstance(sidx, (int, np.integer)) else val[:, 0] + + def rcv(self, s): + """ + Get coordinates of receivers responding to a specific source. + + Parameters + ---------- + s : int + Global source index (0-indexed). + + Returns + ------- + ndarray + Coordinates of receivers in the same group as source s. + """ + gid, _ = self.src_query(s) + return self.rcv_prv[gid] + + @property + def all_src(self): + """ndarray : Coordinates of all sources concatenated.""" + return np.hstack(self.src_prv) + + @property + def all_rcv(self): + """ndarray : Coordinates of all receivers concatenated.""" + return np.hstack(self.rcv_prv) + + def plot(self, *args, **kwargs): + """ + Plot sources and receivers. + """ + for g in range(self.Ng): + src, rcv = self.group(g) + plt.plot(src[0, :], src[1, :], 'x', label=f'Group {g} sources') + plt.plot(rcv[0, :], rcv[1, :], 'o', label=f'Group {g} receivers') + plt.plot(self.center[0], self.center[1], 'r*') + plt.legend() + + +def src_rcv_circle(Na, N0, R0, Z, theta, aov=2*np.pi): + """ + Generate sources/receivers placed on concentric arcs. + + Parameters + ---------- + Na : int + Number of arcs. + N0 : int + Number of sources/receivers per arc. + R0 : float + Radius of measurement circle. + Z : array_like + Center of measurement circle. + theta : float + Angular aperture of each arc. + aov : float, optional + Total angle of view coverage. + + Returns + ------- + tuple + (Xs, Thetas, Xscell) coordinates, angles, and grouped list of points. + """ + import numpy as np + Xs = np.zeros((2, N0 * Na)) + Thetas = np.zeros(N0 * Na) + Xscell = [] + + Z = np.asarray(Z).flatten() + + for n in range(Na): + tt0 = n / Na * aov + tt = tt0 + np.arange(N0) / N0 * theta + + rr = R0 * np.vstack([np.cos(tt), np.sin(tt)]) + start_idx = n * N0 + end_idx = (n + 1) * N0 + Xs[:, start_idx:end_idx] = rr + Thetas[start_idx:end_idx] = tt + Xscell.append(rr + Z[:, None]) + + Xs = Xs + Z[:, None] + return Xs, Thetas, Xscell + +class Concentric(AcquisitionConfig): + """ + Concentric configuration for sources and receivers. + + Parameters + ---------- + Z : array_like + Center of the measurement circle. + Rs : float + Radius of source circle. + Ns : int + Number of sources per arc. + Rr : float + Radius of receiver circle. + Nr : int + Number of receivers per arc. + viewmode : tuple, optional + (Na, theta, aov) specifying number of arcs, aperture, and total coverage. + grouped : bool, optional + If True, each arc is a separate group. + neutCoeff : list or ndarray, optional + Coefficients for neutral source condition. + neutRad : float, optional + Relative distance between Diracs for neutral source. + """ + + def __init__(self, Z, Rs, Ns, Rr, Nr, viewmode=(1, 2*np.pi, 2*np.pi), grouped=False, neutCoeff=None, neutRad=0.01): + Na, theta, aov = viewmode + + Xs, _, Xscell = src_rcv_circle(Na, Ns, Rs, Z, theta, aov) + Xr, _, Xrcell = src_rcv_circle(Na, Nr, Rr, Z, theta, aov) + + if grouped: + src_prv = Xscell + rcv_prv = Xrcell + else: + src_prv = [Xs] + rcv_prv = [Xr] + + super().__init__(src_prv, rcv_prv, center=Z) + + self.radius_src = Rs + self.radius_rcv = Rr + self.neutRad = neutRad + + if neutCoeff is None or len(np.atleast_1d(neutCoeff)) <= 1: + self.neutCoeff = np.array([1.0]) + else: + self.neutCoeff = np.asarray(neutCoeff) + if not np.isclose(np.sum(self.neutCoeff), 0) or np.all(self.neutCoeff == 0): + raise ValueError("Coefficients of Diracs must satisfy neutrality condition (sum=0) and be non-zero!") + + @property + def nbDirac(self): + """int : Number of Diracs for the neutral source.""" + return len(self.neutCoeff) + + def neutSrc(self, s): + """ + Get the positions of Diracs for the s-th source. + + Parameters + ---------- + s : int + Global source index (0-indexed). + + Returns + ------- + ndarray + Positions of Diracs (2 x nbDirac). + """ + psrc = self.src(s) + if self.nbDirac == 1: + return psrc[:, None] + else: + val = np.zeros((2, self.nbDirac)) + L = self.radius_src * self.neutRad + toto = psrc - self.center + q = np.array([toto[1], -toto[0]]) # tangent direction + q = q / np.linalg.norm(q) * L + + for n in range(self.nbDirac): + val[:, n] = psrc + (n / self.nbDirac) * q + return val + +class Coincided(Concentric): + """ + Coincided sources and receivers on a circle. + """ + def __init__(self, Z, Rs, Ns, viewmode=(1, 2*np.pi, 2*np.pi), grouped=False, neutCoeff=None, neutRad=0.01): + super().__init__(Z, Rs, Ns, Rs, Ns, viewmode, grouped, neutCoeff, neutRad) diff --git a/sies/asymp/__init__.py b/sies/asymp/__init__.py new file mode 100644 index 0000000..df4c496 --- /dev/null +++ b/sies/asymp/__init__.py @@ -0,0 +1,158 @@ +import numpy as np +import scipy.linalg +from ..tools.laplacian import Green2D_Grad + +def lambda_contrast(cnd, pmtt=None, freq=0): + """ + Compute the contrast lambda from the conductivity, the permittivity and the frequency. + """ + cnd = np.atleast_1d(cnd) + if pmtt is None: + pmtt = np.zeros_like(cnd) + else: + pmtt = np.atleast_1d(pmtt) + + if freq < 0: + raise ValueError("Frequency must be a positive scalar.") + + for c in cnd: + if np.isclose(c, 1.0) or c < 0: + raise ValueError("Invalid value of conductivity.") + + # Using convention (1) from Matlab code: f^(w) = \int f(x) exp(-2*pi*1i*x*w) dx + toto = cnd + 2 * np.pi * 1j * pmtt * freq + return (toto + 1) / (toto - 1) / 2 + +def make_kernel_matrix_Kstar(points, tvec, normal, avec, sigma): + M = points.shape[1] + Ks = np.zeros((M, M)) + tvec_norm_square = np.sum(tvec**2, axis=0) + + for j in range(M): + # xdoty = + diff = points[:, j, None] - points + xdoty = diff[0, :] * normal[0, j] + diff[1, :] * normal[1, j] + norm_xy_square = np.sum(diff**2, axis=0) + + # Avoid division by zero at diagonal + mask = np.ones(M, dtype=bool) + mask[j] = False + + Ks[j, mask] = (1 / (2 * np.pi)) * xdoty[mask] * (sigma[mask] / norm_xy_square[mask]) + + # Diagonal term + # avec[:,j] is (2,), normal[:,j] is (2,) + diag_val = (1 / (2 * np.pi)) * (-0.5) * np.dot(avec[:, j], normal[:, j]) / tvec_norm_square[j] * sigma[j] + Ks[j, j] = diag_val + + return Ks + +def make_kernel_matrix_dSLdn(D_points, D_sigma, E_points, E_normal): + Gx, Gy = Green2D_Grad(E_points, D_points) # (N_E, N_D) + # K = (diag(normal_E(1,:))*Gx + diag(normal_E(2,:))*Gy) * diag(sigma_D) + K = (E_normal[0, :, None] * Gx + E_normal[1, :, None] * Gy) * D_sigma[None, :] + return K + +def make_block_matrix(inclusions): + nb_incls = len(inclusions) + KsdS = [[None for _ in range(nb_incls)] for _ in range(nb_incls)] + + for m in range(nb_incls): + for n in range(nb_incls): + if m == n: + KsdS[m][n] = -1.0 * make_kernel_matrix_Kstar( + inclusions[n].points, inclusions[n].tvec, inclusions[n].normal, + inclusions[n].avec, inclusions[n].sigma + ) + else: + KsdS[m][n] = -1.0 * make_kernel_matrix_dSLdn( + inclusions[n].points, inclusions[n].sigma, + inclusions[m].points, inclusions[m].normal + ) + return KsdS + +def make_system_matrix_fast(KsdS, lambda_val): + nb_incls = len(KsdS) + lambda_val = np.atleast_1d(lambda_val) + + Acell = [[None for _ in range(nb_incls)] for _ in range(nb_incls)] + for m in range(nb_incls): + for n in range(nb_incls): + if m == n: + Acell[m][n] = lambda_val[m] * np.eye(KsdS[m][n].shape[0]) + KsdS[m][n] + else: + Acell[m][n] = KsdS[m][n] + + return np.block(Acell) + +def cell2mat(CC, CS, SC, SS): + ord_val = CC.shape[0] + M = np.zeros((2 * ord_val, 2 * ord_val)) + M[0::2, 0::2] = CC + M[0::2, 1::2] = CS + M[1::2, 0::2] = SC + M[1::2, 1::2] = SS + return M + +def theoretical_CGPT(inclusions, lambda_val, ord_val): + KsdS = make_block_matrix(inclusions) + return theoretical_CGPT_fast(inclusions, KsdS, lambda_val, ord_val) + +def theoretical_CGPT_fast(inclusions, KsdS, lambda_val, ord_val): + nb_points = inclusions[0].nb_points + nb_incls = len(inclusions) + lambda_val = np.atleast_1d(lambda_val) + + Amat = make_system_matrix_fast(KsdS, lambda_val) + + epsilon = 1e-8 + if np.min(np.abs(lambda_val - 0.5)) < epsilon: + # Extra condition of L^2_0 function + # [Amat; kron(eye(nbIncls), ones(1, size(Amat0,2)/nbIncls))] + extra = np.zeros((nb_incls, Amat.shape[1])) + for i in range(nb_incls): + extra[i, i*nb_points:(i+1)*nb_points] = 1.0 + Amat = np.vstack([Amat, extra]) + + CC = np.zeros((ord_val, ord_val)) + CS = np.zeros((ord_val, ord_val)) + SC = np.zeros((ord_val, ord_val)) + SS = np.zeros((ord_val, ord_val)) + + for m in range(1, ord_val + 1): + B = np.zeros((nb_points, nb_incls), dtype=complex) + for i in range(nb_incls): + # cpoints = points[0] + 1j * points[1] + cpoints = inclusions[i].points[0, :] + 1j * inclusions[i].points[1, :] + dm = m * (cpoints**(m-1)) + toto = inclusions[i].normal[0, :] * dm + inclusions[i].normal[1, :] * dm * 1j + B[:, i] = toto + + b = B.flatten(order='F') # Column-major flattening + if np.min(np.abs(lambda_val - 0.5)) < epsilon: + b_real = np.append(b.real, np.zeros(nb_incls)) + b_imag = np.append(b.imag, np.zeros(nb_incls)) + else: + b_real = b.real + b_imag = b.imag + + # Solve Amat \ b + phi_m_real_vec = np.linalg.solve(Amat, b_real) if Amat.shape[0] == Amat.shape[1] else np.linalg.lstsq(Amat, b_real, rcond=None)[0] + phi_m_imag_vec = np.linalg.solve(Amat, b_imag) if Amat.shape[0] == Amat.shape[1] else np.linalg.lstsq(Amat, b_imag, rcond=None)[0] + + # reshape back to (nb_points, nb_incls) + # We only take the first nb_points * nb_incls elements if we added extra rows + realphim = phi_m_real_vec[:nb_points*nb_incls].reshape((nb_points, nb_incls), order='F') + imagphim = phi_m_imag_vec[:nb_points*nb_incls].reshape((nb_points, nb_incls), order='F') + + for n in range(1, ord_val + 1): + for i in range(nb_incls): + cpoints = inclusions[i].points[0, :] + 1j * inclusions[i].points[1, :] + zn = (cpoints**n) * inclusions[i].sigma + + CC[m-1, n-1] += np.sum(zn.real * realphim[:, i].real) + CS[m-1, n-1] += np.sum(zn.imag * realphim[:, i].real) + SC[m-1, n-1] += np.sum(zn.real * imagphim[:, i].real) + SS[m-1, n-1] += np.sum(zn.imag * imagphim[:, i].real) + + return cell2mat(CC, CS, SC, SS) diff --git a/sies/pde/__init__.py b/sies/pde/__init__.py new file mode 100644 index 0000000..b6a40a1 --- /dev/null +++ b/sies/pde/__init__.py @@ -0,0 +1,2 @@ +from .base import SmallInclusions, Green2D, SingleLayer_eval +from .conductivity import Conductivity_R2 diff --git a/sies/pde/base.py b/sies/pde/base.py new file mode 100644 index 0000000..70acc46 --- /dev/null +++ b/sies/pde/base.py @@ -0,0 +1,64 @@ +import numpy as np +from ..tools.laplacian import Green2D_Dn, Green2D_Grad + +class SmallInclusions: + """Abstract class for several PDEs involving small inclusions in a homogeneous medium.""" + + def __init__(self, inclusions, cfg): + if not isinstance(inclusions, list): + inclusions = [inclusions] + + self.D = [] + self.nb_incls = 0 + self.cfg = cfg + + for d in inclusions: + self.add_inclusion(d) + + def add_inclusion(self, inclusion): + if self.nb_incls >= 1: + if self.D[self.nb_incls - 1].nb_points != inclusion.nb_points: + raise ValueError("Number of boundary discretization points must be the same for all inclusions.") + if not self.check_inclusions(inclusion): + raise ValueError("Inclusions must be separated from each other.") + + self.D.append(inclusion) + self.nb_incls += 1 + + def check_inclusions(self, inclusion): + for d in self.D: + if not inclusion.isdisjoint(d): + return False + return True + + def plot(self, *args, **kwargs): + import matplotlib.pyplot as plt + for d in self.D: + d.plot(*args, **kwargs) + self.cfg.plot(*args, **kwargs) + + def data_simulation(self, *args, **kwargs): + raise NotImplementedError("data_simulation must be implemented by subclasses") + +def Green2D(X, Y): + """ + 2D Green function G(x) = 1/2/pi * log|x| evaluated at the points X-Y. + """ + X = np.asarray(X) + Y = np.asarray(Y) + XY1 = X[0, :, None] - Y[0, None, :] + XY2 = X[1, :, None] - Y[1, None, :] + DN = XY1**2 + XY2**2 + # G = 1/2/pi * log|x| = 1/4/pi * log|x|^2 + G = (1 / (4 * np.pi)) * np.log(DN) + return G + +def SingleLayer_eval(D, F, X): + """ + Evaluate the single layer potential S_D[F](X) + """ + # G is (N_D, N_X) if we use Green2D(D.points, X) + G = Green2D(D.points, X) + # val = (F * sigma) * G + val = (F.flatten() * D.sigma) @ G + return val diff --git a/sies/pde/conductivity.py b/sies/pde/conductivity.py new file mode 100644 index 0000000..bdef671 --- /dev/null +++ b/sies/pde/conductivity.py @@ -0,0 +1,246 @@ +import numpy as np +from scipy.sparse.linalg import lsqr as sparse_lsqr +from .base import SmallInclusions, SingleLayer_eval +from ..tools.laplacian import Green2D_Dn +from ..asymp import lambda_contrast, make_block_matrix, make_system_matrix_fast +from ..tools import add_white_noise + +class Conductivity_R2(SmallInclusions): + r""" + Class for the conductivity problem in free space. + + Solves: div(rho_D * grad(u)) = 0 in R^2 \ {x_s}. + + Parameters + ---------- + inclusions : list of Shape + List of inclusions. + cnd : list or ndarray + Conductivity constants for each inclusion. + pmtt : list or ndarray + Permittivity constants for each inclusion. + cfg : AcquisitionConfig + Acquisition configuration. + """ + + def __init__(self, inclusions, cnd, pmtt, cfg): + super().__init__(inclusions, cfg) + + self.cnd = np.atleast_1d(cnd) + self.pmtt = np.atleast_1d(pmtt) + + if len(self.cnd) < self.nb_incls or len(self.pmtt) < self.nb_incls: + raise ValueError("Conductivity and permittivity must be specified for each inclusion!") + + for c in self.cnd: + if np.isclose(c, 1.0) or c < 0: + raise ValueError("Conductivity must be positive and different from 1!") + for p in self.pmtt: + if p < 0: + raise ValueError("Permittivity must be positive!") + + self.KsdS = make_block_matrix(self.D) + self.dGdn = self._compute_dGdn() + + def _compute_dGdn(self, sidx=None): + """Construct the right hand vector of the given source(s).""" + if sidx is None: + sidx = np.arange(self.cfg.Ns_total) + else: + sidx = np.atleast_1d(sidx) + + nb_points = self.D[0].nb_points + val = np.zeros((nb_points, self.nb_incls, len(sidx))) + + src = self.cfg.src(sidx) + for i in range(self.nb_incls): + toto = Green2D_Dn(src, self.D[i].points, self.D[i].normal) + val[:, i, :] = toto.T + + return val.reshape(nb_points * self.nb_incls, len(sidx), order='F') + + def _compute_phi(self, freq, sidx=None): + """Compute the boundary function phi given frequency and source(s).""" + nb_points = self.D[0].nb_points + lamb = lambda_contrast(self.cnd, self.pmtt, freq) + + Amat = make_system_matrix_fast(self.KsdS, lamb) + + if sidx is None: + dGdn = self.dGdn + else: + dGdn = self._compute_dGdn(sidx) + + phi = np.linalg.solve(Amat, dGdn) if Amat.shape[0] == Amat.shape[1] else np.linalg.lstsq(Amat, dGdn, rcond=None)[0] + + Phi = [] + idx = 0 + for i in range(self.nb_incls): + Phi.append(phi[idx : idx + nb_points, :]) + idx += nb_points + return Phi + + def data_simulation(self, freq_list): + """ + Simulate MSR data for different frequencies. + + Parameters + ---------- + freq_list : list or ndarray + List of working frequencies. + + Returns + ------- + dict + Contains "MSR" (list of matrices) and "freq". + """ + freq_list = np.atleast_1d(freq_list) + msr_list = [] + + for f in freq_list: + Phi = self._compute_phi(f) + msr = np.zeros((self.cfg.Ns_total, self.cfg.Nr), dtype=complex) + + for i in range(self.nb_incls): + toto = np.zeros((self.cfg.Ns_total, self.cfg.Nr), dtype=complex) + for s in range(self.cfg.Ns_total): + rcv = self.cfg.rcv(s) + toto[s, :] = SingleLayer_eval(self.D[i], Phi[i][:, s], rcv) + msr += toto + msr_list.append(msr) + + return {"MSR": msr_list, "freq": freq_list} + + @staticmethod + def add_white_noise(data, nlvl): + """ + Add white noise to simulated data. + + Parameters + ---------- + data : dict + Simulated data dictionary. + nlvl : float + Noise level. + + Returns + ------- + dict + Updated dictionary with "MSR_noisy" and "sigma". + """ + out = data.copy() + out["MSR_noisy"] = [] + out["sigma"] = np.zeros(len(data["MSR"])) + + for i, msr in enumerate(data["MSR"]): + msr_real_noisy, sigma_real = add_white_noise(msr.real, nlvl, mode=1, rowmajor=True) + msr_imag_noisy, sigma_imag = add_white_noise(msr.imag, nlvl, mode=1, rowmajor=True) + + out["MSR_noisy"].append(msr_real_noisy + 1j * msr_imag_noisy) + out["sigma"][i] = np.abs(sigma_real + 1j * sigma_imag) + + return out + + @staticmethod + def make_matrix_A(Xs, Z, order): + """Construct the linear operator matrix A for sources/receivers.""" + N = Xs.shape[1] + A = np.zeros((N, 2 * order)) + for n in range(N): + toto = Xs[:, n] - Z + R = np.linalg.norm(toto) + T = np.arctan2(toto[1], toto[0]) + for m in range(1, order + 1): + A[n, 2*(m-1):2*m] = [np.cos(m*T), np.sin(m*T)] / (2 * np.pi * m * R**m) + return A + + def make_linop_CGPT(self, ord_val, symmode=False): + """Make the linear operator for reconstruction.""" + if self.cfg.Ng != 1: + raise NotImplementedError("make_linop_CGPT for Ng > 1 is not implemented yet.") + + As = self.make_matrix_A(self.cfg.all_src, self.cfg.center, ord_val) + Ar = self.make_matrix_A(self.cfg.all_rcv, self.cfg.center, ord_val) + + Ms, Ns = As.shape + Mr, Nr = Ar.shape + + def L(X, transp_flag='notransp'): + if transp_flag == 'notransp': + Xm = X.reshape((Ns, Nr), order='F') + if symmode: + Y = As @ (Xm + Xm.T) @ Ar.T + else: + Y = As @ Xm @ Ar.T + return Y.flatten(order='F') + elif transp_flag == 'transp': + Xm = X.reshape((Ms, Mr), order='F') + if symmode: + Y = As.T @ Xm @ Ar + Ar.T @ Xm.T @ As + else: + Y = As.T @ Xm @ Ar + return Y.flatten(order='F') + + class Op: + def __init__(self, L_func, As, Ar): + self.L = L_func + self.As = As + self.Ar = Ar + + return Op(L, As, Ar) + + def reconstruct_CGPT(self, MSR_list, ord_val, maxiter=100000, tol=1e-10, symmode=True, method='lsqr'): + """ + Reconstruct contracted GPT from MSR data. + + Parameters + ---------- + MSR_list : list of ndarray + MSR data matrices. + ord_val : int + Maximum order of reconstruction. + maxiter : int, optional + Maximum iterations for LSQR. + tol : float, optional + Tolerance for solver. + symmode : bool, optional + If True, assume symmetric CGPT matrix. + method : {'lsqr', 'pinv'}, optional + Solver method. + + Returns + ------- + dict + Contains "CGPT", "res", and "rres". + """ + op = self.make_linop_CGPT(ord_val, symmode) + out = {"CGPT": [], "res": [], "rres": []} + + Ns, Nr = self.cfg.Ns_total, self.cfg.Nr + + for msr in MSR_list: + if method == 'lsqr': + from scipy.sparse.linalg import LinearOperator + lin_op = LinearOperator((Ns * Nr, 4 * ord_val**2), matvec=lambda x: op.L(x, 'notransp'), rmatvec=lambda x: op.L(x, 'transp')) + b = msr.flatten(order='F') + + res_real = sparse_lsqr(lin_op, b.real, atol=tol, btol=tol, iter_lim=maxiter) + res_imag = sparse_lsqr(lin_op, b.imag, atol=tol, btol=tol, iter_lim=maxiter) + + X = res_real[0] + 1j * res_imag[0] + CGPT = X.reshape((2 * ord_val, 2 * ord_val), order='F') + + if symmode: + CGPT = (CGPT + CGPT.T) + + out["CGPT"].append(CGPT) + residual = msr - op.L(CGPT, 'notransp').reshape((Ns, Nr), order='F') + out["res"].append(np.linalg.norm(residual, 'fro')) + out["rres"].append(out["res"][-1] / np.linalg.norm(msr, 'fro')) + else: + CGPT = np.linalg.pinv(op.As) @ msr @ np.linalg.pinv(op.Ar.T) + out["CGPT"].append(CGPT) + residual = msr - op.As @ CGPT @ op.Ar.T + out["res"].append(np.linalg.norm(residual, 'fro')) + out["rres"].append(out["res"][-1] / np.linalg.norm(msr, 'fro')) + return out diff --git a/sies/shape/__init__.py b/sies/shape/__init__.py new file mode 100644 index 0000000..69a2c58 --- /dev/null +++ b/sies/shape/__init__.py @@ -0,0 +1,4 @@ +from .base import Shape +from .triangle import Triangle +from .ellipse import Ellipse +from .flower import Flower diff --git a/sies/shape/base.py b/sies/shape/base.py new file mode 100644 index 0000000..7b06e3a --- /dev/null +++ b/sies/shape/base.py @@ -0,0 +1,313 @@ +import numpy as np +from scipy.interpolate import CubicSpline +import matplotlib.pyplot as plt + +class Shape: + """ + Base class for C2-smooth closed boundary. + + Parameters + ---------- + points : ndarray + Coordinates of boundary points, an array of dimension 2 x nb_points. + tvec : ndarray + Tangent vectors of boundary points. + avec : ndarray + Acceleration vectors of boundary points. + normal : ndarray + Outward normal vectors of boundary points. + center_of_mass : array_like, optional + Center of mass of the boundary. If None, it is calculated from points. + name : str, optional + Name of the shape. + """ + + def __init__(self, points, tvec, avec, normal, center_of_mass=None, name=""): + self.points = np.asarray(points, dtype=float) + self.tvec = np.asarray(tvec, dtype=float) + self.avec = np.asarray(avec, dtype=float) + self.normal = np.asarray(normal, dtype=float) + self.name = name + + if center_of_mass is not None: + self.center_of_mass = np.asarray(center_of_mass, dtype=float).flatten() + else: + self.center_of_mass = self.get_com(self.points, self.tvec, self.normal) + + @property + def nb_points(self): + """int : Number of discrete boundary points.""" + return self.points.shape[1] + + @property + def theta(self): + """ndarray : Non tied-off parameterization between [0, 2pi).""" + return 2 * np.pi * np.arange(self.nb_points) / self.nb_points + + @property + def tvec_norm(self): + """ndarray : Norm of the tangent vectors.""" + return np.sqrt(np.sum(self.tvec**2, axis=0)) + + @property + def sigma(self): + """ndarray : Element of curve integration (integration weights).""" + return (2 * np.pi / self.nb_points) * self.tvec_norm + + @property + def box(self): + """ndarray : A minimal rectangular box [width, height] containing the shape.""" + d = self.points - self.center_of_mass[:, None] + w = np.max(d[0, :]) - np.min(d[0, :]) + h = np.max(d[1, :]) - np.min(d[1, :]) + return np.array([w, h]) + + @property + def diameter(self): + """float : Upper bound of the diameter of the shape.""" + d = self.points - self.center_of_mass[:, None] + return 2 * np.max(np.sqrt(np.sum(d**2, axis=0))) + + def __add__(self, z0): + """ + Translate the shape. + + Parameters + ---------- + z0 : array_like + Translation vector [dx, dy]. + + Returns + ------- + Shape + The translated shape. + """ + z0 = np.asarray(z0).flatten() + if z0.size != 2: + raise ValueError("Translation vector must have size 2") + new_points = self.points + z0[:, None] + new_com = self.center_of_mass + z0 + return Shape(new_points, self.tvec, self.avec, self.normal, new_com, self.name) + + def __sub__(self, z0): + """ + Translate the shape (subtraction). + + Parameters + ---------- + z0 : array_like + Translation vector [dx, dy]. + + Returns + ------- + Shape + The translated shape. + """ + return self.__add__(-np.asarray(z0)) + + def __mul__(self, s): + """ + Scale the shape. + + Parameters + ---------- + s : float + Positive scaling factor. + + Returns + ------- + Shape + The scaled shape. + """ + if not isinstance(s, (int, float)) or s <= 0: + raise ValueError("Scaling factor must be a positive scalar") + new_points = self.points * s + new_com = self.center_of_mass * s + new_tvec = self.tvec * s + new_avec = self.avec * s + return Shape(new_points, new_tvec, new_avec, self.normal, new_com, self.name) + + def __matmul__(self, phi): + """ + Rotate the shape. + + Parameters + ---------- + phi : float + Rotation angle in radians. + + Returns + ------- + Shape + The rotated shape. + """ + if not isinstance(phi, (int, float)): + raise ValueError("Rotation angle must be a scalar") + + rot = np.array([[np.cos(phi), -np.sin(phi)], + [np.sin(phi), np.cos(phi)]]) + + new_points = rot @ self.points + new_com = rot @ self.center_of_mass + new_tvec = rot @ self.tvec + new_normal = rot @ self.normal + new_avec = rot @ self.avec + + return Shape(new_points, new_tvec, new_avec, new_normal, new_com, self.name) + + def plot(self, *args, **kwargs): + """ + Plot the boundary points. + + Parameters + ---------- + *args, **kwargs + Arguments passed to matplotlib.plt.plot. + """ + plt.plot(self.points[0, :], self.points[1, :], *args, **kwargs) + + def isinside(self, x): + """ + Check if a point is inside the bounding ball of the shape. + + Parameters + ---------- + x : array_like + Point coordinates [x, y]. + + Returns + ------- + bool + True if inside, False otherwise. + """ + x = np.asarray(x).flatten() + return np.linalg.norm(x - self.center_of_mass) < self.diameter / 2 + + def isdisjoint(self, other): + """ + Check if two shapes are disjoint based on their bounding balls. + + Parameters + ---------- + other : Shape + The other shape to check. + + Returns + ------- + bool + True if disjoint, False otherwise. + """ + dist = np.linalg.norm(self.center_of_mass - other.center_of_mass) + return dist > (self.diameter + other.diameter) / 2 + + @staticmethod + def get_com(points, tvec, normal): + """ + Calculate the center of mass of a shape using the Stokes formula. + + Parameters + ---------- + points : ndarray + Boundary points. + tvec : ndarray + Tangent vectors. + normal : ndarray + Normal vectors. + + Returns + ------- + ndarray + Center of mass coordinates [cx, cy]. + """ + nb_points = points.shape[1] + tvec_norm = np.sqrt(np.sum(tvec**2, axis=0)) + sigma = (2 * np.pi / nb_points) * tvec_norm + + mass = (np.sum(points[0, :] * normal[0, :] * sigma) + + np.sum(points[1, :] * normal[1, :] * sigma)) / 2 + + cx = np.sum(0.5 * (points[0, :]**2) * normal[0, :] * sigma) + cy = np.sum(0.5 * (points[1, :]**2) * normal[1, :] * sigma) + + return np.array([cx, cy]) / mass + + @staticmethod + def boundary_vec_interpl(points0, theta0, theta): + """ + Interpolate boundary vectors using cubic splines. + + Parameters + ---------- + points0 : ndarray + Original boundary points. + theta0 : ndarray + Original parameterization. + theta : ndarray + New parameterization for resampling. + + Returns + ------- + tuple + (points, tvec, avec, normal) interpolated at theta. + """ + p0_periodic = np.hstack([points0, points0[:, :1]]) + t0_periodic = np.append(theta0, theta0[-1] + (theta0[1]-theta0[0] if len(theta0)>1 else 2*np.pi)) + + cs_x = CubicSpline(t0_periodic, p0_periodic[0, :], bc_type='periodic') + cs_y = CubicSpline(t0_periodic, p0_periodic[1, :], bc_type='periodic') + + px = cs_x(theta) + py = cs_y(theta) + points = np.vstack([px, py]) + + tx = cs_x(theta, 1) + ty = cs_y(theta, 1) + tvec = np.vstack([tx, ty]) + + normal = np.vstack([tvec[1, :], -tvec[0, :]]) + normal /= np.sqrt(np.sum(normal**2, axis=0)) + + ax = cs_x(theta, 2) + ay = cs_y(theta, 2) + avec = np.vstack([ax, ay]) + + return points, tvec, avec, normal + + @staticmethod + def rescale(D0, theta0, nb_points, nsize=None, dspl=1): + """ + Rescale and resample the boundary. + + Parameters + ---------- + D0 : ndarray + Initial boundary points. + theta0 : ndarray + Initial parameterization. + nb_points : int + Number of points for resampling. + nsize : tuple, optional + Target [width, height] for rescaling. + dspl : int, optional + Down-sampling factor for smoothing. + + Returns + ------- + tuple + (points, tvec, avec, normal) after rescaling and resampling. + """ + if nsize is not None: + minx, maxx = np.min(D0[0, :]), np.max(D0[0, :]) + miny, maxy = np.min(D0[1, :]), np.max(D0[1, :]) + z0 = np.array([(minx + maxx) / 2, (miny + maxy) / 2]) + D0 = np.vstack([ + (D0[0, :] - z0[0]) * (nsize[0] / (maxx - minx)), + (D0[1, :] - z0[1]) * (nsize[1] / (maxy - miny)) + ]) + + dspl = int(np.ceil(dspl)) + idx = np.arange(0, D0.shape[1], dspl) + + theta = np.linspace(theta0[0], theta0[-1], nb_points, endpoint=False) + + return Shape.boundary_vec_interpl(D0[:, idx], theta0[idx], theta) diff --git a/sies/shape/ellipse.py b/sies/shape/ellipse.py new file mode 100644 index 0000000..65e1854 --- /dev/null +++ b/sies/shape/ellipse.py @@ -0,0 +1,44 @@ +import numpy as np +from .base import Shape + +class Ellipse(Shape): + def __init__(self, a, b, nb_points): + """ + Ellipse shape. + a, b: length of semi-major and semi-minor axis + nb_points: number of discretization points + """ + if a < b: + raise ValueError("Semi-major axis a must be >= semi-minor axis b") + + self.axis_a = a + self.axis_b = b + self.phi_rot = 0 + + com = np.array([0, 0]) + theta = 2 * np.pi * np.arange(nb_points) / nb_points + + points = np.vstack([a * np.cos(theta), b * np.sin(theta)]) + tvec = np.vstack([-a * np.sin(theta), b * np.cos(theta)]) + avec = np.vstack([-a * np.cos(theta), -b * np.sin(theta)]) + + # rotation = [[0 1];[-1 0]] in Matlab + normal = np.vstack([tvec[1, :], -tvec[0, :]]) + normal /= np.sqrt(np.sum(normal**2, axis=0)) + + name = "Circle" if a == b else "Ellipse" + super().__init__(points, tvec, avec, normal, com, name) + + def __mul__(self, s): + new_obj = super().__mul__(s) + new_obj.axis_a = self.axis_a * s + new_obj.axis_b = self.axis_b * s + new_obj.phi_rot = self.phi_rot + return new_obj + + def __matmul__(self, phi): + new_obj = super().__matmul__(phi) + new_obj.axis_a = self.axis_a + new_obj.axis_b = self.axis_b + new_obj.phi_rot = self.phi_rot + phi + return new_obj diff --git a/sies/shape/flower.py b/sies/shape/flower.py new file mode 100644 index 0000000..a91721e --- /dev/null +++ b/sies/shape/flower.py @@ -0,0 +1,72 @@ +import numpy as np +from .base import Shape + +class Flower(Shape): + def __init__(self, a, b, nb_points, nb_petals=5, epsilon=0.3, tau=0): + """ + Flower shape (rotational symmetric shapes). + a, b: typical size of the flower + nb_points: number of discretization points + nb_petals: number of petals + epsilon: size of perturbation + tau: percentage of damage (0 to 1) - Currently only tau=0 is fully supported + """ + if tau != 0: + # damaged flower not implemented for now to keep it simple, + # as it's not used in the basic demo. + raise NotImplementedError("Damaged flower (tau > 0) is not yet implemented.") + + self.axis_a = a + self.axis_b = b + self.nb_petals = nb_petals + self.epsilon = epsilon + self.tau = tau + self.phi_rot = 0 + + theta = 2 * np.pi * np.arange(nb_points) / nb_points + k = 1 # fixed in Flower.m constructor as pertb = 1 + + # Position + radial_pert = (1 + epsilon * np.cos(nb_petals * theta)**k) + points = np.vstack([ + a * np.cos(theta) * radial_pert, + b * np.sin(theta) * radial_pert + ]) + + # Velocity (tangent) + # Using the formula from make_flower.m + tx = -a * (np.sin(theta) * radial_pert + k * epsilon * nb_petals * np.sin(nb_petals * theta) * np.cos(nb_petals * theta)**(k-1) * np.cos(theta)) + ty = b * (np.cos(theta) * radial_pert - k * epsilon * nb_petals * np.sin(nb_petals * theta) * np.cos(nb_petals * theta)**(k-1) * np.sin(theta)) + tvec = np.vstack([tx, ty]) + + # Acceleration + # Using the formula from make_flower.m for k=1 + ax = -a * (np.cos(theta) * (1 + epsilon * np.cos(nb_petals * theta)) - 2 * epsilon * nb_petals * np.sin(theta) * np.sin(nb_petals * theta) + epsilon * nb_petals**2 * np.cos(nb_petals * theta) * np.cos(theta)) + ay = -b * (np.sin(theta) * (1 + epsilon * np.cos(nb_petals * theta)) + 2 * epsilon * nb_petals * np.cos(theta) * np.sin(nb_petals * theta) + epsilon * nb_petals**2 * np.cos(nb_petals * theta) * np.sin(theta)) + avec = np.vstack([ax, ay]) + + # Normal + normal = np.vstack([tvec[1, :], -tvec[0, :]]) + normal /= np.sqrt(np.sum(normal**2, axis=0)) + + super().__init__(points, tvec, avec, normal, center_of_mass=[0,0], name='Flower') + + def __mul__(self, s): + new_obj = super().__mul__(s) + new_obj.axis_a = self.axis_a * s + new_obj.axis_b = self.axis_b * s + new_obj.nb_petals = self.nb_petals + new_obj.epsilon = self.epsilon + new_obj.tau = self.tau + new_obj.phi_rot = self.phi_rot + return new_obj + + def __matmul__(self, phi): + new_obj = super().__matmul__(phi) + new_obj.axis_a = self.axis_a + new_obj.axis_b = self.axis_b + new_obj.nb_petals = self.nb_petals + new_obj.epsilon = self.epsilon + new_obj.tau = self.tau + new_obj.phi_rot = self.phi_rot + phi + return new_obj diff --git a/sies/shape/triangle.py b/sies/shape/triangle.py new file mode 100644 index 0000000..8464975 --- /dev/null +++ b/sies/shape/triangle.py @@ -0,0 +1,70 @@ +import numpy as np +from .base import Shape + +class Triangle(Shape): + def __init__(self, a, angl, nb_points, dspl=10): + """ + Isosceles triangle. + a: length of equal sides + angl: angle between equal sides (radians) + nb_points: number of discretization points + dspl: down-sampling factor for smoothing + """ + self.lside = a + self.angl = angl + + h = a * np.cos(angl / 2) + b = a * np.sin(angl / 2) + + t1 = a / (a + b) / 2 + t2 = b / (a + b) + t3 = t1 + + n1 = int(np.floor(t1 * nb_points)) + n2 = int(np.floor(t2 * nb_points)) + n3 = nb_points - n1 - n2 + + A = np.array([[0], [2/3*h]]) + B = np.array([[-b], [-h/3]]) + C = np.array([[b], [-h/3]]) + + t = np.linspace(0, 1, n1, endpoint=False) + AB = A + (B - A) @ t[None, :] + + t = np.linspace(0, 1, n2, endpoint=False) + BC = B + (C - B) @ t[None, :] + + t = np.linspace(0, 1, n3, endpoint=False) + CA = C + (A - C) @ t[None, :] + + points0 = np.hstack([AB, BC, CA]) + theta0 = np.linspace(0, 2 * np.pi, nb_points, endpoint=False) + + if dspl >= 1: + t0 = n3 // 2 + points = np.roll(points0, t0, axis=1) + points, tvec, avec, normal = Shape.rescale(points, theta0, nb_points, dspl=dspl) + else: + # Fallback for dspl < 1 (not recommended) + tvec = np.hstack([ + np.tile((B - A) / t1, (1, n1)), + np.tile((C - B) / t2, (1, n2)), + np.tile((A - C) / t3, (1, n3)) + ]) / (2 * np.pi) + + normal = np.vstack([tvec[1, :], -tvec[0, :]]) + normal /= np.sqrt(np.sum(normal**2, axis=0)) + avec = np.zeros((2, nb_points)) + + t0 = n3 // 2 + points = np.roll(points0, t0, axis=1) + tvec = np.roll(tvec, t0, axis=1) + normal = np.roll(normal, t0, axis=1) + + super().__init__(points, tvec, avec, normal, center_of_mass=[0,0], name='Triangle') + + def __mul__(self, s): + new_obj = super().__mul__(s) + new_obj.lside = self.lside * s + new_obj.angl = self.angl + return new_obj diff --git a/sies/tools/__init__.py b/sies/tools/__init__.py new file mode 100644 index 0000000..1608d9a --- /dev/null +++ b/sies/tools/__init__.py @@ -0,0 +1,30 @@ +import numpy as np + +def add_white_noise(X, nlvl, mode=0, rowmajor=False): + """ + Add white noise to data. + X: input data matrix + nlvl: noise level + mode: if not 0, each column/row (depending on rowmajor) is treated independently + rowmajor: if True, noise is added row by row + """ + X = np.asarray(X) + if rowmajor: + Y_T, sigma = add_white_noise_mat(X.T, nlvl, mode) + return Y_T.T, sigma + else: + return add_white_noise_mat(X, nlvl, mode) + +def add_white_noise_mat(X, nlvl, mode): + M, N = X.shape + if mode: + Y = np.zeros_like(X) + for c in range(N): + t0 = np.linalg.norm(X[:, c]) / np.sqrt(M) + Y[:, c] = X[:, c] + np.random.randn(M) * t0 * nlvl + else: + t0 = np.linalg.norm(X, 'fro') / np.sqrt(M * N) + Y = X + np.random.randn(M, N) * t0 * nlvl + + sigma = np.linalg.norm(X, 'fro') / np.sqrt(M * N) * nlvl + return Y, sigma diff --git a/sies/tools/laplacian/__init__.py b/sies/tools/laplacian/__init__.py new file mode 100644 index 0000000..64632ef --- /dev/null +++ b/sies/tools/laplacian/__init__.py @@ -0,0 +1,40 @@ +import numpy as np + +def Green2D_Grad(X, Y): + """ + Gradient of the 2D Green function G(x) = 1/2/pi * log|x| evaluated at the points X-Y. + X: (2, M), Y: (2, N) + Outputs: Gx, Gy (M, N) matrices + """ + X = np.asarray(X) + Y = np.asarray(Y) + + # X[0, :] - Y[0, :] using broadcasting + # X[0, :, None] is (M, 1), Y[0, None, :] is (1, N) + # Result is (M, N) + XY1 = X[0, :, None] - Y[0, None, :] + XY2 = X[1, :, None] - Y[1, None, :] + + DN = XY1**2 + XY2**2 + # Handle zero division if any (e.g. X == Y) + DN[DN == 0] = np.inf + + Gx = (1 / (2 * np.pi)) * XY1 / DN + Gy = (1 / (2 * np.pi)) * XY2 / DN + return Gx, Gy + +def Green2D_Dn(X, Y, normal): + """ + Normal derivative of the 2D Green function G(y,x) = 1/2/pi * log|y-x| + with respect to y on a boundary. + X: (2, M), Y: (2, N), normal: (2, N) + Output: Gn (M, N) matrix + """ + # Gx, Gy evaluated at (Y - X) + # So Y is "X" in Green2D_Grad and X is "Y" + Gx, Gy = Green2D_Grad(Y, X) # (N, M) + + # Gn = normal[0,:] * Gx + normal[1,:] * Gy + # normal[0,:] is (N,), Gx is (N, M) -> broadcasting works if we do normal[0,:,None] * Gx + Gn = normal[0, :, None] * Gx + normal[1, :, None] * Gy + return Gn.T # (M, N) diff --git a/tests/test_sies.py b/tests/test_sies.py new file mode 100644 index 0000000..a28816b --- /dev/null +++ b/tests/test_sies.py @@ -0,0 +1,50 @@ +import numpy as np +import pytest +from sies import Shape, Triangle, Ellipse, Coincided, Conductivity_R2, theoretical_CGPT, lambda_contrast + +def test_shape_dsl(): + # Test translation + e = Ellipse(1, 0.5, 100) + e2 = e + [1, 2] + assert np.allclose(e2.center_of_mass, [1, 2]) + + # Test scaling + e3 = e * 2 + assert np.isclose(e3.diameter, 2 * e.diameter) + + # Test rotation + e4 = e @ (np.pi/2) + # After 90 deg rotation, semi-major axis (x) becomes semi-minor (y) and vice versa + assert np.allclose(e4.box, [e.box[1], e.box[0]], atol=1e-7) + +def test_conductivity_simulation(): + D = [Ellipse(0.5, 0.25, 256)] + cfg = Coincided([0, 0], 10, 32) + P = Conductivity_R2(D, [10.0], [0.0], cfg) + + freqs = [0, 10] + data = P.data_simulation(freqs) + + assert len(data["MSR"]) == 2 + assert data["MSR"][0].shape == (32, 32) + +def test_cgpt_reconstruction(): + # Use a simple circle for very stable reconstruction + D = [Ellipse(0.5, 0.5, 256)] + cfg = Coincided([0, 0], 10, 64) + P = Conductivity_R2(D, [5.0], [0.0], cfg) + + freq = 0 + data = P.data_simulation([freq]) + + ord_val = 1 + # Theoretical CGPT + lamb = lambda_contrast([5.0], [0.0], freq) + M_theo = theoretical_CGPT(D, lamb, ord_val) + + # Reconstruct + out = P.reconstruct_CGPT(data["MSR"], ord_val, method='pinv') + M_rec = out["CGPT"][0] + + # For a circle, M should be very close to theoretical + assert np.linalg.norm(M_theo - M_rec) / np.linalg.norm(M_theo) < 1e-2