.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "gallery/time_fourier_demo.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. .. rst-class:: sphx-glr-example-title .. _sphx_glr_gallery_time_fourier_demo.py: Learning a time-dependent force field — time-Fourier basis ========================================================== Recover an **unknown, time-varying** force law from trajectories alone, by expanding time in a Fourier dictionary. The companion :doc:`time_dependent_forcing_demo` infers a *known, recorded* drive. Here the drive is **not** recorded: a population of overdamped particles sits in a common harmonic trap whose **centre** :math:`a(t)` *and* **stiffness** :math:`k(t)` both wander in time, .. math:: \mathrm{d}x = -k(t)\,\bigl[x - a(t)\bigr]\,\mathrm{d}t + \sqrt{2D}\,\mathrm{d}W . SFI never sees :math:`k(t)` or :math:`a(t)`. It is handed only the trajectories and a :func:`~SFI.bases.time_fourier` dictionary, which reads the **auto-injected** ``time`` clock (no protocol has to be recorded). The linear estimator reconstructs both time-dependent fields, and :term:`PASTIS` keeps only the harmonics the data support. .. rubric:: Tags synthetic · overdamped · linear · 1D · time-dependent .. GENERATED FROM PYTHON SOURCE LINES 28-30 .. code-block:: Python :dedent: 1 .. GENERATED FROM PYTHON SOURCE LINES 50-58 A trap that drifts and stiffens ------------------------------- The centre :math:`a(t)` and stiffness :math:`k(t)` are smooth but *not* single sinusoids — each is a short sum of harmonics of the fundamental :math:`\omega = 2\pi / T_\text{tot}` (one period over the whole run). They are the ground truth we will try to recover; the estimator is never told them. .. GENERATED FROM PYTHON SOURCE LINES 58-78 .. code-block:: Python import SFI from SFI.bases import X, time_fourier, unit_vector_basis, extra_scalar from SFI.langevin import OverdampedProcess from SFI.trajectory import time_series_extra from SFI.utils.plotting import plot_time_profile_comparison, timeseries dt = 0.02 Nframes = 2500 Nparticles = 48 D_true = 0.08 t = np.arange(Nframes) * dt T_tot = t[-1] w = 2 * np.pi / T_tot k_t = 2.0 + 0.8 * np.cos(w * t) + 0.4 * np.cos(2 * w * t) # stiffness > 0 b_t = 0.9 * np.sin(w * t) + 0.3 * np.cos(2 * w * t) # additive drift k·a a_t = b_t / k_t # implied trap centre .. GENERATED FROM PYTHON SOURCE LINES 79-89 Simulate an ensemble -------------------- A single trapped particle cannot separate *"the trap is stiffer"* from *"the trap has moved"* — at any instant it just sits near the minimum. A **population** breaks that degeneracy: at each time the spread of particles around :math:`a(t)` fixes the slope :math:`-k(t)` and the offset :math:`k(t)\,a(t)`. The particles are independent, so we simulate ``Nparticles`` single-particle runs and merge them into one ensemble collection. .. GENERATED FROM PYTHON SOURCE LINES 89-108 .. code-block:: Python F_true = extra_scalar("neg_k") * X(dim=1) & extra_scalar("drift") * unit_vector_basis(1) proc = OverdampedProcess(F=F_true, D=D_true, theta_F=jnp.array([1.0, 1.0])) proc.set_extras(extras_global={ "neg_k": time_series_extra(-k_t), "drift": time_series_extra(b_t), }) rng = np.random.default_rng(0) runs = [] for i in range(Nparticles): proc.initialize(jnp.array([float(rng.normal() * 0.4)])) ci = proc.simulate(dt=dt, Nsteps=Nframes, key=random.PRNGKey(i), oversampling=3, compute_observables=False) runs.append(ci) coll = runs[0].merge(runs[1:]) # ensemble = many single-particle datasets print(f"ensemble: {coll.datasets[0].T} frames x {Nparticles} particles") .. rst-class:: sphx-glr-script-out .. code-block:: none ensemble: 2500 frames x 48 particles .. GENERATED FROM PYTHON SOURCE LINES 109-115 Trajectories and the hidden protocol ------------------------------------ The cloud of particles tracks the moving centre; the bottom panel shows the stiffness they actually feel. Neither curve is given to the estimator. .. GENERATED FROM PYTHON SOURCE LINES 115-117 .. code-block:: Python :dedent: 1 .. image-sg:: /gallery/images/sphx_glr_time_fourier_demo_001.png :alt: Ensemble in a drifting, stiffening trap :srcset: /gallery/images/sphx_glr_time_fourier_demo_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none examples/gallery/time_fourier_demo.py:129: UserWarning: The figure layout has changed to tight fig.tight_layout() .. GENERATED FROM PYTHON SOURCE LINES 134-143 Inference with a time-Fourier dictionary ---------------------------------------- Tensor a Fourier-in-time dictionary with the spatial terms: ``time_fourier(n) * X`` carries the **stiffness**, and ``time_fourier(n) * unit_vector_basis(1)`` the **moving centre**. :func:`~SFI.bases.time_fourier` reads the auto-injected ``time`` extra; with ``period=None`` its fundamental is the full trajectory duration. Nothing about :math:`a(t)` or :math:`k(t)` is supplied. .. GENERATED FROM PYTHON SOURCE LINES 143-152 .. code-block:: Python n_modes = 6 B = time_fourier(n_modes) * X(dim=1) & time_fourier(n_modes) * unit_vector_basis(1) inf = SFI.OverdampedLangevinInference(coll) inf.compute_diffusion_constant(method="MSD") inf.infer_force_linear(B) inf.compute_force_error() .. GENERATED FROM PYTHON SOURCE LINES 153-159 Reconstructed stiffness and centre ---------------------------------- The two coefficient blocks give the time functions :math:`-k(t)` (the ``x`` block) and :math:`k(t)\,a(t)` (the constant block); dividing recovers the centre. Both match the hidden ground truth. .. GENERATED FROM PYTHON SOURCE LINES 159-161 .. code-block:: Python :dedent: 1 .. image-sg:: /gallery/images/sphx_glr_time_fourier_demo_002.png :alt: Recovered time-dependent force field (never shown the protocol), stiffness $k(t)$, trap centre $a(t)$ :srcset: /gallery/images/sphx_glr_time_fourier_demo_002.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none examples/gallery/time_fourier_demo.py:196: UserWarning: The figure layout has changed to tight fig.tight_layout() k(t) reconstruction NMSE = 0.0024 a(t) reconstruction NMSE = 0.0012 .. GENERATED FROM PYTHON SOURCE LINES 204-211 Sparse selection of harmonics ----------------------------- The dictionary offers ``n_modes`` harmonics per spatial term, but the data support only a few. :term:`PASTIS` prunes the rest — the surviving terms are exactly the harmonics built into :math:`k(t)` and :math:`k(t)\,a(t)`. .. GENERATED FROM PYTHON SOURCE LINES 211-215 .. code-block:: Python inf.sparsify_force(criterion="PASTIS") inf.print_report() .. rst-class:: sphx-glr-script-out .. code-block:: none --- StochasticForceInference Report --- Average diffusion tensor: [[0.0788167]] Measurement noise tensor: [[6.256998e-05]] Force estimated information: 1186.5589599609375 Force: estimated normalized mean squared error (sampling only): 0.010956048988142434 Force Coefficient Table ──────────────────────────────────────────────────────────────────── # Label Coefficient Std.Err SNR Sig ──────────────────────────────────────────────────────────────────── 0 1·x -1.97601e+00 4.13784e-02 47.8 ** 1 cos(1wt)·x -8.41564e-01 6.12111e-02 13.7 ** 3 cos(2wt)·x -4.20467e-01 5.81531e-02 7.2 * 15 sin(1wt)·e0 8.69886e-01 2.77417e-02 31.4 ** 16 cos(2wt)·e0 3.03681e-01 2.46001e-02 12.3 ** ──────────────────────────────────────────────────────────────────── 5/26 basis functions in support, sig: 5* / 4** / 0*** (|SNR| ≥ 2 / 10 / 100) (Std.err. reflects sampling error only; discretization bias is not included.) Zeroed (21): sin(1wt)·x, sin(2wt)·x, cos(3wt)·x, sin(3wt)·x, cos(4wt)·x, sin(4wt)·x, cos(5wt)·x, sin(5wt)·x, cos(6wt)·x, sin(6wt)·x, 1·e0, cos(1wt)·e0, sin(2wt)·e0, cos(3wt)·e0, sin(3wt)·e0, cos(4wt)·e0, sin(4wt)·e0, cos(5wt)·e0, sin(5wt)·e0, cos(6wt)·e0, sin(6wt)·e0 .. GENERATED FROM PYTHON SOURCE LINES 216-233 Notes ----- - **Time as an extra.** The reserved ``time`` extra is injected automatically, per frame, by the trajectory layer (see :meth:`~SFI.trajectory.TrajectoryDataset.build_extras`), so :func:`~SFI.bases.time_fourier` needs no recorded protocol. With ``period=None`` the fundamental is the **full trajectory duration**; pass ``period=`` to fit a known repeat time. - **Why an ensemble.** The stiffness and the centre are jointly identifiable only because many particles sample a spread of positions at each instant; a single tightly-trapped trajectory cannot tell a stiffer trap from a displaced one. - **Beyond traps.** Any time-dependent force field — ramps, oscillatory drives, slow aging — can be learned the same way: tensor :func:`~SFI.bases.time_fourier` with whatever spatial basis the problem needs, and let PASTIS keep the active modes. .. GENERATED FROM PYTHON SOURCE LINES 233-235 .. code-block:: Python stamp_output() .. rst-class:: sphx-glr-script-out .. code-block:: none [Generated: 2026-06-30 10:13] .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 31.602 seconds) .. rst-class:: sphx-glr-example-tags 🏷 Tags: synthetic, overdamped, linear, 1D, time-dependent .. _sphx_glr_download_gallery_time_fourier_demo.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: time_fourier_demo.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: time_fourier_demo.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: time_fourier_demo.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_