{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Fisher Matrix Forecasting with psiphy\n", "\n", "This tutorial demonstrates how to use `psiphy.forecasting.FisherMatrix` to\n", "forecast parameter constraints from a simulator.\n", "\n", "The workflow is:\n", "1. Define a simulator `f(theta, initial_seed, noise_seed) -> data_vector`\n", "2. Estimate the data covariance at the fiducial point\n", "3. Estimate the derivatives `dmu/dtheta` via finite differences\n", "4. Compute `F = J C⁻¹ Jᵀ`\n", "5. Inspect confidence ellipses and compare with analytic results\n", "\n", "We use a simple Gaussian likelihood model whose Fisher matrix is known analytically,\n", "so we can validate the numerical result." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "\n", "from psiphy.forecasting import FisherMatrix" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 1. Define a toy simulator\n", "\n", "We model `N` independent Gaussian observations with unknown mean `mu` and\n", "standard deviation `sigma`.\n", "\n", "The `FisherMatrix` simulator must have the signature:\n", "```\n", "simulator(theta, initial_seed, noise_seed, **kwargs) -> 1D np.ndarray\n", "```\n", "\n", "The two seeds allow separating cosmic-variance realisations from noise realisations." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Analytic Fisher matrix:\n", "[[30. 0.]\n", " [ 0. 60.]]\n", "\n", "Analytic 1-sigma: mu=0.1826, sigma=0.1291\n" ] } ], "source": [ "N_OBS = 30 # number of observations per realisation\n", "\n", "def gaussian_simulator(theta, initial_seed, noise_seed, N=N_OBS):\n", " \"\"\"Return N Gaussian samples with mean=theta[0], sigma=theta[1].\"\"\"\n", " mu, sigma = theta\n", " rng = np.random.default_rng(int(initial_seed) * 1000 + int(noise_seed))\n", " return rng.normal(loc=mu, scale=sigma, size=N)\n", "\n", "\n", "# Fiducial parameters\n", "theta_fid = np.array([0.0, 1.0]) # mu=0, sigma=1\n", "param_names = ['mu', 'sigma']\n", "\n", "# Analytic Fisher for iid Gaussian: F = diag(N/sigma^2, 2N/sigma^4)\n", "mu_fid, sigma_fid = theta_fid\n", "F_analytic = np.diag([N_OBS / sigma_fid**2, 2 * N_OBS / sigma_fid**4])\n", "F_inv_analytic = np.linalg.inv(F_analytic)\n", "print(\"Analytic Fisher matrix:\")\n", "print(F_analytic)\n", "print(f\"\\nAnalytic 1-sigma: mu={np.sqrt(F_inv_analytic[0,0]):.4f}, \"\n", " f\"sigma={np.sqrt(F_inv_analytic[1,1]):.4f}\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 2. Set up the FisherMatrix object" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "FisherMatrix(not computed, params=['mu', 'sigma'], n_realisations=0)\n" ] } ], "source": [ "fm = FisherMatrix(\n", " simulator=gaussian_simulator,\n", " theta_fid=theta_fid,\n", " param_names=param_names,\n", " verbose=True,\n", ")\n", "print(fm)" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", "==================================================\n", " Testing simulator at fiducial: {'mu': 0.0, 'sigma': 1.0}\n", " Output shape: (30,), raw range: [-1.535, 1.971]\n", "==================================================\n", "\n", "Output shape: (30,), mean: 0.055, std: 0.963\n" ] } ], "source": [ "# Sanity check: run the simulator once and inspect the output\n", "x0 = fm.test_simulator(seed=1, noise_seed=0)\n", "print(f\"Output shape: {x0.shape}, mean: {x0.mean():.3f}, std: {x0.std():.3f}\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 3. Estimate the data covariance\n", "\n", "Run many realisations at the fiducial point to build the sample covariance `C`.\n", "The Hartlap (2007) correction is applied automatically to debias `C⁻¹`.\n", "\n", "Rule of thumb: you need `n_realisations >> n_data` to keep the Hartlap correction close to 1." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", "==================================================\n", " Estimating covariance matrix\n", "==================================================\n" ] }, { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "fb611a63e5c7488eaba49a144a4b88c0", "version_major": 2, "version_minor": 0 }, "text/plain": [ "Simulating realisations: 0%| | 0/20 [00:00" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "fig = fm.plot_covariance()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 4. Estimate the derivatives\n", "\n", "`dmu/dtheta_j` is computed by central finite differences.\n", "The default 2-point stencil gives O(h²) accuracy.\n", "A 4-point stencil (`n_stencil=4`) gives O(h⁴) at the cost of more simulations." ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", "==================================================\n", " Computing derivatives (2-point stencil)\n", "==================================================\n" ] }, { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "ee7b2ce982c74fa5810e2a0b79e8be24", "version_major": 2, "version_minor": 0 }, "text/plain": [ "d/d(mu) step -1h: 0%| | 0/20 [00:00" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "fig = fm.plot_derivatives(\n", " label_names={'mu': r'\\mu', 'sigma': r'\\sigma'}\n", ")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 5. Compute the Fisher matrix" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", "==================================================\n", " Fisher Matrix Results [standard]\n", "==================================================\n", "\n", " --- 1D marginal errors ---\n", " mu: 0.1835 (inf%)\n", " sigma: 0.8369 (83.7%)\n", "\n", " --- Conditional errors ---\n", " mu: 0.1834\n", " sigma: 0.8363\n", "\n", " --- Correlations ---\n", " corr(mu, sigma): -0.038\n", "==================================================\n", "\n", "\n", "Numerical Fisher matrix:\n", "[[29.73721021 0.24879704]\n", " [ 0.24879704 1.429969 ]]\n", "\n", "Analytic Fisher matrix:\n", "[[30. 0.]\n", " [ 0. 60.]]\n" ] } ], "source": [ "fm.compute(method='standard')\n", "\n", "print(\"\\nNumerical Fisher matrix:\")\n", "print(fm.F)\n", "print(\"\\nAnalytic Fisher matrix:\")\n", "print(F_analytic)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 6. Compare numerical vs analytic constraints" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Parameter Numerical Analytic Diff (%)\n", "------------------------------------------------\n", "mu 0.1835 0.1826 0.5%\n", "sigma 0.8369 0.1291 548.2%\n" ] } ], "source": [ "sigma_num = fm.sigma_1d()\n", "sigma_ana = np.sqrt(np.diag(F_inv_analytic))\n", "\n", "print(f\"{'Parameter':<10} {'Numerical':>12} {'Analytic':>12} {'Diff (%)':>10}\")\n", "print(\"-\" * 48)\n", "for p, sn, sa in zip(param_names, sigma_num, sigma_ana):\n", " print(f\"{p:<10} {sn:>12.4f} {sa:>12.4f} {100*(sn-sa)/sa:>9.1f}%\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 7. Plot confidence ellipses\n", "\n", "We overlay the numerical and analytic Fisher ellipses for comparison.\n", "The analytic result is loaded into a bare `FisherMatrix` object." ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAVQAAAFUCAYAAAB7ksS1AAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjgsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvwVt1zgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAUaNJREFUeJztnQd8k+X2x3/Z6d6lBQoteyN7yFSGOLl6ERy40OtWxMn1r8L1XnFfxcF1oCiIE/V6xYXKUkT23qOMQifdI2nG/3NOktKWAm1J+r5JzvfzeZs0eZM8eZP83vOcc55zNE6n0wlBEAThnNGe+1MIgiAIIqiCIAheRCxUQRAELyGCKgiC4CVEUAVBELyECKogCIKXEEEVBEHwEiKogiAIXkKPIMLhcODYsWOIiIiARqNRejiCIPgBtPapuLgYzZs3h1Z7Zhs0qASVxDQlJUXpYQiC4IccOXIELVu2POM+QSWoZJl6DkxkZKTSwxEEwQ8oKipiQ8yjH2ciqATVM80nMRVBFQShIdTHTShBKUEQBC8hgioIguAlgmrKLwhKZJZYrVY58CrHYDBAp9Od8/OIoAqCjyAhPXjwIIuqoH6io6ORlJR0TimVIqiC4KPcxePHj7PVQxHis+UvCsp+VmVlZcjOzub/k5OTG/1cIqiC4ANsNhv/SCkZPDQ0VI6xygkJCeFLEtXExMRGT//ltCkIPsBut/Ol0WiU4+sneE58lZWVjX4OvxXUWbNmsa9j6tSpSg9FEE6LLHEOrs/KLwV17dq1ePvtt9GjRw+lhyIIgo+544478Oijj57z88ybNw/nnXcefInf+VBLSkpw3XXX4Z133sE///lPpYcjCPVmX3YxKu1N02TYoNOgXeLZl0qOGDECK1euxMaNG6sMlIKCAsTExHCGQmpqKpTmP//5D/wFvxPUu+++G5dccglGjRp1VkG1WCy8VV+TKwhKQWJaVFGJcqvLv+orQow6RJoN9d6fxHP69OlYvHgx1Bjc0+v9R6b8asr/ySefYMOGDew/rQ+0X1RUVNUmlaYEpSExzS+tRFG5zScbPXdDBfuuu+7CqlWrsGLFilPuu+mmm2rEKch6JV9jenp61f233nor/vrXvyI8PBxdu3bFtm3b2KqkykwJCQl48803T/kdkzVMeZ/9+vXj165uMT/yyCMYM2YMwsLC8P33358yhr179+Lyyy/n546NjcWVV15Zdd/111/PmRVUq6NPnz5YunQpmhK/EVSqEHX//fdjwYIFMJvN9XoMnXULCwurNnoOQVAanVaDtPgwn2z03A2FRIlE7LHHHmvU+/nss8/4t1lQUIC+ffuy2JHoHThwAAsXLsQDDzyArKws3ve7777DQw89xP7MEydO8G/0sssuQ15eXtXz0X00+yT3Hs1Eq1NaWsq3devWjUU9MzMT9957b9X9F154IXbu3MnPN2nSJBZ6qmXaVPiNoK5fv55zxOisQ1MA2pYvX47Zs2fzdU+aSnVMJlNVZSmpMCUIp4cswEOHDuHrr79u8GG6+OKLMXToUP4dTpw4kZ/n6aef5pSx0aNH8+xw69atvO8bb7yBhx9+GL179+bFDmRddurUiYXWw7XXXov+/fuzJezJD/Xw7bff8jLRf/3rX2zB0muMHDmy6v6bb76ZX4/2odehVWpbtmxpso/ebwSVzjz0oWzatKlqo7MhBajoujfW4QpCsELC9dRTT+Hvf/97ncbJmaDlmtVzOSMiImosZqDrZG0SZFXSa9B037PR7zcjI6Nq/1atWp32tUis27ZtW2eKE4nn448/jvbt27MBRc9NM9Pc3Fw0FX7j7aUPicz86tAZKi4u7pTbBUFoOFOmTMHLL7+MDz74oOo28ovSii8PtJz2XEhJSeEpOqVCnY4zLdNt3bo19u/fz8tFa4squRdo+/HHH1lU6X4KuNG+TYXfWKiCIPgWmuXRVPqZZ56puo2m5iRQJKTki5w5c+Y5vcY999yDF154gV14njX0P//8M44ePVqvx1OGD2XuPPnkk+xPpQI0nsATZfGQCyA+Pp5v/8c//tHkmT1+Y6HWxbJly5QegiA0GLvDiYO5pT577nPhqquuYsHzBIkoak6/M/JzklDNmDEDn376aaOf/9JLL0V5eTluu+02DlpRnIP8peRbrQ9kMZMAU6DL4xogHyptN954I99HVixN+ckv3NSZPRpnU9rDCkNnK3JYk19FWqAIvqSiooIT49PS0qqyUnYeL2rSPNTOydI37Vw/s4bqhl9bqILgT9DqJRK6hiTdn8trCU2PCKogNBH1WQoq+DcSlBIEQfASIqiCIAheQgRVEATBS4igCoIgeAkRVEEQBC8hgioIguAlRFAFQRC8hOShCkJTkbMbsDe+o2aD0BmAhI5QkhEjRmD8+PGNbqRJy0z/+OMPdO/eHf6CCKogNBUkphWFQOXJ6k0+wRAKmKMa9JBbbrkF77//Pnbs2IHOnTujqUlNTcUrr7zCAuzBU/LPnxBBFYSmhMS0LA/Q+uin57ABVIq0AYJKwkVV96ly/9y5c/Hiiy/6ZmxBgPhQBaHJf3V6IK6tb7ZGCDX1eKLaws899xw+/PBDVFZW1mi7TNX3ExMT0axZM7YiPVCn1CFDhrAQU3+na665pkYrk+r06tWrRp1VYuzYsXj++ecxYcIEHD58mB9P03xPrVSqZ0rFpz18/PHH6NmzJxcooYpSND61IYIqCEEOWaXU+YJ6MFF90v/9739V923fvp0rL1FFfSrbR/2gqMCzpxD0s88+y/2iqDFfRkbGaftSUfFqcil4oH2pLODkyZPx+eefcyk+EkyylutqG01jolqq//73v7l31dq1a1lc1YYIqiAEMeQzXb16NdcSJevwL3/5CwusB+qIQb2ZqEcTBZmotJ3HaiRBIwuV7iPrddq0aaetUUyCvWbNGi6PR5AlTP2mkpOT6zVO6pxKjQAvuOACFnKymMnqVRsiqIIQxJB4kjB6rD0SVqrQ7+nxVL1fFEGuAU8X0X379uGKK66oatt8/fXXn7Z/E7UioX090366pIZ69YV6SVFbE7UjgioIQQr5SufPn489e/awcNJGliQ16auPf5J8nS1atGArl4owL1iw4Iz9m2jaT5bpqlWr2NdK7aPr00eKIJ8pCbjaEUEVhCDlm2++YSHcsGFDVSfhzZs344knnsB777131uZ29FhqnknW6ZEjR7h1ytk6F9Nz3nXXXSzc1P/JA7kMPL7Zurj99tvx6quvcut46m5KLeUpKKY2RFAFoamh1Ka8/b7Z6LkbMN2nyDr1i/JYqLTdd999OHbs2FkFlTqkfvvttyyoNJ2/6qqrzrg/Re1pmk+iXXu6T62lX3/9dXYNkODWhvJT6fXuvvtubkfSr18/biuvNqSnlCA0VX+izG1Nm9ifpL726h9++CGnXpFVrDakp5Qg+BO0HJSEroGrmBr9WiqjpKQEs2fPrtMCDRRkpZQgNBUKr61Xkvnz53MQa9SoUZxJEKiIoAqC4HMmT57MW6AjQSlBEAQvIYIqCD7kbJFyIbA+KxFUQfABOp2OL61WqxxfP4HqGBC0lLaxiA9VEHyAXq9HaGgocnJy+Ad6tpVAgrKWKYkpLRaIjo6uOhk2BhFUQfABlMROhT8oF5XWoQvqh8S0du2ChiKCKgg+gpZWUkEPmfarH5pFnItl6neCOmfOHN7S09P5/65du+LJJ5/EuHHjlB6aoFLKrDaY9DrotBqsPpCH7ceKUGl3wO5w8qXDCWjoR6DVYHC7OPRpHYucYgt2HC9CYoSJt5hQI7Ra2qtx0FS/aqWUEPD4jaC2bNmSi9m2a9euqvwXrR+mAgkkrkJwYrHZWTRtdgce+3IrsooqkFlYgcyiChRX2PDj1GHomBSB/246hq82HoVRp4VBp4Vep4FWowEFdm0OJyLMehbUNQdP4O6FJ5dFGnQa9E+LxUe3DuT/X/tlL9ISwtC9RRRaxYby1F4QAmItP7VeoAo3VBasPlB1HCqsUFhYyAUdBP/CanNg+7FCbD5SgC1HC7H5aAEsNgd+e/QCvv/OBetB+pYYYUZSlBnNIk0Y2TER0aFGFly9Tluv18gurkB2sQXZRa5Ls16Hq/uloNRiw5h/r0BGQTnvG2nWo1uLKLx+bW/EhhlRWFaJyBC9iGyA0RDd8BsLtTpUr5HaJpSWlmLQoEGn3c9isfBW/cAI/gWJFwlbr1YxSM8rxV/eXMVWZufkCAxuG4+eKdEcpSVLcc71fU77PPURU8Ko16JlTChvtQkz6fH7Yxcgr8SCrRmF2JZRiJ3HixEV4kqzuXneGhzJL8eozokY1bkZzm8XD7Ph3P1ygv/gVxYqlesiAaWqMNSuYeHChbj44otPu/+MGTMwc+bMU24XC1W90Ndxb3YJftyWiR93ZGJbRhF6tIzCN/cMYd8niVin5Aie5quN3/bmYsXeHCzZkYWDuaUwG7T49G+DWPTrayEL/m2h+pWgUrSUuiNSk65Fixbh3Xff5YKzXbp0qbeFmpKSIoKqQihIRL7NPw/kYeLbqxFm1GFEp0SM7ZqEkR0TEGFWX/WkM7E/pwQ/78jCjYNT2Uq9ff46dh9c0j0ZE/qkICrUv95PMFMUqIJaG6pc07ZtW7z11lv12l98qOqCvnor9uZi3u8HOeL+wS39WVhX7s3h6XwgTZf/uykDi7ccx9Ld2RwMu7xnc0wb0wHJUSFKD00Idh9q9R9kdQtU8A8qKu34fN0RzFuVjv05peiSHIlbhqTx50lW6gWdmiHQuOK8FrxRWtZn645g0fqjVW6LLUcL0KFZRECdQIIVvxFUapFAOac0Zaeui5988gm3rP3hhx+UHprQACEl0bDaHXj+x90Y0i4es67sgX6pMUETGU+IMOHuke1w14i2/J7JIr/xvTWgaeLVfVNwy/lpnKEg+Cd+I6hZWVlcT/H48eNsfvfo0YPFlHp7C+qmqKISbyzdh0XrM/DTA8M4xWjVYxf4nV/Um3hOIGSRf3nX+Vj45yF8uvYIPliVjpsGp+KhsR35PsG/8GsfakMRH2rTQtbXx2sO45Wf9/Kqpb8Na4s7hrdBqNFvzuNNSnFFJd5ZeRDbMwrx7o19+TbKsxVXgLIEjQ9VUDePLdqKLzcexV97t8SDYzrKVPYskMU+bXSHqrxaCs498Olm3H9hO0zs14pzZAV1Ixaq4FV2Hnetl+/RMhq7M4thczjQtXkTNKULQI7ml+HlJXvw1cYMpMSE4tGLOuHi7klB42/2RwtVTnmCV6Cke/KTXvbab3hrxQG+jdbQi5g2Hlqt9fLV5+GH+4ehfWI41xj4fP1R+caqGJnyC+fMkRNlmPbZJqw7lI87h7fF1FEd5Kh6EToxzb2pHy8UGNI+nm/bl12CtglhYq2qDBFU4Zwt0xveW8NFRWiZJVVmEnzDqC6u/FyqbXDJ7JWcdvbMld3RLFLSrNSC+FCFRpFfaoXD6URcuAlbjxaidXwoIoM4Daqp+XF7Jh7/ahuXL3zi0i6Y0KelWKs+Qnyogk9Zm34CF726Ak9/u4P/794ySsS0iaEaBz9PG4bRnZvhkS+24MWfdjf1EIQ6kCm/0OA16Q9/vgXntYrGY+M6y9FTEKrz+vLE83BJj2Suy0o4HM5z6jAgnBsS5RfqzZvL9uH+Tzbh0p7JWDBlgOSVqoQLOzdjPyrVab149koOXgnKIIIq1Bsq7Dx1VHu8NKGnJJmrEJNBx21Zbpu/jlu1kLUqNC0SlBLOuhzy113ZXClJUD8korN/3cvLfS/qmoSXru7JnQaExiNLTwWvcKygHLfMW8uXlKJDEX1B3ZD/lPKAqSTiw19swa7MIm4+KDQNcuoS6oRE9Oq3/uCuoIvuHCxi6meM6ZqEQW3juD4AWa3UBbZ5tBSz9jXiQxVOgRLHr31nNV///I5BaN8sQo6SH+Ipjzhn+X5eEkwdYwXfIoIqnILZqEPXFlH4+LaBYtUEANf2b4UWMSG45u3V2Hg4X+nhBDQiqEKNABSty6cVT29c2xspsae2Uhb8j5gwIz66dQC3Wblh7hpsOlKg9JACFhFUgaEljLfPX4+b563l9flC4E3/593SHx2SIrgzgOAbJCglsIBO+3QzV4uaf0t/6GSlTUASbtLjw1v6w+QuVE2fu3zW3kUsVAGv/rIX3287jtmTemFAmzg5IgEM5aTqdVqsOXgCY19ZgczCCqWHFFCIoAY51Nb4nRUHuPXGRd2SlB6O0ESkxoeizGLDlA/Wcr8vwTuIoAY51NZ48X1DcNeIdkoPRWhCEiPMXLQ6PbcUUz/ZJMtUvYQIapBCfZ/eXrGfg1FtEsKlQlEQ0jk5ErOv6YUlO7Pwyi97lR5OQCBBqSDlue93Yd6qdAxuG19V+i1QTxzUirmSNocDNruTb7M5nLwKDKA/rnJ3FIvTa7XQ6zS8UTEY8jea9a7LQK1U9dxVPdCndYzSQwkIRFCDkB+2Hce7vx3Ek5d2CSgxpSWWpVYbyivtKLfaUWa1s4B6sJOgOpy80b6OWo8nydTpNNBpSFC1fOnBoNcg1KCH2ahFiEGHMKM+YKz6q/um8GWpxYb0vFJprHgOiKAGGZS4TwWix3VLws3np8LfsdkdKK6woaiiki/J6rQ7nezKqKi0w1Lp4H5XJy3S+kN6SpaqSaeDyaCFyVAJk17HQkv3RZj1vAiCLgPBgn3uh134ZvMx7rKaFCV9qhqDCGoQ9iIiAXjurz38tgeR0+lEUbkNJ8qsbFWRUFbYbCipsKHUakelzaWcZECSAHKqkFZTNZ2n3EvatLXeP/XIotxMsmr50uFApd0lzqWlrtchjHoNW6gUHadxEKEmHWJDjYgKMfit5UqZHvT9ePiLzfjg5v5++z6UROqhBiGFZZWICjX4pTVKInqi1MqiWV5pQ7HFhtIKG+wOgIzEMKMeZoOON6M7gd1bIk7iSlYvbeRaoNfUul8z3Kx3uwGA2DAjbyTm/sbyPTm48b01mHl5V9w42P9nMN5A6qEKp0CW3M87s3B5z+Z+J6YURMoutnCnVbIcSUQLy6yw2pww6DS8rDKMhVTrM6ubnpcsUxLpyBADCywFu+i4llhs7G6gsdB9FADLLbbyTIBak4QY/UdYh3dIwA2DWmPW9zu5V1W81MBtEDLlDxJe+XkP5q8+hP5psUiO8o+6mBQ4yi2xsJiSdVpQVonCiko4HLTiR4f4cANCjcp8hUlgPZYwFd4mq7WwvBJ5JS4LmpZ5WmwGFlpyAzSLMvmNxTp9XGeM7tJMxLQRiKAGAVS1/b3f09lH5i9iStZoVnEFB5RISGmqTxlOFASKDjXAoLIgkEdc48PJv1vJ4kqWK7kCyMKmoBlVfUqMMKlu7LUhi3po+wS2wrdlFHGbcKF+qPuTrcasWbPQr18/REREIDExEePHj8fu3dKLvD5W3v99tQ2pcaG4bWgbqB0Sn4O5pTiaX478kkoczitjqy/CpEdqXBiv7FKzIFGwi4SzdVwo4sKMvLzz8Iky5BRXIKeoAnuzSlBAJwc/4KuNGRj/5u9S7q8BqPebWYvly5fj7rvvxurVq7FkyRLYbDaMGTMGpaWlSg9N1fxvyzGuIvX0+G6q71RKVp1HcDIKyrltB4kndfKkpZL+VBmJXALRoUa0ig1DTKiRswEOnSjj93bkRDkO5ZXyyUPNXNazObo1j8SDn21S/VjVgt9G+XNycthSJaEdNmyY16N1gQL9EKiy0Pnt4qFmK5oElKb2lPqUXVIBDTQ8PQ6Ujp0UTKNCNOQGIP8vWdp0gmsRHcLCq1Z2Hi/CJbNX4olLu+Dm89MQjBQ1QDfUbbKcAXpzRGysdHQ8HZQ/SRaemsWUBP9AbgkHcrKKKtgqpZVIZJUGipgSZF1TsjxtFMAiNwClr5G1Sg0R1WrX0Hr/if1SuC21v7gqlMQvv7H05Zs2bRqGDBmCbt26nXY/i8XCW/UzTbBAFhE1ZpvYrxWmDFGnZUHCQksdyyx2HC8s5zxPskop9ShQoeh/iCEMOSUWZBVZUG6zwwlXfiudRNS44uqB0R24fYpSGRX+hPo+vXpwzz33YMuWLfj444/PGsgiU92zpaS41iwHAz9sy8SerBLVFr2g/lX7c0p4in80v4xToVJiQgNaTGtYq5FmnvYXldnY3UGZAftySlhY1Qb5r2m6Ty4Kcs8IAeRDvffee/H1119jxYoVSEs7s+VVl4VKohroPlT6SMe9upJ/sPOnDIAag09UU4D8iVQxnqo6JUeZVWmd+Roq4kJuDo3GySltlHqVGh/GlqzaePXnvdiTXcwNHIOJokD0oZJIkGX65Zdf4tdffz2rmBImk4kPQPUtGPhlZzZ2ZRbjnpHqKxpN+ZgkpsXlLjENNeq4xXEwiqkn57NlTAh0Wi2nitFJhoo+kwWvxir/i7ccx58H8pQeimrxm28xpUwtWLAACxcu5FzUzMxM3srLy5UemurYm12CQW3iVNcfipZpUl4prR6ipH2ywmjqW7tISbBBgUOK9lMwjk4ydHwO5ZXxyUdNXNajOXqmROOfi3fK1N/fp/ynW6P9/vvv46abbqrXcwRT2pTaOlqSb5B9phYbjhdU8Nr7ZpEmv6145Qvop0iBKiq8QjUAqBYALRCgWgVq4Y/9ebjmndWYd3M/jOiYiGCgqAG6oT5HzWnwE91XnK1HC9G+WTj74tQCrcOnaD75C48XutKiRExPhU4udFyyisApZBqNmS3Vdonq+TwHtonFeSnR2H6sKGgENSCn/EL98k6vn/snXvtVXf2ByDdYwWJazgEoysUUy/TMokq+ZZr+e1LL1LJSicb3+R2DcLcK/fNqQAQ1gFi6K4cj6OPPawG1kF3s8gnSVJYmGeIzra+omjlNiZL+ybInS1UtKUvk8yWB33A4X+mhqA4R1ADiyw1H0b1FFNo3i4BaglDZRRacKLNwfycSCTUXNlETFKhLjjTzElyy7OlYUr6qWliw+hAmvbWal9MKJ5Fvd4BASzeX7s7Glb3VYZ2SNXUkv4wDLPmllYgJVa52qb+id7tHqJA2nZi4jGGpOpZ/XtmrJbeToRq7wklEUAMEKsTcq1UMVwhSA1QUmhrk5RRZYNbruCWI0HAoGEWWPWVHUHFtcgGQr1xpqOsDdUslS1WNq7uUQgQ1QKC11p/dPkgVVdbpB0YCn19m5fX5tGJLglCNh/J1I8165NJJymbnIJ8asl6oay59xl9vzFB6KKpBBDUAoKpFa9NPqCYSTL4+ElWaotJUX+11WP2B+AgTd24ly5/8qVRcRWlax4VhyvlpXFBbcCHf9ABg+d4cTPjPH9w2RGmoxBtVjyILlQRAfmzeC1JRkZIKq4OtQvKpqmGq/X+XdsHYrklKD0M1iKAGAKv25XLyd2KkWemhuCwoqw3lVgdP9YN9Wam31/2TxU9BPuq1Rcn/aoDSp2iGJIigBgSr9ufh/LZxqnA9UCCKItHU0lmi+t6HLH46SVGfLWqrQsEqNVShevFH6e9GiIXq51DlJqr+PqhtvCqS+Mk6JVGVqL5vIDGl5n8kpBU2qtilfIDqkh7JWJN+QjUWs5KIoPo5eaVWblNBa6yVhCojVVRSy2crTHqxTn0JFeGmY5xTbGXXCgX/lGRslyT2l3+/9TiCHRFUP4cKVXx//1DFG73RFJQsJvqBSyDK98SFG3kmwE0NFV6tRDmpQ9snYLEIqgiqP0NTPSpGrPSUj9K16IddWG6DQadRZbX5QIP801S1K7/cygEqquGgJNf0b4X+abGKfxeVRixUP4bWUY94cRl+3ZWt6DgojcfhdLJfT8S06YgONbCVWl5pQ57CeamjuzTDw2M7Bf0CDhFUP4aauhFp8WGKR/cp2dzpcPn3hKaB2mzTjKCwzIZSi52rUinJvuxirA7y9igiqH7M/pxS/kGlxIYqNgb6EVMwikr0UaqUVJNqWsh3TpkVVrtrMYWSvLPiIGZ8sx3BjAiqH7M/u4SX/ykpYtRMzu50orzSrqpWHcECtUmhpRPU9JD8qErWTO2bGoPdWcWK+3OVRATVj6G8v7YJyk73yW9KVirFIqjKvND0eankt6bPgT4DmikoRb9UCkoBGw4Fb+FpCcf6MXOu76NoKTeyhqhwNAVFyPUg031lCDfrUVRh4+8CWYeUxqQEreNCudoZLUMd2Sk4+02JoPo5Jr1yViH57sgiIVGlFB7V43RC46DpqJOvu2+EUx/C1zS2cmicDjhpEq1xTd6cWgOg1YEjbny7+moT0LEnrw+lrpkNrmm/VoGOtxqNBpf1TGaBD1aC9537OduPFWLqJ5vwzg19kapQlJ8iyzaHA5U2J2JDz1FQSeCcdkCrh9ZaDPOJHdBZCl2btZA1MK/7rbxr0ppnoC/PgcZmgdZugcZuQWb/6aiI64q4be8ibucCvs11fwUK2o7HsSHPwFSwFx0Wjarxsg6tAdtv2c/X2/7vKoTkbatx/+EL3kBhm8sQv20ukv98Gg6dCU6diS+LU0YiY9iL0FqLkPbDDTXuo8ujQ5+H0xCKmF0LYSw+DIc+FA5jJOymKJQl9oI1MhXaylLe6DZ6TGOFjCL+xZZKxIWbeNqvlJX61GVdEcyIoPop1LRtb3YJohRMU6JAFEX4ibosVH1ZFgylmbzpyzJhKMtGXufJsIUlIW7bXETv/4bFyCOaeV1vRuaA/4P5xC60/XZC1fM4dGZYIlOrBNVUsB86Sz7fzgJmDIdT43p9a1QbFLUaXUPYKmI78X2Voc1wZNhLbitTy5eexxGZ/R6FzlLgsmTZinWgLOE8vq+4xTDYhzwLrd3qEmu7BdaIVlWPrYhuXyXutOkqi6us3NDsDQg//ge0lWX8frWOShwdMosFNfLgd0hZ8WDV+7SbIlGaPBhHRs7m12+5fBrsJMLmGB4/HbuS5MFw6mtWFiM/KhVLoWl/iVU5QXU6nbxyi4JlwVgcJ/jecYBQ5I6kKpL3SdakRgNbcRai9/2MyMLjaLbvBAxlWXBqtDhy4Rzerf2iMdBbXAEKp0aPytAEFKaOY1EgkaiIaQ+7MYqtM9rK47vzvuXx3bDnr7+674s8xXI7NGbuaYdWnHIBb3XhMEWhoMNJoa5NScvhp73PEtuRtzqf1xiJjGEvnPaxZMVWQWJtpyIiril5SYshSB89l08oHovcFuryP2psFTAWH3HdV5EPfXkuC/2O6zbArjcj5dd7EJqzEZVhybCGJCFWG4vS9pehLGQAUF4AVBQCkc0BnaFJF5sMeOYXvD25D8YEYZ1UEVQ/haZ1ZJXofO0rO74ZSP8dyNvn3vYD3a+C48KZ0OQfQuc/HkalIRL28CRUhibBGpVW9dDDF/4HdmM4bKFJsIXEVVlsBAnb6cSNfJqW6ADt+05WsdtnS9jCklEcllznruQuOHDZopM3OCphKMuB3ewq1ViUOhaV4c2hL82CoSwTzYs342B8V1S07AfHjv9C+7/72IWCmFQgrh2QNhwYdJfrhFicCUQked0nnBBh4tkKVUALRkRQ/RTK/6RplVcoyQaObQSydwI5u4DsHcC454FWA4Fdi4HfZwNxbYHYNkCPq4E2w2GxOdii/PGKDdCZwrmRXG1Kmw9Ck+D2v2rIB+twX9L/nus1glBVD4KmWmCqJuQK8AhNLcFxuwucWh1b49Do+TqdLMgK5wCWr9AaWEA9kG+XturWYZnVhtZ07FNHIeL6RUD+IddJkE6GhUdPft4vdwLMUUBiFyChk+uyz02A/tyK7Gg0GrSKDWWXVDAiguqnjO/VAv3TGlFUmkTkxAHg6DqXOJJALLgSyNwKGCOAxE5AUg/A6A50DZkGDH8M0NZMWbaUWfkHXq4xI9abCwucDmjIT+mgzQE4bW6BdEDjqHmdL93iWfdz0eYSVE2VZtYtrjWpJqK1hNXJF9pTjkfNh5PY6l0iy8JLgqtzC6/r/6rrOiOcWqNXLEWyDCltigKFpfo4RLSrGYCrwhQBTFroOoHSdnQtsPMboP9trvu/vtu1T+tBQKtBQHjDUqBaxYXikFiogj/RJiGct3phrwTWzgUO/+HaSrJcAkE/mOhWwOWvA6GxQFTKqT9sQ91tVax2B+wOB2cTNTj/lEXTFbxxBXlcAqq1WQBnZR2iaHOJK1wCqqHVQHQbSHAdbqvUfb/DLbju/30LWaVu61SrhROu62yl0m1uS9ZlvZI167Fq6bbqz6NhUWVxdQusg6+7gmr1FVuTwfU50OzhjP2mjKFAp0tcmwc6vvQ6dLKhk8Xu74A/Xb5wxLYFrvkYSOgIWEpcJ9szjKl1bCh+25eLYEQsVD/lv5syWMgu7l63/w1Z24FDq1xWB/nR/njDFaA471qg1WAgpT8QEu3at7krkt0Q6PdHS06JMxlrJG5aW7krPchWxqlMGke1NecOEkErtHayPq0s/pQrWmWNnmI9VntqTW2x0sGpN5xyu+vH77Yw6U8d0/kde/bhXy/NweMP3okuHdqdfF33ezy5p7Oai8FR7dJtLdPY7dVF/cxi7NQZ3LmuBji05D4gITXUGCNnLOjNrrQrQ/hpRdZzYqNyinTCaxCeD5Ge9/LXXNeLjrm+Q3QSjmzhuu3LvwH5B4EeE10zHPpO1eLhizricV1nBCMiqH7Kp2uPcJuRGoJadgLYtBDY/AmQtRUIjXd96clXdv8mr/r3SEw9M+cajficTmgrS6ryK0lEWYQcDs4JJYsUdiu0NKWn69WsSBJAEhfaHBS4qfJTei5PWoTVA1zeYG/6MazZsJUvO3fpdsr9jV4h73FRVLekq4sxnTjouFD6WI1X0cNBJwetyX08yFo1u5Rdo4NDHwaHIRR2Y0SNIBetWGNBtXnBOiex7P5X1+ah3y3Axo+AZbOAn2ewPx3jXgASOqhisYnfCWpRUREiIyOhBCtWrMALL7yA9evX4/jx4/jqq68wfvx4BCO0dpuW+lVhswJvDgLKTwAdxwEj/w6QD80TZPBysIRW43gsVMo0oBQfyuHUV+S7pu0OG7S2CvdWXs0q1fB01kGWmCHMLaB6t1UWgD9Ej/+UrtOCqzPtS4LLvmOycl2WutZW6l7dxU/myq3Vh/AJh/JW9eRq4PzVGNjMMWyl0kIL+mhIWL2+HJi+U7RRStaO/wJbPwdC3b78zG0c4MouteGm99fi6fFd0ae1sq15VC+oMTEx+Oyzz3DVVVehqSktLUXPnj1x8803K/L66kPjis5TlNYQ4pqqNe8FhCf4/JXtDidP+8niDCnKgt5eyoKgsxbz5hFQEkuHLgQOyifVh7CPUDgNZJFrQ+DESYvzZO6qa9UXnbi0nLPqyu8lNwDl9NLttHgiwhaKEwb6/EPYSvVZfQWa9fS+wbUR1lJg3sVAXHsYLn0XO44XcVucYEPfmJUQc+bMwfPPP88pEn379sXkyZMxYMAA+Jpx48bxJrj8Ze0LVgLvPgmMeAwY9hDQYUyTHRqq0K+1FCKi5AD0Gg30FSegs7kKXjsMYWwtkd+P/beCF3JXzZzMD/caB3KXsG/aWsw5qGT+0vTfjAiElBVBE9UZNnsT1smlQNV1XwCf34zoBaPRWjNdFS2um5pGnb42b96M/v37Y8SIEdi9ezeGDx+OBx54wPujE07LqDQzJmW9BLS9ABh8X5MfKZqGmkuOwlBZDEPxEWgdFbCZ42CJTOOVO7QcVMTUd5ClT6vLKiNa8hJYmymC/bAhZUd5hkC1AygLo0lJ6Q/8bRk0hlA8a5iLkorgq4vaKPNh4cKFGD16dNX/W7duZV9my5Yt8eCDrnXJasBisfBW3f8bKNzZ+hiwOhe4aNY5J2M3Bp2liKPa5oosOMzRsPOqG/X7QHfs2oN9B9JPuf23P9bWuKxNuzap6NLpZOBFdeIaEg+HMQq6E4dhtuTCaYuGnabhHpO2qQhPAC58Es0X/R+2lRYg2GiwoMbFxSElJaXGbd27d8fs2bMxdepUVQnqrFmzMHPmTAQixQ4jIjypLbSKqYnRarWcYeOgBHbOJa2E0w+iu/968TWsXb/5tPd//e2PvNWmX5+e+OhddzqRSqFsAQ3scGqM7JrjzAgl6P5X7DcMxchYZYuf+4WgUlBo7ty5HG2vTrt27XDkyBGoienTp2PatGk1LNTaJwN/Zca2OJwXMQWTk3soM4DQGGgMJpSHJMOhKYCx5ChHnilAQj5Ub6c1eYvHH7r3tBYqCen4S8diyKB+dVqoqoQCgZXF7EulalflhnBUaCNhM1I5wLoXZfiUA8uA9N9wwdCHTrsoJJBpsKD+85//xMiRI5GRkYG77roLPXr0QHl5OZ555hmkpZ0sjKEGTCYTb4FIuNmEj3TjMZmirem/AYUZQPcJZ8my9x46nY4LoVQW2VEeEgqtkwIkha6KU9CyqDoMFNWnTT29pmjafrqpOwkqienlFzddcK/BVIv4exZMEJSXWhkWBytMsFpCYQ1zJ+I3FeReWPUasPw5LsLy5fp0NIuPx/nt4hFMNFhQBw4ciNWrV+P+++/noBRNLQiz2YzPP/8cvqSkpAT79u2r+v/gwYPYtGkTYmNj0arVydqUwQBVRa/qH7T9K2Dtu8DvrwDDHwE6XuJzvyplGegMISgJT0O5oRgmFENjioDGZnVZTJWlrpqg7hVNZL1SpJrzT2mljy+LiAQSlIvqrsHK+by00sy9GILrp5rjYaMaDJRypQ9FuSYGZToTtHp9zQUXvmTjAuDnmUB5PjD0QWDEdLw9+3cMSLOJoNZ32r9s2TJkZ2dzkr3D4eC0qfh4356N1q1bx9axB890/sYbb8S8efMQTMSEGpFfZuUTmuaSl4Ce1wC//AP4/CYgJAa44RvAh+4Ao17LP1i9XoMyYxxCwptDay3h/EinxQxQuT5eHVXOyf3UXkRXfnJ9t2vZpbFqZVTAJ/ifCVotxcttTyb0V12vWkmmdSXzczqaK6nftWpKzylqlNhPJywr1cm1WKDTaHxX2pFqQ+z7GWjeG4ho5qpkRSl7wx911YZwL3/VK9iNVynOKUkwMTGxSfNCq1vEwU7LmFD+0uaXVfISVLTsC9z4jat60NYvXIUsiB/+DuhNQOvzXWktZu+scjO6fyy8MsfduphSpWirDG/hno7SEtQy1/JTXoLpyp90rd23nLS+OBptP2UJavWlpqeu23cvRXUXH1GFCPMaf8cpa/xPt+6fywt6elxVQQsh3Gv6DeGuIim0/FTn/qlqDLAbQt0ulTD3clRNjfxgj/vaq3panAUc+t3lXtrxNVCWB1w2G+hzI3DhU6fUFiiusHmvvKQfEXzvOEAY1TkRu54ed6oVktgZuPAJ13Vef1gKbP0M+O1ll/A06wqMnwMkdXdZGo30b5r0rl8tWSGnrBvXeHyoYTWq1LO4cpUp15JUqhhVs0iKe7mlg9b7u4qj8Ebj9IjRGSpI1V7vX7Xuv1pzvZPS5T5u7ttpbK7LUugqTtRo4ld1hKvf5l6j7yol6F6vf8bqVjVPAlSZykluGW3kaYqiwH0bFUVxt3oxhJ617xQ36HNfb3SjPnqfVEP18CqgwzhXKtSv/3BN7WPSXAV2ekwCktw1D+pwLRSzoKrHd95UiKD6KfWaTtEX/bJXgUtfcf9A3OX7wpu57v/mXuDwaqD1YJeFS0WGSXCpFuZZn1oDg55aR2tQWuGoV5V6O211rV2vVcbPc3mmGqc1Kjy5rcLqIucpRMJi7HrQycdXu+Kpk9qpRTQG9OyITi1ioa8octc95T2qXXiuuytFkWCT6OlPWs21y/Tx5WkDhTXL9rl6YLnK9/ES3UZkSlB9Bc9JVl8fQXW3s2E2zAf2/uj6TpTmuF7/mqST0/kLnnBV+T/rUzpxQadEtEusZ3nJAELjDKI5NKVNRUVFobCwULECL97k1g/WYkBaHG4b1qZxT7D3Z2DvTy6RpXJ/JEiXvAz0m+Ka2u36Dohv79qohQYJcTVr5FBeKTKLKnAsv4KrtJNftcmq8FcVmXYLZ42C00255FFz+uLR1QpL11Vk2heBuWMF5fwRJUeFoHNyRM0TL51cDix3+TypMwNvu4EHtrmWjn56vatiGRWVplq5Lb3nIgoW3RAL1Y8h/ykVoWg07Ue5NoKKO+fuPWm9UhoWFxk+5O5JTzlHVwBXfwhUFAErX0KsMRFObSwqrFGw61KA2DoKVDcWeh5qL0Jf0bNVaaqLaiJ8sqZp9Wr9tS49t59SK7W6u+Ckheqauruq8qsCt6WpKc5Ey6L1iD1yAvpt+UB+uitIOf4N1z4LKbXO4FoMQu6hdhdyUIyZuMArQ7HaHMgrtSA+3OS74iwqRQTVj6HK6Ol5rjzEc4YCVx6fGNFzomsjoeW+RPtOugKowdv2LxFenIkIuxWpnv72N7tS2qj1Ma0rp+Z8lWHUvK8ZSloM5W6nVBWJXQCN7EFfb9x+yur43VSMfM+UHWEt5LKIlui2fNyi9n2N8OOrYCjNcrXnLs1ETo87kNX9DkTkb0WXNdNc/muqZ0pR92j3Yhby2U7dCkQk+zRtbW92MS6Z/Rv+e/f56JniLmIeJIig+jHUu2fF3hzfvggJLRUPrlZAmK/TD9PhwM4D6SjOOQxbyQlEuC25ypAEhJTnITRrLf/YqZX0wYsWoCQsCQmb30Szja+4e9C72kcXpl6M7D7ToCvPQ+LGV6tu97SYLm41mkVYX3qcn9/TGoRTh1Scz0pi6ArEVbj8xJVlsIXE84nFWJiOyEM/ulpHu1tIU2fY44NmspB2+rg/B8e0VbVQgT1XLYElpiPMJ3bCnLedT1hlzfrxSas0eSBnfWQnnI/1125Gs4REzgQ5haiWPn/fxe78aInyC35FalwYckusKKqoRKQCEVWNVgtTVCIKNFHIMlXA7C5onNV/es39bBUu3yF5EtIuhjUiBTqqUO8WE/qfoNvI8qLCK3Q7pV5RoGb7zXtd7/fHmxFyYkeN5z504RwUpV2C2B3zkbD59ZNiqzejpPn5yOr3GFvLLVc+wiu4XNN1V0uUo0Of59zN+M1zYCrY53IJuN0C+R2uRmnzwQg7vhqxOxecTHWyW2CJaovMgU9w8KzdV+PcQTVPjywLdk9Yxu2hyVKPPri4xngz+z6CnPPugalwPxI3znafOCJdl+YY9wHTILfbrfwePCcVyjO1RlA/U/Dxpa5gtbGUUz0FE/QhYYpWzT9Raq3KlQ42xEL1Ywa3i8OHt/SHWcEfj2fFFmlUSYUNMZQTWwsSLQ+W2E681QUtZd171c9V/7NAWV01VomMIc9yVwCPeNFWHu9avFAR15lFsPp9ZMHx81D3ACp6zb5gT0+okw4AqidqKjzgflHymWo5h5b/pcLNFbnu1tK04stc1XKEcmXJlVGV1uQRc3e6WG73v6Gg3V+qiXwIrOEuC7G41YXYceP20x7X3B63o6GUV9ph1Gt4wUWoUbnvxIGcEkSHGngLNiTKL5wTFpsdezJLkFlYwVPOlNgmLGos1OBgbikizXrER5jQtXkkp7YpwfQvt2JPVjEW3Tk4ID4hifIHEUt2ZGHL0QI8OMa9MqqJoalliFGLcJOeU6h80sdIqFdkndrSmI06/iyUElNi1pXdUW49QxvrAEa++X5Oem4p3ll5gIVMKSJDDAg16aqm/ULTQ9N9Ov4hBh3CTMp78kIUdDkoiQiqn9M3NQYVlQ5syyhUbAxRIQb225FlRAEyoekptdhgNmirPgelyCyswOiXl/OsKRgRQfVzujaP4h/SunRXF0ylpv1hJh2iQvWotDv5xy00HTQ7KbPaEW428Go1Ja3D/Tkl2JtdokjWiRoQQfVz6AfUKyUGa9NPKDoOCoSY9XqYDFoUlImV2pSQm4Wm+2SZxigcWd+fU8KVyFrG1GqFHSQo72wRzplbhqQ1fYfLWpBFQmIaE2Lk4FRFpR1mQ3D60ZoacrOQmFIN1CiFBXVXZjFS40ODshYqIYIaAIzu4l5/rzBxYUYWUn0pUFBeiSQRVJ9Dx5vcLPEReg4MKpnQT/x5IA8D2sQhWBFBDRCW7c5GTrEFE/oq14SQVsZkFVn4MqfYiooQsVJ9DblXqIQiJfKrYWXS7Gt6KS7qShKcdnkAsmJPLp7/cTfnIioFFTROiDBxGhXVSs0rOU1NU8FruaclFhuiw4yc+xsdYlBFkLRdENZB9SCCGiBc0iOZLVTFg1PhRrZQqHQb5UbSD17wDSfKrCBjkFdHhZsaX6HfS7y5bB/m/3Fqi+5gQgQ1QOjdKhotokPw7ZZjio6DVugkRZkRZtRz+k5eiUX6gPkoVYqi+9GhRui1WvZfKwnVqf9o9WEcyPVSOUk/RQQ1QCAhu7h7En7YlgmbgqumPIn+FCAha9XmcHIhbMG75JZYQIF0yq6g46y0dXrkRDkyCspxflvfdj5WOxKUCiAoIEU1MEnElI4LNI8K4fXcVHGIyrlR0ETSqLwDLZwotdjRLNJVET8u3MfFuuvBqv253GW1f5tYBDMiqAFEh2YRvKkBmu5TgIraGpdZbcgqquBKVLQ0Umg8dDzJOqXjS11Fk6PMp3a+VQAqdN69ZXTQrpDyIFP+AIOE6/++3oqCMuUj7IkRJrZMEyPMbDVL1P/cyS+18rFMiDDyct+66s8qwS3np2Ha6GpdHYIUEdQAg6yVz9YdxcI1h1Xh1yUXBE3148KNKCyvlKj/OSbxkz+a3CiUSdE8Wj3LO/umxmJ4hwQEOyKoAQalz/zlvBb4YFU65ykqDYkpRf2jQ4y8PJIsaBIGoWFQfjFVcqLlvbFhRv6c1eKTnvHNdizdna30MFSBCGoAMmVoGq9YWrxV2RQqD/TjJ6sqMdLExVyOF1Yonongb9CJyAEnkiLNXO+UAlJq4MiJMsxbla4KF5MaEEENQCgwNaxDAt77TT1J1pQjS0JAQRSNxsmiSgEWoX5+UyrP1yzSzFZpSkyoohX5q/PVxgz2k4/t6urfFexIlD9AefzizrzGWy1QnmRqXCj25ZQgOSoER/PLeQrrElj1jFNt0EqzvFIrYsIMvFgiJTaErXw1QMn8X244inHdkhFqFCkh1PHJCF6nY1IE2iSEq2qVEpV0o9bXlPKTHG3mpalkqappjGqC0s3opEO+57gwE7tMKFVKLWw8UoD0vDJc2buF0kNRDSKoAczerGKMeHEZdmUWQS3QlJVENcyo52CVR1Rl+l+TCvdxoZNPsyiXD5qm/Gqic1Ik3ri2NwYGcbk+vxfUN998E2lpaTCbzejTpw9Wrlyp9JBUS+u4MNBk+l+Ld6rKCiRfamp8GFteHkv1WEG5opWy1CamdDwook8uEUqWV1sFfPo+kdhTUR41LCxQC34lqJ9++immTp2Kxx9/HBs3bsTQoUMxbtw4HD6sfM6lGiFf2/SLO2Pl3lws25MDNUFimuYW1RYxZljtDhzNL1NFqpfSPlNaE29kMQ1BhFmPVrHqCUJ5+L+vt+Hpb3coPQzV4VeC+vLLL2PKlCm49dZb0blzZ7zyyitISUnBnDlzlB6aahnTpRkGpMWylaq2VCWyVNvEhyPM6LLASDOO5JcFbfJ/fpmVfaZhRj1nRVBZPnKPKF34pDZkPX+27ohqUrfUhN8IqtVqxfr16zFmzJgat9P/q1atqvMxFosFRUVFNbZggyybJy7twlM08smpDZo2UkFiqlBFq6rCjHoWFSqooiY3hS+h90m1bGlpLjXZI98yVd8nC15tYkq8veIAnwyvG9Ba6aGoDr8R1NzcXNjtdjRrVrN/Ev2fmZlZ52NmzZqFqKioqo2s2WCkW4soLHlgOBcnUatrgixVEhESk9gwV4UqmvoGuguA3h+lkFGjvYRIE1eOomh+igqn+QQJ/8drDvPafRJVwU8F1UPtLxmd3U/3xZs+fToKCwurtiNHjiBYIUuHov6L1h+FWsfXKi6UI9qxYSb2q1KQilwANBUORGuVahvQ+6MVUOTyoBYmdKm2aH51/jyYx7OKGwelKj0UVeI3p5j4+HjodLpTrNHs7OxTrFYPJpOJN8HFt1uOY86y/eibGsMZAGqEKlNRknhGfjkXADlR6poKk181QUXr18+12j5ZerT6KTJUj/gwE4sUrYCiSzVzaY/mGNkxUaxTf7dQjUYjp0ktWbKkxu30/+DBgxUblz9x+/A2bAE+9PlmVacoUeS/fWI411ONDzejZayZ81R5dVVRBQuSP0LHnGqZHsorg9Vu55SxxHAz4iNMaJcQrnox/Wl7JrsoZKofAIJKTJs2De+++y7ee+897Ny5Ew888ACnTN1xxx1KD80vIMvvpQnnYd2hfA4sqBlyAVB5ujYJYYgKMXLqEPkYKyptLEjZxf5TYIXcFVQ85PCJUp7mk4+4VSy9LwNax4dyRF+Nwafq/LY3F3+bvx4/7ag7XiH42ZSfmDhxIvLy8vCPf/wDx48fR7du3fDdd9+hdWuJNtaX/mmx+NuwNnhn5QFMHtSarUE1E+a2Vmk9u77IggiTnkWJ/KpF5TaEm/WIMhtUad2R4NNYKeBEEwIaKzXTo7YlLuvb5BdJ8RabHU/+dxv6p8biku7JSg9H1WicgejtPw2UNkXRfgpQRUZGIlihH0huiZUtI3+bMlMX1exiC1etL3KLVaXNyZkCZPHRCUJpkaJeWoUVldyVlOKlkSF6Xu1EPmHPElK1FDipD6/9shev/rIX390/VDUtdtSqG+o2TwSfQD9sElMSoxV7cjjQ4A+QUCZGUlqVkU8IRp2WU61KrTYWV/JPUrAnxKDj9iBk3ZI16GvIJqHlsyXu5nkk/FTpKz7CyMVM9FoNCypZpWq0pM8E5QS/vnQf19gNRjFtKCKoQczXGzPw1DfbuZLRoLb+U+CCqlZRvir1rKKTAolrmFEPm8NR1RGUXAQsunotzHotTAYdTHSp155zficJJln5FZUOXndPG03pqdMsLRUlIQ8xuCxlEn/a/MkirQ6thnp10nlcX1c4OzLlD2JIGK55ZzWnKH177xDVNHxr7DSb/KrFFTaORLu6rdrZerVU2lFpd8Lj3CJR1es0LHh6rZYvaavuKaB9aSORpuNELgabnS4d/FyERguY9TqEGLQsomT5k1aT24HcD7SpPdh0tmr8al0MotYpvwhqkEM/mstf/43rp354ywC/taSqQxYjWa4UtCKhJUhgKVXJQlalzQG73Qm708EieaZkARJI8hqQVazXkBC7XCZUCcqoc03ftVrwlJ42slD9WUQ97MsuxuWv/44Zl3XF1f2Cc4WhB/GhCvWGLJC3JvfF9XP/xO/7czlp29+h5H/aEiMAh8Pl3+TN6rokUa2OE07ej+xOjxVLQkquAV0t9wD9S35Q8tNS6w/PawVaYes7F2xgP/ulPSWq3xDEhypwKtXyh0dwubhAg6xFmo5XT0anIJJnCl/pcFmpnpxWT8oLy6gGMGhd7gEKbtGmdAaBr6Fj8/hX23gRxf/uPV9amzQQEVSBITGlH9Oby/Zz3ueYAG66RpYnReHJsAxBYFmX58r32zK58d4rE89DOzLxhQbh/w4zwWvQdHfr0UJM/XQTth8rlCMbhIzq3AyvX9sL43tJn6jGIIIqnPwyaDV4eWJPXu556wfrkF2kvvqpgm+gFV3bMgo5KOkveclqRARVOGW9/7s39OOo+A3vreGe8EJgQ1kRf/twHW6Zt5avC41HBFU4BUqanz9lAK+MCbQItlATyrGd+skmbD5agDnX95bP+xwRQRXqhMR09jW9OEWIchILyyrlSAUYFISkoidLdmbh9Wt6o0/rWKWH5PeIoApntWDuWLAB1767mguTCIFDTokFS3ZkYdZfumNUl7qLtAsNQwRVOCOUd0lR36wiCya+vRpZEqgKCGh5LnVH+OXB4UG/EsqbiKAKZ6VTUiQ+u30gyiw2TPjPH7xcVfDfaf7sX/bi+nf/5M4HVA1L8B4iqEK9aJMQjs/uGMTLLQvEn+q37psZ32zHy0v2YFiH+CYpbRhsyEopod60jAnF4vuGshuApoy/7srCRd1krbe/rM+/7+ONWLo7B//6SzdcN0C6XPgCOUUJDcKzlv3H7ZkcrHr8q60sroK6oeDTH/vz8O6NfUVMfYhYqEKjuKxnc7Z6nvh6O/ZkFePN6/pwRXpBXZwotXKB6yvOa4GBbeK4/YrgO8RCFRrNxH6t8PHfBiI9z1VTlX68gnr4fV8uhr+wFN9tPc7/i5j6HrFQhXOiT+sYrvb/v83H2BKiKDIVWQmEIsv+HHz6z/L9+PeSPRjcLh5D28crPaSgQSxU4Zwhy+fWoW34+qINGZj0zmqk55bKkVUAalQ46e0/8OJPu3HbsDaYe2NfSY1qQkRQBa/SMiYExwvLcdGrK/DuygNsLQlNB7VhoU6wn/5tEB69qJOkRjUx0lNK8DoUrHr+h92YtyqdXQJvTe6D+HAJWPmKgjIrd6+9bWgbdGsR5bPXCVaKGtCkTyxUwSclAGdc3hWf3T4IyVFmRIe4VuNQ3ybBuyzfk4Oxr6zAst05vDZfUBYJSgk+7VVFG7HlaAEe+HQTHhrTERd1S+I2JELjOVZQjv/7eht+3ZWNIe3i8eKEnlx2UVAWEVShSaAuoS1iQnHnRxvQs2UU+/coAi00DGomSC2tqU5tdnEF3ri2Ny7uLicotSA+VKFJWbU/l/2rm44U4OkrumLyoFT5BOoBVdJ/7/eD+Gj1YU5Ti3GnqImlry4fqlioQpMyuG08vrorDj9uz0K/1Bi+7ZedWWgdFypdNuuAKkJRju9LP+3h0onXD2wNrdtdImKqPkRQhSaHhID8qARZWSQWO44X4fx2cZg8sDV33qRprQA88sUWbus8ukszzJ/Sn6t+CerFb6b8//rXv7B48WJs2rQJRqMRBQUFPjXdhabDYrPj+62ZmL/6ENYfykdSpBnf3Hs+F0AONqh997zf03Fxj2SM7JiInceL2CLtmBSh9NCClqJAnPJbrVZMmDABgwYNwty5c5UejuBFTHod94GnjQTl5x3ZSHDnrT77/S4M75CAgW1iA3aK6ymF+N7v6Vhz8ASaR5kxomMi39c5WU78/oTfCOrMmTP5ct68eUoPRfAhXZtH8UZQC+slOzJ5XTrls17YOZHdAUPbJ1SVEfRXSi02lFntXKHri/VH8fevtrJPmaL2Y7uKy8Nf8RtBbQwWi4W36qa74D9QJPvnacPx58ETXH/1551Z+GFbJv78+yi+f+mubPRMieaiLP4AVeOi9/DT9kys2JuLiX1T8PT4bpz21Lt1NLeaEfybgBbUWbNmVVm2gn9C03yq40nbk5d24dVAZJ0Wllfi1g/XcVCrb+tY9EuLQfcWURjZKZFdCEpD4zpWWAGDTsO+4M/WHcFji7aAAhZ9WsXgYfcCByI61Mib4P8oGpSaMWPGWQVv7dq16Nu3b9X/NOWfOnVqvYJSdVmoKSkpEpQKECix/ded2bxaiPJaqdfVtpljYdRrMev7ndBrNejeIhrdW0axX9LXPtgNh/OxYk8Othwt5JVhuSVWTB3VHlNHdeDqW6v252FUl8SgDLb5M34TlLrnnnswadKkM+6Tmtr4xG+TycSbEJiQME3q34o3z5SaxJTILKxgAXtj6X7+n27/6NYB6Jcai282H8O69BNcdpB8mIkRJrRNCEdKbCjKrXYcyS8DuWgr7U7Y7JQ8j6qiIx+sSsehvDLOCaUts8i1WolcD99tOY4vNhxFj5bRuKZ/K77s1SqaH5caH8abENgoKqjx8fG8CYI3qO5LfXVSL552ZxVZsC2jEBkF5bx4gMguquBoOglivruD6y3np+HJy7pgZ2YRrnxzVY3njQk1YOOTY/j6og1HOZjULNKEVrGh6JcWi+hQV/GXh8Z2xOOXdA7YbAQhgHyohw8fxokTJ/jSbrdzPirRrl07hIdLsrNwKiRsVDCkdtEQKobtKYhNKUvklzW4swY6NIvAl3cN5spYtLiA3AZmw8lFBt/cM+S0h5rW1wvBjd8k9t9000344IMPTrl96dKlGDFiRL2eQxL7BUFoKA3RDb8RVG8ggioIQkORAtOCIAgKIBUoBEEQvIQIqiAIgpcQQRUEQfASIqiCIAheQgRVEATBS4igCoIgeAkRVEEQBC8hgioIguAlRFAFQRC8hAiqIAiClxBBFQRB8BIiqIIgCF5CBFUQBMFLiKAKgiB4CRFUQRAELyGCKgiC4CVEUAVBELyECKogCIKXEEEVBEHwEiKogiAIXkIEVRAEwUuIoAqCIHgJEVRBEAQvIYIqCILgJURQBUEQvIQIqiAIgpcQQRUEQfASIqiCIAheQgRVEAQhmAQ1PT0dU6ZMQVpaGkJCQtC2bVs89dRTsFqtSg9NEAShCj38gF27dsHhcOCtt95Cu3btsG3bNtx2220oLS3Fiy++qPTwBEEQGI3T6XTCD3nhhRcwZ84cHDhwoN6PKSoqQlRUFAoLCxEZGenT8QmCEBg0RDf8wkKtC3pzsbGxZ9zHYrHwVv3ACIIgBLUPtTb79+/Ha6+9hjvuuOOM+82aNYvPLJ4tJSWlycYoCELwoaigzpgxAxqN5ozbunXrajzm2LFjuOiiizBhwgTceuutZ3z+6dOnsyXr2Y4cOeLjdyQIQjCjqA81NzeXtzORmpoKs9lcJaYjR47EgAEDMG/ePGi1DTsfiA9VEISA9aHGx8fzVh8yMjJYTPv06YP333+/wWIqCILga/wiKEWW6YgRI9CqVStOk8rJyam6LykpSdGxCYIg+JWg/vTTT9i3bx9vLVu2rHGfn2Z9CYIQgPjFvPmmm25i4axrEwRBUAt+IaiCIAj+gAiqIAiClxBBFQRB8BIiqIIgCF5CBFUQBMFLiKAKgiB4CRFUQRCEYErs9xaevFUp4ycIQn3x6EV98t6DSlCLi4v5Usr4CYLQGP2gIikBWbG/MVAbFaoLEBERwaUB1X5WJOGnkoP+1F1Axi3HO9C+JySRJKbNmzc/a1GmoLJQ6WDUrgWgdujLpvYvXF3IuOV4B9L35GyWqQcJSgmCIHgJEVRBEAQvIYKqUkwmE5566im+9Cdk3HK8A/l7cjaCKiglCILgS8RCFQRB8BIiqIIgCF5CBFUQBMFLiKAqRH5+PiZPnsz5bbTR9YKCgrO2gqEFCdW3gQMH1tjHYrHg3nvv5W6yYWFhuPzyy3H06FHFxl1ZWYlHH30U3bt35/FQcvQNN9zACyyqQ00Ya7+3SZMmNXqcb775JtLS0rgFOXXKXbly5Rn3X758Oe9H+7dp0wb/+c9/Ttln0aJF6NKlCwdS6PKrr75q9Pi8Me4vv/wSo0ePRkJCAudyDho0CD/++GONfajdeu3jSltFRYVi4162bFmdY9q1a1eTH2+vQ0Epoem56KKLnN26dXOuWrWKN7p+6aWXnvExN954Iz/u+PHjVVteXl6Nfe644w5nixYtnEuWLHFu2LDBOXLkSGfPnj2dNptNkXEXFBQ4R40a5fz000+du3btcv7xxx/OAQMGOPv06VNjv+HDhztvu+22Gu+NHtsYPvnkE6fBYHC+8847zh07djjvv/9+Z1hYmPPQoUN17n/gwAFnaGgo70f70+Po8V988UXVPvRedTqd85lnnnHu3LmTL/V6vXP16tWNGqM3xk33P/fcc841a9Y49+zZ45w+fTo/nj53D++//74zMjKyxnGlzZs0dNxLly6lQLhz9+7dNcZU/TvaFMfbF4igKgB96egLVf3LQUJDt5HonElQr7jiitPeTwJEX2z6gnvIyMhwarVa5w8//KDYuGtDAkCPqf6DI0GlH6I36N+/P59YqtOpUyfnY489Vuf+jzzyCN9fndtvv905cODAqv+vvvpqPplUZ+zYsc5JkyZ5ZcyNGXdddOnSxTlz5swaghoVFeX0JQ0d91K3oObn55/2OZviePsCmfIrwB9//MHT5QEDBlTdRlN3um3VqlVnfCxNlxITE9GhQwfcdtttyM7Orrpv/fr1PMUeM2ZM1W00xe7WrdtZn9fX465OYWEhT/Gio6Nr3P7RRx+xq6Jr16546KGHqorZNASr1crHofoxIOj/042R3lft/ceOHYt169bx8TzTPt44ro0dd121KuiYxcbG1ri9pKQErVu35mXXl156KTZu3OiVMZ/ruHv16oXk5GRceOGFWLp0aY37fH28fUVQreVXC5mZmSyKtaHb6L7TMW7cOEyYMIF/HAcPHsQTTzyBCy64gL/Q5GeixxqNRsTExNR4XLNmzc74vL4ed3XId/fYY4/h2muvrbGG+7rrrmMfXFJSErZt24bp06dj8+bNWLJkSYPGmJubC7vdzu+5vseAbq9rf5vNxs9HP/rT7eON49rYcdfmpZdeQmlpKa6++uqq2zp16sR+VPJhU0GSV199Feeffz4f2/bt2ysy7uTkZLz99tvsayWf//z581lUyVgYNmwY7+Pr4+0rRFC9yIwZMzBz5swz7rN27Vq+rKvaFblgzlQFa+LEiVXXyers27cvi+vixYtx5ZVXnvZxZ3teX4/bA1l7FGgiS4qCGNUha7v6e6MfO72/DRs2oHfv3mgotcdztjHWtX/t2xv6nI2hsa/x8ccf8+f43//+t8ZJj2YQ1QOXJKZ0PF977TXMnj1bkXF37NiRNw8UTKOqUy+++GKVoDb0OdWCCKoXueeee84amU5NTcWWLVuQlZV1yn05OTmnnJXPBJ3pSVD37t3L/5N1R1MwisRXt1LJLTB48GBFx01iSpYTWda//vrrWSsM0Y/eYDDwe2uIoJLLQKfTnWLJ0DE43RjpuNW1v16vR1xc3Bn3acjn5e1xe/j0008xZcoUfP755xg1atRZK67169ev6juj5LirQ6K/YMGCqv99fbx9htJO3GDEE9z5888/q26jQE9Dgzu5ublOk8nk/OCDD2oEpSii7uHYsWNeD0o1dNxWq9U5fvx4Z9euXZ3Z2dn1eq2tW7fy8y5fvrxRQZI777yzxm2dO3c+Y1CK7q8OBVlqB6XGjRtXYx8Kmng7KNWQcRMLFy50ms1m51dffVWv13A4HM6+ffs6b775ZqeS467NVVddxRkpTXm8fYEIqkLQl6NHjx4cJaete/fup6QfdezY0fnll1/y9eLiYueDDz7I6SQHDx7kSOmgQYM4RaqoqKiGELRs2dL5888/c/rMBRdc4PW0qYaMu7Ky0nn55ZfzmDZt2lQjTcZisfA++/bt48j02rVr+b0tXryYo8S9evVq1Lg9aTxz587lk8DUqVM5jSc9PZ3vpx/65MmTT0mbeuCBB3h/elzttKnff/+d03ieffZZTuOhS1+lTdV33CSmNIY33njjtOlmM2bM4JPp/v37nRs3bmQhpcdUPyk29bj//e9/8wmAUr22bdvG99PJc9GiRU16vH2BCKpCUP7odddd54yIiOCNrtdOI6EvGaW9EGVlZc4xY8Y4ExIS+MvbqlUrTqM6fPhwjceUl5c777nnHmdsbKwzJCSExa72Pk05bhJI+r+ujU4KBI1v2LBhPGaj0ehs27at87777jslx7YhkMi0bt2an6937941LF06bpSmVZ1ly5axgNP+qampzjlz5pzynJ9//jmfLOj4k+BXFwBv0ZBx0/W6jivt54HEjb4r9Hz03aHvEJ2UlRz3c889x58xWdYxMTHOIUOG8ElUiePtbaTalCAIgpeQPFRBEAQvIYIqCILgJURQBUEQvIQIqiAIgpcQQRUEQfASIqiCIAgiqIIgCOpCLFRBEAQvIYIqCILgJURQBUEQvIQIqhDU/Pbbb1wmkAode6ASg1R389ChQ4qOTfA/RFCFoGbTpk3o3Lkzdzyofhu1Z6Fas4LQEERQhaCGWoFQb6PqkKD27NlTsTEJ/osIqhDUkHied955NW6jJnYiqEJjEEEVghZqLrd9+/ZTLFTqY1VbZAWhPoigCkHL7t27UV5ezq22q7cvzsjIEAtVaBQiqEJQT/cJ6gBKTeu+//573HDDDXxb9ai/INQXEVQhqAV19OjRnCZFrav//ve/49lnn+WOrG+88YbSwxP8EGmBIgQtY8eO5RbVs2bNUnooQoAgFqoQ1ClTPXr0UHoYQgAhgioEJZmZmcjKyhJBFbyKTPkFQRC8hFiogiAIXkIEVRAEwUuIoAqCIHgJEVRBEAQvIYIqCILgJURQBUEQvIQIqiAIgpcQQRUEQRBBFQRBUBdioQqCIHgJEVRBEAR4h/8HcajgJIlYgwsAAAAASUVORK5CYII=", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Wrap the analytic Fisher in a bare FisherMatrix for plotting\n", "fm_analytic = FisherMatrix(\n", " simulator=None,\n", " theta_fid=theta_fid,\n", " param_names=param_names,\n", " verbose=False,\n", ")\n", "fm_analytic.F = F_analytic\n", "fm_analytic.F_inv = F_inv_analytic\n", "\n", "fig = fm.plot_ellipses(\n", " labels='Numerical',\n", " others={'Analytic': fm_analytic},\n", " label_names={'mu': r'\\mu', 'sigma': r'\\sigma'},\n", " n_sigma_list=[1, 2],\n", " filled=True,\n", ")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 8. Using the NoisyLine toy model\n", "\n", "The same workflow applies to any simulator. Here we use `psiphy.toy_models.NoisyLine`\n", "to forecast constraints on the slope and intercept of a noisy linear model.\n", "This example shows **correlated** parameter constraints." ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", "==================================================\n", " Estimating covariance matrix\n", "==================================================\n" ] }, { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "da7be5fb88e74e71bae8b99ee909e3b1", "version_major": 2, "version_minor": 0 }, "text/plain": [ "Simulating realisations: 0%| | 0/20 [00:00" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "Correlation (slope, intercept): -0.861\n" ] } ], "source": [ "fig = fm_line.plot_ellipses(\n", " labels='NoisyLine',\n", " label_names={'slope': r'\\alpha', 'intercept': r'\\beta'},\n", " filled=True,\n", ")\n", "\n", "print(f\"Correlation (slope, intercept): {fm_line.correlation_matrix()[0,1]:.3f}\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Summary\n", "\n", "| Step | Method |\n", "|------|--------|\n", "| Define simulator | `f(theta, seed, noise_seed) -> array` |\n", "| Create object | `FisherMatrix(simulator, theta_fid, param_names)` |\n", "| Covariance | `fm.estimate_covariance(seeds, noise_seeds)` |\n", "| Derivatives | `fm.estimate_derivatives(delta_frac, seeds)` |\n", "| Compute | `fm.compute(method='standard')` |\n", "| Inspect | `fm.sigma_1d()`, `fm.plot_ellipses()` |\n", "| Save / load | `fm.save('results.npz')` / `FisherMatrix.load('results.npz')` |\n", "\n", "Advanced options:\n", "- `method='CoultonWandelt2023'` — bias-corrected estimator for large data vectors\n", "- `method='lfim'` — score-trick estimator\n", "- `use_crn=True` in `estimate_derivatives` — Common Random Numbers for low-noise derivatives\n", "- `fm.convergence_test_delta_frac()` — choose the optimal step size" ] } ], "metadata": { "kernelspec": { "display_name": "beorn", "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.11.14" } }, "nbformat": 4, "nbformat_minor": 5 }