BOTTS.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "44cd7798",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.style\n",
    "\n",
    "matplotlib.style.use([\"seaborn-v0_8-ticks\"])\n",
    "from matplotlib.ticker import EngFormatter, ScalarFormatter"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0e6a8b61",
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib widget"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "eb63415c",
   "metadata": {},
   "outputs": [],
   "source": [
    "from numpy.polynomial.polynomial import polyfit, polyval\n",
    "\n",
    "\n",
    "def detrend(x, w):\n",
    "    t = np.arange(-len(x) // 2, len(x) // 2)\n",
    "    return x - polyval(t, polyfit(t, x, 1, w=w)).T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5e0c5600",
   "metadata": {},
   "outputs": [],
   "source": [
    "opts = dict(alpha=1 / 25, ls=\"None\", marker=\".\", ms=2)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "83ddca5e-b791-4be7-86a0-fcb4f66e5a8e",
   "metadata": {},
   "source": [
    "# Data loading"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3698065c",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pathlib\n",
    "from ipywidgets import interact, Layout, Select\n",
    "\n",
    "folder = pathlib.Path(\"data\").resolve(strict=True)\n",
    "\n",
    "print(\"Selected folder: \\n\", folder)\n",
    "\n",
    "filename = None\n",
    "\n",
    "\n",
    "# This bit of vodoo shows a selection box in the page with all the files in the above folder\n",
    "# To pick the file, click.\n",
    "@interact\n",
    "def populate_filenames(\n",
    "    x=Select(\n",
    "        options=((x.name, x) for x in folder.iterdir() if x.suffix == \".parquetzstd\"),\n",
    "        description=\"Available files\",\n",
    "        layout=Layout(width=\"auto\"),\n",
    "        rows=20,  # number of rows visible in the box\n",
    "    )\n",
    "):\n",
    "    global filename\n",
    "    filename = x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "509574d9",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "print(filename)\n",
    "waves_df = pd.read_parquet(filename)\n",
    "display(waves_df.describe())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "39a89fd4",
   "metadata": {},
   "outputs": [],
   "source": [
    "import dataclasses\n",
    "\n",
    "\n",
    "@dataclasses.dataclass\n",
    "class TTSParms:\n",
    "    f_hi: float  # High cutoff frequency for TTS; Hz\n",
    "    f_lo: float  # Low cutoff frequency for TTS; Hz\n",
    "    t_ref: float  # Reference temperature; °C\n",
    "    f_ref: float  # Reference frequency; Hz\n",
    "    fs: float  # sampling frequency; Hz"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6e86e579",
   "metadata": {},
   "outputs": [],
   "source": [
    "import json\n",
    "\n",
    "tmp = waves_df.xs(key=\"Time\", axis=1, level=1)\n",
    "fs = (tmp.count() / (tmp.max() - tmp.min())).median()\n",
    "# parms = TTSParms(f_hi=14, f_lo=0.014, t_ref=55, f_ref=1, fs=fs)  # 154.56603773584905)\n",
    "# parms = TTSParms(f_hi=fs/2, f_lo=0, t_ref=30, f_ref=1, tukey=False, fs=fs)\n",
    "\n",
    "# json.dump(dataclasses.asdict(parms), filename.with_suffix(\".json\").open(mode=\"w\"))\n",
    "parms = TTSParms(**json.load(filename.with_suffix(\".json\").open()))\n",
    "# parms.f_hi=10\n",
    "# parms.f_lo=.01\n",
    "# parms.t_ref=175\n",
    "\n",
    "display(fs, parms)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eda79b74-35ca-45b5-986f-8385b3340f39",
   "metadata": {},
   "source": [
    "# Diagnosing individual chirps"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4d2352a1",
   "metadata": {},
   "outputs": [],
   "source": [
    "inx = 2\n",
    "eps = waves_df[(str(inx), \"Strain\")].values / 100\n",
    "sig = waves_df[(str(inx), \"Stress\")].values\n",
    "display(waves_df[str(inx), \"Temperature\"].median())\n",
    "\n",
    "t = np.arange(len(eps)) / parms.fs\n",
    "print(\"t\", t)\n",
    "trend_slice = (0.5 < t) & (t < 1) | (t[-1] - 0.5 < t) & (t < t[-1])\n",
    "fft_slice = (0.5 < t) & (t < t[-1])\n",
    "\n",
    "eps = detrend(eps, trend_slice)\n",
    "sig = detrend(sig, trend_slice)\n",
    "\n",
    "t = t[fft_slice]\n",
    "eps = eps[fft_slice]\n",
    "sig = sig[fft_slice]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5035c225",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1, 1, figsize=(6.25, 3), layout=\"constrained\")\n",
    "ax.plot(t, eps, color=\"tab:blue\", alpha=1)\n",
    "ax.yaxis.get_major_formatter().set_powerlimits([-3, 3])\n",
    "ax.set_ylabel(\"Strain (blue)\")\n",
    "ax.set_xlabel(\"Time (s)\")\n",
    "tax = ax.twinx()\n",
    "tax.yaxis.get_major_formatter().set_powerlimits([-3, 3])\n",
    "tax.plot(t, sig, color=\"tab:orange\", alpha=0.4)\n",
    "tax.set_ylabel(\"Stress (Pa, orange)\")\n",
    "\n",
    "None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "365965d3-218b-49a2-82e5-5850640c1499",
   "metadata": {},
   "outputs": [],
   "source": [
    "epsc = np.fft.rfft(eps)\n",
    "sigc = np.fft.rfft(sig)\n",
    "youngsc = sigc / epsc\n",
    "\n",
    "freq = np.fft.rfftfreq(len(t), 1 / parms.fs)\n",
    "keep = (parms.f_lo < freq) & (freq < parms.f_hi)\n",
    "# keep = (0<freq)\n",
    "display(keep)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "828b5f18",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axs = plt.subplots(2, 1, layout=\"constrained\", figsize=(4, 5), sharex=True)\n",
    "axs[0].loglog(freq[keep], np.abs(epsc)[keep])  # ,**opts)\n",
    "axs[1].loglog(freq[keep], np.abs(sigc)[keep], color=\"tab:orange\")  # , **opts)\n",
    "axs[1].set_xlabel(\"Frequency (Hz)\")\n",
    "axs[1].set_ylabel(\"Complex\\nstress magnitude (Pa)\")\n",
    "axs[0].set_ylabel(\"Complex\\nstrain magnitude\")\n",
    "None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7d2d917a",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "plt.loglog(freq[keep], youngsc.real[keep])  # ,**opts)\n",
    "plt.loglog(freq[keep], youngsc.imag[keep])  # ,**opts)\n",
    "plt.xlabel(\"Frequency (Hz)\")\n",
    "plt.ylabel(\"Complex modulus (Pa)\")\n",
    "plt.legend([\"Storage\", \"Loss\"])\n",
    "None"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2cf6e9a6-c929-441a-bc8c-f74b3743ae96",
   "metadata": {},
   "source": [
    "# Analyze whole chirp series"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7429c05d",
   "metadata": {},
   "outputs": [],
   "source": [
    "freq = np.fft.rfftfreq(len(t), 1 / parms.fs)\n",
    "epsc = np.fft.rfft(\n",
    "    0.01\n",
    "    * detrend(waves_df.xs(key=\"Strain\", level=1, axis=\"columns\").values, trend_slice)[\n",
    "        fft_slice\n",
    "    ],\n",
    "    axis=0,\n",
    ")\n",
    "sigc = np.fft.rfft(\n",
    "    detrend(waves_df.xs(key=\"Stress\", level=1, axis=\"columns\").values, trend_slice)[\n",
    "        fft_slice\n",
    "    ],\n",
    "    axis=0,\n",
    ")\n",
    "\n",
    "print(sigc.shape)\n",
    "noise_keep = (parms.f_hi < freq) & (freq < parms.fs / 2)\n",
    "# noise_keep = (30<freq) & (freq<60)\n",
    "eps_std = np.nanmedian(np.abs(epsc)[noise_keep], axis=0)\n",
    "sig_std = np.nanmedian(np.abs(sigc)[noise_keep], axis=0)\n",
    "print(eps_std)\n",
    "print(sig_std)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d088c42e",
   "metadata": {},
   "outputs": [],
   "source": [
    "with plt.style.context(\"fast\"):\n",
    "    plt.figure()\n",
    "    plt.loglog(freq, np.abs(epsc) / eps_std, **opts, color=\"blue\")\n",
    "    plt.loglog(freq, np.abs(sigc) / sig_std, **opts, color=\"orange\")\n",
    "None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b2060f9e",
   "metadata": {},
   "outputs": [],
   "source": [
    "epsc = epsc[keep]\n",
    "sigc = sigc[keep]\n",
    "freq = freq[keep]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4cf81997",
   "metadata": {},
   "outputs": [],
   "source": [
    "display(freq.max(), freq.min())\n",
    "curve_slice = np.s_[0 : sigc.shape[-1]]\n",
    "curve_space = np.arange(curve_slice.start, curve_slice.stop)\n",
    "temp_space = (\n",
    "    np.median(\n",
    "        waves_df.xs(\"Temperature\", level=1, axis=\"columns\").values[:, curve_slice],\n",
    "        axis=0,\n",
    "    )\n",
    "    + 273.15\n",
    ")\n",
    "display(temp_space)\n",
    "assert temp_space.min() - 1 <= parms.t_ref + 273.15 <= temp_space.max() + 1, (\n",
    "    temp_space.min(),\n",
    "    parms.t_ref + 273.15,\n",
    "    temp_space.max(),\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "76dea726",
   "metadata": {},
   "outputs": [],
   "source": [
    "youngsc = sigc / epsc\n",
    "youngs_real = youngsc.real\n",
    "youngs_imag = youngsc.imag\n",
    "youngs_mod = np.abs(youngsc)\n",
    "youngsc_std = np.sqrt(\n",
    "    (sig_std / np.abs(epsc)) ** 2 + (eps_std * np.abs(youngsc) / np.abs(epsc)) ** 2\n",
    ")\n",
    "youngs_std = youngsc_std\n",
    "display(youngs_std.shape)\n",
    "phase = np.angle(youngsc)\n",
    "tand = youngsc.imag / youngsc.real\n",
    "tand_std = np.sqrt(youngsc.real**-2 + (tand / youngsc.real) ** 2) * youngsc_std"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6302b0c5",
   "metadata": {},
   "outputs": [],
   "source": [
    "viridis_colors = plt.cm.viridis(np.linspace(0, 1, sigc.shape[1], endpoint=True))\n",
    "my_cmap = plt.cm.colors.ListedColormap(viridis_colors)\n",
    "my_sm = plt.cm.ScalarMappable(\n",
    "    cmap=my_cmap,\n",
    "    norm=plt.Normalize(vmin=min(temp_space) - 273.16, vmax=max(temp_space) - 273.15),\n",
    ")\n",
    "with plt.style.context(\"fast\"):\n",
    "    fig, axs = plt.subplots(\n",
    "        1, 2, figsize=(6.25, 3), layout=\"constrained\", sharex=True, sharey=True\n",
    "    )\n",
    "    for ax, thing, name in zip(axs, (youngs_real, youngs_imag), (\"Storage\", \"Loss\")):\n",
    "        ax.set_prop_cycle(\n",
    "            \"color\",\n",
    "            viridis_colors,\n",
    "        )\n",
    "        ax.loglog(\n",
    "            freq,\n",
    "            thing[:, curve_slice],\n",
    "            alpha=1,\n",
    "            marker=\".\",\n",
    "            markersize=7.5,\n",
    "            fillstyle=\"full\",\n",
    "            markeredgewidth=0,\n",
    "            linestyle=\"none\",\n",
    "            label=(\n",
    "                np.round(\n",
    "                    np.median(\n",
    "                        waves_df.xs(\"Temperature\", level=1, axis=\"columns\").values[\n",
    "                            :, curve_slice\n",
    "                        ],\n",
    "                        axis=0,\n",
    "                    )\n",
    "                ).astype(int)\n",
    "                if name == \"Storage\"\n",
    "                else None\n",
    "            ),\n",
    "        )\n",
    "        ax.set_title(name)\n",
    "    # axs[0].set_ylim(ymin=10**5,)\n",
    "    axs[0].set_ylabel(\"Young's Modulus (Pa)\")\n",
    "    axs[0].set_xlabel(\"Frequency (Hz)\")\n",
    "    axs[1].set_xlabel(\"Frequency (Hz)\")\n",
    "    # fig.legend(title=\"Temp. (°C)\", loc=\"outside right\")\n",
    "    plt.colorbar(my_sm, ax=axs[1], label=\"Temp. (°C)\")\n",
    "None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9dd7c558-0af7-4ef1-a9e7-986ed723d29f",
   "metadata": {},
   "outputs": [],
   "source": [
    "with plt.style.context(\"fast\"):\n",
    "    fig, axs = plt.subplots(\n",
    "        1, 2, figsize=(6.25, 3), layout=\"constrained\", sharex=True, sharey=False\n",
    "    )\n",
    "    for ax, thing, name, attr in zip(\n",
    "        axs,\n",
    "        (youngs_mod, phase * 180 / np.pi),\n",
    "        (\"mag\", \"phase\"),\n",
    "        (\"semilogx\", \"semilogx\"),\n",
    "    ):\n",
    "        ax.set_prop_cycle(\n",
    "            \"color\",\n",
    "            viridis_colors,\n",
    "        )\n",
    "        attr = getattr(ax, attr)\n",
    "        attr(\n",
    "            freq,\n",
    "            thing,\n",
    "            \".\",\n",
    "            label=(\n",
    "                np.round(\n",
    "                    np.median(\n",
    "                        waves_df.xs(\"Temperature\", level=1, axis=\"columns\").values[\n",
    "                            :, curve_slice\n",
    "                        ],\n",
    "                        axis=0,\n",
    "                    )\n",
    "                ).astype(int)\n",
    "                if name == \"mag\"\n",
    "                else None\n",
    "            ),\n",
    "        )\n",
    "        ax.set_title(name)\n",
    "    axs[0].set_yscale(\"log\")\n",
    "    axs[0].set_ylabel(\"Young's Modulus (Pa)\")\n",
    "    axs[0].set_xlabel(\"Frequency (Hz)\")\n",
    "    axs[1].set_xlabel(\"Frequency (Hz)\")\n",
    "    # fig.legend(title=\"Temp. (°C)\", loc=\"outside right\")\n",
    "    plt.colorbar(my_sm, ax=axs[1], label=\"Temp. (°C)\")\n",
    "None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fa7bf57c",
   "metadata": {},
   "outputs": [],
   "source": [
    "with plt.style.context(\"fast\"):\n",
    "    stuff = [youngs_real / youngs_std, youngs_imag / youngs_std]\n",
    "    stuffnames = [\"Storage\", \"Loss\"]\n",
    "    fig, axs = plt.subplots(\n",
    "        1, len(stuff), figsize=(6.5, 3), layout=\"constrained\", sharex=True, sharey=True\n",
    "    )\n",
    "    for i, (ax, thing, axname) in enumerate(zip(axs, stuff, stuffnames)):\n",
    "        ax.set_prop_cycle(\n",
    "            \"color\",\n",
    "            viridis_colors,\n",
    "        )\n",
    "        ax.loglog(\n",
    "            freq,\n",
    "            thing[:, curve_slice],\n",
    "            # alpha=0.25,\n",
    "            marker=\".\",\n",
    "            markersize=5,\n",
    "            fillstyle=\"full\",\n",
    "            markeredgewidth=0,\n",
    "            linestyle=\"none\",\n",
    "        )\n",
    "        ax.set_title(axname)\n",
    "        ax.set_xlabel(\"Frequency (Hz)\")\n",
    "    # axs[0].set_ylim(ymin=.6, ymax=6000)\n",
    "    axs[0].set_ylabel(\"Young's Modulus SNR\")\n",
    "    plt.colorbar(my_sm, ax=axs[-1], label=\"Temperature (°C)\")\n",
    "None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e8f11cb5-06d0-49f0-ad3b-ca3a8c9d03a4",
   "metadata": {},
   "outputs": [],
   "source": [
    "with plt.style.context(\"fast\"):\n",
    "    stuff = [np.abs(epsc), youngs_std]\n",
    "    stuffnames = [\"Magnitude\", \"Error\"]\n",
    "    fig, axs = plt.subplots(\n",
    "        1, len(stuff), figsize=(6.5, 3), layout=\"constrained\", sharex=True, sharey=False\n",
    "    )\n",
    "    for i, (ax, foo, axname) in enumerate(zip(axs, stuff, stuffnames)):\n",
    "        ax.set_prop_cycle(\n",
    "            \"color\",\n",
    "            [\n",
    "                plt.cm.viridis(i)\n",
    "                for i in np.linspace(0, 1, sigc.shape[1], endpoint=True)\n",
    "            ],\n",
    "        )\n",
    "        ax.loglog(\n",
    "            freq,\n",
    "            foo[:, curve_slice],\n",
    "            # alpha=0.25,\n",
    "            marker=\".\",\n",
    "            markersize=5,\n",
    "            fillstyle=\"full\",\n",
    "            markeredgewidth=0,\n",
    "            linestyle=\"none\",\n",
    "            label=(\n",
    "                np.round(\n",
    "                    np.median(\n",
    "                        waves_df.xs(\"Temperature\", level=1, axis=\"columns\").values[\n",
    "                            :, curve_slice\n",
    "                        ],\n",
    "                        axis=0,\n",
    "                    )\n",
    "                ).astype(int)\n",
    "                if not i\n",
    "                else None\n",
    "            ),\n",
    "        )\n",
    "        ax.set_title(axname)\n",
    "        ax.set_xlabel(\"Frequency (Hz)\")\n",
    "    axs[0].set_ylabel(\"Complex Strain\")\n",
    "    axs[1].set_ylabel(\"Young's Modulus (Pa)\")\n",
    "    plt.colorbar(my_sm, ax=axs[1], label=\"Temp. (°C)\")\n",
    "None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f193ac8b",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(4 * 0.8, 3 * 0.8), layout=\"constrained\")\n",
    "i = 5\n",
    "plt.loglog(\n",
    "    freq,\n",
    "    youngs_imag[:, i],\n",
    "    marker=\"o\",\n",
    "    markersize=3,\n",
    "    fillstyle=\"full\",\n",
    "    markeredgewidth=0,\n",
    "    linestyle=\"none\",\n",
    "    label=\"55 °C\",\n",
    ")\n",
    "plt.loglog(\n",
    "    freq,\n",
    "    youngs_imag[:, i + 1],\n",
    "    marker=\"s\",\n",
    "    markersize=3,\n",
    "    fillstyle=\"full\",\n",
    "    markeredgewidth=0,\n",
    "    linestyle=\"none\",\n",
    "    label=\"60 °C\",\n",
    ")\n",
    "plt.loglog(\n",
    "    freq * 0.11,\n",
    "    youngs_imag[:, i + 1],\n",
    "    alpha=0.25,\n",
    "    marker=\"s\",\n",
    "    markersize=3,\n",
    "    fillstyle=\"full\",\n",
    "    markeredgewidth=0,\n",
    "    linestyle=\"none\",\n",
    "    label=\"60 °C(shifted)\",\n",
    ")\n",
    "plt.xlabel(\"Frequency (Hz)\")\n",
    "plt.ylabel(\"Storage Young's Modulus (Pa)\")\n",
    "plt.legend(loc=\"upper left\")\n",
    "None"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7a1887c8-569b-4ec3-93da-aa052fcae4d2",
   "metadata": {},
   "source": [
    "# Pairwise polynomial shift optimization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1e117ab4-c83f-42b0-aa4a-5608a2a1e3a7",
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.optimize import minimize_scalar\n",
    "from sklearn.pipeline import Pipeline\n",
    "from sklearn.linear_model import LinearRegression\n",
    "from sklearn.preprocessing import FunctionTransformer, PolynomialFeatures\n",
    "from itertools import repeat, starmap, tee, islice, chain\n",
    "from tqdm.notebook import tqdm\n",
    "import time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0d6a5120-ec69-4a33-aa40-c9f4f7e9f5f2",
   "metadata": {},
   "outputs": [],
   "source": [
    "POLY_REGRESSOR = lambda: Pipeline(\n",
    "    [\n",
    "        (\"logscale\", FunctionTransformer(func=np.log, inverse_func=np.exp)),\n",
    "        (\"polyfeat\", PolynomialFeatures(8)),\n",
    "        (\"linearfit\", LinearRegression(fit_intercept=False)),\n",
    "    ]\n",
    ")\n",
    "\n",
    "\n",
    "def poly_fit(freq, y, y_std):\n",
    "    lr_real = POLY_REGRESSOR()\n",
    "    lr_imag = POLY_REGRESSOR()\n",
    "    freq = freq.reshape(-1, 1)\n",
    "    y = y.reshape(2, -1)\n",
    "    y_std = y_std.reshape(2, -1)\n",
    "    lr_real.fit(freq, y[0], linearfit__sample_weight=y_std[0] ** -2)\n",
    "    lr_imag.fit(freq, y[1], linearfit__sample_weight=y_std[1] ** -2)\n",
    "\n",
    "    class D:\n",
    "        @staticmethod\n",
    "        def predict(basis):\n",
    "            basis = basis.reshape(2, -1)\n",
    "            m1 = lr_real.predict(basis[0][:, None])\n",
    "            m2 = lr_imag.predict(basis[1][:, None])\n",
    "            return np.concatenate((m1, m2))\n",
    "\n",
    "    return D, np.concatenate((freq, freq))\n",
    "\n",
    "\n",
    "def weighted_loss(m, y, ystd):\n",
    "    return (((y - m) / ystd) ** 2).mean()\n",
    "\n",
    "\n",
    "def pairwise_shift(freq1, yc1real, yc1imag, ystd1, freq2, yc2real, yc2imag, ystd2):\n",
    "    def hobjective(shift):\n",
    "        freq = np.concatenate((freq1, freq2 * shift))\n",
    "        lr, basis = poly_fit(freq, y, ystd)\n",
    "        m = lr.predict(basis)\n",
    "        return weighted_loss(m, y, ystd)\n",
    "\n",
    "    def vobjective(shift):\n",
    "        y = np.concatenate((yc1real, shift * yc2real, yc1imag, shift * yc2imag))\n",
    "        ystd = np.concatenate((ystd1, shift * ystd2, ystd1, shift * ystd2))\n",
    "        freq = np.concatenate((freq1, freq2 * best_hshift))\n",
    "        lr, basis = poly_fit(freq, y, ystd * shift)\n",
    "        m = lr.predict(basis)\n",
    "\n",
    "        return weighted_loss(m, y, ystd)\n",
    "\n",
    "    y = np.concatenate((yc1real, yc2real, yc1imag, yc2imag))\n",
    "    ystd = np.concatenate((ystd1, ystd2, ystd1, ystd2))\n",
    "    best_hshift = minimize_scalar(\n",
    "        hobjective,\n",
    "        bounds=(freq.min() / freq.max() / 2, 1),\n",
    "        # bounds=(0.001, 0.003),\n",
    "        options=dict(disp=0, xatol=0),\n",
    "    ).x\n",
    "    best_vshift = 1\n",
    "    # Maybe you want to try optimizing vertical shifts?\n",
    "    # best_vshift = minimize_scalar(vobjective, bracket=(0.5, 2), options=dict(disp=0)).x\n",
    "    # y *= best_vshift\n",
    "    # best_hshift = minimize_scalar(\n",
    "    #     hobjective,\n",
    "    #     bracket=(freq.min() / freq.max(),best_hshift, 0.75),\n",
    "    #     options=dict(disp=0),\n",
    "    # ).x\n",
    "    return best_hshift, best_vshift"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bcc5efb1",
   "metadata": {},
   "outputs": [],
   "source": [
    "# %%time\n",
    "inx = 2\n",
    "display(\n",
    "    pairwise_shift(\n",
    "        freq,\n",
    "        youngs_real[:, inx],\n",
    "        youngs_imag[:, inx],\n",
    "        youngs_std[:, inx],\n",
    "        freq,\n",
    "        youngs_real[:, inx + 1],\n",
    "        youngs_imag[:, inx + 1],\n",
    "        youngs_std[:, inx + 1],\n",
    "    )\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7793f465",
   "metadata": {},
   "outputs": [],
   "source": [
    "ratios = list(\n",
    "    tqdm(\n",
    "        map(\n",
    "            pairwise_shift,\n",
    "            repeat(freq),\n",
    "            youngs_real.T[:-1],\n",
    "            youngs_imag.T[:-1],\n",
    "            youngs_std.T[:-1],\n",
    "            repeat(freq),\n",
    "            youngs_real.T[1:],\n",
    "            youngs_imag.T[1:],\n",
    "            youngs_std.T[1:],\n",
    "        ),\n",
    "        total=len(temp_space) - 1,\n",
    "    )\n",
    ")\n",
    "ratios.insert(0, (1, 1))\n",
    "ratios, vratios = _ = np.array(ratios).T\n",
    "shifts, vshifts = np.cumprod(_, axis=1)\n",
    "ref_shift = np.interp(273.15 + parms.t_ref, temp_space, shifts)\n",
    "print(ratios.tolist())  # ,vratios)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c0438b6e-5d12-467d-b8ae-5cdc1a3412bb",
   "metadata": {},
   "outputs": [],
   "source": [
    "# PT_BOTTS_pairwise_shifts = \"o\",\"BOTTS-Pair\",np.array([1.0, 0.7667005144691353, 0.04331814368918157, 0.0024594764431351353, 0.00011887261323194619, 7.592370043280243e-06, 7.882068570016968e-07, 1.3090807730023502e-07, 3.091946428940733e-08, 9.354690070567252e-09, 3.3875285494582026e-09, 1.6798476171517379e-09, 1.2549269210204607e-09])\n",
    "PT_BOTTS_pairwise_shifts = \"o\", \"BOTTS-Pair\", shifts\n",
    "PT_DFS_pairwise_shifts = (\n",
    "    \"+\",\n",
    "    \"DFS-TTS-Pair\",\n",
    "    np.array(\n",
    "        [\n",
    "            1.0,\n",
    "            0.1675308142817378,\n",
    "            0.005828201539792645,\n",
    "            0.00017069244815959546,\n",
    "            7.245146331674599e-06,\n",
    "            4.845450761236291e-07,\n",
    "            4.9780265356208074e-08,\n",
    "            7.859503721864011e-09,\n",
    "            1.7282219624909232e-09,\n",
    "            4.796186278514707e-10,\n",
    "            1.5745211187321228e-10,\n",
    "            6.087406477337344e-11,\n",
    "            2.8442653780748484e-11,\n",
    "        ]\n",
    "    ),\n",
    ")\n",
    "PT_trios_shifts = (\n",
    "    \"x\",\n",
    "    \"DFS-TTS-TRIOS\",\n",
    "    np.array(\n",
    "        [\n",
    "            1.00000,\n",
    "            0.455918,\n",
    "            0.0497035,\n",
    "            1.97562e-3,\n",
    "            9.04780e-5,\n",
    "            6.10305e-6,\n",
    "            6.63254e-7,\n",
    "            1.11122e-7,\n",
    "            2.41008e-8,\n",
    "            6.06065e-9,\n",
    "            1.67503e-9,\n",
    "            2.62684e-10,\n",
    "            1.91396e-11,\n",
    "        ]\n",
    "    ),\n",
    ")\n",
    "\n",
    "fig, ax = plt.subplots(1, 1, layout=\"constrained\", figsize=(4 * 0.8, 3 * 0.8))\n",
    "\n",
    "for symbol, label, this_shift in (\n",
    "    PT_BOTTS_pairwise_shifts,\n",
    "    PT_DFS_pairwise_shifts,\n",
    "    PT_trios_shifts,\n",
    "):\n",
    "    this_ref_shift = np.interp(273.15 + parms.t_ref, temp_space, this_shift)\n",
    "    ax.semilogy(\n",
    "        temp_space - 273.15,\n",
    "        this_shift / this_ref_shift,\n",
    "        linestyle=\"none\",\n",
    "        marker=symbol,\n",
    "        label=label,\n",
    "        markerfacecolor=\"none\",\n",
    "    )\n",
    "ax.plot([parms.t_ref] * 2, [0, 1], linestyle=\":\", color=\"black\")\n",
    "ax.text(\n",
    "    parms.t_ref - 4,\n",
    "    5e-6,\n",
    "    f\"$T_{{ref}}$ = {parms.t_ref:.0f} °C\",\n",
    "    rotation=\"vertical\",\n",
    ")  # horizontalalignment=\"right\")\n",
    "ax.set_ylabel(\"$a_T$\")\n",
    "ax.set_xlabel(\"Temperature (°C)\")\n",
    "ax.legend(loc=\"upper right\")\n",
    "None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "294c49ad",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axs = plt.subplots(2, 1, layout=\"constrained\", figsize=(4, 6))\n",
    "axs[0].semilogy(temp_space - 273.15, shifts, \"-+\")\n",
    "axs[0].set_ylabel(\"$a_T$\")\n",
    "axs[1].plot(temp_space - 273.15, vshifts, \"-+\")\n",
    "axs[1].set_ylabel(\"$b_T$\")\n",
    "axs[1].set_xlabel(\"Temperature (°C)\")\n",
    "None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f20dcb72",
   "metadata": {},
   "outputs": [],
   "source": [
    "# arrhenius fit\n",
    "p = np.polynomial.polynomial.polyfit(1 / temp_space, np.log10(shifts), deg=1)\n",
    "display(p)\n",
    "fit = np.polynomial.polynomial.polyval(1 / temp_space, p)\n",
    "plt.figure()\n",
    "plt.plot(\n",
    "    1 / temp_space,\n",
    "    np.log10(shifts),\n",
    ")\n",
    "plt.plot(\n",
    "    1 / temp_space,\n",
    "    fit,\n",
    ")\n",
    "None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e944e356",
   "metadata": {},
   "outputs": [],
   "source": [
    "# wlf fit\n",
    "from scipy.optimize import curve_fit\n",
    "\n",
    "\n",
    "def wlf(t, c1, c2, t0):\n",
    "    return -c1 * (t - t0) / (c2 + t - t0)\n",
    "\n",
    "\n",
    "def wlfinv(loga, c1, c2, t0):\n",
    "    return t0 - loga * c2 / (loga + c1)\n",
    "\n",
    "\n",
    "fit_slice = np.s_[:]\n",
    "beta, cov = curve_fit(\n",
    "    wlf,\n",
    "    temp_space[fit_slice],\n",
    "    np.log10(shifts[fit_slice]),\n",
    "    p0=[10, 100, 273],\n",
    "    bounds=(0, np.inf),\n",
    ")\n",
    "plt.figure()\n",
    "plt.plot(\n",
    "    temp_space,\n",
    "    np.log10(shifts),\n",
    ")\n",
    "plt.plot(\n",
    "    temp_space,\n",
    "    wlf(temp_space, *beta),\n",
    ")\n",
    "print(beta)\n",
    "None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3a1c7740",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axs = plt.subplots(1, 2, figsize=(9.5, 4), layout=\"constrained\")\n",
    "for ax, thing in zip(\n",
    "    axs, (youngs_real * vshifts[None, :], youngs_imag * vshifts[None, :])\n",
    "):\n",
    "    ax.set_prop_cycle(\n",
    "        \"color\",\n",
    "        [\n",
    "            plt.cm.viridis(i)\n",
    "            for i in np.linspace(1, 0, len(curve_space) + 1, endpoint=False)[1:]\n",
    "        ],\n",
    "    )\n",
    "    ax.loglog(np.outer(freq, shifts / ref_shift), thing[:, curve_space], **opts)\n",
    "None"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2b38abaf-0ec0-469c-bf5d-ac895faaac80",
   "metadata": {},
   "source": [
    "# Cross-validated prony fit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "597a8d37-498a-4d3a-8bd4-4e2d4e7e9aa0",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.linear_model import (\n",
    "    LinearRegression,\n",
    "    ElasticNet,\n",
    "    ElasticNetCV,\n",
    "    Ridge,\n",
    "    RidgeCV,\n",
    "    Lasso,\n",
    "    LassoCV,\n",
    ")\n",
    "from sklearn.model_selection import RepeatedKFold"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b641f716-c2c5-4e2b-b0d2-37a1b211b0bb",
   "metadata": {},
   "outputs": [],
   "source": [
    "# %%time\n",
    "SOLID = True\n",
    "RELAXATIONS_PER_DECADE = 2\n",
    "SEED = 6\n",
    "# PRONY_REGRESSOR = lambda: LinearRegression(positive=True, fit_intercept=False)\n",
    "# PRONY_REGRESSOR = lambda: ElasticNet(\n",
    "#     l1_ratio=0.1, alpha=1.2e-6, positive=True, fit_intercept=False, max_iter=10000\n",
    "# )\n",
    "PRONY_REGRESSOR = lambda: ElasticNetCV(\n",
    "    cv=RepeatedKFold(n_repeats=3, random_state=SEED),\n",
    "    l1_ratio=[0.02, 0.05, 0.1, 0.3, 0.5, 0.7, 0.9, 0.95, 0.98][0],\n",
    "    n_alphas=20,\n",
    "    eps=1e-30,\n",
    "    positive=True,\n",
    "    fit_intercept=False,\n",
    "    max_iter=10000,\n",
    "    selection=\"random\",\n",
    "    random_state=4294967295 - SEED,\n",
    ")\n",
    "# REGRESSOR = lambda: LassoCV(n_alphas=100, eps=1e-6, positive=True, fit_intercept=False)\n",
    "# REGRESSOR = lambda: Lasso(alpha=1e1, positive=True, fit_intercept=True)\n",
    "\n",
    "\n",
    "def prony_relaxation_space(freq):\n",
    "    log_fmax, log_fmin = np.log10(np.max(freq)) + 2, np.log10(np.min(freq)) + (\n",
    "        2 if SOLID else -2\n",
    "    )\n",
    "    return np.logspace(\n",
    "        -log_fmin, -log_fmax, int((log_fmax - log_fmin) * RELAXATIONS_PER_DECADE)\n",
    "    )\n",
    "\n",
    "\n",
    "def prony_basis(freq, relaxations):\n",
    "    dt = dimensionless_time = np.outer(freq, relaxations)\n",
    "\n",
    "    dt2 = dt * dt\n",
    "    dt2p1 = dt2 + 1\n",
    "    ep_basis = dt2 / (dt2p1)\n",
    "    epp_basis = dt / (dt2p1)\n",
    "\n",
    "    if SOLID:\n",
    "        ep_basis = np.concatenate(\n",
    "            (ep_basis, np.ones_like(ep_basis, shape=(len(dt), 1))), axis=1\n",
    "        )\n",
    "        epp_basis = np.concatenate(\n",
    "            (epp_basis, np.zeros_like(epp_basis, shape=(len(dt), 1))), axis=1\n",
    "        )\n",
    "    return np.concatenate((ep_basis, epp_basis), axis=0)\n",
    "\n",
    "\n",
    "def prony_fit(freq, y, y_std):\n",
    "    relaxations = prony_relaxation_space(freq)\n",
    "    basis = prony_basis(freq, relaxations)\n",
    "    lr = PRONY_REGRESSOR()\n",
    "    lr.fit(\n",
    "        basis,\n",
    "        y,  # assume already concatenated # np.concatenate((y.real, y.imag))\n",
    "        y_std**-2,  # assume already concatenated # np.concatenate((y_std, y_std))\n",
    "    )\n",
    "    return lr, basis\n",
    "\n",
    "\n",
    "reduced_freq = np.outer(freq, shifts / ref_shift)\n",
    "relaxations = prony_relaxation_space(reduced_freq)\n",
    "\n",
    "args = (\n",
    "    np.concatenate(((youngs_real * vshifts).ravel(), (youngs_imag * vshifts).ravel())),\n",
    "    np.concatenate(\n",
    "        ((youngs_std * vshifts).ravel(), (youngs_std * vshifts).ravel()), axis=-1\n",
    "    ),\n",
    ")\n",
    "\n",
    "import warnings\n",
    "\n",
    "with warnings.catch_warnings():\n",
    "    from sklearn.exceptions import ConvergenceWarning\n",
    "\n",
    "    warnings.filterwarnings(\"ignore\", category=ConvergenceWarning)\n",
    "    lr, basis = prony_fit(reduced_freq, *args)\n",
    "    display(lr)\n",
    "    display(lr.alpha_, lr.l1_ratio_, lr.n_iter_, lr.mse_path_.mean(axis=-1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "be9673bc-989b-4a8a-bf6a-dfa550ad5093",
   "metadata": {},
   "outputs": [],
   "source": [
    "# prony coefficients\n",
    "plt.figure()\n",
    "\n",
    "g_inf = lr.coef_[-1]\n",
    "g = lr.coef_[:-1]\n",
    "display(np.sum(g > 0))\n",
    "plt.loglog(relaxations[g > 0], g[g > 0], \".\")\n",
    "None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8a366d67-9d38-401d-9cd9-59c18f561d77",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Alfrey's rule\n",
    "dense_relaxation_space = np.logspace(\n",
    "    -np.log10(np.max(reduced_freq)) + 2 * 0,\n",
    "    -np.log10(np.min(reduced_freq)) - 2 * 0,\n",
    "    500,\n",
    ")\n",
    "dimensionless_time = np.outer(dense_relaxation_space, 1 / relaxations[g > 0])\n",
    "H = np.sum(g[g > 0] * dimensionless_time * np.exp(-dimensionless_time), axis=1)\n",
    "plt.figure()\n",
    "plt.loglog(dense_relaxation_space, H)\n",
    "plt.ylim(bottom=1)\n",
    "plt.xlabel(\"relaxation time\")\n",
    "plt.ylabel(\"relaxation spectrum magnitude\")\n",
    "None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c4d40780-3645-413d-a59f-7adf152d0ca1",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(layout=\"constrained\", figsize=(4 * 0.8, 3 * 0.8))\n",
    "dlines = plt.loglog(\n",
    "    reduced_freq.ravel(),\n",
    "    (args[0]).reshape((2, -1)).T,\n",
    "    \".\",\n",
    "    label=[\"Data Storage\", \"Data Loss\"],\n",
    "    zorder=0,\n",
    "    **opts,\n",
    ")\n",
    "\n",
    "sortind = np.argsort(reduced_freq.ravel())\n",
    "mle = lr.predict(basis).reshape(2, -1)\n",
    "plines = plt.loglog(\n",
    "    reduced_freq.ravel()[sortind],\n",
    "    mle.T[sortind],\n",
    "    \"k\",\n",
    "    linewidth=1,\n",
    "    label=[\"Prony Storage\", \"Prony Loss\"],\n",
    ")\n",
    "\n",
    "plt.ylim(bottom=1e6, top=5e9)\n",
    "\n",
    "plt.xlabel(f\"Frequency at $T_{{ref}}$ = {parms.t_ref} °C (Hz)\")\n",
    "plt.ylabel(\"Young's modulus (Pa)\")\n",
    "plt.legend([\"Storage\", \"Loss\", \"Prony\"], loc=\"upper left\")\n",
    "None"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b358abbf-4c4c-46cd-a209-c5e5ecb5fe15",
   "metadata": {},
   "source": [
    "# Bootstrap prony fit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "25aba912-f0e5-4e39-af79-06143b35752f",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "def boot(*args):\n",
    "    for items in zip(*args):\n",
    "        itemstack = np.stack(items, axis=-1)\n",
    "        yield rng.choice(itemstack, len(itemstack)).T\n",
    "\n",
    "\n",
    "def recombine(*iters):\n",
    "    for item in zip(*iters):\n",
    "        yield (chain.from_iterable(item))\n",
    "\n",
    "\n",
    "NBOOTS = 500  # you can still plot all the completed samples if you interrupt the cell\n",
    "# SEED = 5\n",
    "PRONY_REGRESSOR = lambda: ElasticNet(\n",
    "    l1_ratio=lr.l1_ratio_,\n",
    "    alpha=lr.alpha_,\n",
    "    positive=True,\n",
    "    fit_intercept=False,\n",
    "    max_iter=10000,\n",
    "    selection=\"random\",\n",
    "    random_state=SEED,\n",
    ")\n",
    "# PRONY_REGRESSOR = lambda: LinearRegression(positive=True, fit_intercept=False)\n",
    "\n",
    "rlr = PRONY_REGRESSOR()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cbb04a97-b72d-4d39-8404-aeb1544ad857",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "boot_coefs = np.full((NBOOTS, len(relaxations) + 1), np.nan)\n",
    "rng = np.random.default_rng(SEED)\n",
    "\n",
    "for r in tqdm(range(NBOOTS)):\n",
    "    SOLID = False\n",
    "    boot_example = list(\n",
    "        boot(\n",
    "            repeat(freq),\n",
    "            youngs_real.T,\n",
    "            youngs_imag.T,\n",
    "            youngs_std.T,\n",
    "        )\n",
    "    )\n",
    "    rratios = list(\n",
    "        starmap(\n",
    "            pairwise_shift,\n",
    "            recombine(boot_example, boot_example[1:]),\n",
    "        )\n",
    "    )\n",
    "\n",
    "    rratios.insert(0, (1, 1))\n",
    "    rratios, rvratios = _ = np.array(rratios).T\n",
    "    rshifts, rvshifts = np.cumprod(_, axis=1)\n",
    "\n",
    "    ### WARNING: NO VSHIFTS\n",
    "    ry = np.concatenate([x[1] for x in boot_example] + [x[2] for x in boot_example])\n",
    "    ry_std = np.concatenate([x[3] for x in boot_example] + [x[3] for x in boot_example])\n",
    "\n",
    "    SOLID = True\n",
    "    rref_shift = np.interp(273.15 + parms.t_ref, temp_space, rshifts)\n",
    "    rreduced_freq = np.multiply(\n",
    "        np.array([x[0] for x in boot_example]), rshifts[:, None] / rref_shift\n",
    "    )\n",
    "\n",
    "    # lr, rbasis = prony_fit(rreduced_freq, *args)  # can't use, unwrap below\n",
    "\n",
    "    # must ensure all coefs correspond to global relaxations\n",
    "    # relaxations = prony_relaxation_space(reduced_freq) # precalculated in cross-validation cell\n",
    "\n",
    "    rbasis = prony_basis(rreduced_freq, relaxations)\n",
    "    rlr.fit(rbasis, ry, ry_std**-2)\n",
    "    boot_coefs[r] = rlr.coef_\n",
    "\n",
    "boot_coefs = np.array(boot_coefs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "293a77b9-5fd4-4e64-b6f2-7f2c84017c4a",
   "metadata": {},
   "outputs": [],
   "source": [
    "display(rreduced_freq.shape, relaxations.shape, relaxations)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6e887ba9-afbc-40c8-a105-6c00472d5d76",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "\n",
    "plt.loglog(\n",
    "    rreduced_freq.ravel(),\n",
    "    ry.reshape((2, -1)).T,\n",
    "    label=[\"Data Storage\", \"Data Loss\"],\n",
    "    **opts\n",
    ")\n",
    "plt.loglog(\n",
    "    rreduced_freq.ravel(),\n",
    "    rlr.predict(rbasis).reshape((2, -1)).T,\n",
    "    \"k\",\n",
    "    label=[\"Prony Storage\", \"Prony Loss\"],\n",
    "    **opts\n",
    ")\n",
    "plt.ylim(bottom=1e5, top=5e9)\n",
    "plt.xlim(left=1e-6, right=1e7)\n",
    "plt.ylabel(\"Young's modulus (Pa)\")\n",
    "plt.legend()\n",
    "None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9ed775da-935e-4475-8e1d-eb8755dafc69",
   "metadata": {},
   "outputs": [],
   "source": [
    "display(\n",
    "    reduced_freq.shape,\n",
    "    rreduced_freq.shape,\n",
    "    relaxations.shape,\n",
    "    np.shape(boot_coefs),\n",
    "    basis.shape,\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d9c6d87a-8185-4f23-8c08-dc217973fdcf",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "plt.loglog(\n",
    "    1 / relaxations,\n",
    "    np.transpose(boot_coefs)[:-1],\n",
    "    \"k\",\n",
    "    alpha=0.01,\n",
    ")\n",
    "plt.loglog(1 / relaxations, g, \".\")\n",
    "# plt.xlim(left=.01,right=1e10)\n",
    "None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "40c61c1c-54f0-47ba-8c61-7bddd6e44ae6",
   "metadata": {},
   "outputs": [],
   "source": [
    "dense_red_freq_space = 1 / dense_relaxation_space\n",
    "mle_basis = prony_basis(dense_red_freq_space, relaxations).reshape(\n",
    "    2, -1, len(relaxations) + 1\n",
    ")\n",
    "mle_05, mle_50, mle_95 = mle_quantiles = np.nanpercentile(\n",
    "    mle_basis @ boot_coefs.T, [5, 50, 95], axis=-1\n",
    ")\n",
    "plt.figure(layout=\"constrained\")\n",
    "plt.loglog(\n",
    "    dense_red_freq_space, mle_50.T, label=[\"Prony Storage\", \"Prony Loss\"], linewidth=1\n",
    ")\n",
    "plt.fill_between(\n",
    "    dense_red_freq_space, mle_05[0], mle_95[0], color=\"black\", linewidth=0, alpha=0.5\n",
    ")\n",
    "plt.fill_between(\n",
    "    dense_red_freq_space, mle_05[1], mle_95[1], color=\"black\", linewidth=0, alpha=0.5\n",
    ")\n",
    "\n",
    "# plt.loglog(dense_red_freq_space,(mle_basis@boot_coefs.T)[1],'k',alpha=1/510)\n",
    "\n",
    "plt.ylabel(\"Young's modulus (Pa)\")\n",
    "plt.xlabel(f\"Frequency at $T_{{ref}}$ = {parms.t_ref} °C (Hz)\")\n",
    "plt.legend()\n",
    "None"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c7d142ae-ca72-489f-a414-70e6a1d175ac",
   "metadata": {},
   "source": [
    "## Compare DFS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "77991f47-d349-42a0-9012-2a7462ef0b23",
   "metadata": {},
   "outputs": [],
   "source": [
    "len(temp_space)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ab0ab8f9-7718-4c2c-b6e0-7e9345ed5078",
   "metadata": {},
   "outputs": [],
   "source": [
    "dfs_xydata = np.array(\n",
    "    [\n",
    "        [2.88930816e04, 2.35809997e09],\n",
    "        [4.84048149e03, 2.05624000e09],\n",
    "        [1.68394703e02, 1.33180006e09],\n",
    "        [4.93183084e00, 4.12248992e08],\n",
    "        [2.09334604e-01, 7.52838000e07],\n",
    "        [1.40000004e-02, 3.64805000e07],\n",
    "        [1.43830527e-03, 3.04376000e07],\n",
    "        [2.27085283e-04, 2.90306000e07],\n",
    "        [4.99336582e-05, 2.85859000e07],\n",
    "        [1.38576602e-05, 2.85666000e07],\n",
    "        [4.54927672e-06, 2.88316000e07],\n",
    "        [1.75883932e-06, 2.91687000e07],\n",
    "        [8.21795917e-07, 2.97732000e07],\n",
    "        [3.63743255e04, 2.40706995e09],\n",
    "        [6.09382037e03, 2.08962995e09],\n",
    "        [2.11996900e02, 1.39641997e09],\n",
    "        [6.20882267e00, 4.60131008e08],\n",
    "        [2.63537311e-01, 8.30694000e07],\n",
    "        [1.76250003e-02, 3.77599000e07],\n",
    "        [1.81072358e-03, 3.06503000e07],\n",
    "        [2.85884147e-04, 2.90460000e07],\n",
    "        [6.28629082e-05, 2.84966000e07],\n",
    "        [1.74458041e-05, 2.87342000e07],\n",
    "        [5.72721437e-06, 2.88671000e07],\n",
    "        [2.21425305e-06, 2.91446000e07],\n",
    "        [1.03458235e-06, 2.97469000e07],\n",
    "        [4.57924364e04, 2.41417011e09],\n",
    "        [7.67164415e03, 2.16678989e09],\n",
    "        [2.66887548e02, 1.47785997e09],\n",
    "        [7.81642307e00, 5.06838016e08],\n",
    "        [3.31772902e-01, 9.32335040e07],\n",
    "        [2.21884996e-02, 3.89044000e07],\n",
    "        [2.27955963e-03, 3.10149000e07],\n",
    "        [3.59905824e-04, 2.92139000e07],\n",
    "        [7.91394943e-05, 2.86680000e07],\n",
    "        [2.19629055e-05, 2.86720000e07],\n",
    "        [7.21011582e-06, 2.87352000e07],\n",
    "        [2.78757174e-06, 2.90755000e07],\n",
    "        [1.30245841e-06, 2.95340000e07],\n",
    "        [5.76493321e04, 2.44448998e09],\n",
    "        [9.65803955e03, 2.18273997e09],\n",
    "        [3.35991926e02, 1.54586995e09],\n",
    "        [9.84030563e00, 5.71596992e08],\n",
    "        [4.17677847e-01, 1.04231000e08],\n",
    "        [2.79337000e-02, 4.05646000e07],\n",
    "        [2.86979905e-03, 3.13223000e07],\n",
    "        [4.53095140e-04, 2.92872000e07],\n",
    "        [9.96308419e-05, 2.86889000e07],\n",
    "        [2.76496936e-05, 2.86458000e07],\n",
    "        [9.07700909e-06, 2.87357000e07],\n",
    "        [3.50934918e-06, 2.91006000e07],\n",
    "        [1.63969999e-06, 2.95546000e07],\n",
    "        [7.25761196e04, 2.45549005e09],\n",
    "        [1.21587364e04, 2.22003994e09],\n",
    "        [4.22988252e02, 1.59059994e09],\n",
    "        [1.23881955e01, 6.19676032e08],\n",
    "        [5.25824607e-01, 1.16966000e08],\n",
    "        [3.51664014e-02, 4.22704000e07],\n",
    "        [3.61285849e-03, 3.16327000e07],\n",
    "        [5.70412282e-04, 2.93565000e07],\n",
    "        [1.25427644e-04, 2.87650000e07],\n",
    "        [3.48088589e-05, 2.86650000e07],\n",
    "        [1.14272633e-05, 2.87294000e07],\n",
    "        [4.41800341e-06, 2.90332000e07],\n",
    "        [2.06425744e-06, 2.93258000e07],\n",
    "        [9.13679726e04, 2.52575002e09],\n",
    "        [1.53069509e04, 2.26403994e09],\n",
    "        [5.32510959e02, 1.64974003e09],\n",
    "        [1.55958229e01, 6.89894976e08],\n",
    "        [6.61974332e-01, 1.33571000e08],\n",
    "        [4.42719012e-02, 4.44351000e07],\n",
    "        [4.54832192e-03, 3.21030000e07],\n",
    "        [7.18106921e-04, 2.95735000e07],\n",
    "        [1.57904137e-04, 2.88840000e07],\n",
    "        [4.38217817e-05, 2.86639000e07],\n",
    "        [1.43860802e-05, 2.87072000e07],\n",
    "        [5.56193988e-06, 2.90200000e07],\n",
    "        [2.59874761e-06, 2.93409000e07],\n",
    "        [1.15025417e05, 2.43572992e09],\n",
    "        [1.92703018e04, 2.29150003e09],\n",
    "        [6.70391313e02, 1.71342003e09],\n",
    "        [1.96339700e01, 7.47174016e08],\n",
    "        [8.33375979e-01, 1.52975008e08],\n",
    "        [5.57349995e-02, 4.71127000e07],\n",
    "        [5.72599579e-03, 3.27821000e07],\n",
    "        [9.04042694e-04, 2.97644000e07],\n",
    "        [1.98789452e-04, 2.88360000e07],\n",
    "        [5.51683327e-05, 2.87289000e07],\n",
    "        [1.81109948e-05, 2.87352000e07],\n",
    "        [7.00206469e-06, 2.89488000e07],\n",
    "        [3.27162811e-06, 2.93574000e07],\n",
    "        [1.44808407e05, 2.72194995e09],\n",
    "        [2.42598703e04, 2.35040998e09],\n",
    "        [8.43972578e02, 1.77074995e09],\n",
    "        [2.47177014e01, 8.14913984e08],\n",
    "        [1.04915810e00, 1.74464992e08],\n",
    "        [7.01662004e-02, 5.02472000e07],\n",
    "        [7.20860091e-03, 3.32486000e07],\n",
    "        [1.13812221e-03, 2.98733000e07],\n",
    "        [2.50261069e-04, 2.89960000e07],\n",
    "        [6.94528093e-05, 2.87352000e07],\n",
    "        [2.28003894e-05, 2.87705000e07],\n",
    "        [8.81507632e-06, 2.90111000e07],\n",
    "        [4.11873537e-06, 2.93060000e07],\n",
    "        [1.82302960e05, 2.09378995e09],\n",
    "        [3.05413633e04, 2.37644006e09],\n",
    "        [1.06249839e03, 1.83948006e09],\n",
    "        [3.11177385e01, 8.80388992e08],\n",
    "        [1.32081162e00, 1.99146000e08],\n",
    "        [8.83340016e-02, 5.35809000e07],\n",
    "        [9.07508972e-03, 3.41104000e07],\n",
    "        [1.43281079e-03, 3.00090000e07],\n",
    "        [3.15059979e-04, 2.91825000e07],\n",
    "        [8.74358955e-05, 2.86868000e07],\n",
    "        [2.87039860e-05, 2.87105000e07],\n",
    "        [1.10975222e-05, 2.89036000e07],\n",
    "        [5.18517997e-06, 2.93551000e07],\n",
    "        [2.29506001e05, 2.44329011e09],\n",
    "        [3.84493272e04, 2.39518003e09],\n",
    "        [1.33760723e03, 1.89703002e09],\n",
    "        [3.91749411e01, 9.51798016e08],\n",
    "        [1.66280456e00, 2.27535008e08],\n",
    "        [1.11206003e-01, 5.79391000e07],\n",
    "        [1.14248696e-02, 3.48349000e07],\n",
    "        [1.80380327e-03, 3.04442000e07],\n",
    "        [3.96637311e-04, 2.92606000e07],\n",
    "        [1.10075353e-04, 2.88735000e07],\n",
    "        [3.61362045e-05, 2.87367000e07],\n",
    "        [1.39709631e-05, 2.88924000e07],\n",
    "        [6.52775972e-06, 2.92419000e07],\n",
    "        [2.88930808e05, 2.25561011e09],\n",
    "        [4.84048136e04, 2.44013005e09],\n",
    "        [1.68394698e03, 1.94643994e09],\n",
    "        [4.93183071e01, 1.02092000e09],\n",
    "        [2.09334599e00, 2.61351008e08],\n",
    "        [1.40000001e-01, 6.31465000e07],\n",
    "        [1.43830523e-02, 3.58598000e07],\n",
    "        [2.27085276e-03, 3.05164000e07],\n",
    "        [4.99336569e-04, 2.94678000e07],\n",
    "        [1.38576598e-04, 2.88805000e07],\n",
    "        [4.54927660e-05, 2.88535000e07],\n",
    "        [1.75883928e-05, 2.89546000e07],\n",
    "        [8.21795895e-06, 2.92113000e07],\n",
    "        [3.63743240e05, 2.83956992e09],\n",
    "        [6.09382012e04, 2.46500992e09],\n",
    "        [2.11996891e03, 1.98321997e09],\n",
    "        [6.20882241e01, 1.09905997e09],\n",
    "        [2.63537300e00, 2.97115008e08],\n",
    "        [1.76249996e-01, 6.90476960e07],\n",
    "        [1.81072350e-02, 3.73628000e07],\n",
    "        [2.85884135e-03, 3.09941000e07],\n",
    "        [6.28629056e-04, 2.96065000e07],\n",
    "        [1.74458034e-04, 2.91616000e07],\n",
    "        [5.72721413e-05, 2.89999000e07],\n",
    "        [2.21425295e-05, 2.89601000e07],\n",
    "        [1.03458230e-05, 2.93052000e07],\n",
    "        [4.57924364e05, 2.48066995e09],\n",
    "        [7.67164415e04, 2.48574003e09],\n",
    "        [2.66887548e03, 2.03843994e09],\n",
    "        [7.81642307e01, 1.16831002e09],\n",
    "        [3.31772902e00, 3.37899008e08],\n",
    "        [2.21884996e-01, 7.67076000e07],\n",
    "        [2.27955963e-02, 3.86129000e07],\n",
    "        [3.59905824e-03, 3.14767000e07],\n",
    "        [7.91394943e-04, 2.98455000e07],\n",
    "        [2.19629055e-04, 2.91606000e07],\n",
    "        [7.21011582e-05, 2.89534000e07],\n",
    "        [2.78757174e-05, 2.89992000e07],\n",
    "        [1.30245841e-05, 2.92633000e07],\n",
    "        [5.76493298e05, 2.56333005e09],\n",
    "        [9.65803916e04, 2.51097011e09],\n",
    "        [3.35991913e03, 2.08479002e09],\n",
    "        [9.84030524e01, 1.23843994e09],\n",
    "        [4.17677830e00, 3.83539008e08],\n",
    "        [2.79336989e-01, 8.53048000e07],\n",
    "        [2.86979893e-02, 4.02024000e07],\n",
    "        [4.53095122e-03, 3.19398000e07],\n",
    "        [9.96308379e-04, 3.01823000e07],\n",
    "        [2.76496925e-04, 2.94161000e07],\n",
    "        [9.07700872e-05, 2.92175000e07],\n",
    "        [3.50934904e-05, 2.90737000e07],\n",
    "        [1.63969993e-05, 2.93416000e07],\n",
    "        [7.25761181e05, 2.68987008e09],\n",
    "        [1.21587362e05, 2.53516006e09],\n",
    "        [4.22988243e03, 2.14198003e09],\n",
    "        [1.23881953e02, 1.30705997e09],\n",
    "        [5.25824596e00, 4.31500992e08],\n",
    "        [3.51664007e-01, 9.54642000e07],\n",
    "        [3.61285842e-02, 4.26585000e07],\n",
    "        [5.70412270e-03, 3.26713000e07],\n",
    "        [1.25427641e-03, 3.00984000e07],\n",
    "        [3.48088582e-04, 2.96128000e07],\n",
    "        [1.14272631e-04, 2.91861000e07],\n",
    "        [4.41800331e-05, 2.92413000e07],\n",
    "        [2.06425740e-05, 2.93807000e07],\n",
    "        [9.13679726e05, 2.66642995e09],\n",
    "        [1.53069509e05, 2.55257011e09],\n",
    "        [5.32510959e03, 2.16378010e09],\n",
    "        [1.55958229e02, 1.37087002e09],\n",
    "        [6.61974332e00, 4.79663008e08],\n",
    "        [4.42719012e-01, 1.08774000e08],\n",
    "        [4.54832192e-02, 4.51292000e07],\n",
    "        [7.18106921e-03, 3.32548000e07],\n",
    "        [1.57904137e-03, 3.03186000e07],\n",
    "        [4.38217817e-04, 2.97139000e07],\n",
    "        [1.43860802e-04, 2.93701000e07],\n",
    "        [5.56193988e-05, 2.93283000e07],\n",
    "        [2.59874761e-05, 2.93249000e07],\n",
    "        [1.15025414e06, 2.76546995e09],\n",
    "        [1.92703013e05, 2.57286989e09],\n",
    "        [6.70391295e03, 2.22704000e09],\n",
    "        [1.96339695e02, 1.43663002e09],\n",
    "        [8.33375956e00, 5.38400000e08],\n",
    "        [5.57349980e-01, 1.23358000e08],\n",
    "        [5.72599563e-02, 4.80782000e07],\n",
    "        [9.04042670e-03, 3.40167000e07],\n",
    "        [1.98789447e-03, 3.09972000e07],\n",
    "        [5.51683312e-04, 2.99472000e07],\n",
    "        [1.81109944e-04, 2.96384000e07],\n",
    "        [7.00206450e-05, 2.96321000e07],\n",
    "        [3.27162803e-05, 2.99870000e07],\n",
    "        [1.44808407e06, 2.71963008e09],\n",
    "        [2.42598703e05, 2.60763008e09],\n",
    "        [8.43972578e03, 2.23856000e09],\n",
    "        [2.47177014e02, 1.51048000e09],\n",
    "        [1.04915810e01, 6.00414016e08],\n",
    "        [7.01662004e-01, 1.41479008e08],\n",
    "        [7.20860091e-02, 5.10889000e07],\n",
    "        [1.13812221e-02, 3.53071000e07],\n",
    "        [2.50261069e-03, 3.11212000e07],\n",
    "        [6.94528093e-04, 3.01021000e07],\n",
    "        [2.28003894e-04, 2.96693000e07],\n",
    "        [8.81507632e-05, 2.95531000e07],\n",
    "        [4.11873537e-05, 2.96955000e07],\n",
    "        [1.82302957e06, 2.74002995e09],\n",
    "        [3.05413628e05, 2.61007002e09],\n",
    "        [1.06249837e04, 2.27563008e09],\n",
    "        [3.11177380e02, 1.55470003e09],\n",
    "        [1.32081160e01, 6.62268992e08],\n",
    "        [8.83340001e-01, 1.61538000e08],\n",
    "        [9.07508957e-02, 5.50373000e07],\n",
    "        [1.43281077e-02, 3.63304000e07],\n",
    "        [3.15059974e-03, 3.15691000e07],\n",
    "        [8.74358940e-04, 3.02066000e07],\n",
    "        [2.87039856e-04, 2.98257000e07],\n",
    "        [1.10975220e-04, 2.98221000e07],\n",
    "        [5.18517989e-05, 2.97059000e07],\n",
    "        [2.29505985e06, 2.76734003e09],\n",
    "        [3.84493246e05, 2.63878989e09],\n",
    "        [1.33760714e04, 2.31211008e09],\n",
    "        [3.91749385e02, 1.62355994e09],\n",
    "        [1.66280445e01, 7.24009984e08],\n",
    "        [1.11205995e00, 1.84403008e08],\n",
    "        [1.14248688e-01, 5.96131000e07],\n",
    "        [1.80380315e-02, 3.73292000e07],\n",
    "        [3.96637284e-03, 3.18553000e07],\n",
    "        [1.10075346e-03, 3.02486000e07],\n",
    "        [3.61362021e-04, 2.97697000e07],\n",
    "        [1.39709622e-04, 2.96359000e07],\n",
    "        [6.52775928e-05, 2.97162000e07],\n",
    "        [2.88930802e06, 2.76901990e09],\n",
    "        [4.84048126e05, 2.65683994e09],\n",
    "        [1.68394695e04, 2.34917990e09],\n",
    "        [4.93183060e02, 1.67359002e09],\n",
    "        [2.09334594e01, 7.88316032e08],\n",
    "        [1.39999998e00, 2.10956000e08],\n",
    "        [1.43830520e-01, 6.51287000e07],\n",
    "        [2.27085272e-02, 3.88538000e07],\n",
    "        [4.99336558e-03, 3.23718000e07],\n",
    "        [1.38576595e-03, 3.04960000e07],\n",
    "        [4.54927650e-04, 2.98401000e07],\n",
    "        [1.75883924e-04, 2.97680000e07],\n",
    "        [8.21795878e-05, 2.97440000e07],\n",
    "        [3.63743258e06, 2.77038003e09],\n",
    "        [6.09382042e05, 2.67361997e09],\n",
    "        [2.11996902e04, 2.37694003e09],\n",
    "        [6.20882273e02, 1.73298995e09],\n",
    "        [2.63537313e01, 8.57160000e08],\n",
    "        [1.76250005e00, 2.43412000e08],\n",
    "        [1.81072359e-01, 7.20784960e07],\n",
    "        [2.85884149e-02, 4.06288000e07],\n",
    "        [6.28629088e-03, 3.28805000e07],\n",
    "        [1.74458042e-03, 3.07442000e07],\n",
    "        [5.72721442e-04, 3.00508000e07],\n",
    "        [2.21425307e-04, 2.99080000e07],\n",
    "        [1.03458236e-04, 2.98362000e07],\n",
    "        [4.57924351e06, 2.79011994e09],\n",
    "        [7.67164395e05, 2.69050010e09],\n",
    "        [2.66887541e04, 2.40341990e09],\n",
    "        [7.81642286e02, 1.78670003e09],\n",
    "        [3.31772894e01, 9.23057024e08],\n",
    "        [2.21884990e00, 2.79875008e08],\n",
    "        [2.27955957e-01, 8.00124960e07],\n",
    "        [3.59905814e-02, 4.27340000e07],\n",
    "        [7.91394921e-03, 3.35888000e07],\n",
    "        [2.19629049e-03, 3.09295000e07],\n",
    "        [7.21011562e-04, 3.01660000e07],\n",
    "        [2.78757166e-04, 2.99776000e07],\n",
    "        [1.30245838e-04, 2.99523000e07],\n",
    "        [5.76493323e06, 2.79023002e09],\n",
    "        [9.65803958e05, 2.70944000e09],\n",
    "        [3.35991927e04, 2.43164006e09],\n",
    "        [9.84030566e02, 1.83625997e09],\n",
    "        [4.17677848e01, 9.88563968e08],\n",
    "        [2.79337001e00, 3.15992992e08],\n",
    "        [2.86979906e-01, 8.95311040e07],\n",
    "        [4.53095141e-02, 4.51368000e07],\n",
    "        [9.96308421e-03, 3.44789000e07],\n",
    "        [2.76496936e-03, 3.12890000e07],\n",
    "        [9.07700911e-04, 3.03802000e07],\n",
    "        [3.50934919e-04, 3.00593000e07],\n",
    "        [1.63970000e-04, 3.00112000e07],\n",
    "        [7.25761156e06, 2.81025997e09],\n",
    "        [1.21587358e06, 2.72059008e09],\n",
    "        [4.22988229e04, 2.45210010e09],\n",
    "        [1.23881949e03, 1.88201997e09],\n",
    "        [5.25824578e01, 1.05824000e09],\n",
    "        [3.51663995e00, 3.58104992e08],\n",
    "        [3.61285830e-01, 1.00493000e08],\n",
    "        [5.70412251e-02, 4.80442000e07],\n",
    "        [1.25427637e-02, 3.54388000e07],\n",
    "        [3.48088570e-03, 3.16938000e07],\n",
    "        [1.14272627e-03, 3.04438000e07],\n",
    "        [4.41800316e-04, 3.01391000e07],\n",
    "        [2.06425733e-04, 3.00762000e07],\n",
    "        [9.13679665e06, 2.81035008e09],\n",
    "        [1.53069498e06, 2.72902989e09],\n",
    "        [5.32510923e04, 2.47672013e09],\n",
    "        [1.55958219e03, 1.93237005e09],\n",
    "        [6.61974287e01, 1.12208998e09],\n",
    "        [4.42718983e00, 4.03097984e08],\n",
    "        [4.54832162e-01, 1.13787000e08],\n",
    "        [7.18106873e-02, 5.14073000e07],\n",
    "        [1.57904126e-02, 3.65964000e07],\n",
    "        [4.38217787e-03, 3.20795000e07],\n",
    "        [1.43860793e-03, 3.05249000e07],\n",
    "        [5.56193951e-04, 3.02228000e07],\n",
    "        [2.59874744e-04, 3.01463000e07],\n",
    "        [1.15025421e07, 2.82073011e09],\n",
    "        [1.92703025e06, 2.71596006e09],\n",
    "        [6.70391338e04, 2.50260992e09],\n",
    "        [1.96339708e03, 1.97767002e09],\n",
    "        [8.33376010e01, 1.18937997e09],\n",
    "        [5.57350016e00, 4.43440992e08],\n",
    "        [5.72599600e-01, 1.28419000e08],\n",
    "        [9.04042728e-02, 5.58330000e07],\n",
    "        [1.98789459e-02, 3.81279000e07],\n",
    "        [5.51683348e-03, 3.26839000e07],\n",
    "        [1.81109955e-03, 3.08939000e07],\n",
    "        [7.00206495e-04, 3.04504000e07],\n",
    "        [3.27162824e-04, 3.02862000e07],\n",
    "        [1.44808409e07, 2.82982989e09],\n",
    "        [2.42598707e06, 2.76435994e09],\n",
    "        [8.43972593e04, 2.52434995e09],\n",
    "        [2.47177019e03, 2.01423002e09],\n",
    "        [1.04915811e02, 1.24910003e09],\n",
    "        [7.01662016e00, 4.80158016e08],\n",
    "        [7.20860103e-01, 1.41368000e08],\n",
    "        [1.13812223e-01, 6.01602000e07],\n",
    "        [2.50261073e-02, 3.94112000e07],\n",
    "        [6.94528105e-03, 3.30617000e07],\n",
    "        [2.28003898e-03, 3.09508000e07],\n",
    "        [8.81507647e-04, 3.03818000e07],\n",
    "        [4.11873544e-04, 3.02330000e07],\n",
    "        [1.82302952e07, 2.83787008e09],\n",
    "        [3.05413620e06, 2.76451994e09],\n",
    "        [1.06249835e05, 2.54237005e09],\n",
    "        [3.11177372e03, 2.05954995e09],\n",
    "        [1.32081156e02, 1.30521997e09],\n",
    "        [8.83339977e00, 5.14884000e08],\n",
    "        [9.07508932e-01, 1.54779008e08],\n",
    "        [1.43281073e-01, 6.43009000e07],\n",
    "        [3.15059965e-02, 4.12170000e07],\n",
    "        [8.74358917e-03, 3.38257000e07],\n",
    "        [2.87039848e-03, 3.12180000e07],\n",
    "        [1.10975217e-03, 3.05112000e07],\n",
    "        [5.18517975e-04, 3.03063000e07],\n",
    "        [2.29505990e07, 2.84626995e09],\n",
    "        [3.84493254e06, 2.77990989e09],\n",
    "        [1.33760717e05, 2.56048998e09],\n",
    "        [3.91749393e03, 2.09536998e09],\n",
    "        [1.66280448e02, 1.35128000e09],\n",
    "        [1.11205997e01, 5.40694976e08],\n",
    "        [1.14248691e00, 1.69226000e08],\n",
    "        [1.80380318e-01, 6.84076960e07],\n",
    "        [3.96637293e-02, 4.31653000e07],\n",
    "        [1.10075348e-02, 3.46014000e07],\n",
    "        [3.61362028e-03, 3.15614000e07],\n",
    "        [1.39709625e-03, 3.06835000e07],\n",
    "        [6.52775942e-04, 3.03826000e07],\n",
    "        [2.88930807e07, 2.84086989e09],\n",
    "        [4.84048134e06, 2.79079987e09],\n",
    "        [1.68394698e05, 2.58159002e09],\n",
    "        [4.93183068e03, 2.12445005e09],\n",
    "        [2.09334598e02, 1.39553997e09],\n",
    "        [1.40000000e01, 5.68672000e08],\n",
    "        [1.43830523e00, 1.85922000e08],\n",
    "        [2.27085276e-01, 7.24932000e07],\n",
    "        [4.99336567e-02, 4.46208000e07],\n",
    "        [1.38576597e-02, 3.54468000e07],\n",
    "        [4.54927658e-03, 3.20297000e07],\n",
    "        [1.75883927e-03, 3.08481000e07],\n",
    "        [8.21795892e-04, 3.05157000e07],\n",
    "    ]\n",
    ").T\n",
    "trios_xydata = np.array(\n",
    "    [\n",
    "        [2.29393507e03, 2.35809997e09],\n",
    "        [1.04584629e03, 2.05624000e09],\n",
    "        [1.14016602e02, 1.33180006e09],\n",
    "        [4.53194400e00, 4.12248992e08],\n",
    "        [2.07550657e-01, 7.52838000e07],\n",
    "        [1.40000004e-02, 3.64805000e07],\n",
    "        [1.52146161e-03, 3.04376000e07],\n",
    "        [2.54906653e-04, 2.90306000e07],\n",
    "        [5.52856703e-05, 2.85859000e07],\n",
    "        [1.39027376e-05, 2.85666000e07],\n",
    "        [3.84241006e-06, 2.88316000e07],\n",
    "        [6.02580040e-07, 2.91687000e07],\n",
    "        [4.39049997e-08, 2.97732000e07],\n",
    "        [2.88790037e03, 2.40706995e09],\n",
    "        [1.31664576e03, 2.08962995e09],\n",
    "        [1.43538756e02, 1.39641997e09],\n",
    "        [5.70539372e00, 4.60131008e08],\n",
    "        [2.61291449e-01, 8.30694000e07],\n",
    "        [1.76250003e-02, 3.77599000e07],\n",
    "        [1.91541147e-03, 3.06503000e07],\n",
    "        [3.20909264e-04, 2.90460000e07],\n",
    "        [6.96007091e-05, 2.84966000e07],\n",
    "        [1.75025534e-05, 2.87342000e07],\n",
    "        [4.83731975e-06, 2.88671000e07],\n",
    "        [7.58605220e-07, 2.91446000e07],\n",
    "        [5.52732578e-08, 2.97469000e07],\n",
    "        [3.63564113e03, 2.41417011e09],\n",
    "        [1.65755423e03, 2.16678989e09],\n",
    "        [1.80704089e02, 1.47785997e09],\n",
    "        [7.18264532e00, 5.06838016e08],\n",
    "        [3.28945538e-01, 9.32335040e07],\n",
    "        [2.21884996e-02, 3.89044000e07],\n",
    "        [2.41135352e-03, 3.10149000e07],\n",
    "        [4.03999713e-04, 2.92139000e07],\n",
    "        [8.76218596e-05, 2.86680000e07],\n",
    "        [2.20343484e-05, 2.86720000e07],\n",
    "        [6.08980795e-06, 2.87352000e07],\n",
    "        [9.55024753e-07, 2.90755000e07],\n",
    "        [6.95847169e-08, 2.95340000e07],\n",
    "        [4.57700658e03, 2.44448998e09],\n",
    "        [2.08673969e03, 2.18273997e09],\n",
    "        [2.27493247e02, 1.54586995e09],\n",
    "        [9.04242574e00, 5.71596992e08],\n",
    "        [4.14118401e-01, 1.04231000e08],\n",
    "        [2.79337000e-02, 4.05646000e07],\n",
    "        [3.03571792e-03, 3.13223000e07],\n",
    "        [5.08606125e-04, 2.92872000e07],\n",
    "        [1.10309520e-04, 2.86889000e07],\n",
    "        [2.77396349e-05, 2.86458000e07],\n",
    "        [7.66662333e-06, 2.87357000e07],\n",
    "        [1.20230640e-06, 2.91006000e07],\n",
    "        [8.76020751e-08, 2.95546000e07],\n",
    "        [5.76210279e03, 2.45549005e09],\n",
    "        [2.62704638e03, 2.22003994e09],\n",
    "        [2.86396676e02, 1.59059994e09],\n",
    "        [1.13837255e01, 6.19676032e08],\n",
    "        [5.21343536e-01, 1.16966000e08],\n",
    "        [3.51664014e-02, 4.22704000e07],\n",
    "        [3.82173772e-03, 3.16327000e07],\n",
    "        [6.40296386e-04, 2.93565000e07],\n",
    "        [1.38871287e-04, 2.87650000e07],\n",
    "        [3.49220883e-05, 2.86650000e07],\n",
    "        [9.65169503e-06, 2.87294000e07],\n",
    "        [1.51361221e-06, 2.90332000e07],\n",
    "        [1.10284343e-07, 2.93258000e07],\n",
    "        [7.25406170e03, 2.52575002e09],\n",
    "        [3.30725730e03, 2.26403994e09],\n",
    "        [3.60552256e02, 1.64974003e09],\n",
    "        [1.43312694e01, 6.89894976e08],\n",
    "        [6.56332994e-01, 1.33571000e08],\n",
    "        [4.42719012e-02, 4.44351000e07],\n",
    "        [4.81128544e-03, 3.21030000e07],\n",
    "        [8.06085844e-04, 2.95735000e07],\n",
    "        [1.74828690e-04, 2.88840000e07],\n",
    "        [4.39643290e-05, 2.86639000e07],\n",
    "        [1.21507710e-05, 2.87072000e07],\n",
    "        [1.90552594e-06, 2.90200000e07],\n",
    "        [1.38839839e-07, 2.93409000e07],\n",
    "        [9.13231900e03, 2.43572992e09],\n",
    "        [4.16358861e03, 2.29150003e09],\n",
    "        [4.53908217e02, 1.71342003e09],\n",
    "        [1.80419921e01, 7.47174016e08],\n",
    "        [8.26273959e-01, 1.52975008e08],\n",
    "        [5.57349995e-02, 4.71127000e07],\n",
    "        [6.05704711e-03, 3.27821000e07],\n",
    "        [1.01480155e-03, 2.97644000e07],\n",
    "        [2.20096194e-04, 2.88360000e07],\n",
    "        [5.53477892e-05, 2.87289000e07],\n",
    "        [1.52969083e-05, 2.87352000e07],\n",
    "        [2.39891408e-06, 2.89488000e07],\n",
    "        [1.74788933e-07, 2.93574000e07],\n",
    "        [1.14969073e04, 2.72194995e09],\n",
    "        [5.24164700e03, 2.35040998e09],\n",
    "        [5.71436534e02, 1.77074995e09],\n",
    "        [2.27135201e01, 8.14913984e08],\n",
    "        [1.04021718e00, 1.74464992e08],\n",
    "        [7.01662004e-02, 5.02472000e07],\n",
    "        [7.62536979e-03, 3.32486000e07],\n",
    "        [1.27755934e-03, 2.98733000e07],\n",
    "        [2.77084665e-04, 2.89960000e07],\n",
    "        [6.96787315e-05, 2.87352000e07],\n",
    "        [1.92576647e-05, 2.87705000e07],\n",
    "        [3.02005361e-06, 2.90111000e07],\n",
    "        [2.20046208e-07, 2.93060000e07],\n",
    "        [1.44737470e04, 2.09378995e09],\n",
    "        [6.59884178e03, 2.37644006e09],\n",
    "        [7.19395884e02, 1.83948006e09],\n",
    "        [2.85946240e01, 8.80388992e08],\n",
    "        [1.30955568e00, 1.99146000e08],\n",
    "        [8.83340016e-02, 5.35809000e07],\n",
    "        [9.59977059e-03, 3.41104000e07],\n",
    "        [1.60835171e-03, 3.00090000e07],\n",
    "        [3.48828882e-04, 2.91825000e07],\n",
    "        [8.77203147e-05, 2.86868000e07],\n",
    "        [2.42439604e-05, 2.87105000e07],\n",
    "        [3.80202176e-06, 2.89036000e07],\n",
    "        [2.77021728e-07, 2.93551000e07],\n",
    "        [1.82213815e04, 2.44329011e09],\n",
    "        [8.30745582e03, 2.39518003e09],\n",
    "        [9.05666437e02, 1.89703002e09],\n",
    "        [3.59985258e01, 9.51798016e08],\n",
    "        [1.64863416e00, 2.27535008e08],\n",
    "        [1.11206003e-01, 5.79391000e07],\n",
    "        [1.20854042e-02, 3.48349000e07],\n",
    "        [2.02479636e-03, 3.04442000e07],\n",
    "        [4.39149872e-04, 2.92606000e07],\n",
    "        [1.10433416e-04, 2.88735000e07],\n",
    "        [3.05213607e-05, 2.87367000e07],\n",
    "        [4.78646539e-06, 2.88924000e07],\n",
    "        [3.48749954e-07, 2.92419000e07],\n",
    "        [2.29393501e04, 2.25561011e09],\n",
    "        [1.04584626e04, 2.44013005e09],\n",
    "        [1.14016599e03, 1.94643994e09],\n",
    "        [4.53194388e01, 1.02092000e09],\n",
    "        [2.07550652e00, 2.61351008e08],\n",
    "        [1.40000001e-01, 6.31465000e07],\n",
    "        [1.52146157e-02, 3.58598000e07],\n",
    "        [2.54906646e-03, 3.05164000e07],\n",
    "        [5.52856689e-04, 2.94678000e07],\n",
    "        [1.39027372e-04, 2.88805000e07],\n",
    "        [3.84240996e-05, 2.88535000e07],\n",
    "        [6.02580024e-06, 2.89546000e07],\n",
    "        [4.39049985e-07, 2.92113000e07],\n",
    "        [2.88790024e04, 2.83956992e09],\n",
    "        [1.31664570e04, 2.46500992e09],\n",
    "        [1.43538750e03, 1.98321997e09],\n",
    "        [5.70539348e01, 1.09905997e09],\n",
    "        [2.61291438e00, 2.97115008e08],\n",
    "        [1.76249996e-01, 6.90476960e07],\n",
    "        [1.91541139e-02, 3.73628000e07],\n",
    "        [3.20909251e-03, 3.09941000e07],\n",
    "        [6.96007062e-04, 2.96065000e07],\n",
    "        [1.75025526e-04, 2.91616000e07],\n",
    "        [4.83731955e-05, 2.89999000e07],\n",
    "        [7.58605188e-06, 2.89601000e07],\n",
    "        [5.52732555e-07, 2.93052000e07],\n",
    "        [3.63564113e04, 2.48066995e09],\n",
    "        [1.65755423e04, 2.48574003e09],\n",
    "        [1.80704089e03, 2.03843994e09],\n",
    "        [7.18264532e01, 1.16831002e09],\n",
    "        [3.28945538e00, 3.37899008e08],\n",
    "        [2.21884996e-01, 7.67076000e07],\n",
    "        [2.41135352e-02, 3.86129000e07],\n",
    "        [4.03999713e-03, 3.14767000e07],\n",
    "        [8.76218596e-04, 2.98455000e07],\n",
    "        [2.20343484e-04, 2.91606000e07],\n",
    "        [6.08980795e-05, 2.89534000e07],\n",
    "        [9.55024753e-06, 2.89992000e07],\n",
    "        [6.95847169e-07, 2.92633000e07],\n",
    "        [4.57700640e04, 2.56333005e09],\n",
    "        [2.08673960e04, 2.51097011e09],\n",
    "        [2.27493237e03, 2.08479002e09],\n",
    "        [9.04242538e01, 1.23843994e09],\n",
    "        [4.14118385e00, 3.83539008e08],\n",
    "        [2.79336989e-01, 8.53048000e07],\n",
    "        [3.03571780e-02, 4.02024000e07],\n",
    "        [5.08606105e-03, 3.19398000e07],\n",
    "        [1.10309516e-03, 3.01823000e07],\n",
    "        [2.77396338e-04, 2.94161000e07],\n",
    "        [7.66662303e-05, 2.92175000e07],\n",
    "        [1.20230635e-05, 2.90737000e07],\n",
    "        [8.76020716e-07, 2.93416000e07],\n",
    "        [5.76210267e04, 2.68987008e09],\n",
    "        [2.62704632e04, 2.53516006e09],\n",
    "        [2.86396670e03, 2.14198003e09],\n",
    "        [1.13837253e02, 1.30705997e09],\n",
    "        [5.21343525e00, 4.31500992e08],\n",
    "        [3.51664007e-01, 9.54642000e07],\n",
    "        [3.82173764e-02, 4.26585000e07],\n",
    "        [6.40296372e-03, 3.26713000e07],\n",
    "        [1.38871284e-03, 3.00984000e07],\n",
    "        [3.49220875e-04, 2.96128000e07],\n",
    "        [9.65169483e-05, 2.91861000e07],\n",
    "        [1.51361218e-05, 2.92413000e07],\n",
    "        [1.10284340e-06, 2.93807000e07],\n",
    "        [7.25406170e04, 2.66642995e09],\n",
    "        [3.30725730e04, 2.55257011e09],\n",
    "        [3.60552256e03, 2.16378010e09],\n",
    "        [1.43312694e02, 1.37087002e09],\n",
    "        [6.56332994e00, 4.79663008e08],\n",
    "        [4.42719012e-01, 1.08774000e08],\n",
    "        [4.81128544e-02, 4.51292000e07],\n",
    "        [8.06085844e-03, 3.32548000e07],\n",
    "        [1.74828690e-03, 3.03186000e07],\n",
    "        [4.39643290e-04, 2.97139000e07],\n",
    "        [1.21507710e-04, 2.93701000e07],\n",
    "        [1.90552594e-05, 2.93283000e07],\n",
    "        [1.38839839e-06, 2.93249000e07],\n",
    "        [9.13231876e04, 2.76546995e09],\n",
    "        [4.16358850e04, 2.57286989e09],\n",
    "        [4.53908205e03, 2.22704000e09],\n",
    "        [1.80419916e02, 1.43663002e09],\n",
    "        [8.26273936e00, 5.38400000e08],\n",
    "        [5.57349980e-01, 1.23358000e08],\n",
    "        [6.05704694e-02, 4.80782000e07],\n",
    "        [1.01480152e-02, 3.40167000e07],\n",
    "        [2.20096188e-03, 3.09972000e07],\n",
    "        [5.53477877e-04, 2.99472000e07],\n",
    "        [1.52969079e-04, 2.96384000e07],\n",
    "        [2.39891402e-05, 2.96321000e07],\n",
    "        [1.74788928e-06, 2.99870000e07],\n",
    "        [1.14969073e05, 2.71963008e09],\n",
    "        [5.24164700e04, 2.60763008e09],\n",
    "        [5.71436534e03, 2.23856000e09],\n",
    "        [2.27135201e02, 1.51048000e09],\n",
    "        [1.04021718e01, 6.00414016e08],\n",
    "        [7.01662004e-01, 1.41479008e08],\n",
    "        [7.62536979e-02, 5.10889000e07],\n",
    "        [1.27755934e-02, 3.53071000e07],\n",
    "        [2.77084665e-03, 3.11212000e07],\n",
    "        [6.96787315e-04, 3.01021000e07],\n",
    "        [1.92576647e-04, 2.96693000e07],\n",
    "        [3.02005361e-05, 2.95531000e07],\n",
    "        [2.20046208e-06, 2.96955000e07],\n",
    "        [1.44737468e05, 2.74002995e09],\n",
    "        [6.59884167e04, 2.61007002e09],\n",
    "        [7.19395872e03, 2.27563008e09],\n",
    "        [2.85946236e02, 1.55470003e09],\n",
    "        [1.30955566e01, 6.62268992e08],\n",
    "        [8.83340001e-01, 1.61538000e08],\n",
    "        [9.59977043e-02, 5.50373000e07],\n",
    "        [1.60835169e-02, 3.63304000e07],\n",
    "        [3.48828876e-03, 3.15691000e07],\n",
    "        [8.77203132e-04, 3.02066000e07],\n",
    "        [2.42439600e-04, 2.98257000e07],\n",
    "        [3.80202169e-05, 2.98221000e07],\n",
    "        [2.77021723e-06, 2.97059000e07],\n",
    "        [1.82213803e05, 2.76734003e09],\n",
    "        [8.30745527e04, 2.63878989e09],\n",
    "        [9.05666376e03, 2.31211008e09],\n",
    "        [3.59985234e02, 1.62355994e09],\n",
    "        [1.64863405e01, 7.24009984e08],\n",
    "        [1.11205995e00, 1.84403008e08],\n",
    "        [1.20854034e-01, 5.96131000e07],\n",
    "        [2.02479622e-02, 3.73292000e07],\n",
    "        [4.39149843e-03, 3.18553000e07],\n",
    "        [1.10433409e-03, 3.02486000e07],\n",
    "        [3.05213587e-04, 2.97697000e07],\n",
    "        [4.78646506e-05, 2.96359000e07],\n",
    "        [3.48749931e-06, 2.97162000e07],\n",
    "        [2.29393496e05, 2.76901990e09],\n",
    "        [1.04584624e05, 2.65683994e09],\n",
    "        [1.14016596e04, 2.34917990e09],\n",
    "        [4.53194379e02, 1.67359002e09],\n",
    "        [2.07550647e01, 7.88316032e08],\n",
    "        [1.39999998e00, 2.10956000e08],\n",
    "        [1.52146154e-01, 6.51287000e07],\n",
    "        [2.54906641e-02, 3.88538000e07],\n",
    "        [5.52856677e-03, 3.23718000e07],\n",
    "        [1.39027369e-03, 3.04960000e07],\n",
    "        [3.84240988e-04, 2.98401000e07],\n",
    "        [6.02580011e-05, 2.97680000e07],\n",
    "        [4.39049976e-06, 2.97440000e07],\n",
    "        [2.88790039e05, 2.77038003e09],\n",
    "        [1.31664577e05, 2.67361997e09],\n",
    "        [1.43538757e04, 2.37694003e09],\n",
    "        [5.70539377e02, 1.73298995e09],\n",
    "        [2.61291452e01, 8.57160000e08],\n",
    "        [1.76250005e00, 2.43412000e08],\n",
    "        [1.91541149e-01, 7.20784960e07],\n",
    "        [3.20909267e-02, 4.06288000e07],\n",
    "        [6.96007097e-03, 3.28805000e07],\n",
    "        [1.75025535e-03, 3.07442000e07],\n",
    "        [4.83731979e-04, 3.00508000e07],\n",
    "        [7.58605226e-05, 2.99080000e07],\n",
    "        [5.52732583e-06, 2.98362000e07],\n",
    "        [3.63564103e05, 2.79011994e09],\n",
    "        [1.65755419e05, 2.69050010e09],\n",
    "        [1.80704084e04, 2.40341990e09],\n",
    "        [7.18264513e02, 1.78670003e09],\n",
    "        [3.28945529e01, 9.23057024e08],\n",
    "        [2.21884990e00, 2.79875008e08],\n",
    "        [2.41135345e-01, 8.00124960e07],\n",
    "        [4.03999702e-02, 4.27340000e07],\n",
    "        [8.76218573e-03, 3.35888000e07],\n",
    "        [2.20343478e-03, 3.09295000e07],\n",
    "        [6.08980779e-04, 3.01660000e07],\n",
    "        [9.55024728e-05, 2.99776000e07],\n",
    "        [6.95847150e-06, 2.99523000e07],\n",
    "        [4.57700659e05, 2.79023002e09],\n",
    "        [2.08673969e05, 2.70944000e09],\n",
    "        [2.27493247e04, 2.43164006e09],\n",
    "        [9.04242576e02, 1.83625997e09],\n",
    "        [4.14118402e01, 9.88563968e08],\n",
    "        [2.79337001e00, 3.15992992e08],\n",
    "        [3.03571793e-01, 8.95311040e07],\n",
    "        [5.08606127e-02, 4.51368000e07],\n",
    "        [1.10309520e-02, 3.44789000e07],\n",
    "        [2.77396350e-03, 3.12890000e07],\n",
    "        [7.66662335e-04, 3.03802000e07],\n",
    "        [1.20230640e-04, 3.00593000e07],\n",
    "        [8.76020754e-06, 3.00112000e07],\n",
    "        [5.76210247e05, 2.81025997e09],\n",
    "        [2.62704623e05, 2.72059008e09],\n",
    "        [2.86396660e04, 2.45210010e09],\n",
    "        [1.13837249e03, 1.88201997e09],\n",
    "        [5.21343507e01, 1.05824000e09],\n",
    "        [3.51663995e00, 3.58104992e08],\n",
    "        [3.82173751e-01, 1.00493000e08],\n",
    "        [6.40296351e-02, 4.80442000e07],\n",
    "        [1.38871279e-02, 3.54388000e07],\n",
    "        [3.49220863e-03, 3.16938000e07],\n",
    "        [9.65169450e-04, 3.04438000e07],\n",
    "        [1.51361213e-04, 3.01391000e07],\n",
    "        [1.10284336e-05, 3.00762000e07],\n",
    "        [7.25406121e05, 2.81035008e09],\n",
    "        [3.30725708e05, 2.72902989e09],\n",
    "        [3.60552231e04, 2.47672013e09],\n",
    "        [1.43312684e03, 1.93237005e09],\n",
    "        [6.56332950e01, 1.12208998e09],\n",
    "        [4.42718983e00, 4.03097984e08],\n",
    "        [4.81128511e-01, 1.13787000e08],\n",
    "        [8.06085790e-02, 5.14073000e07],\n",
    "        [1.74828678e-02, 3.65964000e07],\n",
    "        [4.39643261e-03, 3.20795000e07],\n",
    "        [1.21507701e-03, 3.05249000e07],\n",
    "        [1.90552581e-04, 3.02228000e07],\n",
    "        [1.38839830e-05, 3.01463000e07],\n",
    "        [9.13231934e05, 2.82073011e09],\n",
    "        [4.16358877e05, 2.71596006e09],\n",
    "        [4.53908234e04, 2.50260992e09],\n",
    "        [1.80419927e03, 1.97767002e09],\n",
    "        [8.26273989e01, 1.18937997e09],\n",
    "        [5.57350016e00, 4.43440992e08],\n",
    "        [6.05704733e-01, 1.28419000e08],\n",
    "        [1.01480159e-01, 5.58330000e07],\n",
    "        [2.20096202e-02, 3.81279000e07],\n",
    "        [5.53477912e-03, 3.26839000e07],\n",
    "        [1.52969089e-03, 3.08939000e07],\n",
    "        [2.39891417e-04, 3.04504000e07],\n",
    "        [1.74788939e-05, 3.02862000e07],\n",
    "        [1.14969075e06, 2.82982989e09],\n",
    "        [5.24164709e05, 2.76435994e09],\n",
    "        [5.71436544e04, 2.52434995e09],\n",
    "        [2.27135205e03, 2.01423002e09],\n",
    "        [1.04021720e02, 1.24910003e09],\n",
    "        [7.01662016e00, 4.80158016e08],\n",
    "        [7.62536992e-01, 1.41368000e08],\n",
    "        [1.27755936e-01, 6.01602000e07],\n",
    "        [2.77084669e-02, 3.94112000e07],\n",
    "        [6.96787327e-03, 3.30617000e07],\n",
    "        [1.92576650e-03, 3.09508000e07],\n",
    "        [3.02005366e-04, 3.03818000e07],\n",
    "        [2.20046212e-05, 3.02330000e07],\n",
    "        [1.44737464e06, 2.83787008e09],\n",
    "        [6.59884149e05, 2.76451994e09],\n",
    "        [7.19395852e04, 2.54237005e09],\n",
    "        [2.85946228e03, 2.05954995e09],\n",
    "        [1.30955562e02, 1.30521997e09],\n",
    "        [8.83339977e00, 5.14884000e08],\n",
    "        [9.59977017e-01, 1.54779008e08],\n",
    "        [1.60835164e-01, 6.43009000e07],\n",
    "        [3.48828866e-02, 4.12170000e07],\n",
    "        [8.77203109e-03, 3.38257000e07],\n",
    "        [2.42439594e-03, 3.12180000e07],\n",
    "        [3.80202159e-04, 3.05112000e07],\n",
    "        [2.77021716e-05, 3.03063000e07],\n",
    "        [1.82213807e06, 2.84626995e09],\n",
    "        [8.30745544e05, 2.77990989e09],\n",
    "        [9.05666396e04, 2.56048998e09],\n",
    "        [3.59985241e03, 2.09536998e09],\n",
    "        [1.64863408e02, 1.35128000e09],\n",
    "        [1.11205997e01, 5.40694976e08],\n",
    "        [1.20854036e00, 1.69226000e08],\n",
    "        [2.02479627e-01, 6.84076960e07],\n",
    "        [4.39149852e-02, 4.31653000e07],\n",
    "        [1.10433411e-02, 3.46014000e07],\n",
    "        [3.05213593e-03, 3.15614000e07],\n",
    "        [4.78646517e-04, 3.06835000e07],\n",
    "        [3.48749938e-05, 3.03826000e07],\n",
    "        [2.29393500e06, 2.84086989e09],\n",
    "        [1.04584626e06, 2.79079987e09],\n",
    "        [1.14016598e05, 2.58159002e09],\n",
    "        [4.53194386e03, 2.12445005e09],\n",
    "        [2.07550651e02, 1.39553997e09],\n",
    "        [1.40000000e01, 5.68672000e08],\n",
    "        [1.52146156e00, 1.85922000e08],\n",
    "        [2.54906645e-01, 7.24932000e07],\n",
    "        [5.52856686e-02, 4.46208000e07],\n",
    "        [1.39027372e-02, 3.54468000e07],\n",
    "        [3.84240994e-03, 3.20297000e07],\n",
    "        [6.02580021e-04, 3.08481000e07],\n",
    "        [4.39049983e-05, 3.05157000e07],\n",
    "    ]\n",
    ").T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "502121af-e002-4187-b32c-f0d786c6cb6c",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(layout=\"constrained\", figsize=(4 * 0.8, 3 * 0.8))\n",
    "plt.loglog(\n",
    "    reduced_freq[::600].reshape(-1),\n",
    "    args[0].reshape((2,) + reduced_freq.shape)[0, ::600].reshape(-1),\n",
    "    \"o\",\n",
    "    markerfacecolor=\"none\",\n",
    ")\n",
    "plt.plot(dfs_xydata[0, ::6], dfs_xydata[1, ::6], \"+\", markerfacecolor=\"none\")\n",
    "plt.plot(trios_xydata[0, ::6], trios_xydata[1, ::6], \"x\", markerfacecolor=\"none\")\n",
    "plt.xlabel(f\"Frequency at $T_{{ref}}$ = {parms.t_ref} °C (Hz)\")\n",
    "plt.ylabel(\"Young's modulus (Pa)\")\n",
    "plt.legend([\"BOTTS-Pair\", \"DFS-TTS-Pair\", \"DFS-TTS-TRIOS\"], loc=(-0.03, 0.65))\n",
    "plt.xlim(left=1.1e-9)\n",
    "None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e282950b-b413-44a5-afae-32e90d834e90",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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",
   "version": "3.12.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}