From 10529f7a1f0cd2f68366b7bb57df1fff34244923 Mon Sep 17 00:00:00 2001
From: Wolfgang Kastaun <wolfgang.kastaun@aei.mpg.de>
Date: Wed, 12 Feb 2025 15:38:27 +0100
Subject: [PATCH] Added notebook comparing original shift inversion to new one
---
.../shift_inversion_comparison.ipynb | 363 ++++++++++++++++++
1 file changed, 363 insertions(+)
create mode 100644 docs/notebooks/shift_inversion_comparison.ipynb
diff --git a/docs/notebooks/shift_inversion_comparison.ipynb b/docs/notebooks/shift_inversion_comparison.ipynb
new file mode 100644
index 0000000..8d466c3
--- /dev/null
+++ b/docs/notebooks/shift_inversion_comparison.ipynb
@@ -0,0 +1,363 @@
+{
+ "cells": [
+ {
+ "cell_type": "markdown",
+ "id": "3efcb1ab-d4fb-4d84-aaa6-9222b0feb9a2",
+ "metadata": {},
+ "source": [
+ "# Compare original shift inversion algorithm to new one"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "7a3b9115-11a1-4f19-aaba-36d246860aa9",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "%matplotlib ipympl\n",
+ "from matplotlib import pyplot as plt"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "6cca7a2a-1853-4555-9f12-c27fbb30158e",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import numpy as np"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "c7043bd6-dcde-4840-bb1a-fa819d1c3936",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "from lisainstrument import dsp"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "3d80234a-1eb8-4d8c-b76e-c4fb6fa0f3d8",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "from lisainstrument.shift_inversion_numpy import make_shift_inverse_lagrange_numpy"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "d8e70bed-5bb0-496c-a00a-730613971b2f",
+ "metadata": {},
+ "source": [
+ "Below is the original shift inversion method from the Instrument class. Only modifications are to make it self-contained, and removed logging."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "556511dc-17b4-42bd-867e-0fc9da431fd3",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "def legacy_invert_scet_wrt_tps(\n",
+ " scet_wrt_tps: np.ndarray,\n",
+ " clockinv_tolerance: float,\n",
+ " physics_fs: float,\n",
+ " interpolation_order: int,\n",
+ " clockinv_maxiter: int = 5,\n",
+ "):\n",
+ " edge = min(100, len(scet_wrt_tps) // 2 - 1)\n",
+ " error = 0\n",
+ "\n",
+ " niter = 0\n",
+ " next_inverse = scet_wrt_tps\n",
+ " while not niter or error > clockinv_tolerance:\n",
+ " if niter >= clockinv_maxiter:\n",
+ " msg = \"Legacy fixed point iter did not converge\"\n",
+ " raise RuntimeError(msg)\n",
+ " inverse = next_inverse\n",
+ "\n",
+ " next_inverse = dsp.timeshift(\n",
+ " scet_wrt_tps, -inverse * physics_fs, interpolation_order\n",
+ " )\n",
+ " \n",
+ " error = np.max(np.abs((inverse - next_inverse)[edge:-edge]))\n",
+ " niter += 1\n",
+ "\n",
+ " return inverse"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "127215c1-5082-4994-8ef4-f6eb75dcbe11",
+ "metadata": {},
+ "source": [
+ "Setting up some example data"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "ed10795d-2d67-4715-9c18-a6e3993b2605",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "order = 31\n",
+ "nsamp = 3000\n",
+ "fsample = 16.0\n",
+ "dt = 1 / fsample\n",
+ "\n",
+ "f_mod = 0.007\n",
+ "a_mod = 1e-2 / (2 * np.pi * f_mod)\n",
+ "max_it = 5\n",
+ "tol = 1e-10"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "d28a7126-215b-49d5-83fd-b6d6783160fb",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "\n",
+ "def dx_from_x(x):\n",
+ " return np.cos(2 * np.pi * f_mod * x) * a_mod\n",
+ "\n",
+ "xi = np.arange(nsamp) * dt\n",
+ "dxi = dx_from_x(xi)\n"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "ccc111a8-3492-4569-a667-8731f60219ac",
+ "metadata": {},
+ "source": [
+ "Set up new shift inversion operator (unchunked pure numpy version)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "932daa39-b57d-4d97-b3c3-81471560c46a",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "op_np = make_shift_inverse_lagrange_numpy(\n",
+ " order=order,\n",
+ " fsample=fsample,\n",
+ " max_abs_shift=a_mod * 1.01,\n",
+ " max_iter=max_it,\n",
+ " tolerance=tol,\n",
+ " )\n",
+ "\n"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "4b15c610-c7f2-44ea-b8e6-7e7fc62c5bf5",
+ "metadata": {},
+ "source": [
+ "Compute shift inversion with both"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "5c3e3fea-5fdf-49ad-98d3-78062a04d801",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "ai_np = op_np(dxi)\n",
+ "ai_leg = legacy_invert_scet_wrt_tps(dxi, tol/100, fsample, order, max_it+2)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "fba2473e-3a8c-48c5-9cb5-13175098b70f",
+ "metadata": {},
+ "source": [
+ "Expected result"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "eb55600f-0cd2-4702-bd95-217d4722e248",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "ai_ex = dx_from_x(xi - ai_np)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "fc7a324c-4532-4c66-8317-3620d55e6d3c",
+ "metadata": {},
+ "source": [
+ "The original version seems to have more boundary artifacts"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "f8b72f4b-00bb-4f57-ae42-5a4293a03d65",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "plt.close()\n",
+ "plt.plot(xi, ai_np, \"k+\")\n",
+ "plt.plot(xi, ai_leg, \"rx\")\n",
+ "plt.plot(xi, ai_ex, \"g-\")\n"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "c665ea55-9d2b-4e8e-8ad9-e0113541633a",
+ "metadata": {},
+ "source": [
+ "The error within the valid points is within the tolerance for both."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "74d29668-57f2-4800-8e5b-bc82253a2185",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "plt.close()\n",
+ "plt.semilogy(xi, np.abs(ai_np-ai_ex), \"k+\")\n",
+ "plt.semilogy(xi, np.abs(ai_leg-ai_ex), \"rx\")\n",
+ "plt.axhline(y=tol)\n",
+ "plt.axvline(x=xi[op_np.margin_left])\n",
+ "plt.axvline(x=xi[-op_np.margin_right])\n"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "912c544e-e926-41cb-ba93-d414c538c712",
+ "metadata": {},
+ "source": [
+ "To investigate the differences inside the valid region, fix a bug in the legacy method. It used to return the second-best guess instead of the latest iteration result. Also, we make the margin size for the error measure a parameter."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "45775517-6b5b-4b78-a345-f51e5bba6439",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "def modified_invert_scet_wrt_tps(\n",
+ " scet_wrt_tps: np.ndarray,\n",
+ " clockinv_tolerance: float,\n",
+ " physics_fs: float,\n",
+ " edge : int,\n",
+ " interpolation_order: int,\n",
+ " clockinv_maxiter: int = 15,\n",
+ "):\n",
+ " \n",
+ " error = 0\n",
+ "\n",
+ " niter = 0\n",
+ " next_inverse = scet_wrt_tps\n",
+ " while not niter or error > clockinv_tolerance:\n",
+ " if niter >= clockinv_maxiter:\n",
+ " msg = \"Legacy fixed point iter did not converge\"\n",
+ " raise RuntimeError(msg)\n",
+ " inverse = next_inverse\n",
+ "\n",
+ " next_inverse = dsp.timeshift(\n",
+ " scet_wrt_tps, -inverse * physics_fs, interpolation_order\n",
+ " )\n",
+ " \n",
+ " error = np.max(np.abs((inverse - next_inverse)[edge:-edge]))\n",
+ " niter += 1\n",
+ " print( niter )\n",
+ " return next_inverse"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "f711b48f-fffa-4732-b609-7c2d500f59ff",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "#ai_mod = modified_invert_scet_wrt_tps(dxi, tol, fsample, op_np.margin_left, order)\n",
+ "ai_mod = modified_invert_scet_wrt_tps(dxi, tol, fsample, 100, order)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "5566a405-1994-4f25-905c-6ed6d2ee9f6a",
+ "metadata": {},
+ "source": [
+ "Now the error inside the valid region agrees perfectly (the margin size was not the issue)."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "cba3a269-ae9b-4236-87a3-be7966af0dc8",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "plt.close()\n",
+ "plt.semilogy(xi, np.abs(ai_np-ai_ex), \"k+\")\n",
+ "plt.semilogy(xi, np.abs(ai_mod-ai_ex), \"rx\")\n",
+ "plt.axhline(y=tol)\n",
+ "plt.axvline(x=xi[op_np.margin_left])\n",
+ "plt.axvline(x=xi[-op_np.margin_right])\n"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "9153fb52-c7de-4352-bf61-c1cda7cf7944",
+ "metadata": {},
+ "source": [
+ "The boundary artifacts are still larger for the legacy algorithm. Reason yet unkown."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "63689287-a4bf-4f23-957a-51f4046fd0ad",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "plt.close()\n",
+ "plt.plot(xi, ai_np, \"k+\")\n",
+ "plt.plot(xi, ai_mod, \"rx\")\n",
+ "plt.plot(xi, ai_ex, \"g-\")\n",
+ "plt.xlim(xi[-20], xi[-1]+dt)\n"
+ ]
+ }
+ ],
+ "metadata": {
+ "kernelspec": {
+ "display_name": "Python 3 (ipykernel)",
+ "language": "python",
+ "name": "python3"
+ },
+ "language_info": {
+ "codemirror_mode": {
+ "name": "ipython",
+ "version": 3
+ },
+ "file_extension": ".py",
+ "mimetype": "text/x-python",
+ "name": "python",
+ "nbconvert_exporter": "python",
+ "pygments_lexer": "ipython3"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
--
GitLab