From 58ec08dea127f4b7c03dbc4bc62afe9a7ecfa87e Mon Sep 17 00:00:00 2001 From: Jonathan Brodrick Date: Mon, 8 Jun 2026 15:53:08 +0200 Subject: [PATCH] Add `project` helper to compute projection operator using Golub-Pereyra gradients --- docs/api/functions.md | 10 + docs/api/linear_solve.md | 4 +- docs/examples/variable_projection.ipynb | 262 ++++++++++++++ lineax/__init__.py | 4 + lineax/_misc.py | 11 + lineax/_project.py | 353 ++++++++++++++++++ lineax/_solve.py | 124 ++++++- lineax/_solver/normal.py | 13 +- lineax/_solver/qr.py | 17 +- lineax/_solver/svd.py | 33 +- lineax/_tags.py | 47 +++ mkdocs.yml | 1 + pyproject.toml | 2 + tests/test_project.py | 463 ++++++++++++++++++++++++ 14 files changed, 1311 insertions(+), 33 deletions(-) create mode 100644 docs/examples/variable_projection.ipynb create mode 100644 lineax/_project.py create mode 100644 tests/test_project.py diff --git a/docs/api/functions.md b/docs/api/functions.md index 82c8f489..03b25de5 100644 --- a/docs/api/functions.md +++ b/docs/api/functions.md @@ -2,6 +2,16 @@ We define a number of functions on [linear operators](./operators.md). +## Transformations + +These functions transform an operator to a new one (e.g. representing its inverse or column-space projection). + +::: lineax.invert + +--- + +::: lineax.project + ## Computational changes These do not change the mathematical meaning of the operator; they simply change how it is stored computationally. (E.g. to materialise the whole operator.) diff --git a/docs/api/linear_solve.md b/docs/api/linear_solve.md index 0eb4af57..65b9502a 100644 --- a/docs/api/linear_solve.md +++ b/docs/api/linear_solve.md @@ -6,6 +6,4 @@ This is the main entry point. ## invert -A convenience function for obtaining the inverse of an operator as a [`lineax.FunctionLinearOperator`][]. - -::: lineax.invert \ No newline at end of file +A convenience function for obtaining the inverse of an operator as a [`lineax.FunctionLinearOperator`][]; see [`lineax.invert`][]. \ No newline at end of file diff --git a/docs/examples/variable_projection.ipynb b/docs/examples/variable_projection.ipynb new file mode 100644 index 00000000..fa4c598c --- /dev/null +++ b/docs/examples/variable_projection.ipynb @@ -0,0 +1,262 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "a12cb0ab", + "metadata": {}, + "source": [ + "# Nonlinear least squares with variable projection\n", + "\n", + "Many curve-fitting problems are *separable*: the model is linear in some parameters and nonlinear in others. A classic example is a sum of exponentials,\n", + "\n", + "$$y(t) \\approx \\sum_{j=1}^n c_j\\, e^{-k_j t},$$\n", + "\n", + "which is linear in the amplitudes $c$ but nonlinear in the decay rates $k$. Collecting the basis functions into a matrix $\\Phi(k)_{ij} = e^{-k_j t_i}$, the model is $\\Phi(k)\\,c$ and the least-squares problem is\n", + "\n", + "$$\\min_{k,\\,c}\\; \\tfrac{1}{2}\\,\\lVert \\Phi(k)\\,c - y \\rVert^2.$$\n", + "\n", + "**Variable projection** ([Golub & Pereyra, 1973](https://doi.org/10.1137/0710036)) eliminates the linear parameters. For any fixed $k$ the optimal amplitudes are $c^\\star(k) = \\Phi(k)^\\dagger y$, so substituting them back gives a residual in the nonlinear parameters alone:\n", + "\n", + "$$r(k) = y - \\Phi(k)\\,\\Phi(k)^\\dagger\\, y = (I - P(k))\\,y, \\qquad P(k) = \\Phi(k)\\,\\Phi(k)^\\dagger.$$\n", + "\n", + "Here $P(k)$ is the orthogonal projector onto $\\mathrm{range}\\,\\Phi(k)$, which is exactly what [`lineax.project`][] returns — together with its efficient gradient. So $r(k)$ is a one-liner that we can hand straight to a nonlinear least-squares solver." + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "id": "227191a3", + "metadata": { + "execution": { + "iopub.execute_input": "2026-06-07T17:55:50.629117Z", + "iopub.status.busy": "2026-06-07T17:55:50.629029Z", + "iopub.status.idle": "2026-06-07T17:55:52.721614Z", + "shell.execute_reply": "2026-06-07T17:55:52.721001Z" + } + }, + "outputs": [], + "source": [ + "import equinox as eqx\n", + "import jax\n", + "import jax.numpy as jnp\n", + "import jax.random as jr\n", + "import lineax as lx\n", + "import matplotlib.pyplot as plt\n", + "import optimistix as optx\n", + "\n", + "\n", + "jax.config.update(\"jax_enable_x64\", True)\n", + "\n", + "# Synthetic data: a sum of two decaying exponentials, plus a little noise.\n", + "t = jnp.linspace(0.0, 4.0, 200)\n", + "true_rates = jnp.array([0.5, 3.0]) # nonlinear parameters, k\n", + "true_amplitudes = jnp.array([2.0, 1.0]) # linear parameters, c\n", + "\n", + "\n", + "def features(rates):\n", + " # Phi(k): each column is a basis function exp(-k_j t); shape (len(t), len(rates)).\n", + " return jnp.exp(-t[:, None] * rates[None, :])\n", + "\n", + "\n", + "y = features(true_rates) @ true_amplitudes + 0.01 * jr.normal(jr.PRNGKey(0), t.shape)" + ] + }, + { + "cell_type": "markdown", + "id": "41fc0a76", + "metadata": {}, + "source": [ + "The projector is built from the operator $\\Phi(k)$. Here $\\Phi(k)$ is dense so we wrap it in a [`lineax.MatrixLinearOperator`][], but [`lineax.project`][] accepts any [`lineax.AbstractLinearOperator`][] — so for a large-scale problem $\\Phi$ could be matrix-free.\n", + "\n", + "We solve with [`lineax.SVD`][]: during the fit the optimiser may try rates that nearly coincide, making $\\Phi$ rank-deficient, and the SVD pseudoinverse handles that gracefully (the default full-rank solver would error). Eliminating the amplitudes then takes a single line, and we fit the rates with [Optimistix](https://github.com/patrick-kidger/optimistix)." + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "32ef0ef3", + "metadata": { + "execution": { + "iopub.execute_input": "2026-06-07T17:55:52.726772Z", + "iopub.status.busy": "2026-06-07T17:55:52.726597Z", + "iopub.status.idle": "2026-06-07T17:55:53.484156Z", + "shell.execute_reply": "2026-06-07T17:55:53.483484Z" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "true rates: [0.5 3. ]\n", + "fitted rates: [0.50192609 3.03565001]\n", + "fitted amplitudes: [2.00523518 0.99997341]\n" + ] + } + ], + "source": [ + "def varpro_residuals(rates, args):\n", + " operator = lx.MatrixLinearOperator(features(rates))\n", + " return y - lx.project(operator, lx.SVD()).mv(y) # (I - P) y, amplitudes eliminated\n", + "\n", + "\n", + "solver = optx.LevenbergMarquardt(rtol=1e-8, atol=1e-8)\n", + "init_rates = jnp.array([1.0, 2.0])\n", + "solution = optx.least_squares(varpro_residuals, solver, init_rates, max_steps=100)\n", + "\n", + "fitted_rates = solution.value\n", + "# Recover the amplitudes for the fitted rates with a single linear solve.\n", + "fitted_operator = lx.MatrixLinearOperator(features(fitted_rates))\n", + "fitted_amplitudes = lx.linear_solve(fitted_operator, y, lx.SVD()).value\n", + "print(\"true rates: \", true_rates)\n", + "print(\"fitted rates: \", fitted_rates)\n", + "print(\"fitted amplitudes:\", fitted_amplitudes)" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "0399dddc", + "metadata": { + "execution": { + "iopub.execute_input": "2026-06-07T17:55:53.485596Z", + "iopub.status.busy": "2026-06-07T17:55:53.485444Z", + "iopub.status.idle": "2026-06-07T17:55:53.604845Z", + "shell.execute_reply": "2026-06-07T17:55:53.604336Z" + } + }, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAhgAAAFzCAYAAAB8X3AUAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjksIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvJkbTWQAAAAlwSFlzAAAPYQAAD2EBqD+naQAAS9pJREFUeJzt3QdYFHf+P/D37FLsKCoCgliwgAW7orG3xJIYU4yJJYkxTXMxpnr3T0zMXUwuzRQvppyJJvFnjRpLVOyJ2CsqYkNQAbuiiJTd+T+fr+4eICDowO7C+/U882yb3Z2ZZZn3fqum67oOIiIiIgOZjHwxIiIiIgYMIiIiKhIswSAiIiLDMWAQERGR4RgwiIiIyHAMGERERGQ4BgwiIiIyHAMGERERGc4NpYzVakVCQgIqVqwITdMcvTlEREQuQ8bmvHLlCvz9/WEy5V9GUeoChoSLwMBAR28GERGRyzpx4gQCAgLyXafUBQwpubAdnEqVKjl6c4iIiFxGcnKy+pFuO5fmp9QFDFu1iIQLBgwiIqLCK0gTAzbyJCIiIsMxYBAREZHhGDCIiIjIcKWuDQYRUVYWiwUZGRk8KEQ3ubu7w2w2424xYBBRqXX16lWcPHlS9e0nov814JQuqBUqVMDdYMAgolJbciHholy5cqhevToH3iPCjYG0zp49q74b9evXv6uSDAYMIiqVpFpE/plKuChbtqyjN4fIach34vjx4+o7cjcBw6GNPL/55hs0a9bMPiZFeHg4/vjjj3yfM3fuXDRq1AhlypRB06ZNsWzZsmLbXiIqeThlAFHRfCccGjCkjufDDz/Ejh07sH37dnTv3h0PPPAA9u/fn+v6kZGRGDJkCEaOHIldu3Zh4MCBatm3b1+xbzsRERHlTdOdrHWTt7c3Pv74YxUicho8eDBSUlKwZMkS+33t27dH8+bNMXXq1AIPc+rl5YXLly8bNpLnhu27cSlqBao06YVObVoa8ppEVLSuX7+O2NhY1KlTR5WIEtHtvxuFOYeanKnB1axZs1SAkKqS3GzatAk9e/bMdl+fPn3U/XlJS0tTByTrYqSIA6ehLRqN++M+wJqF09RtIiJn9u6776ofZsWpdu3amDx5crG+JzmWwwNGVFSU6grj6emJ559/HgsWLEBoaGiu6yYlJaFGjRrZ7pPbcn9eJk2apNKWbTF6JtVNR88jUm+qrncy7cPmY+cNfX0iIqO99tprWL16tVMf2J9++gmVK1d29GaQKweMhg0bYvfu3diyZQteeOEFjBgxAgcOHDDs9cePH6+KcmyLzKJqpPB6VbHeciNgtDMdQHjQ7WeYIyJyJPlRV7VqVX4IVLIDhoeHB4KDg9GqVStV2hAWFoYvvvgi13V9fX1x+nT2Kgi5LffnRUpGbL1UimIG1V6hNfDK0EFIcauC8loaelaMN/T1icj5SdXoxMUHiqWKtGvXrvjb3/6GN954Q7VZk/9/UuWRVXx8vGowL0FC/uc9+uij2f535qwiWbduHdq2bYvy5curUoOOHTsiLi5OdVU0mUyqEX5WUtURFBQEq9Wa6zaeOXMGAwYMUN1/pR7/119/vWWdzz77TPUElPeUkuUXX3xRDXxm256nnnpK/SiUHg2y2Pbx559/RuvWrdV04bLvjz/+uHo/cj4ODxg5yR+stJvIjbTNyFmsFxERkWebjeLSq7EfyjfqcePGsbUO3RYiKl4SKkbN2I7pkcfVZXGEjOnTp6sTs5T8/vvf/8bEiRPV/0Lb/1AJFxcuXMD69evV/ceOHVON5HOTmZmpeuN16dIFe/fuVW3ann32WXVSl3YT0u7txx9/zPYcuf3kk0+q8JEbeUxKi9euXYt58+bhP//5zy0hQJ775Zdfql6Dsj9r1qxRoUl06NBBhRgJR4mJiWqRah0hYzO8//772LNnDxYuXKhCkLwfOSHdgd566y19/fr1emxsrL537151W9M0feXKlerxYcOGqftsNm7cqLu5uemffPKJHh0drU+YMEF3d3fXo6KiCvyely9fll4z6tJQO3/W9QmVdP277sa+LhEVidTUVP3AgQPq8m689/t+ve5bS/WgN5eoy4mL9+tFqUuXLvo999yT7b42bdrob775prou/z/NZrMeHx9vf3z//v3q/97WrVvVbfnfGRYWpq6fP39ePbZu3bpc32/27Nl6lSpV9OvXr6vbO3bsUP+n5f92bmJiYrK9l5D/13Lf559/nud+zZ07V69atar99o8//qh7eXnd9nhs27ZNvfaVK1duuy7d/XejMOdQh5ZgSKIdPny4aofRo0cPbNu2DStWrECvXr3sxXySXG0k1c6cORPfffedqkqRZCwJtkmTJnC4ut1uXCbsBFIvOnpriKiYSDssi67DrGnqsn3dom/bIAMUZuXn52cvIYiOjlZVDlkbtEvDean6kMdykmoWKQGQHnlSrSFV1Fn/70rphozmKA3wbY0vu3Xrpko3ciPv4ebmpqq9bWRwxJwNNletWqX+79esWVNVdwwbNgznz5/HtWvX8t13GTdJtrNWrVrqeVLyYjtfkHNxaMD473//q4q3pEpEvhzyB2cLF7Z6OPljzuqRRx5BTEyMeo4MsNW3b184Ba+aQLUGgG4FYv909NYQUTGRdljfD2+NJzvWVpdyuzhmu8xKqjPyag9REFLlIVUj8iNu9uzZaNCgATZv3mxvJyc/BGWd9PR09SPv6aefvqvtl//7/fv3V0Fp/vz5KjRMmTJFPSbvkRcZxkCCkFSdSLsO+VFqCz75PY8cw+naYLi0m6UY29f8xvEwiEoRCRVv9w8tlnBxOyEhIar9Q9Yec9Iz79KlS3kOASBatGihet3JiMlSKixBwuaZZ55RPwClLYW02Rg0aFCeryOlFbKOhAYb+VEo728jj0kg+vTTT9VgiRJoEhISsr2OBBsZHymrgwcPqlIOGQG6U6dO6r3YwNN5MWAYaJdHC3VZ/WxksTX2IiLKShplSu+MJ554Ajt37sTWrVtVCYRUJUjvi5xkxEYJFlKCIT1HVq5cicOHD6ugYiPXJQi8+eabarqG/CaHkyrve++9F88995xqhCphQgJK1udIz0FprPnVV1+pBqjSMyTnaMxSBSO9SqRh/7lz51TViVSLSPCwPe/3339XDT7JOTFgGGjF1WBk6GYEaacRpJ3loFtEVOykumTRokWoUqUKOnfurAJH3bp1VdVHbmS6eikZeOihh1RJgvQgGT16tAoIWcn0DVINUZDqEalO8ff3V6FGSjvkNX18fOyPSxs66ab60UcfqdISqe6QYQqykuoaGXxRer/I7J7SW0YupdpcJr2U0hgpyfjkk0/u+FhRKZuLpKgVxVwkNlJi4TVrANqaYjA+YyS6P/GmUxSZEtGtOBdJ4UhJgZzYpSsrlWzXS9pcJCWBhImqTfuo62PrnGK4ICKXJ9UU0qD+66+/xksvveTozSEXwoBhsHrtB6hL7zOb8P7vUWyHQUQubcyYMarLqYwgere9R6h0YcAwmn9LZLhXhHv6ZezYvI6NPYnIpUmbBxkWQNpwyHgYRAXFgGE0sxuOlr/Rm6SjtlcNvsMZVomIqLRhwCgCet3u6rKLeW+xjexHRETkTBgwikBIpxuD0LQ2HcaPg4PZ2JOIiEodBoyiUCUI8AmFCRZ0c4sqkrcgIiJyZgwYRaXBje6qOLS8yN6CiIjIWTFgFJUG9964PBwBWDKL7G2IiIwgE5DJKKC7d+8u8HNkFlaZbTU/0r117NixcFbvvvsumjdvXuTvo2mamv27qMls4zKTrslkwuTJk4tt/3Lj5pB3LQ0C2gBlq9yYuv3kViCog6O3iIgoT3JSkmnaq1WrVqqO0muvvWboAGJyQpcgsTtHUJNjK8O3FyUZZVPGLZFh2GXodxlxUyaVy7p/Egpl4rniCDsswSiyI2sG6vdWVzcu+5UDbhGR05I5RmSMC19fX7i5lYzfnQWdvr1ChQqoWrXoe/r5+vrC09OzSN8jPj5eTSLXr18/+Pn5qXlmimv/csOAUYT2lmuvLqsnruWAW0RkSPG3TCImv0qzeuCBB+yjbB49elTdrlGjhjq5tGnTRk21nnOmUplbRGZZlfkkZDKynFUkMlW6THAm81HITKgyS+oXX3yR63a99957aiIyeS2ZoCy/k7sM2iWlBjVr1kT58uXRrl07rFu3Lt/9thXzf/vtt6qkRU6cjz76qJoPI2d1zb/+9S91jGR7RVRUFLp37672QU60sq8y/HnO187qhx9+UDPIyjwcMiW8TFOf1cmTJ9Wsst7e3mofZJZamTlWBiWTY7Fnzx51LGWR+3KrIrnddtn2RyZzk7Ag68gkdBIgciPvI7PoCpncTt5PPtOs+yfXp0+fribDs23f7Y793SgZUdVJ/ZHaGKG6CQ1Mp1BbO6MG3OLkZ0ROSuZ9zLjmmPd2LydnoNuu9sgjj6ji7rVr16JHjx7qvgsXLmD58uVYtmyZui0nqb59+6oTrfxinjFjBgYMGICYmBg13bmNnLjeeecdTJgwIdf3khATEBCgJjiTk1tkZKQ6CcrJTk7uNjKdupyI5UQlJ7SnnnpKrS/vnxspwj9w4ABmzZqlgsCCBQvU9O5ywq1fv36e+37kyBHMmTMHixcvVlUBEn5efPFFNRNr1m2RkBMREaFup6SkoE+fPggPD8e2bdtw5swZNXW8bIPtxJ+TvJ4cF5l7pUWLFti1axdGjRqlgsSIESPU8ZVZYiUgyXTxUjKxc+dOdbxk5leZt0U+j1U3Q51UU+RU0O2Sz1mOt1zK/svrS1iQ7clJHpPwJbPnbt26VV2X0JeVBLvo6Gh1/GTGWyEhqagwYBShlg1rY9vORgg3H0AX0y60r9u3KN+OiO6GhIsP/B1zDP+eAHiUv+1qUod/3333YebMmfaAMW/ePNVuolu3bvap0GWxkZIKOYnLyVBOYDby6/nVV1+135ZwkJW7u7v6NW4jJRmbNm1SJ/msAcPDwwPTpk1TpQqNGzfGxIkT8frrr6v3lYaGOYvw5cQmlxIubCc9OSHL/R988EG+M3xKWJITu/jqq69UVcCnn36qTvJCQoCUPsg2ie+//97+PHlMSHCQwCVTxUspT04SuOQ1ZZp5235LIJLSEwkYcuzPnj2rgoHt5BwcHGx/vpQaSTWT781tyo28RkG2Sz5vuV+qr6QkRfZXQlRuAcNWEiIkWOT2/rJtsp6UIuW3fUZhFUlRz67a8sbkZy8FHGXpBRHdtSeeeALz589XJwnbL+7HHnvMfjKXX9hy0pYi/sqVK6uTivxqlZN6VlKsfztTpkxRE53JCUteR6pocr6OhBkJFzbyq1y24cSJE7e8npRSSNVLgwYN1OvZlvXr16uqHZH1fqlusZHSF1u4sL2PlBpIyYyNVBHYwoWQ/Zbts53ERceOHW95XtaSBdkOKR3Juh3//Oc/7dsnVUhSsnE3v/yjC7hdEtiyzv8ipRlS2uEqWIJRxBp0egTY8xGqndsKpF0BPCsW9VsS0Z1WU0hJgqPeu4DkV66u61i6dKlqX/Hnn3/i888/tz8u4UKqCKQKRH5Zyy/Whx9++JZ2EVlPbrmRKgx5Lfk1LyfzihUr4uOPP1ZtDe6UBA85Ye7YseOWidPkRC6y9r6Q6o7CuN0+FWT7bCUf0jYkK9v2yvEsLu7u7tluS5uJnO1vnBkDRlGrGgx41wUuHAOOrQNCbpRoEJGTkTYQBaimcDRp7yDF91JyIfXy0pixZcuW9sc3btyoGgg++OCD9pNmzuqPgpDX6dChg2rnYGP7FZ+VNGhMTU21n3g3b96swoK0AchJfvlLCYb8Cu/UqVOu75u1uiErKTlJSEiwV63I+0ipja0xZ26kFEfaNEjJhC18yH7l9TypmpDXP3bsmCopyk2zZs1UNYy0fcmtFENKUGQf81PY7TJSQbbPKKwiKY5/WjcH3dq1aha7qxLRXZOTn5RgSNuHnCdCaSj522+/qZIAOfk//vjjd/SrV15n+/btWLFiBQ4dOoS3335btTvISUpGpEpB2ilIQ1NpwyBtPXK2vxBSNSLbK71XZBtjY2NVg8RJkyap/bldsJI2ELJPUmrzt7/9TbUFya8tgbyX7XnS+FIaS0oj2WHDhuXa/kJIuxPZni+//FLtt1TrSPsQGVtCSO8ReU/p4SGhQMKIVFlJ+xRbDx3Zr927d+PcuXP2qqy73S6jyPbt3btXVcXI9uXVK8UIDBjFYIdnG3UZcO5PPDdjK0MGEd0VaaApv57lJCEBIis5EUrjQCl9kOoU6a2QtYSjoJ577jlVUiK9E6S64Pz589lKM2yksamEkc6dO6t177//ftUdMi9yspaAIQ1M5de6nKgluGTt4ZJXyYZsj/SQ6d27typJyNl9NCdpGyIBSUobpDpJqopke6XhZF6kN4eUUMh2SpsO6TEipQ3S2NNWArBy5Ur4+PiobZF1PvzwQ3sVigxwJb1iunXrptqu/N///Z8h22UUaSAqx13a4Mj2SUgqKpoulXmliHTPkW5D0n+6sPV7d+qfi/bgpZ33wUtLwaPpE9C0w714u39osbw3EeVOWvHLL005ccivSXJeeY2OaYTx48erEpG//vrL8Ncuid+NwpxDWYJRDNrV98Uqawt1/V7TFrSv65hR1YiI6Ab5bS1tSqTbp/TWIOMxYBRTd9V6nW8UYw6ptAe9QnyK422JiCgP8gs8NDRUVXn8/e9/53EqAgwYxaR510GARwWUvZYInNpZXG9LRFQiqkiMrh6RMUKkAaZUjQQFBRn62nQDA0ZxcS9rn/wM0YuK7W2JiIgcgQGjOIXery4ubJuHiP1JxfrWRERExYkBoxityQzDdd0d3umn8Pkv89ldlcgJlLKOdETF9p1gwChGf8Vfx3rrjWlz+5q3qdlVicgxbOMW5De1OFFplH7zO5FzOPfC4lDhxSi8XlUs2dQGfczbVHfVWHZXJXIYmfFSBjySmTFlzofcRp4kKm2sVqv6Tsh3Q74jd4MBo5i7q7o9OhKWRd8h2JSA4GoXZfT74twEIsoycZTMTikDCsXFxfG4EN0kYVtGVpXvyN1gwChm3ZoHA/u7A4dXYu3CH5B5z+ucxp3IQWQMBBnmmtUkRNm/F0aU6DFgOMD+yt3QGCtR4+RK9J3RDd8Pb82QQeQg8o+UQ4UTGY+Vjg6wLL0FMnUTQk1xqGdKYmNPIiIqcRgwHKB5w7rYaG2irvfXNnJuEiIiKnEYMBzU2LNq+FB1/dkqOzg3CRERlTgMGA7SpMfjgFtZlL96HEjg3CRERFSyMGA4imdFoFHfG9f3znXYZhAREZW4gDFp0iS0adMGFStWhI+PDwYOHIiYmJh8n/PTTz+pvrlZF5dtAd5ssLq4umM2Vu075eitISIiKhkBY/369Rg9ejQ2b96MiIgIZGRkoHfv3khJScn3eZUqVUJiYqJ9cdVBclalN8Z5vSIqZF7A9JkzODcJERGVGA4dB2P58uW3lE5IScaOHTvQuXPnPJ8npRa+vr5wdZGxyUi0tscwcwQGmTdi87EHOB4GERGVCE7VBuPy5cvq0tvbO9/1rl69iqCgIAQGBuKBBx7A/v3781w3LS0NycnJ2RZnmptkQWZHdb2XaRs61Crn6E0iIiIqWQFDJlgZO3YsOnbsiCZNbowRkZuGDRti2rRpWLRoEX755Rf1vA4dOuDkyZN5tvPw8vKyLxJKnKm76gtDh+Cihz8qaNfRQ9vu6E0iIiIyhKYbNfH7XXrhhRfwxx9/4K+//kJAQECBnyftNkJCQjBkyBC8//77uZZgyGIjJRgSMqS0RNpyOIU1/wQ2fAw0uBd4fLajt4aIiChXcg6VH+sFOYc6RQnGmDFjsGTJEqxdu7ZQ4ULINMstWrTAkSNHcn3c09NTHYSsi9Np+qi6sByKwLpd0Y7eGiIiorvm0IAhhScSLhYsWIA1a9agTp06hX4Ni8WCqKgoNe2yq4o464Uoa22YYcHqeVPZm4SIiFyeQwOGdFGVdhQzZ85UY2EkJSWpJTU11b7O8OHDMX78ePvtiRMnYuXKlTh27Bh27tyJoUOHqm6qzzzzDFzVpqPnscDaRV0fbF7Pyc+IiMjlOTRgfPPNN6oep2vXrqoEwrbMnv2/dgjx8fFqrAubixcvYtSoUardRd++fVV9UGRkJEJDQ+GqpDfJb5kdkKa7oYkpFr2qJDl6k4iIiEpGI09nbKBSnCIOnIbfqtFociECaPMM0O9TR28SERGRazfypBtdVpv0G/2/uUky/ldNRERE5GoYMJxJnS5A5VpA2mUs+PUbNvYkIiKXxYDhTEwmHKn5oLrqd2wuRs3YzpBBREQuiQHDySwxdYNF19DedAD1TEnsUUJERC6JAcPJNA4JxXprmLr+kGkd2tet6uhNIiIiKjQGDCds7Fm100h1fXjZjdh85DSrSYiIyOUwYDihsO6PIc2zKipknMeJLQvZFoOIiFwOA4YzcvPAjsp91NXBpjUwaxrbYhARkUthwHBSessn1WU30274I4ltMYiIyKUwYDipju3a4ZxvJ5g0HT83jVJtM4iIiFwFA4YTq9b9JXVZO34+kJ7i6M0hIiIqMAYMZxbcC6hSG7h+GYia6+itISIiKjAGDGdmMgFtRqmrp1d9iYj9nGWViIhcAwOGk1tbrjeu6Z6okXoUP/zyC8fEICIil8CA4eT+PJGJRdZ71PURbivZXZWIiFwCA4aTC69XFT9l9lLXe5u2Ifl0LEsxiIjI6TFgODnpnvrasEGIKdMMbpoVtWPncGRPIiJyegwYLhIy9vg/pq4/Zl6NcloGq0qIiMipMWC4CO+WA3HCWh1VtSsYZFqHuPMprCohIiKnxYDhIno2qYmUVs+p68+Yl2H9wdOsKiEiIqfFgOFCGt33Iq6ZK6G26TR6als5CRoRETktBgxX4lEeiQ2HqavPuy2BRbdyEjQiInJKDBgupl6/cbCYPRFmOop591o5CRoRETklBgxXU74azC1vlGK0PjnD0VtDRESUKwYMVxQ+GtBMwJEIIGmfo7eGiIjoFgwYrsi7LhD6gLq6Z8777K5KREROhwHDRW3xG6ouQ89H4N0ZfzBkEBGRU2HAcFErLvpjo7UJ3DULRrv/zpE9iYjIqTBguPAkaF9kPKiuP2xah641rjt6k4iIiOwYMFx4fpJRw4YhtmIreGgWlN3yBatJiIjIaTBguHjIONd6nLre7OxiTJixnCGDiIicAgOGi/sjuS4irY1VKQbbYhARkbNgwCgBbTEmZwxS1x8xrUXXGmmO3iQiIiIGjJLWFqPTaY7uSUREjscSjBISMuo8NFFdt+yYgT+37XL0JhERUSnHgFFCRFyrj0hLKMx6JhIWvcvGnkRE5FAMGCXEpqPn8ZllsLr+sHk9ju7f5uhNIiKiUowBowQ19txurY8VljYwazp6JHyLiYsPsCSDiIgcggGjBLXD+H54axxtNg5WmFH/4gbs27Qco2ZsZ8ggIqJix4BRwkLGi4/0xa5q/dXtN91mwqyB85QQEVHpChiTJk1CmzZtULFiRfj4+GDgwIGIiYm57fPmzp2LRo0aoUyZMmjatCmWLVtWLNvrKlLCX0Oq7oFWpsPooW1D+7pVHb1JRERUyjg0YKxfvx6jR4/G5s2bERERgYyMDPTu3RspKSl5PicyMhJDhgzByJEjsWvXLhVKZNm3b1+xbrsz69yqGRJDR6rrk7wWYPOR06wmISKiYqXpuq7DSZw9e1aVZEjw6Ny5c67rDB48WAWQJUuW2O9r3749mjdvjqlTp972PZKTk+Hl5YXLly+jUqVKKLGuX0b6Z83gkX4Jf894BjMt3VUbDalGISIiuhOFOYc6VRsM2WDh7e2d5zqbNm1Cz549s93Xp08fdX9u0tLS1AHJupQKZbyw1meEuvqK2xx44RomrzrEkgwiIioWThMwrFYrxo4di44dO6JJkyZ5rpeUlIQaNbL/Cpfbcn9e7TwkbdmWwMBAlBbmdqNw1OqH6loyRrstQHRiMnuVEBFR6QoY0hZD2lHMmjXL0NcdP368KhmxLSdOnEBp0bNpIC53uTGE+FPm5aiNBJg1jb1KiIiodASMMWPGqDYVa9euRUBAQL7r+vr64vTp09nuk9tyf248PT1VPVHWpTRp2eNRnPXrCnfNggnuP8OiW9mrhIiISnbAkPalEi4WLFiANWvWoE6dOrd9Tnh4OFavXp3tPumBIvdT7qo//BmsJnd0Me3Bgp5X2dCTiIhKdsCQapFffvkFM2fOVGNhSDsKWVJTU+3rDB8+XFVz2Lz88stYvnw5Pv30Uxw8eBDvvvsutm/froIK5aFqPZjCX1RXG+7+AP9atJuNPYmIqOQGjG+++Ua1i+jatSv8/Pzsy+zZs+3rxMfHIzEx0X67Q4cOKpB89913CAsLw7x587Bw4cJ8G4aSDI7xOtLKVEe5q3Fw2zaVjT2JiKj0jINRHErNOBi5WPjTJxh4/H2k6J7ok/4p+nRshbf7hzp6s4iIyEW47DgYVLTKt3kC26wNUF5LwwS3H9G+Tt7jjRAREd0NBoxSpFdjP2Tc9zksmht6mXfgwo7f2BaDiIiKBANGKdMh/B7ENRqlrnc++jFembGBIYOIiAzHgFEKzSrzKI5bfeGnXcDrbnM58BYRERmOAaMUalO/Jv6R+ZS6Psy8EqaEHSzFICIiQzFglEIyo+qTQ5/CpvI9YdJ0PHjyY7wwYzNDBhERGYYBoxSHjI31xuGiXgGhpji86LaYVSVERGQYBoxSLKxRMN7NGK6ujzH/BkvCXpZiEBGRIRgwSnkpRv8nXsbOch3hoVnwyMkP8CKrSoiIyAAMGKVcr8a+WBv8d1zQK6CxKQ4vuS1kVQkREd01BgxCs0b18XbG0+pIvGheiF5VknhUiIjorjBgkKoqGTh0DPZX6QE3zYr2e/4BZKbxyBAR0R1jwCB7yGj8zPdI86wKnDmA43PH88gQEdEdY8Agu4i4TIy+8qS6Xjvmv9ix5jceHSIiuiMMGGS36eh5rNVbY0ZmL3W73sZX8clvG9l1lYiICo0Bg+zC61WFRdfxoWUoYqwBqGy5gJa7/oFRM7bhmenbGDSIiKjAGDAoWzuM74e3xpCODTAz8B2k6e7obt6FEeaVWB19BqNmbGfIICKiAmHAoFtCxtv9Q3FPx674V+bj6r6/u81EQy0eZk3jGBlERFQgDBiUZ9Do9PjfsadsO3hqGZji/iXK6NfQvm5VHjEiIrotBgzKd5TPsDEzcb2sL+qZErA6eB56hfjwiBER0W0xYFD+yldDmcd/Bkxu8D35B7BlKo8YERHdFgMG3V5gW6D3v25cX/n/gPgtPGpERJQvBgwqmHbPISnwPsCaiev/Nwy4epZHjoiI8sSAQQUSEX0GPQ4/jCNWf5RJPY3orx/Bqn2nePSIiChXDBhU4FE+r2vl8ELGWKTongi5vgvxs8ZxAC4iIsoVAwYVapTPI3oAxmW8oO572m05qh6azQG4iIjoFgwYVKhRPnuE+GCFtS0+y3hY3f++2zS0NR3iAFxERJSNW/abRPmHDFkiDpzGnK3VsPRoPPqZt2KK+2c4WKMrDx0REd15CcaIESOwYcOGwj6NSlppxpPtUOaR75BUtj6qa8kIXvMsVu+NdfSmERGRqwaMy5cvo2fPnqhfvz4++OADnDrFngSlVY9mdXC4+3c4p1eCX+ohWOc+jYh9CY7eLCIicsWAsXDhQhUqXnjhBcyePRu1a9fGfffdh3nz5iEjI6NotpKc1tqkMng+41Vc193Ry7wT1xe/hoj9SY7eLCIicsVGntWrV8e4ceOwZ88ebNmyBcHBwRg2bBj8/f3xyiuv4PDhw8ZvKTlt75Lt1voYlzEaVl3DgLSl2DJzIruvEhGVcnfViyQxMRERERFqMZvN6Nu3L6KiohAaGorPP//cuK0kp+9dEufbE5MsN6Z3/3/uv8IjZjG7rxIRlWKFDhhSDTJ//nz0798fQUFBmDt3LsaOHYuEhARMnz4dq1atwpw5czBx4sSi2WJyypAxtmcDfJ/ZFzMye6n7Pnf/D9qaYth9lYiolCp0N1U/Pz9YrVYMGTIEW7duRfPmzW9Zp1u3bqhcubJR20guU5LRRnVf9Tt6XrXH+MH9Y+z3buXoTSMiIgfQdF3XC/OEn3/+GY888gjKlCkDV5ScnAwvLy/VG6ZSpUqO3pwSafXe46i/chhqXd0LVKgBPL0c8K7r6M0iIqJiPIcWOmC4OgaMYpJ6CfipH3B6H1ClNvD0CqCib3G9OxEROfgcyqHCqWiUrQwMnY9r5QOBi8dx5Yf7gdSLPNpERKUEAwYVmYgTGvpcGIczemVUvByDo5PvxZo9R3nEiYhKAYcGDBlyfMCAAWr8DE3T1CBe+Vm3bp1aL+eSlMSBnZx1ivcE+GJ4+lu4pJdHvbSDqDB/CNbsOeboTSMiopIcMFJSUhAWFoYpU6YU6nkxMTFqDA7b4uPjU2TbSHc/xXuMXgtD08cjWS+nuq5WXTxCNQQlIqKSy6GzqcoQ47IUlgQKdoN1nUG4Zm+Lx6poYET6W5jhMQlhmXuxYe5QrNZmokfTWo7eTCIiKgIu2QZDxt6Q8Th69eqFjRs3Onpz6DYh44cRbVTQSPdriZEZryNF90RncxQq/f4kPvtjLyYuPqCmgCciopLDpQKGhIqpU6eqkURlCQwMRNeuXbFz5848n5OWlqa61WRdyHGjfW61NsIzGa8jVfdAm4wdaBP5AmZFHuSw4kREJYxLBYyGDRviueeeQ6tWrdChQwdMmzZNXeY378mkSZNUn13bIqGEHFtlkuzbHk9nvKFKMjqZ92Ga+79RUUvjsOJERCWISwWM3LRt2xZHjhzJ8/Hx48erAUFsy4kTJ4p1+yj3koxN1lA8mf4Wruhl0d4UjR/dJ+GeAA8eLiKiEsLlA8bu3btV1UlePD091WhjWRdyjpKMZh3vxcKm/0GquSJamw6hwcqhWLfrID8eIqISwKG9SK5evZqt9CE2NlYFBm9vb9SqVUuVPpw6dQozZsxQj0+ePBl16tRB48aNcf36dfzwww9Ys2YNVq5c6cC9oDsNGbIAodi8sYoKFzWvRePagkF4fefH6B3e+ubjRETkihxagrF9+3a0aNFCLWLcuHHq+jvvvKNuyxgX8fHx9vXT09Px6quvomnTpujSpQv27Nmjpofv0aOHw/aB7t7KCzUwJGMCEnRv1DedwivxY/Dhz4vYs4SIyIVxsjNyOOmiOmrGdtTEOUz3+BDBpgRc1Cvgt5DJGPnYI47ePCIiuomTnZFLtskICQnFI+nvYJc1GFW0qxgaMwY7V8919OYREdEdYAkGOV1pxsItMXj02D/QxbwXGboZK+pPwE6vXmrocbbLICJyHJZgkMuSAFGjWjU8l/k6Flk6wF2zoP+Rd+Cx5SuMmrENz0zfxrYZREQuwOW7qVLJIyUV13UzXs0cjWmZ96r73nKbiUluP2B9dAJH/SQicgEMGOS0bTJGdKyL8/e8h/cyhsGiaxjithY/un+EKloKR/0kInJybINBLtEuI3r9HIxMeh/ltTQctfohsf/PuKdtG0dvGhFRqZKcnKym3ZCRsW83cCUDBrmMTZHrELp2FLwyzuCa2QsLG/0bR8o0Y+NPIqJiwkaeVCKFd+iK3ff+hj3WuihnuYyH9o3G5c3T2SaDiMgJsQ0GuZT1CWY8nvEO/rC0gaeWiU/dp2KC2wxsPXLa0ZtGRERZMGCQy/UwSdE9MCbjZXyZOVDd95Tbcjx8YDSmLI7ExMUH2I2ViMgJsA0GuWSjz83HzqOMuxllj/6BEac/REUtFYm6N17IGIvd1mDVC4WDchEROa4NhkNnUyW6u5lYgYmLLRh0sjymun2KeqZEzHafiAmZT2Hyqht/+AwZRESOwSoScvkqk8NWfzyY/j5WWFqrdhkfun+PJ85+hjEzIlldQkTkIKwioZJTZeKmocaeKRh67ReYNB37rLUxucrfMbhPV5ZkEBEZgONgGHRwyDXDxi+//Befuf8HVbUruKKXxd8zRuL+oX9jyCAiukscB4NKLWlzMXToSLxS5StstTZSjT+/8vga3mteBzJSHb15RESlBttgUIkMGcN6d8CQ9H/g68wHYdU1tDq3CKc/64iNmyMdvXlERKUC22BQiW+b0TBlO7od+Aeqa8lI0T0R2/Y9JNZ+EJuOXeAw40REhcAqEqKbJRlv9w/FwXKtMSD9Q0RaQtVkaU22vYX0WcOwKDKKw4wTERURVpFQqejKmqRXxvCMv+PfGY8iQzejn3krlnm8gS6mKE79TkRUBBgwqFSUZMjIng39K2OqdSAeTH9PTfleQ7uE6R6TEH7o31i9N85ercLhxomI7h7bYFCpIeFh1IztMGsa3PXr+LTyfPS7vkQ9dshaE//xfhMLk6qpxy26zuHGiYhyYBsMonxKMp7sWBtfDe+IHY3/gZEZb+Cs7oUGplP498VX8Dfzb9D0DBUypIEoERHdGVaRUKls+CmX0jZjtaU57k37CMstbeChWTDOfR4WebyNhohF+7pVHb25REQuiwGDUNpLNFqEBOP5jLF4OeMlXNAroLEpDovLvA2PPz/EqqgTjt5MIiKXxDYYRFnGzOjsr6PRjndR49RKdVyirbVwuc+XaN+xG48TEZV6yYWYboMBgyiHib/vx9kts/Cu249qPhMLTIhvNBJ1Bk0EPMrxeBFRqZVciIDBKhKiHMKDq2GxpT36pH2MJZZ2MMOKOge/x7Uv2gBHVvN4EREVAEswiPKoMpm86hCiE5PRXduBie4/wl+7oB5LrDUAv1Z+HmGN6nOGViIqVZJZRWLMwaHSzTZuhhTzlUUqXnObixHmFTBpOi7p5fFB5uO4UP9RDG4bxKBBRKVCMgOGMQeHKGtJhlUHwrSj+MD9B9XTRGyzNsQ7GU+iZqM2GNymFoMGEZVoyQwYxhwcopwjgMoIn2ZYMNL8B8a6zUc5LQ0WXcMvll74NPNhtA2pi4a+lZCabuFMrURU4jBgGHRwiHJ2Y7UNvjV7WzwORB/AP9x/RT/zFnXfeb0iPsp8DHMtXaBpJlXiIeNsyHgbREQlAQOGQQeH6HahQ4LGtZg1eM9tOuqbTqn7d1vrqWqTfXo9hPhXwtieDRgyiKhEYMAw6OAQFTRozNt6DAGHf8ZYt99QUUuFVdcwy9IVn1oG47xeiSUZRFQiMGAYdHCIChs09h2MQe+E/6DxuT/UfZf1cvjK8hC0tqPwj/vDeECJyKUxYBh0cIju1Lb1S1F+9XiE3uxtcq1CLRxu9joWXW+lBvISm46eR1kPMxuEEpHLYMAw6OAQ3Y2IfQm4vm0GuiV8hwoZN6Z+325tiH9mPIHderAaX8Mqw+lqYINQInIJHCqcyAn0auKPMu2eQtsrH+OLzEFI1T3Q2hSDhZ7v4Ev3r+CvnVXrSW8T6QIrvVSIiEoKzkVCVISkGiRNK4vPMx9G17TPMDezs2oAer95E1Z7vIq33GbCS0tR42vYusASEZUEDg0YGzZswIABA+Dv7w9N07Bw4cLbPmfdunVo2bIlPD09ERwcjJ9++qlYtpXoToTXq3pjcC5Nw2l4Y0X9CdjSewHO+7SHp5aJ592WYHO5cVjZdhd6BVdUDUUnLj6gLomIXJmbI988JSUFYWFhePrppzFo0KDbrh8bG4t+/frh+eefx6+//orVq1fjmWeegZ+fH/r06VMs20xUGDLIlgy2ZRukyz7oVoeuwOGVQMQ7KHv2IBrs/RiX9k/D+uv3Y561O6ZtjEXPEB8OP05ELstpZlOVEowFCxZg4MCBea7z5ptvYunSpdi3b5/9vsceewyXLl3C8uXLC/Q+bORJTsVqwb7l38Nr8ycINN1okxFn9VFVKoutHWCBiWNoEJHTKLGNPDdt2oSePXtmu09KLuR+IpdkMuM3S2f0zvgU72SMwFndC0GmM5js8R8s9RiPPuYd2Hz0nKO3kojItapICispKQk1amSf10FuS6JKTU1F2bJlb3lOWlqaWmxkXSJna6chVSK/Wu/F3LQuGO+9Hg+kzEUj0wl8a/oUl2JXYeeav2HJ1RCU9XTjuBlE5BJcKmDciUmTJuG9995z9GYQFaKdxkNYu3sM3Dd/hbZn5qDyhT1ouWEkYA3GF5kP4U+9mQoknEiNiJyZS1WR+Pr64vTp7K3r5bbUA+VWeiHGjx+v6opsy4kTJ4ppa4kKFzLe7h9qbwTarXkDpHb+f+h47TP8kHkfruvuaGk6gukeH2G++wR0N+/B7K1x7HFCRE7LpQJGeHi46jmSVUREhLo/L9KdVQJI1oXIVcbQuKBVwT8zh6FT2hf4PrOvGqyrhekIprl/hDHHnkPc5gUYNWMbu7USkdNxaBXJ1atXceTIkWzdUHfv3g1vb2/UqlVLlT6cOnUKM2bMUI9L99Svv/4ab7zxhuraumbNGsyZM0f1LCEqqW0zZAyNs3plbKn/KtK9X0LLkzPQ8sxvaG46hv+a/o295rpYuOwJ7I7vh9QMnfObEJFTcGg3VRk0q1u3brfcP2LECDWA1pNPPonjx4+r9bI+55VXXsGBAwcQEBCAt99+W61XUOymSq5EBty6ZQwNAOt37sfB3/6FYeZVKKfdaMQcYw3A1Mz7sdjaHlbNjfObEJHhONmZQQeHyNnDx48rt6Lzudl4wrwKFbVUdf8Ja3V8Z+mH+dZuGNKxgWrbQURkBAYMgw4OkSuEjFEztsMLKSpkPO32B6ppN7pin9UrYYP3wzgZPBSXrWVZdUJEd40Bw6CDQ+RK1Shl3M3ITEtBm4vL0PzEDFTLvNHjKlkvi18tPfFj5r04p1Vh1QkR3TEGDIMODpGr+ufve3Fxy//hOfMiNDCdUvel62YstoZjWmZfwK8ZxvZsYO+tIg1Ks7bxICLKDQNGPhgwqDRVnZhhRTfTLjzrthRtTQftj0daQ/FDZl+stTaHSTOrGV85cBcR3Q4DhkEHh6ikVJ1cz7CgdloMAmN+xD1pf8JNs6p1jln98F/LfVho7YzHOja8pUGovAZLOIjIhgEjHwwYVJpJYHhnxgo85bYCj5nXoJJ2Td1/Ua+AjZXvR3zwUJxDFVVlIlQpiKaxhIOIFAaMfDBgUGlnK9noGOgJ/9j5qLrvv6iemWRvp7HEGo7pmb1RvVEHrD14VoULCRlPdqzNLq9EpVxyIWoBSvxkZ0SUnTTmtDfoDHsD75vuQ+KWeXjavAytTYcwyPyXWqLjguGl9cQyhCNVd1eDfRERFRQDBlEp1z7YB6Mi22K5pS2aakcwwm0l+pk2I8R6BJ96HMG7ppk412Aw6vg1cfSmEpELcehQ4Y7AKhKi/BuEros5g9OJJ/GoaS2ecFuNmto5tY4ODUeqdMKWaoOQWLW9mveE3VuJSpfkQlSRMGAQUe5dXDUN0DPxSVgSqh2Yjk7mffZ1jll98Yu1F+ZldsKnw7tyDA2iUiKZAcOYg0NUWmWdZE26qU6PPI4gnMIwcwQeMm9ApZvznlzX3bGrQhfE1HwIZ71bIjXDylINohIsmQHDmINDRP8r0TABkNEzyuM6Bpr/wlDzKoSY4u2H6IjVH7Ot3VSpxr+Hd1f3cQwNopKFAcOgg0NEuQ/apS7TM2FK3IngE/PR3xSJ8jenjZeurn+6heO/qZ2xVW+MTF1DzxAfDG5Ti1UpRC6OAcOgg0NEBSvdqIhU9DdH4jHzWoSZjtkfj7P6YLalG+ZZOuMMqjBoELk4BgyDDg4RFb4Hipa0F4NNa/CAeaO9rUambsJ6axjmWzpjlbUlOofUZIkGkQtiwDDo4BDRnfdA8dBT8Xy1vbgneSlamQ7b17mkl8fvlg6qVMOnYXsMbhvEqhMiF8GAYdDBIaK764EiI4bK7XUbN8I/boEaIdRPu2Bf95C1pgoaHQeNRpdWTXm4iZwcA4ZBB4eIjCNBY87WWFw/tFZ1db3XtA1ltAz1mBUmHK3UDqmNB6NZ9yGAexn7c9gThch5MGAYdHCIyHgSGmZvi8eW6OMYYN6MB80b0MZ0yP74Na0cLgbdh6Sg/nhkhTs0zczZXImcBAOGQQeHiIq+KiXufApiY/ZioGkDHjT/hYCbQ5OLM3plLLaE43dLR2T6hqFroxpITbdwMC8iB2HAMOjgEFHxNQzV1C0rWmuH8KB5I/q5bYEXrtrXO2b1wyJLByzWO6jrHFuDqPgxYBh0cIioeKtNVkWfsY8Y6qll4h5tD4aU3YqOmVtQVku3r7/bWleVaiy1tMM/h/dmLxSiYsKAYdDBIaLiDxqTVx1CdGIyrDpUd9dujapjU3Qcept2qLE17jFFwU2TCHLDfrdQmJs8iDMBfbAu0R1lPcysRiEqIgwYBh0cInLsWBoWXcf3w1ur+22DeSWcikPFo0twv3kTWmdpHCq2WxuoUo3l1nZI1L3Vc6WrLBEZgwHDoINDRM4xlkZeJR0XEo+jj2kr+pk3Z+uJIrZbG+J80H0oGzaIJRtEBmHAMOjgEJHrlHTUwAX0NW9F31zDxo2SjRXWtkjQq7Jkg+gOMWAYdHCIyHVKOrJWo5w+eQwVji3LtWRjt7UedpfvCEuDfjhlDkR4cDVWoxAVEAOGQQeHiFzTxMUHMD3yuL1k4z7zjWqUVtphmDTdvl6s1RcrrK3RsvdQtL2nDyIOnuXIoUT5YMAw6OAQUcmoPukZUgMNfStiT3QMAs+uRy/TdnQ07VNdYW1S3L2xKDUMq/U2+MsSik4hAWjoW4k9UoiyYMDIBwMGUeltKGoLHjLWRlmkoqtpD3qZt6O7aTcqadfsz72ql8E6axhWWlpjvd4cl/XyamAvBg4q7ZIL8SNd03X9f+WFpQADBlHpZgse0lZjXcwZNeaGWc9EuCkaPU3b0du8Hb7aRfv6mboJO/QGWGNpgTXWFjiKmrDqGhuKUqmUzIBhzMEhotJVlTK6WzAOJV7C6Zgt6GPerqpSGphOZXvOCWt1rNVbwK3RvXj80ccB97IO236i4saAYdDBIaLSW5ViK+U4fyIGHrGrVDVKuOkAPG9OMS8s5jIw1+uK6ArhWJERhsYhoeyRQiVaMkswjDk4RERZA0dFUzrMcRvgfWodupt3wU+7kO0ARVsDUSb0PtRp/wBWXQlC5PErnPmVShQGDIMODhFR3l1grQjR4tHNtAvdzbvRQjsMc5YusNJQdJO1MdZbm8Fapyu6dQhn6Qa5PAYMgw4OEVF+08tnbSFfBVfQ2bQH3cy7cY9pH6ppydmeF2f1QaxXW1z074yYsi3RqkEtFTjk9TYdPc+SDnIJDBgGHRwiottNL591rA2JHGsPnoVVtyBUi0MX0150MkWhlSkGHpolW8+UnXp9pAV1xcdHAhCNusjI0jOFoYOcFQOGQQeHiKgwDURz9koJC/DCnpOXUQGpaGc6oMJGZ9Ne1DUlZXudi3oFbLaGQq/TGV6hPfHEwvMwayb7bLKcEZacBQOGQQeHiMiI0GEr7ZABvqwAamlncI8pCv3LR6NJWvZBvkSSXgWR1saqDcfZau3wRJ97GDLIKbhcwJgyZQo+/vhjJCUlISwsDF999RXatm2b67o//fQTnnrqqWz3eXp64vr16wV6LwYMInKErF1fY5KS7dUr0DPxUI3TqH5uKzqa96OVdihbV1gRa62BeK/WqNK4J5rdMwAR8VbVbqOsh5lDmVOxKsw51A0ONnv2bIwbNw5Tp05Fu3btMHnyZPTp0wcxMTHw8fHJ9TmyU/K4jSZfUiIiJyalGbaqDumJYqtGMWtu2G9uhIPWmphiGYhyWgaeDjqDgEvb0eDaDjTTjqGO6TTqXFkKbJblFQRYAxFobYzN1kbYoTfCtI2V1FDmg9vcaDhK5AwcHjA+++wzjBo1yl4qIUFj6dKlmDZtGt56661cnyOBwtfXt5i3lIjIGOH1qmLaxlh7yOja0Af7E5LV7Wu6O8I6PwDgAQyasR0VcQ1tTAfRwbQfHUwHEGqKQ4jphFqewnL1eoetNbH1cCMsjmmElfXuQe/wVgwaVLoDRnp6Onbs2IHx48fb7zOZTOjZsyc2bdqU5/OuXr2KoKAgWK1WtGzZEh988AEaN26c67ppaWlqyVq8Q0TkSFLKII03s7bVaB5Y+ZYRRWWdG+03ymGttaXqFuuNZDWiaFtTNNqaDqqgUd90Si1PYDVwYgri46ojoe498A/rCQR1AKrUQUT0GVarULFyaBuMhIQE1KxZE5GRkQgPD7ff/8Ybb2D9+vXYsmXLLc+R4HH48GE0a9ZM1QF98skn2LBhA/bv34+AgIBb1n/33Xfx3nvv3XI/G3kSkat3i5Up6E+eOoFrRzainemgKulooknJSPZ/69fL+mDl1XrYZm2EzdYQHIM/LLqJPVSo5DbyvJOAkVNGRgZCQkIwZMgQvP/++wUqwQgMDGTAIKISMW9KzgBSSbuO5lqMKt2Q0NFMOwpPLTPb61zSy2OXtT5OVGgKjzrhiPVoiNYNAlmtQiWnkWe1atVgNptx+vTpbPfL7YK2sXB3d0eLFi1w5MiRXB+XHiayEBGVpIaiud1vCyBx52vh84PNYcnUURbp6Ol1ArWv7lGBo6XpMCprKWrEUaTuBg78DIuuIXpnENZ4haF6aGc0bd8LqFzLXq0ibUbYeJQKy6EBw8PDA61atcLq1asxcOBAdZ+0q5DbY8aMKdBrWCwWREVFoW/fvkW8tUREzi1r0LCNu5EKDyxNrgerXg/RDWpgi09ZxEdvgfeFXWipHVKBo6Z2Hk2048CV48CWRcAW4Kp7VWRerws3vQG+iawP7fGH0bNpoHofjjRKBeHwcTCkm+qIESPw7bffqrEvpJvqnDlzcPDgQdSoUQPDhw9X1SiTJk1S60+cOBHt27dHcHAwLl26pMbPWLhwoWosGhoaetv34zgYRFQaSAiYvOoQohOTYdWh2m482bE23u4fah9x1Dbwly/Oq6DRSi2HEKodzza0uciAG85UCMF+UwMsPueHKD0Yx/XqN9uCVMo2HgcDSMnlMlUkYvDgwTh79izeeecdNdBW8+bNsXz5chUuRHx8vOpZYnPx4kXVrVXWrVKliioBkTYcBQkXRESlha1KI+vQ5dJ2I2cvlhsDf/lgWXRV/GFtr3qqeCIdYdoxtLgZOCR8yORtNa9GoSai0Nvjxntc0Ctg75F62HO4Hg7q9bBoY118HVBLDY8u7yldcTnUeenl8BKM4sYSDCIqTfJqGFqQniqjuwVjXcwZRCdeRgDOqFFGw0xH0dx0FCFa3C2NR8UJa3Xs0etht7Ue9un10LxdF7z1QGuWapQQLtOLxBEYMIiICj+Bm606xaRBVbm0qlkOGQlRaG46hmamowjTjqKelgBTji6yVpiQ4B6Ev1KDsF+viyhrbbw0ZCB6NAuyvx8bkroOBgyDDg4REWWfR+V6hiVb+Mg6v8qW6OMIM8WiiXYUvbxOoX5mDCqln7nlEFpgwoVydXFIq4uIS744oNfBfmstTB7eST3OwOG8GDAMOjhERHTnpR8y58ofkTtV4GhmOqYGAWtiikV17dYRla26hhMmf+zODFKBY6+1Np59dCC6NW+Q7fUZPhyLAcOgg0NERHfOVr1ia9MR6F0Opy6moLp+UQUN6RrbxCRLLPy0C7m+xrXyAbhSpTEir/ph6ZmqiNGDcEKvqnqv2CZ3swUPzi5b9BgwDDo4RERkXKmGyK09hwSPtIuJCNGOo7Eq5TiOplosAk1nc/8/rpfDQb0Woq2BuFq5ISLO++CwHoAUlLG/JnuvFA0GDIMODhERFX17DnFrScc1VNSv3ijh0GLRUGaQ1eIQrCXAPccYHbYqluN6DRU8DlprIaVKI/jUb43Tmg/Cg6txJFKDMGAYdHCIiMhxJR22wCFs16X3yrWEaIRo8WhoileXIaZ4VNcu5/q6V/SyiNZrwa9+KwQ2agVUbwRUD0FEXCarVe4AA4ZBB4eIiJwjcOTsOptzzI5BDdzRvnwSrhzfhSpXD6ngUU87dcuIpDZn9Uo4bA3AIT0AR/QAxFgD4F+/OQJqBqhRSXNrz8FGpmDAyA8DBhFRyR1ELOu4HWZkItiUgIaIRyNTPBpqJ1HfdBIB2rk8X/OMXhmHrTWzBI+a8K/fAoE1a2LK2iP2QFNa23gkc6AtYw4OERG5fjuPuPMpWHvwrL26pRyuo74KG6fUZYOb1/MLHqcleEjosPrjuO6P4MYtMLRfT6BSTfusszlLPSIK0bvFVUpHGDAMOjhEROT68hqN9MZEbRXVIGFS3VIRqapapYHpJILl8maJh8w2m5frWhkctvgiVvdX4SMWfuqyd6eO+GLDyVveU4Zfz1kFk7PNiTOXjjBgGHRwiIioZI9GmtfjttFJJXh4aamoi5N4NjQD2rlDMF+4MSx6Le003DSJD7k7qVfDMasfjur+arkRQvxwFlVghZYl6PjYS1iyznrrjBgwDDo4RERUuuU2Oun0yOMqDLgjU4WMulqiChzq0pSgrlfWUvJ8zWu6J+L0Gqpbbbzui8vlArDzalWc0H2RoFdGjxBf+yBizlaFwoBh0MEhIiIqfHXLaVTXriIIpzCg5lWEuJ9B5WuxqJYWj0qpp/It9UjVPRB/M3zo3nWgVQ1GVIo3fovzxBmtKjJ1zV7N4oiwwYBh0MEhIiIqTHVL1hIOcy7VHauiTiAm5gBqZJ7CufholEmOQ20tCbW10wg0nYFZxZbcpenuiNN9EKf7qiVWr4GqgQ3RtmUr3NOqOSJiLhR5KQcDhkEHh4iI6G7mX/k+nwabOdd9qUtt/L5hswobQTdDh4SPIO00amlnch3BNOsMtYm6N05IALH6oFXzFkitEIitlyqhTmgb9GhW25APkgHDoINDRERkxPgcBV036yBiWathNN2CR4I1tKp4AXv27rpZ6pGEQO2sCh9ltfQ83+OhtAl4ftgThpRqMGAYdHCIiIicrRomIksI0QDo0FEdlxGonUGQdkZdhpQ5j6oZier6g+n/RL+OLQzpmcKAYdDBISIiclYRuQyZHhbghT0nL9/SCNWosTUKcw51u+t3IyIiomLXK7SGvVrFVs0ijTz3nUpWYUNCRohfJYzt2cAhXVsZMIiIiEpA0LCZtjHWXqLhqHAhGDCIiIhKiF6hNVR1SEEbmRYlBgwiIqISXKLhKFJFQ0RERGQoBgwiIiIyHAMGERERGY4Bg4iIiAzHgEFERESGY8AgIiIiwzFgEBERkeEYMIiIiMhwDBhERERkuFI3kqeu6/YZ4YiIiKjgbOdO27k0P6UuYFy5ckVdBgYGOnpTiIiIXPZcKtO250fTCxJDShCr1YqEhARUrFgRmqYZlugksJw4cQKVKlVCScB9cg38nFxDSfucStr+CO5TwUhkkHDh7+8Pkyn/VhalrgRDDkhAQECRvLZ80UrKl82G++Qa+Dm5hpL2OZW0/RHcp9u7XcmFDRt5EhERkeEYMIiIiMhwDBgG8PT0xIQJE9RlScF9cg38nFxDSfucStr+CO6T8UpdI08iIiIqeizBICIiIsMxYBAREZHhGDCIiIjIcAwYREREZDgGjAKaMmUKateujTJlyqBdu3bYunVrvuvPnTsXjRo1Uus3bdoUy5Ytgyvv008//aRGPs26yPOcyYYNGzBgwAA1wpxs38KFC2/7nHXr1qFly5aqBXlwcLDaT1fdH9mXnJ+RLElJSXAWkyZNQps2bdRIuj4+Phg4cCBiYmJu+zxn/j7dyT45+/fpm2++QbNmzeyDToWHh+OPP/5w2c+osPvj7J9Pbj788EO1nWPHjoWzfE4MGAUwe/ZsjBs3TnXL2rlzJ8LCwtCnTx+cOXMm1/UjIyMxZMgQjBw5Ert27VL/cGTZt28fXHWfhHwxExMT7UtcXBycSUpKitoPCU4FERsbi379+qFbt27YvXu3+mI+88wzWLFiBVxxf2zk5Jb1c5KTnrNYv349Ro8ejc2bNyMiIgIZGRno3bu32te8OPv36U72ydm/TzLasZywduzYge3bt6N79+544IEHsH//fpf8jAq7P87++eS0bds2fPvttypE5afYPyfppkr5a9u2rT569Gj7bYvFovv7++uTJk3Kdf1HH31U79evX7b72rVrpz/33HMuu08//vij7uXlpbsK+dNesGBBvuu88cYbeuPGjbPdN3jwYL1Pnz66K+7P2rVr1XoXL17UXcWZM2fUNq9fvz7PdVzh+1TYfXK175OoUqWK/sMPP5SIz+h2++NKn8+VK1f0+vXr6xEREXqXLl30l19+Oc91i/tzYgnGbaSnp6vU27Nnz2zzmcjtTZs25focuT/r+kJKB/Ja3xX2SVy9ehVBQUFqkqPbpX9X4Oyf051q3rw5/Pz80KtXL2zcuBHO7PLly+rS29u7xHxOBdknV/o+WSwWzJo1S5XISNWCq39GBdkfV/p8Ro8erUpicx5/Z/icGDBu49y5c+oPskaNGtnul9t51W3L/YVZ3xX2qWHDhpg2bRoWLVqEX375Rc1K26FDB5w8eRKuKq/PSWZVTE1NhauRUDF16lTMnz9fLfKPsWvXrqoKzBnJ35BUS3Xs2BFNmjTJcz1n/z7dyT65wvcpKioKFSpUUO2Tnn/+eSxYsAChoaEu+xkVZn9c4fMREpTk+y3tgAqiuD+nUjebKt0ZSfpZ07582UJCQlS93/vvv8/D6gTkn6IsWT+jo0eP4vPPP8fPP/8MZ/zlJXW/f/31F0qKgu6TK3yf5G9J2iZJicy8efMwYsQI1d4kr5OysyvM/rjC53PixAm8/PLLqt2PszZAZcC4jWrVqsFsNuP06dPZ7pfbvr6+uT5H7i/M+q6wTzm5u7ujRYsWOHLkCFxVXp+TNO4qW7YsSoK2bds65Ql8zJgxWLJkieopIw3w8uPs36c72SdX+D55eHionlWiVatWqiHhF198oU6yrvgZFWZ/XOHz2bFjh2qUL73gbKRkWv7+vv76a6Slpan/8478nFhFUoA/SvljXL16tf0+KS6T23nV38n9WdcXkjLzq+9z9n3KSf6QpchRiuVdlbN/TkaQX2zO9BlJe1U5EUvx9Jo1a1CnTh2X/5zuZJ9c8fsk/yPkpOWKn1Fh98cVPp8ePXqobZLvuG1p3bo1nnjiCXU9Z7hwyOdUJE1HS5hZs2bpnp6e+k8//aQfOHBAf/bZZ/XKlSvrSUlJ6vFhw4bpb731ln39jRs36m5ubvonn3yiR0dH6xMmTNDd3d31qKgo3VX36b333tNXrFihHz16VN+xY4f+2GOP6WXKlNH379+vO1Nr6l27dqlF/rQ/++wzdT0uLk49Lvsj+2Vz7NgxvVy5cvrrr7+uPqcpU6boZrNZX758ue6K+/P555/rCxcu1A8fPqz+1qQ1uclk0letWqU7ixdeeEG1zl+3bp2emJhoX65du2Zfx9W+T3eyT87+fZJtlV4wsbGx+t69e9VtTdP0lStXuuRnVNj9cfbPJy85e5E4+nNiwCigr776Sq9Vq5bu4eGhunhu3rw524c6YsSIbOvPmTNHb9CggVpfukIuXbpUd+V9Gjt2rH3dGjVq6H379tV37typOxNbN82ci20/5FL2K+dzmjdvrvarbt26qnuaq+7PRx99pNerV0/9I/T29ta7du2qr1mzRncmue2PLFmPu6t9n+5kn5z9+/T000/rQUFBavuqV6+u9+jRw34ydsXPqLD74+yfT0EDhqM/J07XTkRERIZjGwwiIiIyHAMGERERGY4Bg4iIiAzHgEFERESGY8AgIiIiwzFgEBERkeEYMIiIiMhwDBhERERkOAYMInIaMr28THdORK6PAYOIiIgMx6HCicgpPPnkk5g+fXq2+2JjY1G7dm2HbRMR3TkGDCJyCpcvX8Z9992HJk2aYOLEieq+6tWr5zrtNBE5PzdHbwARkfDy8oKHhwfKlSsHX19fHhQiF8c2GERERGQ4BgwiIiIyHAMGETkNqSKxWCyO3gwiMgADBhE5DekxsmXLFhw/fhznzp2D1Wp19CYR0R1iwCAip/Haa6+pXiOhoaGqB0l8fLyjN4mI7hC7qRIREZHhWIJBREREhmPAICIiIsMxYBAREZHhGDCIiIjIcAwYREREZDgGDCIiIjIcAwYREREZjgGDiIiIDMeAQURERIZjwCAiIiLDMWAQERGR4RgwiIiICEb7//H1Cx4lghdQAAAAAElFTkSuQmCC", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "plt.figure(figsize=(6, 4))\n", + "plt.plot(t, y, \".\", markersize=4, label=\"noisy data\")\n", + "plt.plot(t, features(fitted_rates) @ fitted_amplitudes, label=\"variable-projection fit\")\n", + "plt.xlabel(\"t\")\n", + "plt.ylabel(\"y\")\n", + "plt.legend()\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "b80973e3", + "metadata": {}, + "source": [ + "## Why eliminate the linear parameters?\n", + "\n", + "Variable projection optimises only the $n$ nonlinear rates instead of all $2n$ parameters at once, and the reduced problem is better conditioned — so a second-order solver reaches the solution in far fewer steps. Below we compare against fitting all parameters $(k, c)$ jointly. Both start from the same rates, and to be fair the joint solver is given the best amplitudes for that initial guess, $c_0 = \\Phi(k_0)^\\dagger y$." + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "id": "e48a78d3", + "metadata": { + "execution": { + "iopub.execute_input": "2026-06-07T17:55:53.606059Z", + "iopub.status.busy": "2026-06-07T17:55:53.605961Z", + "iopub.status.idle": "2026-06-07T17:55:54.548793Z", + "shell.execute_reply": "2026-06-07T17:55:54.548303Z" + } + }, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiMAAAFzCAYAAADsTAnbAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjksIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvJkbTWQAAAAlwSFlzAAAPYQAAD2EBqD+naQAARjZJREFUeJzt3Qd4VFXawPE3hSQEAqFD6ALSpYOgYqGrrILr2gELLC4qyrei7CrYUVkVEZTVVdG1oSJYEKQKwiIgiEgHRemdJKQX5nveE2acSSaTSTJ9/j+fMTN35t45c3LJfeec95wTYbFYLAIAAOAnkf56YwAAAEUwAgAA/IpgBAAA+BXBCAAA8CuCEQAA4FcEIwAAwK8IRgAAgF8RjAAAAL+K9u/bB7azZ8/KoUOHJCEhQSIiIvxdHAAAgobOqXrmzBlJSkqSyEjXbR8EIy5oINKwYUNP/34AAAgb+/fvlwYNGrh8DcGIC9oiYq3IKlWqePa3AwBACEtNTTVf6K3XUlcIRlywds1oIEIwAgBA6bmT5kACKwAA8CuCEQAA4FcEIwAAwK/IGQEQkEMC8/LyJD8/399FAeBCVFSUREdHl3v6C4IRJ2bMmGFu/CEEfC8nJ0cOHz4sGRkZVD8QBOLj46VevXoSExNT5mNEWPQrCIodllS1alVJSUlhNA3go4kGd+/ebb5t1apVy/xxY8JBIDBp+KBfHo4fP26+vLdo0cJhcrPSXENpGQEQMPQPmwYkOjeBftsCENgqVqwoFSpUkN9//938+42LiyvTcQhGfCl5v0jGyaLb42uIJDb03r5AkClp6mgAofXvlWDEVzSYmN5FJC/byW8hVuSeDcUHFeXZFwCAAEcw4ivaquEsmFB52bJj6TuSWdl5QFExbb+0crGvOTbBCAAgSBGMBIhWP0/xdxEAAPALgpEAsTe6qWRFVHT6XJwlU5rm7fV5mQAA8AWyxAJE0ztmSet/rnF60+cAhKfLLrtM7r//fo++vrTH9LURI0bItdde69X38GUdnDx5UmrXri2//fabbdvf//53r3/G8n7mG2+8UV544QXxBVpGACCAffbZZ2boZDh5+eWXzRwWnrwId+zYUaZOneqXen366aflmmuukSZNmti2bdq0SS6++GIJZI888oj07t1b7rrrLjNfiDfRMuIrOgRXR744o9v1+bLsGxXjel8gDC3cclgGTl0pLR9ZYH7q42Cjczao6tWrS0JCgoQC62cqiV74EhMTvVoWX9WrziT85ptvyp133umw/aeffjIBUiBr166dNGvWTN577z2vvxfBiK/oaBcdgjtqRdFbSUNz7fYdX+MVuSr7KUmu2rrguWZ9GUmDkKXfjjNy8kp1+3zTQRn93kbZeeSMZOedNT/1sW4vzXFK88389ddfl6SkJDNhmz39NnzHHXeY+wsXLjTfhPUiW6NGDbn66qvll19+cfj2fs8995hm9Jo1a8qAAQOcNq2XdByl6/rosfSirsd69NFHXX4eLffkyZOladOmZhKrDh06yKeffuryM1vL6+p9nH2m7Oxsue+++0y3hU6QpZ9l/fr1LrtpSiqfPv/8889L8+bNJTY2Vho1amRaI6zHWrFihWlt0dl89abdJYXr1Z1y6T76mvHjx5tgpm7duvLYY4+5rKevv/7alOnCCy+0bTtw4ICcOHHCFowkJyfL4MGDzXseOXLE6XFKKl9py/buu++a80ePa0/r/bbbbrM91nJ99NFH4m100/iSBhVlHYJ7bt/c2iJbDx6Uxc0fles33CKy62uR39eINO7p6dICfpeZmy9tJn5Tpn0thX6O/WhTqfbf9sQAiY9x70/k9ddfL/fee68sX75c+vTpY7adOnXKBA56MVLp6ekybtw4ueCCCyQtLU0mTpwoQ4YMMc311kmj3nnnHbn77rtl9erVxb6Xu8fRb+Lr1q2TH374QUaNGmUu0CNHjnR6TL3Q67ffmTNnmim9V65cKbfeequZkv/SSy8ttizuvE/hz6QXyzlz5pjtjRs3NkGEBil79uwxF9GylG/ChAnyxhtvyEsvvWQu0rq20Y4dO8y+GoTs2rXLfMt/4oknzDbdrzB3y6XPa/2vXbtW1qxZY4Kdiy66SPr16+e07N9995106dLFYZv+rjSY1G6bn3/+WYYOHWrOG1ddR+6UrzRl03NWg5cvvvjC3FfHjh2T+fPny6JFi2yv6969uwnsNGjRoMpbCEaCTP3EghE3m/IayfWdh4lsfEdk4cMiI5frNHj+Lh4QlqpVqyaDBg2SDz74wBaM6Dd3bQ24/PLLzePrrrvOYZ+33nrLXBS3bdtmLpRKL7R6kXHFnePodPp6YdZWgJYtW5oLnj52FozoReaZZ56RJUuWSM+eBV9qzjvvPFm1apX8+9//dhmMuPM+9p9JA6nXXntNZs2aZepLaRCxePFi05Xx4IMPlrp8nTt3NgHH9OnTZfjw4eZ57Vqw5mNoq42ucaTLC2hrgTOlKZcGgZMmTbJ9Nn3fpUuXFhuM6DTp2mpWOBjR1h09X7Tl6Lnnnis2UCxN+UpTNm1huvnmm+Xtt9+2BSMa8Gkwqa0sVlp27V7TFhsNgryFYCTIVu1tUK0gGDlwOlPkL4+KbJ0rcniTyE8fiHS61d/FAzyqYoUo00JRGtfOWC27j6bZWkSUrm5+fp0Emfu3XqV679K45ZZbzAXl1VdfNd8g33//fTMawdpaoQsAaiuGfmvVJnprl86+fftsQUThb9DOuHMc7RKwX2BQL+I6KkL/pukihPb0m7XmNRS+YOkFqFOnTi7L4s772H8m7U7Kzc0139attCVAv31v377d6XuUVD7dTwMWaxBYFqUpl17w7elqtdqiUJzMzMwi67VoMLJ582YTiGhLRM9zQVZ5y1fasun52q1bNzl48KDUr1/fBDvammL/O9WgRXl7FW2CESfGjBljbtYVBwNJg2oFi4cdOJ0hUrmWyKXjRRY9IrLkcZHWfxKJc70yIhBM9I+iu10lVuP6nW9yRPTvqaYvWH8+0Pf8Uh+rNLRvXfMl9OKif+C1eV5bCeyf12+W+o3Wml+iwYN9UmelSpXcep+SjlMa2tWjtNx6QbLniWZ5dz5Tecqn+Ra+VLgbRc/RwrlC9rR17PTp00WCEe2a0ZYRT5a/QinLpsGcttBo/kj//v1l69atpp7taXdjcV1bnkS7fpCxbxkxiWLd/ypSvZlI+jGR73wzHhwIZAPb1ZOZt3aWVnUTJDY60vyceWsXGdjOeRO9p+i3X73AaIvIhx9+aLottAvBOs/Ezp07zVBJ/QbfunXrIhcod7h7HG01sff999+bZvvCrSKqTZs25qKuLSuaAGp/024YV0rzPtbuE+0ysc+J0W/8moip5XCmpPLp++m3d+2OKI6+p6uW7rKUy116wdcuNKszZ87Ir7/+ar7wajeKtp5t3brV5TG8WT4dtqstItpd07dv3yK/8y1btkiDBg1MUOVNtIwEmXqJceabno4SOJGWI7USYkUGPCPy4Q0i378qonkkNZr5u5iA3wMSvfmadtXo6Ba9uGiCpX1OiY5c0FE32nSuF9aHH3641Md39zi6XRMZ//rXv8rGjRvllVdeKXbyKh3eqhNwPfDAA+ZbtOZapKSkmAtflSpVbHkYzpTmfaytJJrMqjkOmnSp+QmaT6JdAIWHvpamfA899JBJ8NQLtnZlHD9+3PwOrMfURFENnHQUTeXKlYskypalXO7SJFNNsNWgUX9/OqRXgzUNIjRQ0Yv94MGDTRJwcRd8b5ZP80a0frWlTVtICtMWPm018TaCkSATGx0ldRLi5EhqlhxMziwIRs4fINKsj8gvS0UWPSpy0wf+LiYQlq644gpzsdDWC/0jb6V5Izo8UkcvaJeKtppMmzbNIVHQHe4eZ9iwYSZXQXMK9MI3duxYM9KlOE8++aRphtdRK/qtXUd6aKvOP/7xD5flKe37qGeffdYEFTp8VFsJunbtKt988425UJe1fDqkODo62uTSHDp0yARqo0ePtu2vF1sNWjQA0PLu3bvXI+VyR/v27U1ZP/74YxO0aRdNq1atbF1gU6ZMMXkf2qqmSboaUHmq3tyhqQiaFK3dM4VnhM3KypJ58+aZUWHeFmHx5DR3IcaaM6JRuEbggeLPr/1Pfvj9tEy/uZNcfcG5LO1jO0Re6yViyRe5bZ5Is4IMfiCY6B8/vVDofBKFk/4QWJzNauopN910kwlufDHZli/ohV5bNbQVxJrQHEj69Okjbdu2NYGtPR3BM3fuXIehvqX5d1uaa2jg1QpKN6LGqnYrke7nhoYtnCCSn0dNAggqOlmb5lfoHBl6cQwVV111lWkx0lErgeT06dMm2Pj2229NDouzhFjtevMFummCkMOIGnuXPiSyebbI8e0iG97+IzgBgCCgLQe9evUyc7PYd7OEgkBcmLBTp04mINF5TrTLz1lyq68QjIRKy4iKry5y+T9Fvv67yPKnRdpdV7ANADxMv017mnb7eHs+C/zBfhVhf6ObJqhbRgoFI6rL7SK124hknhb59lnfFw4AgFIiGAlC9W0tIxlFF7+KihYZOLng/vr/FCS2AgAQwAhGglBSYkG2clbuWTmZ7mTWxfMuE2l1dcHImm8mFEw/CQBAgCIYCda5RqrEFt9Vo/o/KRIVI/LLMpFdZVv1FAAAXyAYCbURNVbVzxO58G8F97/5h0he2datAADA2whGgnxEzcHiWkZU77+LVK4jcuoXkXX/9l3hAAAoBYb2htrwXnuxCSJ9Jop8PkZk+WSROu1EKtpNHRxfQyTR9UJYAAB4G8FIqHbTWDW5RGf9F8lNF/mv47oDEh0rcs8GAhIAgF/RTRPKLSNK5xuRYkbT5GWLZJz0QukAAHAfwUgITHzGWoeA/40YMaLIqqclLTTnrSnC9W+CroWiKwhHRESYVW7L+l4nT56U2rVrF5mtU1fCLc3n9Yby1uGNN94oL7zwgkfLhLKhm8aJGTNmmFt+fr4E+lwjmbn5cio9R2pULhjqC8A/Xn755VJ9Mfjss8/MQmTeWClXl3yfNWuWmbL9vPPOMyvFVqxYsdTHUU8//bRcc8010qRJE4ftmzZtkosvvliC2SOPPCK9e/c2a7Do6rLwH1pGnNDVC3XlyPXr10tQzzUChKPk/SKHNhW96XYv0ouZtkC4S1stEhISvFKWX375RerVq2cWnatbt65p2SjLe+k6MW+++abceeedRZ776aefTEATzNq1ayfNmjWT9957z99FCXsEI6G6Rg0QjjTgmN5F5PVLi950uxcDksLdNNnZ2XLfffeZQCAuLs60Ith/wSncxaCP9fXjx483gYoGEY899pjD8VesWGFaYLTrRW/OFjrT1917772yb98+8xpt0bB/L3ePo77++muJjY2VCy+80GH7gQMH5MSJE7ZgJDk5WQYPHmw+45EjR5y21OhzGqzVqFFDrr76ahMwufvZVXp6ugwbNkwqV65sAi13u1e0HoYPHy516tQxrUMdOnSQVatW2Z7Xcn/00UduHQveQzASxOon/rFGTbF0+K6OmnFGt+vzQKDSbo+cdPdvqQcKErOd0e36vLvHKucyCnphnTNnjrzzzjuyceNGad68uQwYMEBOnTpV7D762kqVKsnatWvl+eeflyeeeEIWL15sntPgoWfPnjJy5Eg5fPiwuTVsWHRovr5O92vQoIF5TeEWXnePo7777jvp0qVLke3aRaOBhQY6P//8s3Tr1k3q168vy5cvN4FEYRpIjBs3Tn744QdZunSp6TYaMmSInD171q3Prh588EETRH3++eeyaNEi0wWl9erK77//Lt27d5fMzEz54osvZPPmzXLPPfdIlSpVbK/R59etW2eCR/gPOSOhPqJG5xHR4bvWUTPrXhfZ9L5Ig24if36bYb0IbLkZIs8kee54bw10/7X/OCQSU6lMb6MX39dee83kbQwaNMhse+ONN8zFVbs99MLqzAUXXCCTJk0y91u0aCHTp083F+9+/fqZbqCYmBiJj493esG30tdpl0xUVJTT17l7HOvFPCkpyWkwoi0MH3zwgbm4P/fccya4Kc51113n8Pitt96SWrVqme5w7Sop6bOnpaWZetPulD59+tiCFw24XLn77rtNq87HH39s26bHtqefLycnx7ToNG7c2OXx4D0EI+Ew14gGJNbJzXRWVg1GDvwgYvnjWwkAz9EuiNzcXLnooots2zRZVb+Fb9++vdj99IJsT7sjjh075rdfjbYoaBeTs2DE2sowf/5809Liyu7du2XixImm1UO7d6wtItqFYh+MFPfZtT41YOjRo4ftee3OadmypctAasGCBfLjjz+6LJs1sVfzY+A/BCOhMCV8cilyRnTNmqaXiuxdIfLjeyJX/NN7BQTKq0J8QQuFu45sdt36ccdCkboXuP/ePlZ4dI3mc9h3ZfhazZo15fRpnauoaDAydOhQ0zKi+SIl0bwMbXXQ1iFtidDPpEGIBhje+uxaRm0BKinJ1tptpi018B9yRkKkm6ZUc410GV7wU4OR/DwvlQ7wgIiIgq4Sd2/RfwxfdUqfd/dY+t5lpCM09EK4evVq2zZtKdH8jTZt2pT5uHpMT0w54O5xOnXqZLpS7J05c0Z+/fVXM+pQu1J0ro6tW7e6nKdk586dZhitdrG0bt3aaYBTUn1qsKItK1Z6jF27dhW7j74+Ly+vxBaPLVu2mO4eDbzgPwQjQSzpXAJrRk6+nM7IdX/HVleLVKwucuaQyJ4l3isg4GsBkrCtiZiar6C5ITqSRC/omlOhF0Znw2TdpQmjekHW0S/23R3eOo4m3GqgYR886JBezUfRoOqOO+4wn0dbPvQ4zlSrVs2MoHn99ddlz549smzZMpPMWho6gkbfR+tT99cAQkcFaSJscbRLR/Nj9PegXWP6O5g5c6bpMiqcpNu/f/9SlQeeRzASxOIqREntBOtcI6Xo79Q/yh1vLri/8R0vlQ7wA2vC9qgVRW8+Xofp2WefNYmbt912m3Tu3NlciL/55htzcS4rnfXUGghot4LmXHjzOO3btzdlt08A1e6PVq1amSG/asqUKSZ3Q7tt7LtdrDRg0KGzGzZsMF0zDzzwgNmntHSfSy65xAQ+ffv2NUOFnY30sdIA6MsvvzTBh4720dfriBodam2VlZUl8+bNc5l8C9+IsDCXeLFSU1NNZJ2SkuIwFCyQDH11tWzclyyv3tJZrmxfz/0dj+8UmdFdJCJK5IGtIlVKsS/gJXpx2Lt3rzRt2tRp4mQgu+mmm8wFPtQm0NIEVW2R0NYIVy0RwUhHPM2dO9cMFYbn/92W5hoaWmdWGHJ7RE1htVqKNOopYskX2RRafzwBX9K8BO0CWLNmjbRt2zbkKv+qq64y69wcPHhQQo3mlbzyyiv+LgYIRsJo9V5nOp9LZN34XxE/ZuwDwUxbDLp27WoCkdGjR0so0tlbi5sYLZjpmjSuhgfDdxjaG85Twre5RmTBQyLJv4vs/Vak2RWeLyAQ4nToKHNUAOVDN02Qq29rGSnDhD0x8SIX/KXg/gYSWQEA/kEwEq5zjRSec2THfJF050PzAADwJoKREFksT+caSS7NXCNWdduLJHUSOZsr8tOHni8gAAAlIBgJgblGatnmGilD3oh9Iqt21ZRzpVIAAEqLYCSkumrKuNBT+z+LVKgkcnK3yL41ni0cUAZMfwSE179XgpFwH1GjYhNE2g0tuE8iK/zIulgao1OA4GH991p4scPSYGhvCCh3y4jqMkLkx/+KbJsnMuhZkYpln7IaKCudwTQxMdG2dHx8fLxZvRVAYLaIaCCi/171363++y0rgpFwn/jMqn4XkdptRY5tFdn8iUiPUZ4rIFAKdevWNT+tAQmAwKaBiPXfbVkRjISAcnfTKP32qcN8F4wvWDyv+8hyLaEOlP1UjJB69eqZBc1yc8swQgyAz2jXTHlaRKwIRkKsm0abzcrcrK0ToC2eKHJ0i8jBjSINil8RE/A2/QPniT9yAAIfCawhNNdIelnnGrHSPBGdIl5p6wgAAD5AMBICPDLXSOE5R7bMEclO80DpAABwjWAkxFpHyjWiRjXuJVKjhUhOWkFAAgCAlxGMhFjeyMHkcraMaL5J52EF9+mqAQD4QMgHI1999ZW0bNlSWrRoIf/5z38kVHlkRI1Vh5tEIiuIHNwgcmRL+Y8HAEC4BiN5eXkybtw4WbZsmfz4448yZcoUOXnypIQij0x8ZlW5lkirKwvu0zoCAPCykA5G1q1bJ23btpX69etL5cqVZdCgQbJo0SIJRR6Z+MxZIuvm2SK5HjomAADBFoysXLlSBg8eLElJSWbujHnz5hV5zYwZM6RJkyYSFxcnPXr0MAGI1aFDh0wgYqX3Dx48KKHeTeORRcbOu1wksZFIVorIts/LfzwAAIIxGElPT5cOHTqYgMOZ2bNnm26YSZMmycaNG81rBwwYEJbTSFtbRtKy8yQl0wOzVkZGinQ6l8jK4nkAgHANRrRb5amnnpIhQ4Y4ff7FF1+UkSNHyu233y5t2rSRmTNnmoW13nrrLfO8tqjYt4Tofd1WnOzsbElNTXW4BdNcIzUre2iuEatOt4hERIrs+5/Iid2eOSYAAMEUjLiSk5MjGzZskL59+9q2RUZGmsdr1qwxj7t37y5btmwxQUhaWposWLDAtJwUZ/LkyVK1alXbrWHDhhK2SayqSpJIk94F91dOETm06Y9b8v6S99fX2O9Tmn0BAGEjaNemOXHihOTn50udOnUctuvjHTt2mPvR0dHywgsvyOWXXy5nz56V8ePHS40aNYo95oQJE0y3j5W2jARTQKLByKb9yZ5rGdGg4ffVfySy6s0qOlbkng0iiQ2L33d6F5G87KLPlbQvACCsBG0w4q4//elP5uaO2NhYcwtWHp1rRGWcFDlbTP6JBhn6fHEBhT7nLBBxZ18AQFgJ2mCkZs2aZkXPo0ePOmzXx3Xr1pVwVN/Tw3tL8vFwkQoF71kEw4EBAKGeMxITEyNdunSRpUuX2rZpV4w+7tmzp4Qjj+eMlCT5N5Hj253f9DlXdnwtknHKN+UEAAS0gG4Z0aTTPXv22B7v3btXNm3aJNWrV5dGjRqZ/I7hw4dL165dTbLq1KlTzXBgHV1THjqUWG+akxJMGlrXpzk314jOzeJVV78kUqO58+dO7hH56oHi9135nMiqF0TOu0yk7RCRVleJVKzmmHOiXTmFxdegewcAQkxAByM//PCDST61siaXagAya9YsueGGG+T48eMyceJEOXLkiHTs2FEWLlxYJKm1tMaMGWNumsCqo2qCRf3EgpyRM9l5kpqZJ1XjK3j3DZM6iyR1dP5cbBXX+1ZvJnLqF5E9SwpuX97/R2CS1EHkjSvKnvxankAmGPf153uzb2jXFeAjAR2MXHbZZSXOJnrPPfeYG0QqxuhcIzFyIi1H9p/OkKrx5Qyk9I+VXvyLCwr0+bLuO+xzkbwska3zRLbOFTm2VWTP4oJbRLSIJa9sya/lGcUTjPsGa7nZN/DrCvChgA5GUHr1q8WbYESTWNvVL2cwon+k9I9VWb5VubvvpQ8W3I7vEtmmgcm8gsDEle9eFKnWSCSmskhMJZEK8efux4ucOeJ6FI92H0XHiZgurIhzP885/ZvrfVP2i1RMdP68Pudy3wOO3VAO+x4o+77l3Z99w7uuGNWGABFh8chCJqHJ2k2TkpIiVaqU0O0QIMZ8sFHmbz4sj1zVWu665DwJStu+EPn4Nn+XAgh9o1YU39UK+PAaSstICCWwemX1Xn/QBfpc6TxCJC5BJCfd8ZabIZJ+oiAXBQAQNAhGQiiB1SsTnwWirrcX/21Op5t//dKSvwlqg6BpFDzXMKj3D28S+U+f4ve9Y7FIvQucP3d4s8hb/Vzsu6iEffuXbd/y7s++4V1XQIAgGAkxPp9rxBvKkzjrLs0VKTz0ObKEfw7RMSIV4op/zuW+scVPEKfPlXXf8u7PvuFdV0CAIBgJMfZzjQSt8iTOenMEUCDuG6zlZt/AryvAh0hgDbEE1sycfGk9caG5/9Ok/lK1opfnGglEwTifA/OMBEddh9u+gI+uoQQjHqrIQNLlycVyMj1H5t93sbRNCq6cFwBA+F1Dg3ZtGoT4iBoAQNggGHFCh/W2adNGunXrJsEoLEbUAABCBsGIEzqsd9u2bbJ+/XoJRiExogYAEDYIRkIQ3TQAgGBCMBKC6KYBAAQTgpEQRDcNACCYEIyEoPrnckbOZOVJSmauv4sDAIBLBCMhKD4mWmpUKpienCRWAECgIxgJwaG99l01QT0tPAAgLBCMhODQXkUSKwAgWBCMhHjeCBOfAQACHcFIiGJEDQAgWBCMhCgmPgMABAuCkRD1R84IU8IDAAIbwUiIqp9YkDOSylwjAIAARzASoirFRkv1c3ONMLwXABDICEZCGEmsAIBgQDASopOeKZJYAQDBgGAkRCc9U0x8BgAIBgQjIcw2JXwyI2oAAIGLYCSE0U0DAAgGBCMhjG4aAEAwIBgJg7lGUjJzJTUr19/FAQDAKYKREJ9rpFp8BXOfuUYAAIGKYCTE0VUDAAh0BCMhjonPAACBjmAkxDGiBgAQ6AhGQngGVsXqvQCAQEcwEsIzsCpaRgAAgY5gJMRZW0YOJmf6uygAADhFMBLi6p+bEj45I1fOMNcIACAAEYyEuMr2c43QOgIACEAEI2HAlsR6iq4aAEDgIRgJo2nhD5xm9V4AQOAhGAkDjKgBAAQygpEwQDACAAhkBCPhlDOSTDcNACDwEIyEgQbVrTkjJLACAAIPwUgYJbAy1wgAIBARjISBhLgKkshcIwCAAEUwEmZJrAfpqgEABBiCkRBftdeqQeK5JFaCEQBAgImwWCwWfxciUKWmpkrVqlUlJSVFqlSpIsHszlnrZemOYxIVESEt6lSW+/u2kIHt6vm7WACAEFWaaygtI2Fg4ZbDJhBR+RaL7DxyRka/t9FsBwDA3whGwsDUJbslwu6xNoVFRIi8vHS3H0sFAEABgpEwsPdEuglA7Gnn3K/H0/1UIgAA/kAwEgaa1qzk0DIi51pGzqtVyU8lAgDgDwQjYUCTVe1bRiLOtYyM7XO+H0sFAEABgpEwoKNmZt7aWeIqFPy661erKDNv7SID29X1d9EAACAYCaeApE/rOub+iF5NCEQAAAGDlpEw0qh6wcRn+0+xei8AIHAQjISRhtUKgpF9BCMAgABCMBJGGlYvWJ9mP1PCAwACCMFIGLaMHDidIawCAAAI6mDknXfekfnz59sejx8/XhITE6VXr17y+++/e7J88KCkxIpmfpGs3LNyPC2bugUABG8w8swzz0jFigVN/mvWrDGr3D7//PNSs2ZNeeCBBzxdRnhITHSkJFU911VzKpN6BQAEhOiy7LR//35p3ry5uT9v3jy57rrrZNSoUXLRRRfJZZdd5ukywoMaVKsoB5MzTVdNl8bVqFsAQHC2jFSuXFlOnjxp7i9atEj69etn7sfFxUlmJt+4A1nDc8N7951keC8AIIhbRjT4uOuuu6RTp06ya9cuufLKK832rVu3SpMmTTxdRnghiXX/aYIRAEAQt4xojkjPnj3l+PHjMmfOHKlRo4bZvmHDBrnppps8XUZ4Y3gvOSMAgGBuGdGRM9OnTy+y/fHHH/dEmeCLWVhpGQEABHPLyMKFC2XVqlUOLSUdO3aUm2++WU6fPi3BTj9PmzZtpFu3bhKqOSOHU7IkN/+sv4sDAEDZgpEHH3xQUlNTzf2ff/5Z/u///s/kjezdu1fGjRsX9NU6ZswY2bZtm6xfv15CTa3KsWaIb/5ZixxOzvJ3cQAAKFs3jQYd2nKgNGfk6quvNnOPbNy40ZbMisAUGRlhhvf+ejzddNU0qlHQUgIAQFC1jMTExEhGRsFojCVLlkj//v3N/erVq9taTBAEI2pYMA8AEKwtIxdffLHpjtFJztatWyezZ88223WYb4MGDTxdRnhtwTyG9wIAgrRlREfSREdHy6effiqvvfaa1K9f32xfsGCBDBw40NNlhLdG1DC8FwAQrC0jjRo1kq+++qrI9pdeeskTZYKPumn20U0DAAjWYETl5+ebdWm2b99uHrdt21b+9Kc/SVRUlCfLBy8O79X1aQAACMpgZM+ePWbUzMGDB6Vly5Zm2+TJk6Vhw4Yyf/58adasmafLCS+0jJxIy5GMnDyJjylzTAoAgH9yRu677z4TcOjqvTqcV2/79u2Tpk2bmucQ2KrGV5CEuIIA5MBpFjYEAPhXmb4Sr1ixQr7//nszlNdK16d59tlnzQgbBEcS69ZDqWZ47/l1EvxdHABAGCtTy0hsbKycOXOmyPa0tDQzBwkCH0msAICgDkZ0xtVRo0bJ2rVrxWKxmJu2lIwePdoksSLwsXovACCog5Fp06aZnJGePXtKXFycufXq1UuaN28uU6dO9Xwp4bURNUx8BgAIypyRxMRE+fzzz82oGuvQ3tatW5tgBMGBKeEBAEEXjJS0Gu/y5ctt91988cXylQo+66bR0TTazRYREUGtAwACOxj58ccf3XodF7Xg0ODcXCNp2XmSnJEr1SqReAwACPBgxL7lA8EvrkKU1E6IlWNnss208AQjAICgSmBFaCCJFQAQCAhGwljDagV5I6zeCwDwJ4KRMEbLCAAgEBCMhDFbMHKK1XsBAP5DMBLGmGsEABAICEbCmHWukYPJmZJ/1uLv4gAAwhTBSBirV7WiREdGSG6+RY6mZvm7OACAMEUwEsaiIiMkKdE6ooa8EQCAfxCMhLlGtgXzMv1dFABAmAqLYGTIkCFSrVo1+fOf/+zvogRs3ojOwgoAgD+ERTAyduxYeffdd/1djIBeo+YAwQgAwE/CIhi57LLLJCEhwd/FCEhMfAYAkHAPRlauXCmDBw+WpKQks+LvvHnzirxmxowZ0qRJE4mLi5MePXrIunXr/FLWUMSU8AAACfdgJD09XTp06GACDmdmz54t48aNk0mTJsnGjRvNawcMGCDHjh2zvaZjx47Srl27IrdDhw758JMEd8vI0TNZkpWb7+/iAADCULS/CzBo0CBzK86LL74oI0eOlNtvv908njlzpsyfP1/eeustefjhh822TZs2eaQs2dnZ5maVmpoqoa5GpRiJj4mSjJx8M/lZs1qV/V0kAECY8XvLiCs5OTmyYcMG6du3r21bZGSkebxmzRqPv9/kyZOlatWqtlvDhg0l1GnXGNPCAwD8KaCDkRMnTkh+fr7UqVPHYbs+PnLkiNvH0eDl+uuvl6+//loaNGhQbCAzYcIESUlJsd32798v4TS8l7lGAABh2U3jC0uWLHHrdbGxseYWbhjeCwDwp4BuGalZs6ZERUXJ0aNHHbbr47p16/qtXKGG4b0AAH8K6GAkJiZGunTpIkuXLrVtO3v2rHncs2dPv5YtFKeEZxZWAEBYdtOkpaXJnj17bI/37t1rRsdUr15dGjVqZIb1Dh8+XLp27Srdu3eXqVOnmuHA1tE13qDDjPWm+SphlTNyivVpAAC+F2GxWCziR99++61cfvnlRbZrADJr1ixzf/r06TJlyhSTtKpzikybNs1MfuZtOrRXR9VoMmuVKlUkVKVn50nbSd+Y+5sf6y9V4ir4u0gAgCBXmmuo34ORQBYuwYjq/ORiOZWeI/Pvu1jaJlX1d3EAAGF0DQ3onBH4DtPCAwD8hWAEDiNqDpzOoEYAAD5FMAKHYIQRNQAAXyMYcUJH0rRp00a6desm4YIp4QEA/kIw4sSYMWNk27Ztsn79egkXTAkPAPAXghE4tIxozggDrAAAvkQwAiMpsaJERIhk5Z6V42nZ1AoAwGcIRmDEREdKUlXrTKyMqAEA+A7BCGwaVGNaeACA7xGMOBGOo2kcVu+lZQQA4EMEI06E42gah+G9THwGAPAhghHYsHovAMAfCEZg04hZWAEAfkAwgiI5I4dTMiU3/yw1AwDwCYIR2NSqHGuG+J61iBxOzqJmAAA+QTCCP06GyIg/hveSxAoA8BGCEThgwTwAgK8RjDgRrvOM2I+o2cdcIwAAHyEYcSJc5xmxH1Gz/3Smv4sCAAgTBCNwQDcNAMDXCEbgdHjvARJYAQA+QjACpy0jJ9JyJCMnj9oBAHgdwQgcVI2vIAlx0eb+AfJGAAA+QDCC4qeFP5lB7QAAvI5gBEWwei8AwJcIRpwI53lGFKv3AgB8iWDEiXCeZ8R+RA1TwgMAfIFgBEUw1wgAwJcIRlB8y8ipDLFYLNQQAMCrCEZQhHXl3vScfDmdkUsNAQC8imAERcRViJLaCbG21hEAALyJYAROkcQKAPAVghE41fBcV83+U6zeCwDwLoIRuGwZ2Uc3DQDAywhG4BSr9wIAfIVgBE4x1wgAwFcIRpwI9+ng7aeEP5icKflnmWsEAOA9BCNOhPt08Kpe1YoSHRkhufkWOZqa5e/iAABCGMEInIqKjJCkxILWEZJYAQDeRDCCYjWymxYeAABvIRhBiXkj+08z1wgAwHsIRlCsBtUKWkYO0DICAPAighEUiynhAQC+QDCCEqeEJ4EVAOBNBCMoMYH1aGq2ZOXmU1MAAK8gGEGxqleKkfiYKNvkZwAAeAPBCIoVERHBtPAAAK8jGIFLDO8FAHgbwQhcYngvAMDbCEbg1vBeRtQAALyFYMQJVu11MiX8aaaEBwB4B8GIE6za6yRn5BSjaQAA3kEwApcanpsSPiUzV1KzcqktAIDHEYzApUqx0Wa+EcXqvQAAbyAYgdvTwhOMAAC8gWAE7i+YR94IAMALCEZQotz8s+bn5AXbZeDUlbJwy2FqDQDgMQQjcEkDj2+2HjX3z1pEdh45I6Pf20hAAgDwGIIRuDR1yW6JsHtsMWvWiLy8dDc1BwDwCIIRuLT3RLoJQOxZLCK/Hk+n5gAAHkEwApea1qzk0DJiVSshlpoDAHgEwQhcur9vC1vXjL0DpzNl0udbJDsvnxoEAJQLwQhcGtiunsy8tbO0qpsgsdGR5mf/NnXMc++s+V2ue+1/8tsJumwAAGUXYbFoBgCcSU1NlapVq0pKSopUqVKFSrKzfMcxGffxJjmdkSuVY6PlmaHt5U8dkqgjAECpr6G0jKBMLm9VW74ee4l0b1Jd0rLz5L4Pf5QJn/0sWbl02wAASodgBGVWr2pF+WBkD7n3iuYmp+TDdfvk2hmrZc+xNGoVAOA2ghGUS3RUpPxf/5by3zt6SM3KsbLjyBkZ/MoqmbPhADULAHALwQg84uIWNeXrsRfLRc1rSGZuvvzfJz/J/338k2Tk5FHDAACXSGB1gQTW0ss/a5FXl++Rl5bsMtPH16kSK/ExUXIoOcvMWaJDhXWEDgAgtKWSwFo+M2bMkDZt2ki3bt089CsJH1GREXJvnxbywcgLpWrFaDmami17T2RIdt5Z1rUBADhFN40TY8aMkW3btsn69eud1xpKdOF5NaR2QpzDNjN5mrCuDQDAEcEIvGbfqYwi2zQg0STXHUdSqXkAgEEwAp+va6PT7F358nfy0Keb5WhqFr8BAAhzBCPw2bo21p9dGlUzya2zf9gvl035Vl5avEvSsxl1AwDhimAEPl3XZuatXWTO33rJnLt7SedGiWYY8MtLd8tl//pWPlq3z4zGAQCEF4b2usDQXu/SZZEWbDkizy7YYcsvaVknQSZc2UouPb+WRBReKhgAEJLXUIIRD1Ukyi4n76z89/vfZdrS3ZKSmWu2XdKipvRuUUvmbDwge0+kM0cJAAQZghE/VCTKLyUjV6Yv3y3v/O93yck/6/CctpFoB452+zBpGgAEPiY9Q1CqGl9B/nlVG1ky7lJJiIt2PkfJkt1+Kx8AwDtIYEXAaVQj3nTdFDdHybtrfpMzWQXdOQCA4EcwguCao0REJn6+VXo8s1QmfPazbD2U4ofSAQA8iWAEQTVHyV+6NpBmtSpJRk6+fLhun1w1bZUMfXW1fLbxgGTl5vu1zACAsmE0jQsksPrXwi2HzRwkvx5Pl/NqVZKxfc6Xge3qmiHB3/96St5b+7t8s+WI5J2bm6RafAW5vmtDuaVHI9l+OFWmLtnNSBwA8BNG0/ihIuEfx85kycfr98sHa/fJoZQspyNwGIkDAL7HaBqEDV0Z+J4rWsh3D10hbwzraiZLs7LO5WoNSLSlBAAQeMgZQUiIioyQfm3qyDt3dJeYqMhiR+Lc9c4PMnv9PtOiAgAIDI6TOQAhQPNLdh45Y2sZsbdk+1FzUx0aJkrfVrWlT+s60rpeAtPPA4CfkMDqAjkjwZv4Ovq9jWYEjsUitp//vKqVZOaclaXbj8pPBxyHBNdPrCh9WhcEJqmZOTJj+S8kvwJAOZDA6iEEI6E3EsfqWGqWLN1xzAQmq/ackKzcopOsKZJfAaBsCEY8hGAkPOj8JKv3nJAl24/JJz/stw0VtteqboIsvL+3X8oHAMGI0TRAKcRViDLdM5OHtjeJsM5o8uvrK3+R9Ow86hYAPIzRNIAb09CrZ77eIRc/t0ymL9stqayNAwAeQzACuDEN/bCejU2gcjojV/61aJdc/OwyeXHxLknOyKH+AKCcGE3jAjkj4am45Ne8/LMy/+fDMn3ZHtl9LM28tlJMlNzWs4ncdUlTqVk51t9FB4CAQQKrHyoS4ePsWYt8s/WITFu2x6yBo+IqRMotPRpLi9qVZdb/fmNYMICwl1qKaygtIy4QjMAVXbBv6fZj8sqy3UXmLTH/uM7N/Drz1s4ysF09KhNAWEktRTDCDKxAGUVEREjfNnXMZGnf7T4hf/3vD5JpN1+JdYDwxM+3SqXYaGmbVFWqV4qhvgEg3IKR/fv3y2233SbHjh2T6OhoefTRR+X666/3d7EQYkFJ7/NriZPpSYxjZ7LltjfXmft1q8RJ26Qq5tYmqar52aBaRXMMzVXRxfz2nkg3ybKaTEuLCoBwEPLdNIcPH5ajR49Kx44d5ciRI9KlSxfZtWuXVKpUqcR96aZBaQycutLpmjgJcdFSo1KM/HYyw+l+VStWkDpVYmXX0YKk2LJ28RDMAAgk5Iy40KFDB/nqq6+kYcOGHq1IoLg1cWbe2sWMxjmTlSvbD5+RrYdSZOuhVHPbffSM0xlfrXQOtgbV4qVafAVJjI+x+xkjieZ+BXNfE2knL9hhC2JKG8yUJ5Bh39Cuq2AtN/u28PsXmaAKRlauXClTpkyRDRs2mFaMuXPnyrXXXuvwmhkzZpjXaMuGBhOvvPKKdO/evdTvpe8xfPhw2bJli1uvJxiBp9fEKSw7L192H02Ta2esdhmUlFV0ZIQ0qVlJ4mOipGKFKPMzPibazDpbcD9KDiZnylebDxfZ97YLG0v7BlVNYKPdSPozMlIDnQjb/CubD6TIm6v2FgmCRl7SVDo1quZwvMKTyf2477S8/l3Rff/au+i+4mTff69kX2/XFXUdvuflTA8k3gdVMLJgwQJZvXq16T4ZOnRokWBk9uzZMmzYMJk5c6b06NFDpk6dKp988ons3LlTateubV6jXTB5eUWn6V60aJEkJSWZ+6dOnZJLLrlE3njjDenVq5fTsmRnZ5ubfUVqCwotI/BHF4/+UdCA5rnrLjCTrZ3OyJGUcz/1sU64pveTM3LNdPUA4An6ZUPX41owtnf4BCP29NtX4WBEA5Bu3brJ9OnTzeOzZ8+aAOHee++Vhx9+2K3jaoDRr18/GTlypElmLc5jjz0mjz/+eJHtBCPwdxdPWYOZxjXi5Zmh7SUzJ18ycvLP/cyTjNyC+3p7a/Vep8m32kV06bnEXH3a+qdCf1j0P4vIml9Pmp+Fafm7Nan+xwYnr1n/2ylnm025uzZx/W3uh99OF7tvl8au993wO/u6W1fUV/ieW7HRkbLzqUFSHiEztDcnJ8d0rUyYMMG2LTIyUvr27Str1qxx6xj6B3TEiBFyxRVXuAxElL7PuHHjirSMAN6mzaHaLFqaLh572sfrLJh5eFBr6dWspst9V+05UTSQiRBpWTdB3r69e+mDoHPfqj7+a88y7/vJ6F5l3vfTu9nXU3VFXYfvuXVerZIHeYTN2jQnTpyQ/Px8qVOnjsN2faz5I+7QLiDt6pk3b57pztHbzz//7PS1sbGxJnqzvwG+DEi0WVS/jehPdwMR+2BG//joNxr96W6rirP1eDSQ0WCIfcO7roK13OwrPjs/PCWgu2kOHTok9evXl//973/Ss+cf37LGjx8vK1askLVr13q1PCSwIlyUNvGWfcOnroK13Ox7vs/Oj5DPGdFumvj4ePn0008d8kh0RExycrJ8/vnnXi0PwQgAAN6/hgZ0N01MTIwZZbN06VLbNk1g1cf2LSUAACB4+T2BNS0tTfbs2WN7vHfvXtm0aZNUr15dGjVqZBJKtSWka9euZm4RHdqbnp4ut99+u9fKpPOa6E3zVQAAgHf5vZvm22+/lcsvv7zIdg1AZs2aZe7rsF7rpGeagDpt2jQz5Nfb6KYBACDMckYCDcEIAABhnjMCAABCH8EIAADwK4IRJzR5tU2bNmYaegAA4F3kjLhAzggAAGG+No2/WXN7tUIBAID7rNdOd8bJEIy4cOZMwbLsLJYHAEDZr6XaQuIK3TQu6Gyvuj5OQkKCmareU6yrAe/fv5/F+Kgrj+Lcoq68hXOLuiotbRHRQCQpKUkiI12nqNIy4oJWXoMGDcRbWBmYuuLc8j/+HVJfnFveU1KLiBWjaQAAgF8RjAAAAL8iGPGD2NhYmTRpkvkJ6opzyz/4d0h9cW4FDhJYAQCAX9EyAgAA/IpgBAAA+BXBCAAA8CuCEQAA4FcEI35YEbhJkyYSFxcnPXr0kHXr1vm6CEHhscceM7Pe2t9atWrl72IFhJUrV8rgwYPNrIZaL/PmzSsy6+HEiROlXr16UrFiRenbt6/s3r1bwlVJ9TVixIgi59rAgQMlHE2ePNmsVq6zTteuXVuuvfZa2blzp8NrsrKyZMyYMVKjRg2pXLmyXHfddXL06FEJN+7U1WWXXVbk3Bo9erTfyhzICEZ8aPbs2TJu3DgzrHfjxo3SoUMHGTBggBw7dsyXxQgabdu2lcOHD9tuq1at8neRAkJ6ero5dzSwdeb555+XadOmycyZM2Xt2rVSqVIlc57pRSQclVRfSoMP+3Ptww8/lHC0YsUKE2h8//33snjxYsnNzZX+/fubOrR64IEH5Msvv5RPPvnEvF6XzBg6dKiEG3fqSo0cOdLh3NJ/n3DCAp/p3r27ZcyYMbbH+fn5lqSkJMvkyZP5LRQyadIkS4cOHaiXEug/4blz59oenz171lK3bl3LlClTbNuSk5MtsbGxlg8//DDs67Nwfanhw4dbrrnmmrCvG2eOHTtm6mzFihW2c6lChQqWTz75xPaa7du3m9esWbMmrOuwcF2pSy+91DJ27Fi/litY0DLiIzk5ObJhwwbTZG6/9o0+XrNmja+KEVS0a0Gb1s877zy55ZZbZN++ff4uUsDbu3evHDlyxOE807UhtEuQ86x43377rWlqb9mypdx9991y8uRJn/y+Al1KSor5Wb16dfNT/4ZpC4D9+aXdp40aNQr786twXVm9//77UrNmTWnXrp1MmDBBMjIyfPo7DBYslOcjJ06ckPz8fKlTp47Ddn28Y8cOXxUjaOjFc9asWebioE2bjz/+uFxyySWyZcsW00cL5zQQUc7OM+tzKNpFo90MTZs2lV9++UX+8Y9/yKBBg8zFNSoqKqxXLb///vvloosuMhdSpedQTEyMJCYmOrw23M8vZ3Wlbr75ZmncuLH5UrV582Z56KGHTF7JZ5995tfyBiKCEQQkvRhYXXDBBSY40X/UH3/8sdx5551+LRtCy4033mi73759e3O+NWvWzLSW9OnTR8KV5kNo8E+uVtnratSoUQ7nliaV6zmlQa+eY/gD3TQ+os10+i2rcNa5Pq5bt66vihG09JvY+eefL3v27PF3UQKa9VziPCs77RbUf6/hfK7dc8898tVXX8ny5culQYMGDueXdjknJyc7vD6c/44VV1fO6JcqFc7nVnEIRnxEmza7dOkiS5cudWja08c9e/b0VTGCVlpamvk2od8sUDztatCLgv15lpqaakbVcJ6558CBAyZnJBzPNc2J1ovr3LlzZdmyZeZ8sqd/wypUqOBwfmm3g+Zzhdv5VVJdObNp0ybzMxzPrZLQTeNDOqx3+PDh0rVrV+nevbtMnTrVDAO7/fbbfVmMoPD3v//dzA2hXTM6dFCHQ2vL0k033SThTgMz+29WmrSqf+Q0cU4TCbXv+qmnnpIWLVqYP5CPPvqo6bPWeRDCkav60pvmI+lcGRrEacA7fvx4ad68uRkOHY7dDR988IF8/vnnJjfLmgeiSdA6Z43+1G5S/VumdVelShW59957TSBy4YUXSjgpqa70XNLnr7zySjMni+aM6LDo3r17m65AFOLv4Tzh5pVXXrE0atTIEhMTY4b6fv/99/4uUkC64YYbLPXq1TP1VL9+ffN4z549/i5WQFi+fLkZQlj4pkNUrcN7H330UUudOnXMkN4+ffpYdu7caQlXruorIyPD0r9/f0utWrXMkNXGjRtbRo4caTly5IglHDmrJ729/fbbttdkZmZa/va3v1mqVatmiY+PtwwZMsRy+PBhS7gpqa727dtn6d27t6V69erm32Hz5s0tDz74oCUlJcXfRQ9IEfq/wgEKAACAr5AzAgAA/IpgBAAA+BXBCAAA8CuCEQAA4FcEIwAAwK8IRgAAgF8RjAAAAL8iGAHCiC7+FhERUWRtEQSOJk2amNmZgXBCMAJ4yIgRI8J2ynVv1KUGTaNHj3Y6Dbc+p68JB/pZ582b5/djAN5EMAKg3HQlV09r2LChfPTRR5KZmWnblpWVZdb70DV4yis3N1f8JT8/3yyUCaAAwQjgI1u2bJFBgwZJ5cqVpU6dOnLbbbfJiRMnzHOvv/66Wcyu8AXqmmuukTvuuMP2WBfl6ty5s8TFxZml7nWRt7y8PIdvwP/5z39kyJAhEh8fbxbL++KLL4qUZfXq1WaxLj2OLnCmZbO3atUqueSSS8yCXxoU3HfffWZRR/uuhCeffFKGDRtmFksbNWqU2f7GG2+Y1+t7axlefPFFSUxMLFN96efUY3322We2bXpfA5FOnTo5vHbhwoVy8cUXm/fSRcmuvvpqs1CZ1W+//WbqZvbs2XLppZeaz/3++++boEAXfbPup4vk6WKW9i1czrpNOnbsKI899pjtsX7O9u3bS6VKlUyZ//a3v5kF+qxmzZpl3kN/F23atJHY2Fiz0u2xY8fMgpBaz7qooZbJnr630rrU8lsfOwsGdQVZXQ1WP5suMDl58uQSj+HO+fTaa6+Z81bLqK/59NNP3fjtAaXk78VxgFChC69dc801Tp87ffq0WYxtwoQJlu3bt1s2btxo6devn+Xyyy83z586dcosCrhkyRLbPidPnnTYtnLlSkuVKlUss2bNsvzyyy+WRYsWWZo0aWJ57LHHbPvoP+kGDRpYPvjgA8vu3bst9913n6Vy5crmWPaLxrVu3drsv3nzZsvVV19tjpOTk2NeowsSVqpUyfLSSy9Zdu3aZVm9erWlU6dOlhEjRtjeRxeU07L861//Mq/X26pVqyyRkZGWKVOmmIX5ZsyYYRYJq1q1apnr8sUXXzQL/VnpfS2XPmddGFB9+umnljlz5pjP/OOPP1oGDx5sad++vSU/P988v3fvXvO59XPq63799VfLoUOHLM8995xZ8E23bdu2zXLnnXdaEhISHH6P+ln1Pe116NDBMmnSJNtjfX7ZsmXmfZYuXWpp2bKl5e6777Y9r4un6UJ8vXr1MvW5Y8cOS3p6umXQoEHmWGvWrLH88MMP5vmKFSva3u/YsWO2xdd0MTp97IzWecOGDc058ttvv1m+++47cw64Ooa751ONGjUsb7zxhvmdPvLII5aoqChTV4AnEYwAPghGnnzySbM6rL39+/ebP/bWFXV13zvuuMP2/L///W9LUlKS7YKqF+JnnnnG4Rj//e9/zerGtn/QIuaCYZWWlma2LViwwCEY+eijj2yv0UBFL4CzZ882j/WCPGrUKIf30YubBhq6Yqv1An3ttdc6vEZXVr7qqqsctt1yyy3lCkb0wqkrnuoFVm9xcXGW48ePFwlGCtPX6Of8+eefHYKRqVOnOrxO6+7555+3Pc7NzTXBXGmDkcI++eQTcxG30kBA33/Tpk22bfp7123r1q2zbdNAVbfZv58+njt3rsv6uvfeey1XXHGFWbHZGWfHcPd8Gj16tMNrevTo4RBoAZ5ANw3gAz/99JMsX77cdNFYb61atTLPWbsTbrnlFpkzZ45kZ2ebx9pkf+ONN0pkZKTtGE888YTDMUaOHCmHDx+WjIwM23tp94uVdhtoN4p2B9jr2bOn7X716tWlZcuWsn37dtv7aLeC/fsMGDDAdCHt3bvXtl/Xrl0djrlz507p3r27wzb7x9otYX/MZ555Rr777juHbYW7KWrVqiVXXXWVKc/bb79t7tesWbNI/e7evVtuuukm042gn9faFaHvac++zCkpKabuevToYdsWHR1d5HO5Y8mSJdKnTx+pX7++JCQkmC64kydPOvxeYmJiHH43Wt/6fl26dLFt03OiLN1amsy7adMm83vULrVFixaVuI+755P9uWJ9bD1XAE+J9tiRABRL8wc0N+C5554r8pz28yt9Xr+Mzp8/X7p162Yu1C+99JLDMbRPf+jQoUWOoX3+VhUqVHB4Tvv9S5Msqe/z17/+1VzUCrNPHNVApzQ0J0YvmPZBkOYh2G/TXJrCNGdG8yHUjBkznB5b607zJDRnxZp7065duyKJtaUts9JgsKCRwHnyq+ajaI7K3XffLU8//bT5XJpzc+edd5r31/wZpZ9VfxfeoHkfGiguWLDABEZ/+ctfpG/fvi7zO9w9nwBfIBgBfEAvFtrqod/Y9duwM3oB0AuDtg7s2bPHfMvV/eyPoa0PzZs3L3d5vv/+e1tgcfr0adm1a5e0bt3a9j7btm0r9ftoedevX++wzf6xfm5nxyzpfQYOHGgu6noh1xaawrQFQutFAxFNulUaDJSkatWqJhBcu3at9O7d22zT5M0NGzY41Lu2zmhrgVVqaqpDC5G+XoOfF154wdaK9fHHH5f4/toKYn0/DT6Vfo7Cc8BocKmJtiXRFqEbbrjB3P785z+bejt16pQJjpwdw93zSc8VTVS2f1w4gRgoL4IRwIO06d/+m77SURo6N4ZeLLUrQUds6AVCAw4duqqjX6KiomxdNfote+vWrXLrrbc6HGfixInmOQ0i9GKjFz5tateRME899VSpyqnN81oubYn45z//abo+rCNIHnroITPCRlsj7rrrLtOaoMHJ4sWLZfr06cUe89577zUXdR1Zoi0Vy5YtM9/Uy9saoHVj7Raw1pO9atWqmc+iI5I0uNCumYcfftitY48dO1aeffZZM+pIgwMte+Fg4IorrjDdRPqZtAtFfw/25dCLubaUvPLKK+Y1OlJp5syZbgVvGjBoK5SOWNFg7f777zctKPY0gF26dKlcdNFFZhSOft7CtNz62TVI0PPik08+kbp169q6fJwdw93zSY+lXVc6WkkD5XXr1smbb77pVv0CbvNI5gkAk1Cp/6QK3zQhVOnIlCFDhlgSExNNwmirVq0s999/v0PSoSaragKh7qcjHApbuHChbcSFjoTo3r275fXXX3eZqKgJpJpAaZ/A+uWXX1ratm1rRuvoMX766SeHfTSpUkf76EgcHVlzwQUXWJ5++mmXSZ1Ky1K/fn1TPk1wfeqppyx169b1aDKwKpzAunjxYjNCSJNdtazffvutQ11YE1h1pI09TVgdO3asqUv9vYwbN84ybNgwh/dOSUkxybn6Gh2xoqNPCiew6qgf/b3p5x4wYIDl3XffNe+no6iU1r+zRF4d3aJJv1ruRo0amf0K1+0XX3xhad68uSU6Oto854zWe8eOHc3vSsupyak6YqukY7hzPumoKD0XtIw62saa6Ax4UoT+z/3QBQDcpwmRO3bsMPkvwUKTQbV1hBlLC/KN5s6dy8zC8Dq6aQB4zL/+9S/p16+f6drRLpp33nlHXn31VWoYgEsEIwA8RvMJnn/+eTlz5owZZjtt2jSTdwIArtBNAwAA/IpJzwAAgF8RjAAAAL8iGAEAAH5FMAIAAPyKYAQAAPgVwQgAAPArghEAAOBXBCMAAMCvCEYAAID40/8DMi1KLIU7DQ4AAAAASUVORK5CYII=", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "def joint_residuals(params, args):\n", + " rates, amplitudes = params\n", + " operator = lx.MatrixLinearOperator(features(rates))\n", + " return operator.mv(amplitudes) - y\n", + "\n", + "\n", + "def loss_history(residual_fn, y0, max_steps=40):\n", + " # Step the solver manually (the standard Optimistix pattern) to record the loss\n", + " # after each Levenberg-Marquardt step. Optimistix functions return (output, aux).\n", + " fn = lambda z, args: (residual_fn(z, args), None)\n", + " args, options, tags = None, {}, frozenset()\n", + " f_struct, aux_struct = jax.eval_shape(lambda: fn(y0, args))\n", + " step = eqx.filter_jit(\n", + " eqx.Partial(solver.step, fn=fn, args=args, options=options, tags=tags)\n", + " )\n", + " terminate = eqx.filter_jit(\n", + " eqx.Partial(solver.terminate, fn=fn, args=args, options=options, tags=tags)\n", + " )\n", + "\n", + " y = y0\n", + " state = solver.init(fn, y, args, options, f_struct, aux_struct, tags)\n", + " losses = [0.5 * jnp.sum(fn(y, args)[0] ** 2)]\n", + " done, _ = terminate(y=y, state=state)\n", + " while not done and len(losses) <= max_steps:\n", + " y, state, _ = step(y=y, state=state)\n", + " losses.append(0.5 * jnp.sum(fn(y, args)[0] ** 2))\n", + " done, _ = terminate(y=y, state=state)\n", + " return jnp.array(losses)\n", + "\n", + "\n", + "# The joint fit needs an initial amplitude guess; give it the best amplitudes for\n", + "# `init_rates` (a single linear solve) so the comparison is fair.\n", + "init_operator = lx.MatrixLinearOperator(features(init_rates))\n", + "init_amplitudes = lx.linear_solve(init_operator, y, lx.SVD()).value\n", + "\n", + "varpro_loss = loss_history(varpro_residuals, init_rates)\n", + "joint_loss = loss_history(joint_residuals, (init_rates, init_amplitudes))\n", + "\n", + "plt.figure(figsize=(6, 4))\n", + "plt.semilogy(varpro_loss, \"-o\", markersize=4, label=\"variable projection ($k$ only)\")\n", + "plt.semilogy(joint_loss, \"-s\", markersize=4, label=\"joint fit ($k$ and $c$)\")\n", + "plt.xlabel(\"Levenberg--Marquardt step\")\n", + "plt.ylabel(\"loss\")\n", + "plt.legend()\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "4da5c7ec", + "metadata": {}, + "source": [ + "Variable projection reaches the noise floor in a handful of steps, while the joint fit — with twice as many parameters and a more ill-conditioned landscape — lags well behind. This conditioning advantage also tends to widen the basin of convergence, making the fit less sensitive to the initial guess; see [Golub & Pereyra (1973)](https://doi.org/10.1137/0710036) and [O'Leary & Rust (2013)](https://doi.org/10.1007/s10589-012-9492-9)." + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "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.13.5" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/lineax/__init__.py b/lineax/__init__.py index 4f370dee..cb7dca28 100644 --- a/lineax/__init__.py +++ b/lineax/__init__.py @@ -46,12 +46,14 @@ tridiagonal as tridiagonal, TridiagonalLinearOperator as TridiagonalLinearOperator, ) +from ._project import project as project from ._solution import RESULTS as RESULTS, Solution as Solution from ._solve import ( AbstractLinearSolver as AbstractLinearSolver, AutoLinearSolver as AutoLinearSolver, invert as invert, linear_solve as linear_solve, + projection_mv as projection_mv, ) from ._solver import ( BiCGStab as BiCGStab, @@ -76,6 +78,8 @@ MaxRankTag as MaxRankTag, negative_semidefinite_tag as negative_semidefinite_tag, positive_semidefinite_tag as positive_semidefinite_tag, + project_tags as project_tags, + project_tags_rules as project_tags_rules, symmetric_tag as symmetric_tag, tags_from_checks as tags_from_checks, transpose_tags as transpose_tags, diff --git a/lineax/_misc.py b/lineax/_misc.py index 025547d0..5f0bd23a 100644 --- a/lineax/_misc.py +++ b/lineax/_misc.py @@ -15,6 +15,7 @@ import equinox as eqx import jax +import jax.core import jax.numpy as jnp import jax.tree_util as jtu from jaxtyping import Array, ArrayLike, Bool, PyTree # pyright:ignore @@ -27,6 +28,16 @@ def tree_where( return jtu.tree_map(keep, true, false) +def to_shapedarray(x): + """Convert a `jax.ShapeDtypeStruct` leaf to a `jax.core.ShapedArray` (the abstract + value a primitive's `abstract_eval` rule must return); pass other leaves through. + """ + if isinstance(x, jax.ShapeDtypeStruct): + return jax.core.ShapedArray(x.shape, x.dtype) + else: + return x + + def resolve_rcond(rcond, n, m, dtype): if rcond is None: # This `2 *` is a heuristic: I have seen very rare failures without it, in ways diff --git a/lineax/_project.py b/lineax/_project.py new file mode 100644 index 00000000..445d10eb --- /dev/null +++ b/lineax/_project.py @@ -0,0 +1,353 @@ +# Copyright 2023 Google LLC +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +"""Projection operators for variable projection and separable nonlinear least squares. + +This module provides efficient projection operators P = A A^† with optimized +gradients following Golub & Pereyra (1973) for use in variable projection methods. +""" + +from typing import Any + +import equinox as eqx +import equinox.internal as eqxi +import jax +import jax.core +import jax.interpreters.ad as ad +import jax.lax as lax +import jax.numpy as jnp +import jax.tree_util as jtu +from equinox.internal import ω +from jaxtyping import Array, ArrayLike, PyTree + +from ._custom_types import sentinel +from ._misc import to_shapedarray +from ._operator import ( + AbstractLinearOperator, + conj, + FunctionLinearOperator, + IdentityLinearOperator, + linearise, + TangentLinearOperator, +) +from ._solve import ( + AbstractLinearSolver, + AutoLinearSolver, + check_rank_compat, + linear_solve, + projection_mv, +) +from ._tags import MaxRankTag, project_tags, tags_from_checks + + +# The solver-specific projection fast path `projection_mv` lives in `._solve`, with +# per-solver registrations in the solver modules (`_solver/svd.py`, `qr.py`, +# `normal.py`); it is shared with the `linear_solve` JVP's row-space projection. + + +# ============================================================================ +# Public API +# ============================================================================ + + +def project( + operator: AbstractLinearOperator, + solver: AbstractLinearSolver = AutoLinearSolver(well_posed=None), + *, + options: dict[str, Any] | None = None, + state: PyTree[Any] = sentinel, + throw: bool = True, +) -> AbstractLinearOperator: + r"""Returns a [`lineax.FunctionLinearOperator`][] representing the orthogonal + projection ``P = A A^†`` onto the range (column space) of `operator`. + + `project(A).mv(v)` is equivalent to `A.mv(linear_solve(A, v, solver).value)`: it + projects `v` onto `range(A)`. The result is idempotent (`P.mv(P.mv(v)) == P.mv(v)`) + and Hermitian. + + Like [`lineax.invert`][], the solver state (factorisation) is computed once and + reused across every mv. The returned operator provides the efficient + [Golub-Pereyra](https://doi.org/10.1137/0710036) gradient used to perform + [variable projection](https://www.cs.umd.edu/users/oleary/software/varpro.pdf) in + separable nonlinear least-squares problems; see the + [worked example](../examples/variable_projection.ipynb). + + **Arguments:** + + - `operator`: the linear operator `A` to project onto. May be any shape. (Raw + arrays are not accepted; wrap them as `lineax.MatrixLinearOperator(matrix)`.) + - `solver`: the linear solver to use. Defaults to + `AutoLinearSolver(well_posed=None)`. For rank-deficient `A`, use + `AutoLinearSolver(well_posed=False)` or another appropriate solver. + - `options`: additional options passed to the solver. Defaults to `None`. + - `state`: as [`lineax.invert`][]. Defaults to being initialised with gradients + stopped through the operator. + - `throw`: as [`lineax.linear_solve`][]. Defaults to `True`. + + **Returns:** + + A [`lineax.FunctionLinearOperator`][] whose `mv` applies `P @ v`. When `A` has full + row rank the projector is the identity, returned as a + [`lineax.IdentityLinearOperator`][]. + """ + if eqx.is_array(operator): + raise ValueError( + "`project(operator=...)` should be an `AbstractLinearOperator`, not a " + "raw JAX array. If you are trying to pass a matrix then this should be " + "passed as `lineax.MatrixLinearOperator(matrix)`." + ) + + m = operator.out_size() + n = operator.in_size() + + if m == 0 or n == 0: + raise ValueError( + f"operator cannot have zero dimensions, got out_size={m}, in_size={n}" + ) + + # raise if operator is rank-deficient and solver assumes full rank + check_rank_compat(solver, operator) + + if solver.assume_full_rank() and m <= n: + # square/wide matrix with full row rank: P = A A^† = I + return IdentityLinearOperator(operator.out_structure()) + + if options is None: + options = {} + if state == sentinel: + dynamic_operator, static_operator = eqx.partition(operator, eqx.is_array) + stopped_operator = eqx.combine( + lax.stop_gradient(dynamic_operator), static_operator + ) + state = solver.init(stopped_operator, options) + + dynamic_state, static_state = eqx.partition(state, eqx.is_array) + dynamic_state = lax.stop_gradient(dynamic_state) + state = eqx.combine(dynamic_state, static_state) + options = eqxi.nondifferentiable(options, name="`lineax.project(..., options=...)`") + solver = eqxi.nondifferentiable(solver, name="`lineax.project(..., solver=...)`") + + def mv(vector: PyTree[ArrayLike]) -> PyTree[Array]: + (out,) = eqxi.filter_primitive_bind( + _projection_mv_p, operator, state, vector, options, solver, throw + ) + return out + + # `tags_from_checks` drops MaxRankTag for the full rank case. + # Add the cheap supplemental MaxRankTag(min(m, n)) to compensatee. + # The min-rule in `project_tags` keeps whichever bound is tighter. + in_tags = tags_from_checks(operator) | {MaxRankTag(min(m, n))} + tags = project_tags(in_tags) + return FunctionLinearOperator(mv, operator.out_structure(), tags) + + +# ============================================================================ +# Primitive Implementation +# ============================================================================ + + +@eqxi.filter_primitive_def +def _projection_mv_impl( + Phi: AbstractLinearOperator, + state: Any, + v: PyTree[ArrayLike], + options: dict[str, Any], + solver: AbstractLinearSolver, + throw: bool, +) -> tuple[PyTree[Array]]: + """Implementation rule for the projection_mv primitive (base case): P @ v.""" + # Solver-specific fast path (e.g. O'Leary 1990 for SVD: P @ v = U U^H v) + out = projection_mv(solver, state, v, options) + if out is not NotImplemented: + return (out,) + + c = linear_solve(Phi, v, solver, state=state, options=options, throw=throw).value + return (Phi.mv(c),) + + +@eqxi.filter_primitive_def +def _projection_mv_abstract_eval( + Phi: Any, + state: Any, + v: Any, + options: Any, + solver: Any, + throw: Any, +) -> tuple[PyTree[jax.core.ShapedArray]]: + """Abstract evaluation - returns the output shape/dtype PyTree.""" + # P = Φ Φ^† maps out-space to out-space, so P @ v has the same PyTree structure, + # shapes and dtypes as `v` (which lives in the out-structure). Mirroring + # `linear_solve`: get the shape/dtype PyTree via `filter_eval_shape`, then convert + # to `ShapedArray` (the primitive's abstract-value type). (`Phi.out_structure()` + # is unusable here -- the operator's arrays are abstract.) + del Phi, state, options, solver, throw + struct = eqx.filter_eval_shape(lambda x: x, v) + return (jtu.tree_map(to_shapedarray, struct),) + + +def _is_none(x: Any) -> bool: + return x is None + + +@eqxi.filter_primitive_jvp +def _projection_mv_jvp( + primals: tuple[ + AbstractLinearOperator, + Any, + PyTree[ArrayLike], + dict[str, Any], + AbstractLinearSolver, + bool, + ], + tangents: tuple[Any, Any, Any, Any, Any, Any], +) -> tuple[tuple[PyTree[Array]], tuple[PyTree[Array]]]: + """JVP rule using the Golub-Pereyra efficient formulation. + + For P @ v where P = Φ Φ^†: + d(P @ v) = P @ dv + A + B + + where: + A = (I - P) @ (dΦ @ c) + B = (Φ^†)^T @ (dΦ^H @ r) + + **Tangent handling**: dPhi can be: + - AbstractLinearOperator (when Phi is MatrixLinearOperator) + - ndarray (when materializing gradients) + - Zero tangent (no gradient w.r.t. Phi) + """ + Phi, state, v, options, solver, throw = primals + dPhi, dstate, dv, doptions, dsolver, dthrow = tangents + # dstate, doptions, dsolver, dthrow should be None (non-differentiable) + del dstate, doptions, dsolver, dthrow + + # Forward pass: compute c = Φ^† v explicitly (the tangent terms below need it). + c = linear_solve(Phi, v, solver, state=state, options=options, throw=throw).value + Pv = Phi.mv(c) + + # dv term + if any(t is not None for t in jtu.tree_leaves(dv, is_leaf=_is_none)): + dv = jtu.tree_map(eqxi.materialise_zeros, v, dv, is_leaf=_is_none) + (P_dv,) = eqxi.filter_primitive_bind( + _projection_mv_p, Phi, state, dv, options, solver, throw + ) + else: + P_dv = jtu.tree_map(jnp.zeros_like, Pv) + + if all(t is None for t in jtu.tree_leaves(dPhi, is_leaf=_is_none)): + # No gradient w.r.t. Phi + return (Pv,), (P_dv,) + + dPhi_op = TangentLinearOperator(Phi, dPhi) + dPhi_op = linearise(dPhi_op) # Optimize for matvecs + + # Operator tangent path + dPhi_c = dPhi_op.mv(c) + + # A term: (I - P) @ (dΦ @ c) + (P_dPhi_c,) = eqxi.filter_primitive_bind( + _projection_mv_p, Phi, state, dPhi_c, options, solver, throw + ) + A = (dPhi_c**ω - P_dPhi_c**ω).ω + + # Full Golub-Pereyra: B term = (Φ^†)^H (dΦ^H r), where r = (I-P) v is the residual. + # Under the Kaufman approximation (1975, https://doi.org/10.1007/BF01932995) + # this term is dropped as an optimisation which is valid near an optimum as r -> 0. + # In lineax, we do not offer the Kaufman approximation as option to simplify the API + # and avoid a footgun due to inaccurate gradient when the residual is large. + r = (v**ω - Pv**ω).ω + dPhi_H_r = conj(dPhi_op).T.mv(r) # dΦ^H r + # (Φ^†)^H = (Φ^H)^† + transposed_state, transposed_options = solver.transpose(state, options) + conj_transposed_state, conj_transposed_options = solver.conj( + transposed_state, transposed_options + ) + B = linear_solve( + conj(Phi).T, + dPhi_H_r, + solver, + state=conj_transposed_state, + options=conj_transposed_options, + throw=throw, + ).value + return (Pv,), ((P_dv**ω + A**ω + B**ω).ω,) + + +def _is_undefined(x: Any) -> bool: + return isinstance(x, ad.UndefinedPrimal) + + +def _keep_undefined(v: Any, ct: Any) -> Any: + if _is_undefined(v): + return ct + else: + return None + + +@eqxi.filter_primitive_transpose(materialise_zeros=True) # pyright: ignore +def _projection_mv_transpose( + inputs: tuple[Any, Any, Any, Any, Any, Any], + cts_out: tuple[PyTree[Array]], +) -> tuple[Any, Any, Any, Any, Any, Any]: + """Transpose rule for VJP w.r.t. the vector input. + + The map ``v -> P v`` is linear, so its transpose applies ``P^T``. The + projection ``P = Φ Φ^†`` is Hermitian (P^H = P), hence: + * real Φ: P^T = P, so ct_v = P @ ct_out. + * complex Φ: P^T = conj(P), so ct_v = conj(P @ conj(ct_out)). + The complex branch is the conjugate of the real one; on the real path + `iscomplexobj` is False so no conjugation is emitted. + + The cotangent w.r.t. Phi is left to the JVP rule (Phi enters nonlinearly), + so ct_Phi is None here. + """ + (ct_out,) = cts_out + Phi, state, v, options, solver, throw = inputs + + if any(jnp.iscomplexobj(x) for x in jtu.tree_leaves(ct_out)): + (pc,) = eqxi.filter_primitive_bind( + _projection_mv_p, + Phi, + state, + jtu.tree_map(jnp.conj, ct_out), + options, + solver, + throw, + ) + ct_v = jtu.tree_map(jnp.conj, pc) + else: + (ct_v,) = eqxi.filter_primitive_bind( + _projection_mv_p, Phi, state, ct_out, options, solver, throw + ) + ct_v = jtu.tree_map(_keep_undefined, v, ct_v, is_leaf=_is_undefined) + + # Φ enters nonlinearly, so its cotangent is handled by the JVP rule (None + # here). The remaining inputs are non-differentiable. + ct_Phi = jtu.tree_map(lambda _: None, Phi) + ct_state = jtu.tree_map(lambda _: None, state) + ct_options = jtu.tree_map(lambda _: None, options) + ct_solver = jtu.tree_map(lambda _: None, solver) + ct_throw = None + + return ct_Phi, ct_state, ct_v, ct_options, ct_solver, ct_throw + + +# Create the primitive using create_vprim +# This automatically handles batching by vmapping the rules +_projection_mv_p = eqxi.create_vprim( + "projection_mv", + _projection_mv_impl, + _projection_mv_abstract_eval, + _projection_mv_jvp, + _projection_mv_transpose, +) diff --git a/lineax/_solve.py b/lineax/_solve.py index 817d710a..042cf4e3 100644 --- a/lineax/_solve.py +++ b/lineax/_solve.py @@ -14,6 +14,7 @@ import abc import functools as ft +from types import NotImplementedType from typing import Any, Generic, TypeAlias, TypeVar import equinox as eqx @@ -29,7 +30,7 @@ from jaxtyping import Array, ArrayLike, PyTree from ._custom_types import sentinel -from ._misc import inexact_asarray, strip_weak_dtype +from ._misc import inexact_asarray, strip_weak_dtype, to_shapedarray from ._operator import ( AbstractLinearOperator, conj, @@ -57,13 +58,6 @@ # -def _to_shapedarray(x): - if isinstance(x, jax.ShapeDtypeStruct): - return jax.core.ShapedArray(x.shape, x.dtype) - else: - return x - - def _to_struct(x): if isinstance(x, jax.core.ShapedArray): return jax.ShapeDtypeStruct(x.shape, x.dtype) @@ -138,7 +132,7 @@ def _linear_solve_abstract_eval(operator, state, vector, options, solver, throw) throw, check_closure=False, ) - out = jtu.tree_map(_to_shapedarray, out) + out = jtu.tree_map(to_shapedarray, out) return out @@ -239,11 +233,34 @@ def _linear_solve_jvp(primals, tangents): True, ) tmp2 = t_operator_conj_transpose.mv(tmp1) # pyright: ignore - # tmp2 is the y term - tmp3 = operator.mv(tmp2) - tmp4 = (-(tmp3**ω)).ω - # tmp4 is the Ay term - vecs.append(tmp4) + # tmp2 is the y term. The remaining contribution is -A^†A y, where A^†A + # is the row-space projector = the column-space projector of A^H. + # + # When `solver` provides a projection fast path (e.g. VVᴴ for SVD, QQᴴ + # for QR), route this through `project`: its custom primitive carries + # the efficient Golub-Pereyra gradient of the projector, which (unlike a + # bare `projection_mv` matvec on the frozen factorisation) + # gives correct *higher-order* derivatives of the projector. + projection_fast_path = eqx.filter_eval_shape( + projection_mv, + solver, + state_conj_transpose, # pyright: ignore + tmp2, + options_conj_transpose, # pyright: ignore + ) + if projection_fast_path is not NotImplemented: + from ._project import project + + P = project( + operator_conj_transpose, # pyright: ignore + solver, + state=state_conj_transpose, # pyright: ignore + ) + sols.append((-(P.mv(tmp2) ** ω)).ω) + else: + tmp3 = operator.mv(tmp2) + # tmp3 is the A y term; -A^†A y is recovered by the final A^† solve + vecs.append((-(tmp3**ω)).ω) sols.append(tmp2) vecs = jtu.tree_map(_sum, *vecs) # the A^ term at the very beginning @@ -474,9 +491,7 @@ def assume_full_rank(self) -> bool: """ -def _check_rank_compat( - solver: "AbstractLinearSolver", operator: AbstractLinearOperator -): +def check_rank_compat(solver: "AbstractLinearSolver", operator: AbstractLinearOperator): if solver.assume_full_rank(): dim_bound = min(operator.in_size(), operator.out_size()) if max_rank(operator) < dim_bound: @@ -490,6 +505,77 @@ def _check_rank_compat( ) +@ft.singledispatch +def projection_mv( + solver: "AbstractLinearSolver", + state: Any, + vector: PyTree[ArrayLike], + options: dict[str, Any], +) -> PyTree[Array] | NotImplementedType: + """Solver-specific fast path for the projection matvec ``P @ v = A A^† @ v``, + where `A` is the operator that `state` factorises, used by [`lineax.project`][]. + + A `functools.singledispatch` function dispatching on `type(solver)`. The base + implementation returns `NotImplemented`, signalling "no fast path -- use the + generic ``A (A^† @ v)`` solve". Register fast paths with + `@projection_mv.register(MySolver)`. + + !!! warning + + This is a *primal-value* fast path: it applies `P` using the **frozen** + factorisation in `state`, so it must not be differentiated directly (its + result is marked non-differentiable, so accidental autodiff raises rather than + returning a wrong gradient). To differentiate a projection, use + [`lineax.project`][], which supplies the correct gradient. + + **Arguments:** + + - `solver`: the (concrete) lineax solver. + - `state`: the solver state (factorisation) computed by `solver.init`. + - `vector`: the input vector `v`, a PyTree matching the factorised operator's + output structure. + - `options`: solver options (unused by the built-in fast paths). + + **Returns:** + + ``P @ v`` as a (non-differentiable) PyTree, or `NotImplemented`. + """ + del solver, state, vector, options + return NotImplemented + + +_register_projection_mv = projection_mv.register +_projection_mv_msg = ( + "Cannot autodifferentiate `lineax.projection_mv`: it applies the projector using " + "the frozen factorisation, so its derivative would be incorrect. Use " + "`lineax.project` to differentiate a projection." +) + + +def _guarded_register(cls): + # `@projection_mv.register(Cls)` decorator that wraps the implementation's + # (non-`NotImplemented`) output in `eqxi.nondifferentiable` -- see the + # `projection_mv` warning. Baking the guard in at registration means there is no + # unguarded path to it (registry/dispatch all return the guarded form). + def register(func): + @ft.wraps(func) + def guarded(solver, state, vector, options): + out = func(solver, state, vector, options) + if out is NotImplemented: + return out + return eqxi.nondifferentiable(out, msg=_projection_mv_msg) + + return _register_projection_mv(cls, guarded) + + return register + + +# Overriding the overloaded `singledispatch.register` attribute is not expressible to +# the type checker, so we suppress here (the call sites `@projection_mv.register(Cls)` +# are themselves fully type-checked). +projection_mv.register = _guarded_register # pyright: ignore[reportAttributeAccessIssue] + + _qr_token = eqxi.str2jax("qr_token") _diagonal_token = eqxi.str2jax("diagonal_token") _well_posed_diagonal_token = eqxi.str2jax("well_posed_diagonal_token") @@ -793,7 +879,7 @@ def linear_solve( stats={}, ) if state == sentinel: - _check_rank_compat(solver, operator) + check_rank_compat(solver, operator) dynamic_operator, static_operator = eqx.partition(operator, eqx.is_array) stopped_operator = eqx.combine( lax.stop_gradient(dynamic_operator), static_operator @@ -855,7 +941,7 @@ def invert( if options is None: options = {} - _check_rank_compat(solver, operator) + check_rank_compat(solver, operator) if state == sentinel: dynamic_operator, static_operator = eqx.partition(operator, eqx.is_array) stopped_operator = eqx.combine( diff --git a/lineax/_solver/normal.py b/lineax/_solver/normal.py index e46aabdd..469260aa 100644 --- a/lineax/_solver/normal.py +++ b/lineax/_solver/normal.py @@ -20,7 +20,7 @@ from .._operator import conj, linearise, materialise, TaggedLinearOperator from .._solution import RESULTS -from .._solve import AbstractLinearOperator, AbstractLinearSolver +from .._solve import AbstractLinearOperator, AbstractLinearSolver, projection_mv from .._tags import positive_semidefinite_tag from .cholesky import Cholesky @@ -186,3 +186,14 @@ def assume_full_rank(self): - `inner_solver`: The solver to wrap. It should support solving positive definite systems or positive semidefinite systems """ + + +@projection_mv.register(Normal) +def _(solver, state, vector, options): + del options + inner_state, tall, _op_conj_T, inner_options = state + if tall.value: + return NotImplemented # A A^† = A (A^H A)^† A^H - no fast path + # only reachable for rank-deficient inner solvers + # A A^† = A A^H (A A^H)^† = inner inner^† + return projection_mv(solver.inner_solver, inner_state, vector, inner_options) diff --git a/lineax/_solver/qr.py b/lineax/_solver/qr.py index 69e4b42c..fc50fd51 100644 --- a/lineax/_solver/qr.py +++ b/lineax/_solver/qr.py @@ -21,7 +21,7 @@ from jaxtyping import Array, PyTree from .._solution import RESULTS -from .._solve import AbstractLinearSolver +from .._solve import AbstractLinearSolver, projection_mv from .misc import ( pack_structures, PackedStructures, @@ -122,3 +122,18 @@ def assume_full_rank(self): Nothing. """ + + +@projection_mv.register(QR) +def _(solver, state, vector, options): + # A A^† v = Q Q^H v due to R being non-singular, computed using`ormqr` + del solver, options + (a, taus), transpose, packed = state + if transpose.value: + return NotImplemented # A A^† = I in wide case, should be unreachable + v = ravel_vector(vector, packed) + n_full, n_min = a.shape + qHv = jll.ormqr(a, taus, v[:, None], transpose=True)[:n_min, 0] # Q_thin^H v + padded = jnp.concatenate([qHv, jnp.zeros((n_full - n_min,), dtype=qHv.dtype)]) + pv = jll.ormqr(a, taus, padded[:, None], transpose=False)[:, 0] # Q_thin (.) + return unravel_solution(pv, transpose_packed_structures(packed)) diff --git a/lineax/_solver/svd.py b/lineax/_solver/svd.py index 5b872091..c1d3610a 100644 --- a/lineax/_solver/svd.py +++ b/lineax/_solver/svd.py @@ -23,7 +23,7 @@ from .._misc import resolve_rcond from .._operator import AbstractLinearOperator, max_rank from .._solution import RESULTS -from .._solve import AbstractLinearSolver +from .._solve import AbstractLinearSolver, projection_mv from .misc import ( pack_structures, PackedStructures, @@ -80,6 +80,16 @@ def init(self, operator: AbstractLinearOperator, options: dict[str, Any]): packed_structures = pack_structures(operator) return (u, s, vt), packed_structures + def _singular_mask(self, u, s, vt): + m = u.shape[0] + n = vt.shape[1] + rcond = resolve_rcond(self.rcond, n, m, s.dtype) + rcond = jnp.array(rcond, dtype=s.dtype) + if s.size > 0: + rcond = rcond * s[0] + # Not >=, or this fails with a matrix of all-zeros. + return s > rcond + def compute( self, state: _SVDState, @@ -89,14 +99,7 @@ def compute( del options (u, s, vt), packed_structures = state vector = ravel_vector(vector, packed_structures) - m, _ = u.shape - _, n = vt.shape - rcond = resolve_rcond(self.rcond, n, m, s.dtype) - rcond = jnp.array(rcond, dtype=s.dtype) - if s.size > 0: - rcond = rcond * s[0] - # Not >=, or this fails with a matrix of all-zeros. - mask = s > rcond + mask = self._singular_mask(u, s, vt) rank = mask.sum() safe_s = jnp.where(mask, s, 1) s_inv = jnp.where(mask, jnp.array(1.0) / safe_s, 0).astype(u.dtype) @@ -130,3 +133,15 @@ def assume_full_rank(self): precision times `max(N, M)`, where `(N, M)` is the shape of the operator. (I.e. `N` is the output size and `M` is the input size.) """ + + +@projection_mv.register(SVD) +def _(solver, state, vector, options): + # O'Leary 1990: A A^† v = U U^H v, due to unitarity of V + del options + (u, s, vt), packed = state # u already truncated to max_rank if tagged + v = ravel_vector(vector, packed) # out-space -> flat + mask = solver._singular_mask(u, s, vt).astype(u.dtype) + uHv = jnp.matmul(u.conj().T, v, precision=lax.Precision.HIGHEST) + pv = jnp.matmul(u, mask * uHv, precision=lax.Precision.HIGHEST) + return unravel_solution(pv, transpose_packed_structures(packed)) # flat -> out diff --git a/lineax/_tags.py b/lineax/_tags.py index 83051a8a..735b562b 100644 --- a/lineax/_tags.py +++ b/lineax/_tags.py @@ -261,6 +261,53 @@ def _(tags: frozenset[object]): # tridiagonal_tag intentionally absent: inverse of tridiagonal matrix generally dense. +project_tags_rules = [] + + +@project_tags_rules.append +def _(tags: frozenset[object]): + return positive_semidefinite_tag # P is always Hermitian PSD + + +@project_tags_rules.append +def _(tags: frozenset[object]): + if diagonal_tag in tags: # projection of a diagonal op is diagonal + return diagonal_tag + + +@project_tags_rules.append +def _(tags: frozenset[object]): + rank_tags = [t for t in tags if isinstance(t, MaxRankTag)] + if rank_tags: + # tightest bound, drop redundant tags + return min(rank_tags, key=lambda t: t.r) + + +def project_tags(tags: frozenset[object]) -> frozenset[object]: + """Lineax uses "tags" to declare that a particular linear operator exhibits some + property, e.g. symmetry. + + This function takes in a collection of tags representing a linear operator `A`, and + returns a collection of tags that should be associated with the orthogonal projector + `P = A A^†` onto `range(A)`. Specifically, `P` is always Hermitian positive + semidefinite, has max rank equal to that of `A`, and is diagonal if `A` is diagonal. + + **Arguments:** + + - `tags`: a `frozenset` of tags. + + **Returns:** + + A `frozenset` of tags. + """ + new_tags = [] + for rule in project_tags_rules: + out = rule(tags) + if out is not None: + new_tags.append(out) + return frozenset(new_tags) + + def invert_tags(tags: frozenset[object]) -> frozenset[object]: """Lineax uses "tags" to declare that a particular linear operator exhibits some property, e.g. symmetry. diff --git a/mkdocs.yml b/mkdocs.yml index e6cfc389..f46402c8 100644 --- a/mkdocs.yml +++ b/mkdocs.yml @@ -97,6 +97,7 @@ nav: - Examples: - 'examples/classical_solve.ipynb' - 'examples/least_squares.ipynb' + - 'examples/variable_projection.ipynb' - 'examples/structured_matrices.ipynb' - 'examples/no_materialisation.ipynb' - 'examples/operators.ipynb' diff --git a/pyproject.toml b/pyproject.toml index 3c5b6ccf..f7982a24 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -12,12 +12,14 @@ dev = [ docs = [ "hippogriffe==0.2.2", "griffe==1.7.3", + "matplotlib==3.10.9", "mkdocs==1.6.1", "mkdocs-include-exclude-files==0.1.0", "mkdocs-ipynb==0.1.1", "mkdocs-material==9.6.7", "mkdocstrings==0.28.3", "mkdocstrings-python==1.16.8", + "optimistix==0.1.0", "pygments==2.20.0", "pymdown-extensions==10.21.2" ] diff --git a/tests/test_project.py b/tests/test_project.py new file mode 100644 index 00000000..6eabae48 --- /dev/null +++ b/tests/test_project.py @@ -0,0 +1,463 @@ +# Copyright 2023 Google LLC +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +import functools as ft + +import equinox as eqx +import jax +import jax.numpy as jnp +import jax.random as jr +import lineax as lx +import pytest + +from .helpers import ( + finite_difference_jvp, + make_function_operator, + make_jac_operator, + make_jacfwd_operator, + make_jacrev_operator, + make_matrix_operator, + make_trivial_pytree_operator, + tree_allclose, +) + + +# Projection is only interesting for rectangular / rank-deficient `A`, so we use +# tall operators throughout. Constructors that support non-square matrices: +nonsquare_operators = ( + make_matrix_operator, + make_function_operator, + make_jac_operator, + make_jacfwd_operator, + make_jacrev_operator, + make_trivial_pytree_operator, +) + +# Operators that support forward-mode autodiff (excludes `make_jacrev_operator`, +# which uses a `custom_vjp` with no JVP). +fwd_operators = (make_jac_operator, make_jacfwd_operator) + +# Full-rank-assuming solvers: handle tall full-rank `A` (m > n). +full_rank_solvers = ( + lx.AutoLinearSolver(well_posed=None), + lx.QR(), + lx.Normal(lx.Cholesky()), +) +# Pseudoinverse solvers: handle (tagged) rank-deficient `A`. +pseudo_solvers = (lx.SVD(), lx.AutoLinearSolver(well_posed=False)) + +_M, _N = 5, 3 # tall: m > n +_RANK = 2 # exact rank for the rank-deficient cases + + +def _tall_full_rank(getkey, dtype): + return jr.normal(getkey(), (_M, _N), dtype=dtype) + + +def _tall_rank_deficient(getkey, dtype): + # Exactly rank `_RANK`: (m, r) @ (r, n). Use *exact* (tagged) rank deficiency + # -- the SVD fast path applies `U U^H` over every cached singular direction + # and does not rcond-mask, so a merely near-singular matrix would mismatch. + b = jr.normal(getkey(), (_M, _RANK), dtype=dtype) + c = jr.normal(getkey(), (_RANK, _N), dtype=dtype) + return b @ c + + +# (kind, solver, tags): full-rank with full-rank solvers; rank-deficient (tagged) +# with pseudoinverse solvers. +structural_cases = [ + *[("full", s, ()) for s in full_rank_solvers], + *[("deficient", s, lx.MaxRankTag(_RANK)) for s in pseudo_solvers], +] + + +def _make_matrix(getkey, kind, dtype): + if kind == "full": + return _tall_full_rank(getkey, dtype) + else: + return _tall_rank_deficient(getkey, dtype) + + +@pytest.mark.parametrize("make_operator", nonsquare_operators) +@pytest.mark.parametrize("kind, solver, tags", structural_cases) +@pytest.mark.parametrize("dtype", (jnp.float64, jnp.complex128)) +def test_project(make_operator, kind, solver, tags, dtype, getkey): + matrix = _make_matrix(getkey, kind, dtype) + operator = make_operator(getkey, matrix, tags) + A = operator.as_matrix() + P = lx.project(operator, solver) + v = jr.normal(getkey(), (_M,), dtype=dtype) + Pv = P.mv(v) + + # Idempotent: P P v == P v. + assert tree_allclose(P.mv(Pv), Pv, rtol=1e-5, atol=1e-6) + + # Matches the dense projector A A^†. + dense = (A @ jnp.linalg.pinv(A)) @ v + assert tree_allclose(Pv, dense, rtol=1e-5, atol=1e-6) + + # Hermitian: P == P^H, and conj(P).T == P. + Pmat = P.as_matrix() + assert tree_allclose(Pmat, Pmat.conj().T, rtol=1e-5, atol=1e-6) + assert tree_allclose(lx.conj(P).T.as_matrix(), Pmat, rtol=1e-5, atol=1e-6) + + # Tags. (`project` returns a `FunctionLinearOperator` for these tall operators.) + assert isinstance(P, lx.FunctionLinearOperator) + assert lx.is_positive_semidefinite(P) + assert lx.symmetric_tag not in P.tags + # Symmetry is inferred for real dtypes, not for complex (where P is Hermitian). + if jnp.iscomplexobj(jnp.zeros((), dtype=dtype)): + assert not lx.is_symmetric(P) + else: + assert lx.is_symmetric(P) + assert lx.max_rank(P) == lx.max_rank(operator) + + +def test_project_diagonal(getkey): + # Projection of a diagonal operator is diagonal (0/1 on the diagonal). Use a + # diagonal-with-zeros operator with an *exact* rank tag so the SVD fast path + # truncates to the nonzero directions. + for dtype in (jnp.float64, jnp.complex128): + diag = jnp.array([1.0, 2.0, 0.0, 3.0, 0.0], dtype=dtype) + nnz = 3 + operator = lx.MatrixLinearOperator( + jnp.diag(diag), (lx.diagonal_tag, lx.MaxRankTag(nnz)) + ) + P = lx.project(operator, lx.SVD()) + assert lx.is_diagonal(P) + assert lx.max_rank(P) == nnz + v = jr.normal(getkey(), (5,), dtype=dtype) + dense = (operator.as_matrix() @ jnp.linalg.pinv(operator.as_matrix())) @ v + assert tree_allclose(P.mv(v), dense, rtol=1e-5, atol=1e-6) + + +# --------------------------------------------------------------------------- +# Gradients (Golub-Pereyra). The custom JVP/transpose rules should be the +# efficient equivalent of the naive `A (A^† v)` autodiff. We check this against +# two references: +# * `lineax.invert` (`A @ invert(A) @ v`), which uses the *same* solver and +# adjoint machinery -- it handles rank deficiency (via the pseudoinverse +# solver) and shares `project`'s complex (Wirtinger/adjoint) convention, so +# it is valid for full-rank/rank-deficient and real/complex alike; and +# * finite differences, an entirely lineax-independent oracle, used for the +# real full-rank case (complex FD is non-holomorphic, but `project` matches +# it to FD tolerance anyway since the convention bug is fixed). +# --------------------------------------------------------------------------- + + +@pytest.mark.parametrize("kind, solver, tags", structural_cases) +@pytest.mark.parametrize("dtype", (jnp.float64, jnp.complex128)) +def test_project_jvp_invert(kind, solver, tags, dtype, getkey): + matrix = _make_matrix(getkey, kind, dtype) + t_matrix = _make_matrix(getkey, kind, dtype) + vec = jr.normal(getkey(), (_M,), dtype=dtype) + t_vec = jr.normal(getkey(), (_M,), dtype=dtype) + + def run(M, v): + return lx.project(lx.MatrixLinearOperator(M, tags), solver).mv(v) + + def ref(M, v): + op = lx.MatrixLinearOperator(M, tags) + return op.mv(lx.invert(op, solver).mv(v)) + + out, t_out = eqx.filter_jvp(run, (matrix, vec), (t_matrix, t_vec)) + ref_out, t_ref = eqx.filter_jvp(ref, (matrix, vec), (t_matrix, t_vec)) + assert tree_allclose(out, ref_out, rtol=1e-4, atol=1e-6) + assert tree_allclose(t_out, t_ref, rtol=1e-4, atol=1e-6) + + +@pytest.mark.parametrize("solver", full_rank_solvers) +@pytest.mark.parametrize("dtype", (jnp.float64, jnp.complex128)) +def test_project_jvp_fd(solver, dtype, getkey): + # JVP against finite differences -- a lineax-independent ground truth. (The + # complex case is formally non-holomorphic, but with the Hermitian B-term the + # custom JVP still matches FD to FD tolerance.) + matrix = _tall_full_rank(getkey, dtype) + t_matrix = _tall_full_rank(getkey, dtype) + vec = jr.normal(getkey(), (_M,), dtype=dtype) + t_vec = jr.normal(getkey(), (_M,), dtype=dtype) + + def run(M, v): + return lx.project(lx.MatrixLinearOperator(M), solver).mv(v) + + out, t_out = eqx.filter_jvp(run, (matrix, vec), (t_matrix, t_vec)) + expected, t_expected = finite_difference_jvp(run, (matrix, vec), (t_matrix, t_vec)) + assert tree_allclose(out, expected, rtol=1e-4, atol=1e-6) + assert tree_allclose(t_out, t_expected, rtol=1e-3, atol=1e-4) + + +@pytest.mark.parametrize("solver", full_rank_solvers) +@pytest.mark.parametrize("dtype", (jnp.float64, jnp.complex128)) +def test_project_grad_invert(solver, dtype, getkey): + # Reverse-mode gradient of a real-valued loss matches autodiff through + # `lineax.invert` -- the quantity actually consumed downstream (separable NLLS). + matrix = _tall_full_rank(getkey, dtype) + vec = jr.normal(getkey(), (_M,), dtype=dtype) + + def loss_project(M): + P = lx.project(lx.MatrixLinearOperator(M), solver) + return jnp.sum(jnp.abs(P.mv(vec)) ** 2) + + def loss_invert(M): + op = lx.MatrixLinearOperator(M) + return jnp.sum(jnp.abs(op.mv(lx.invert(op, solver).mv(vec))) ** 2) + + assert tree_allclose( + jax.grad(loss_project)(matrix), + jax.grad(loss_invert)(matrix), + rtol=1e-4, + atol=1e-6, + ) + + +@pytest.mark.parametrize("make_operator", fwd_operators) +@pytest.mark.parametrize("solver", full_rank_solvers) +@pytest.mark.parametrize("dtype", (jnp.float64, jnp.complex128)) +def test_project_jvp_operator(make_operator, solver, dtype, getkey): + # The operator-tangent path (TangentLinearOperator / linearise) goes through + # `JacobianLinearOperator`. Finite-differencing the operator pytree does not + # perturb closure-captured arrays, so we instead check the analytic JVP + # through the operator matches the analytic JVP through the equivalent + # `MatrixLinearOperator` (grounded against finite differences / pinv above). + matrix = _tall_full_rank(getkey, dtype) + t_matrix = _tall_full_rank(getkey, dtype) + vec = jr.normal(getkey(), (_M,), dtype=dtype) + + make_op = ft.partial(make_operator, getkey) + operator, t_operator = eqx.filter_jvp(make_op, (matrix, ()), (t_matrix, ())) + make_mat = ft.partial(make_matrix_operator, getkey) + mat_op, t_mat_op = eqx.filter_jvp(make_mat, (matrix, ()), (t_matrix, ())) + + f = lambda op: lx.project(op, solver).mv(vec) + out, t_out = eqx.filter_jvp(f, (operator,), (t_operator,)) + out_ref, t_out_ref = eqx.filter_jvp(f, (mat_op,), (t_mat_op,)) + assert tree_allclose(out, out_ref, rtol=1e-4, atol=1e-6) + assert tree_allclose(t_out, t_out_ref, rtol=1e-4, atol=1e-6) + + +@pytest.mark.parametrize("solver", full_rank_solvers) +def test_project_jacfwd_jacrev(solver, getkey): + matrix = _tall_full_rank(getkey, jnp.float64) + operator = lx.MatrixLinearOperator(matrix) + vec = jr.normal(getkey(), (_M,), dtype=jnp.float64) + f = lambda op: lx.project(op, solver).mv(vec) + jfwd = eqx.filter_jacfwd(f)(operator) + jrev = eqx.filter_jacrev(f)(operator) + assert tree_allclose(jfwd, jrev, rtol=1e-4, atol=1e-6) + + +@pytest.mark.parametrize("solver", full_rank_solvers) +def test_project_second_order(solver, getkey): + # Second-order finite-difference check (real dtype), guarding the primitive + # re-binding in the JVP/transpose rules (which is what makes higher-order + # autodiff recurse). + matrix = _tall_full_rank(getkey, jnp.float64) + t_matrix = _tall_full_rank(getkey, jnp.float64) + tt_matrix = _tall_full_rank(getkey, jnp.float64) + vec = jr.normal(getkey(), (_M,), dtype=jnp.float64) + + def run(M): + return lx.project(lx.MatrixLinearOperator(M), solver).mv(vec) + + # First-order directional derivative as a function of `M`. + jvp1 = lambda M: eqx.filter_jvp(run, (M,), (t_matrix,))[1] + + out, t_out = eqx.filter_jvp(jvp1, (matrix,), (tt_matrix,)) + expected, t_expected = finite_difference_jvp(jvp1, (matrix,), (tt_matrix,)) + assert tree_allclose(out, expected, rtol=1e-3, atol=1e-4) + assert tree_allclose(t_out, t_expected, rtol=1e-2, atol=1e-3) + + +# --------------------------------------------------------------------------- +# Targeted cases. +# --------------------------------------------------------------------------- + + +@pytest.mark.parametrize("solver", full_rank_solvers) +def test_project_identity_shortcut(solver, getkey): + # Full-rank solver + m <= n => P = I exactly, returned as IdentityLinearOperator. + matrix = jr.normal(getkey(), (_N, _M), dtype=jnp.float64) # wide: m < n + operator = lx.MatrixLinearOperator(matrix) + P = lx.project(operator, solver) + assert isinstance(P, lx.IdentityLinearOperator) + + +@pytest.mark.parametrize("solver_cls", (lx.SVD, lx.QR)) +def test_project_fast_path(solver_cls, getkey): + # The SVD / QR registrations reproduce the dense projector, and report a fast + # path (not NotImplemented) for a tall operator. + for dtype in (jnp.float64, jnp.complex128): + matrix = _tall_full_rank(getkey, dtype) + operator = lx.MatrixLinearOperator(matrix) + solver = solver_cls() + P = lx.project(operator, solver) + v = jr.normal(getkey(), (_M,), dtype=dtype) + dense = (matrix @ jnp.linalg.pinv(matrix)) @ v + assert tree_allclose(P.mv(v), dense, rtol=1e-5, atol=1e-6) + state = solver.init(operator, {}) + out = lx.projection_mv(solver, state, v, {}) + assert out is not NotImplemented + assert tree_allclose(out, dense, rtol=1e-5, atol=1e-6) + + +def test_project_options_threaded(getkey): + # `options` must be threaded into every inner solve (as in `invert`), not just used + # in `solver.init`. Exercise the options path with an iterative inner solver and a + # `preconditioner` option: `project` must accept it and still equal the dense + # projector. (A primitive-signature regression that dropped `options` would error.) + matrix = _tall_full_rank(getkey, jnp.float64) + operator = lx.MatrixLinearOperator(matrix) + v = jr.normal(getkey(), (_M,), dtype=jnp.float64) + preconditioner = lx.IdentityLinearOperator(jax.ShapeDtypeStruct((_N,), jnp.float64)) + solver = lx.Normal(lx.CG(rtol=1e-10, atol=1e-10)) + P = lx.project(operator, solver, options={"preconditioner": preconditioner}) + dense = (matrix @ jnp.linalg.pinv(matrix)) @ v + assert tree_allclose(P.mv(v), dense, rtol=1e-5, atol=1e-6) + + +def test_project_closure_operator(getkey): + # `project` differentiates correctly through an operator that closes over a + # parameter (with no explicit pytree leaf for it) -- matching an explicit + # `MatrixLinearOperator`. So `project` needs no separate `check_closure` guard: + # `filter_primitive_bind` converts the closure to operator leaves (seen by the JVP), + # and the solves delegate to the already-guarded `linear_solve`. + matrix = _tall_full_rank(getkey, jnp.float64) + z = jr.normal(getkey(), (_M,), dtype=jnp.float64) + + def loss_closure(m): + op = lx.FunctionLinearOperator( + lambda c: m @ c, jax.ShapeDtypeStruct((_N,), m.dtype) + ) + return jnp.sum(lx.project(op, lx.SVD()).mv(z) ** 2) + + def loss_matrix(m): + return jnp.sum(lx.project(lx.MatrixLinearOperator(m), lx.SVD()).mv(z) ** 2) + + assert tree_allclose( + jax.grad(loss_closure)(matrix), + jax.grad(loss_matrix)(matrix), + rtol=1e-5, + atol=1e-6, + ) + + +def test_projection_mv_nondifferentiable(getkey): + # `projection_mv` is a primal-only fast path: it applies the projector using the + # *frozen* factorisation, so differentiating it directly would give an incorrect + # projector gradient. Its output is `eqxi.nondifferentiable`, so accidental + # autodiff raises rather than silently returning a wrong gradient. (Differentiable + # projections go through `lineax.project`, which supplies the Golub-Pereyra rule.) + matrix = _tall_full_rank(getkey, jnp.float64) + operator = lx.MatrixLinearOperator(matrix) + solver = lx.SVD() + state = solver.init(operator, {}) + v = jr.normal(getkey(), (_M,), dtype=jnp.float64) + # Primal evaluation works. + assert lx.projection_mv(solver, state, v, {}) is not NotImplemented + # Differentiating it directly raises (both forward and reverse mode). + f = lambda vec: lx.projection_mv(solver, state, vec, {}) + with pytest.raises(Exception): + eqx.filter_jvp(f, (v,), (v,)) + with pytest.raises(Exception): + eqx.filter_grad(lambda vec: jnp.sum(f(vec)))(v) + + +def test_project_normal_tall_generic(getkey): + # Tall Normal(Cholesky) has no projection fast path: P = A (A^H A)^† A^H lives + # in m-space, not the n-space inner B B^†. The registration declines + # (NotImplemented) and the generic path is used. + matrix = _tall_full_rank(getkey, jnp.float64) + operator = lx.MatrixLinearOperator(matrix) + solver = lx.Normal(lx.Cholesky()) + state = solver.init(operator, {}) + v = jr.normal(getkey(), (_M,), dtype=jnp.float64) + assert lx.projection_mv(solver, state, v, {}) is NotImplemented + # The generic path still produces the correct projector. + P = lx.project(operator, solver) + dense = (matrix @ jnp.linalg.pinv(matrix)) @ v + assert tree_allclose(P.mv(v), dense, rtol=1e-5, atol=1e-6) + + +def _multileaf_pytree_operator(getkey, dtype): + # Multi-leaf pytree out-structure: in (3,) -> {"a": (2,), "b": (4,)}, tall full rank + M1 = jr.normal(getkey(), (2, 3), dtype=dtype) + M2 = jr.normal(getkey(), (4, 3), dtype=dtype) + in_struct = jax.ShapeDtypeStruct((3,), dtype) + return lx.FunctionLinearOperator(lambda x: {"a": M1 @ x, "b": M2 @ x}, in_struct) + + +@pytest.mark.parametrize("solver", full_rank_solvers) +@pytest.mark.parametrize("dtype", (jnp.float64, jnp.complex128)) +def test_project_multileaf_pytree(solver, dtype, getkey): + # `project` must support operators whose out-structure is a multi-leaf pytree + # (its primitive's abstract-eval / JVP / transpose are pytree-aware). + operator = _multileaf_pytree_operator(getkey, dtype) + P = lx.project(operator, solver) + vec = { + "a": jr.normal(getkey(), (2,), dtype=dtype), + "b": jr.normal(getkey(), (4,), dtype=dtype), + } + Pv = P.mv(vec) + # Idempotent. + assert tree_allclose(P.mv(Pv), Pv, rtol=1e-5, atol=1e-6) + # Matches the dense projector (leaves flatten in pytree order: "a" then "b"). + A = operator.as_matrix() + flat = jnp.concatenate([vec["a"], vec["b"]]) + dense = (A @ jnp.linalg.pinv(A)) @ flat + assert tree_allclose( + jnp.concatenate([Pv["a"], Pv["b"]]), dense, rtol=1e-5, atol=1e-6 + ) + + +@pytest.mark.parametrize("solver", full_rank_solvers) +def test_project_multileaf_pytree_grad(solver, getkey): + # JVP (vs finite differences) and jacfwd == jacrev through a multi-leaf pytree. + M2 = jr.normal(getkey(), (4, 3), dtype=jnp.float64) + vec = { + "a": jr.normal(getkey(), (2,), dtype=jnp.float64), + "b": jr.normal(getkey(), (4,), dtype=jnp.float64), + } + + def run(M1): + op = lx.FunctionLinearOperator( + lambda x: {"a": M1 @ x, "b": M2 @ x}, + jax.ShapeDtypeStruct((3,), jnp.float64), + ) + Pv = lx.project(op, solver).mv(vec) + return jnp.concatenate([Pv["a"], Pv["b"]]) + + M1 = jr.normal(getkey(), (2, 3), dtype=jnp.float64) + t_M1 = jr.normal(getkey(), (2, 3), dtype=jnp.float64) + _, t_out = eqx.filter_jvp(run, (M1,), (t_M1,)) + _, t_fd = finite_difference_jvp(run, (M1,), (t_M1,)) + assert tree_allclose(t_out, t_fd, rtol=1e-3, atol=1e-4) + assert tree_allclose(jax.jacfwd(run)(M1), jax.jacrev(run)(M1), rtol=1e-4, atol=1e-6) + + +def test_project_errors(getkey): + # Raw array rejected. (Under the jaxtyping/beartype import hook used in the + # test suite this surfaces as a `TypeError`; the explicit `ValueError` in + # `project` fires when the hook is absent.) + with pytest.raises((ValueError, TypeError)): + lx.project(jnp.eye(3)) # pyright: ignore + + # Zero-dimension operator rejected. + with pytest.raises(ValueError): + lx.project(lx.MatrixLinearOperator(jnp.zeros((3, 0)))) + + # Full-rank solver on a declared rank-deficient operator raises (rank check). + matrix = _tall_rank_deficient(getkey, jnp.float64) + operator = lx.MatrixLinearOperator(matrix, lx.MaxRankTag(_RANK)) + with pytest.raises(ValueError): + lx.project(operator, lx.AutoLinearSolver(well_posed=None))