Pade approximant based spectral extrapolation

python/meep.i

@@ -1723,6 +1723,7 @@ PyObject *_get_array_slice_dimensions(meep::fields *f, const meep::volume &where
         FluxRegion,
         ForceRegion,
         Harminv,
+        Pade,
         Identity,
         Mirror,
         ModeRegion,

python/simulation.py

@@ -9,6 +9,7 @@
 import warnings
 from collections import OrderedDict, namedtuple
 from typing import Callable, List, Optional, Tuple, Union
+from scipy.interpolate import pade
 
 try:
     from collections.abc import Sequence, Iterable
@@ -88,7 +89,11 @@ def fix_dft_args(args, i):
 
 
 def get_num_args(func):
-    return 2 if isinstance(func, Harminv) else func.__code__.co_argcount
+    return (
+        2
+        if isinstance(func, Harminv) or isinstance(func, Pade)
+        else func.__code__.co_argcount
+    )
 
 
 def vec(*args):
@@ -859,6 +864,168 @@ def amplitude(self, point, component):
         return mp.eigenmode_amplitude(self.swigobj, swig_point, component)
 
 
+class Pade:
+    """
+    Padé approximant based spectral extrapolation is implemented as a class with a [`__call__`](#Pade.__call__) method,
+    which allows it to be used as a step function that collects field data from a given
+    point and runs [Padé](https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.pade.html)
+    on that data to extract an analytic rational function which approximates the frequency response.
+    For more information about the Padé approximant, see: https://en.wikipedia.org/wiki/Padé_approximant.
+
+    See [`__init__`](#Pade.__init__) for details about constructing a `Pade`.
+
+    In particular, `Padé` stores the discrete time series $\\hat{f}[n]$ corresponding to the given field
+    component as a function of time and expresses it as:
+
+    $$\\hat{f}(\\omega) = \\sum_n \\hat{f}[n] e^{i\\omega n \\Delta t}$$
+
+    The above is a "Taylor-like" polynomial in $$n$$ with a Fourier basis and
+    coefficients which are the sampled field data. We then compute the Padé approximant
+    to be the analytic form of this function as:
+
+    $$R(\\omega) = \\frac{P(\\omega)}{Q(\\omega)}$$
+
+    Where $$P$$ and $$Q$$ are polynomials of degree $$m$$ and $$n$$, and $$m + n + 1$$ is the
+    degree of agreement of the Padé approximant to the analytic function $$f(\\omega)$$. This
+    function $$R$$ is stored in the callable method `pade_instance.freq_response`.
+    Be sure to save a reference to the `Pade` instance if you wish
+    to use the results after the simulation:
+
+    ```py
+    sim = mp.Simulation(...)
+    p = mp.Pade(...)
+    sim.run(p, until=time)
+    # do something with p.freq_response
+    ```
+    """
+
+    def __init__(
+        self,
+        c: int = None,
+        pt: Vector3Type = None,
+        m: Optional[Union[int, float]] = None,
+        n: Optional[Union[int, float]] = None,
+        sampling_interval: int = 1,
+        start_time: Optional[int] = 0,
+        stop_time: Optional[int] = None,
+    ):
+        """
+        Construct a Padé object.
+
+        A `Pade` is a step function that collects data from the field component `c`
+        (e.g. $E_x$, etc.) at the given point `pt` (a `Vector3`). Then, at the end
+        of the run, it uses the scipy Padé algorithm to approximate the analytic
+        frequency response at the specified point.
+
+        + **`c` [`Component`]** — Specifies the field component to use for extrapolation.
+         No default.
+        + **`pt` [`Vector3`]** — Specifies the location to accumulate fields. No default.
+        + **`m` [`Optional[Union[int,float]]`]** — Specifies the order of the numerator $$P$$. Behavior
+         is inferred from the supplied data type. If int is provided, directly specifies the order
+         of $$P$$. If float is specified, the order is given the given fractional length of the sample data.
+         Defaults to half the length of the sampled field data.
+        + **`n` [`Optional[Union[int, float]]`]** — Specifies the order of the denominator $$Q$$ with
+         similar behavior to `n`. Defaults to length of sampled data - `m` - 1.
+        + **`sampling_interval` [`int`]** — Specifies the interval at which to sample the field data.
+         Defaults to 1.
+        + **`start_time` [`Optional[int]`]** — Specifies the time (in increments of dt) at which
+         to start sampling the field data. Default 0 (beginning of simulation).
+        + **`stop_time` [`Optional[int]`]** — Specifies the time (in increments of dt) at which
+         to stop sampling the field data. Default None (end of simulation).
+        """
+        self.c = c
+        self.pt = pt
+        self.m = m
+        self.n = n
+        self.sampling_interval = sampling_interval
+        self.start_time = start_time
+        self.stop_time = stop_time
+        self.data = []
+        self.data_dt = 0
+        self.freq_response: Callable = None
+        self.step_func = self._pade()
+
+    def __call__(self, sim, todo):
+        """
+        Allows a Pade instance to be used as a step function.
+        """
+        self.step_func(sim, todo)
+
+    def _collect_pade(self):
+        def _collect1(c, pt):
+            self.t0 = 0
+
+            def _collect2(sim):
+                self.data_dt = sim.meep_time() - self.t0
+                self.t0 = sim.meep_time()
+                self.data.append(sim.get_field_point(c, pt))
+
+            return _collect2
+
+        return _collect1
+
+    def _analyze_pade(self, sim):
+        # Sample the collected field data and possibly truncate the beginning or
+        # end to avoid transients
+        samples = self.data[self.start_time : self.stop_time : self.sampling_interval]
+
+        # Infer the desired behavior for m and n from the types of the supplied arguments
+        if not self.m:
+            self.m = int(len(samples) / 2)
+        elif type(self.m) is float and (0 <= self.m <= 1):
+            self.m = int(len(samples) * self.m)
+        elif type(self.m) is int and self.m >= 1:
+            pass
+        else:
+            raise TypeError(
+                "Order of numerator m must be positive integer or float between 0 and 1."
+            )
+
+        if not self.n:
+            self.n = int(len(samples) - self.m - 1)
+        elif type(self.n) is float and (0 <= self.n <= 1):
+            self.n = int(len(samples) * self.n)
+        elif type(self.n) is int and self.n >= 1:
+            pass
+        else:
+            raise TypeError(
+                "Order of denominator n must be positive integer or float between 0 and 1."
+            )
+
+        print(self.m, self.n)
+
+        # Compute the Padé approximant
+        p, q = pade(samples, self.m, self.n)
+
+        # TODO: make compatible with any monitor type (e.g. poynting flux, energy spectra)
+        # TODO: add decimation factor (comes with monitor type)
+        # TODO: include option to use Trefethen SVD-based algorithm
+
+        # Construct the numerator and denominator polynomials
+        P = lambda w: np.sum(
+            [
+                pi * np.exp(1j * w * self.sampling_interval * sim.fields.dt * i)
+                for i, pi in enumerate(p.coef)
+            ]
+        )
+        Q = lambda w: np.sum(
+            [
+                qi * np.exp(1j * w * self.sampling_interval * sim.fields.dt * i)
+                for i, qi in enumerate(q.coef)
+            ]
+        )
+
+        return lambda freq: P(2 * np.pi * freq) / Q(2 * np.pi * freq)
+
+    def _pade(self):
+        def _p(sim):
+            self.freq_response = self._analyze_pade(sim)
+
+        f1 = self._collect_pade()
+
+        return _combine_step_funcs(at_end(_p), f1(self.c, self.pt))
+
+
 class Harminv:
     """
     Harminv is implemented as a class with a [`__call__`](#Harminv.__call__) method,

python/tests/test_ring.py

@@ -3,6 +3,8 @@
 import unittest
 
 import meep as mp
+import numpy as np
+from scipy.signal import find_peaks
 
 
 class TestRing(unittest.TestCase):
@@ -44,6 +46,7 @@ def init(self):
 
         self.sim.use_output_directory(self.temp_dir)
         self.h = mp.Harminv(mp.Ez, mp.Vector3(r + 0.1), fcen, df)
+        self.p = mp.Pade(mp.Ez, mp.Vector3(r + 0.1), sampling_interval=4)
 
     def test_harminv(self):
         self.init()
@@ -70,6 +73,27 @@ def test_harminv(self):
         self.assertAlmostEqual(ep, 11.559999999999999, places=places)
         self.assertAlmostEqual(fp, -0.08185972142450348, places=places)
 
+    def test_pade(self):
+        self.init()
+
+        self.sim.run(
+            self.p,
+            mp.after_sources(self.h),
+            until_after_sources=300,
+        )
+
+        freqs = np.linspace(
+            self.h.fcen - self.h.df / 2, self.h.fcen + self.h.df / 2, 1000
+        )
+
+        freq_domain = [self.p.freq_response(freq) for freq in freqs]
+
+        idx = find_peaks(
+            np.abs(freq_domain) ** 2 / max(np.abs(freq_domain) ** 2), prominence=1e-4
+        )[0]
+
+        self.assertAlmostEqual(freqs[idx[0]], self.h.modes[0].freq, places=3)
+
 
 if __name__ == "__main__":
     unittest.main()