diff --git a/docs/p/api/at.latticetools.response_matrix.rst b/docs/p/api/at.latticetools.response_matrix.rst
new file mode 100644
index 000000000..b60b5cb70
--- /dev/null
+++ b/docs/p/api/at.latticetools.response_matrix.rst
@@ -0,0 +1,15 @@
+at.latticetools.response\_matrix
+================================
+
+.. automodule:: at.latticetools.response_matrix
+   :inherited-members:
+
+   
+   .. rubric:: Classes
+
+   .. autosummary::
+   
+      ResponseMatrix
+      OrbitResponseMatrix
+      TrajectoryResponseMatrix
+   
\ No newline at end of file
diff --git a/docs/p/index.rst b/docs/p/index.rst
index 714b6ee7c..099afe969 100644
--- a/docs/p/index.rst
+++ b/docs/p/index.rst
@@ -36,7 +36,8 @@ Sub-packages
    howto/multiprocessing
    howto/CavityControl
    howto/Collective
-   Working with MAD-X files <api/at.load.madx>
+   Work with MAD-X files <api/at.load.madx>
+   Use response matrices <notebooks/response_matrices>
 
 .. autosummary::
    :toctree: api
diff --git a/docs/p/notebooks/observables.ipynb b/docs/p/notebooks/observables.ipynb
index 668d8cef2..33e0c26a1 100644
--- a/docs/p/notebooks/observables.ipynb
+++ b/docs/p/notebooks/observables.ipynb
@@ -936,7 +936,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.9.20"
+   "version": "3.11.11"
   }
  },
  "nbformat": 4,
diff --git a/docs/p/notebooks/response_matrices.ipynb b/docs/p/notebooks/response_matrices.ipynb
new file mode 100644
index 000000000..570c0f514
--- /dev/null
+++ b/docs/p/notebooks/response_matrices.ipynb
@@ -0,0 +1,932 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "addfb8d6-8b83-45b7-b3be-650fa0e5bed3",
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "source": [
+    "# Response matrices"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "id": "ce595103-f8d4-4425-9c16-c4f766cd11c6",
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "outputs": [],
+   "source": [
+    "import at\n",
+    "import numpy as np\n",
+    "import math\n",
+    "from pathlib import Path\n",
+    "from importlib.resources import files, as_file\n",
+    "from timeit import timeit\n",
+    "from at.future import VariableList, RefptsVariable"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "id": "f21dfcdb-1fc4-4ab2-a34f-48d1eb467c23",
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "outputs": [],
+   "source": [
+    "with as_file(files(\"machine_data\") / \"hmba.mat\") as path:\n",
+    "    hmba_lattice = at.load_lattice(path)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 3,
+   "id": "f4b3676d-a5d3-4f16-9fd4-20788606b6b6",
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "outputs": [],
+   "source": [
+    "for sx in hmba_lattice.select(at.Sextupole):\n",
+    "    sx.KickAngle=[0,0]\n",
+    "hmba_lattice.enable_6d()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 4,
+   "id": "620def05-b2fe-4cce-8566-9bd10d2018af",
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "outputs": [],
+   "source": [
+    "ring = hmba_lattice.repeat(8)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "64d73aed-8e8b-4be7-9489-0e0d4c70fd95",
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "source": [
+    "A {py:class}`.ResponseMatrix` object defines a general-purpose response matrix, based\n",
+    "on a {py:class}`.VariableList` of attributes which will be independently varied, and an\n",
+    "{py:class}`.ObservableList` of attributes which will be recorded for each\n",
+    "variable step.\n",
+    "\n",
+    "{py:class}`.ResponseMatrix` objects can be combined with the \"+\" operator to define\n",
+    "combined responses. This concatenates the variables and the observables.\n",
+    "\n",
+    "The module also defines two commonly used response matrices:\n",
+    "{py:class}`.OrbitResponseMatrix` for circular machines and\n",
+    "{py:class}`.TrajectoryResponseMatrix` for beam lines. Other matrices can be easily\n",
+    "defined by providing the desired Observables and Variables to the\n",
+    "{py:class}`.ResponseMatrix` base class."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "08ae06fa-dea3-4542-acbc-4e59a14c66a7",
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "source": [
+    "## General purpose response matrix\n",
+    "\n",
+    "Let's take the horizontal displacements of all quadrupoles as variables:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 5,
+   "id": "68a54c6d-be1d-41a9-90f8-f496b22d076c",
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "outputs": [],
+   "source": [
+    "variables = VariableList(RefptsVariable(ik, \"dx\", name=f\"dx_{ik}\", delta=.0001)\n",
+    "                        for ik in ring.get_uint32_index(at.Quadrupole))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "2f70b1b0-bfc0-4e01-8ada-d7d76b77f406",
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "source": [
+    "Variable names are set to _dx\\_nnnn_ where _nnnn_ is the index of the quadrupole in the ring.\n",
+    "\n",
+    "Let's take the horizontal positions at all beam position monitors as observables:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 6,
+   "id": "8490caa3-14c1-4c7f-ba1e-f4e64078b6a5",
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "outputs": [],
+   "source": [
+    "observables = at.ObservableList([at.OrbitObservable(at.Monitor, axis='x')])"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "fee7a274-979c-4503-90b1-b69d248bc317",
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "source": [
+    "We have a single Observable named _orbit[x]_ by default, with multiple values.\n",
+    "\n",
+    "### Instantiation"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 7,
+   "id": "2d64d87b-ba18-4c22-921f-f6765f74080d",
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "outputs": [],
+   "source": [
+    "resp_dx = at.ResponseMatrix(ring, variables, observables)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "9f283b27-6898-4c00-9ac5-ec9ba48bd9f5",
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "source": [
+    "At that point, the response matrix is empty.\n",
+    "\n",
+    "### Matrix Building\n",
+    "\n",
+    "A general purpose response matrix may be filled by several methods:\n",
+    "\n",
+    "1. Direct assignment of an array to the {py:attr}`~.ResponseMatrix.response` property.\n",
+    "   The shape of the array is checked,\n",
+    "2. {py:meth}`~.ResponseMatrix.load` loads data from a file containing previously\n",
+    "   saved values or experimentally measured values,\n",
+    "3. {py:meth}`~.ResponseMatrix.build_tracking` computes the matrix using tracking,\n",
+    "4. For some specialized response matrices\n",
+    "   {py:meth}`~OrbitResponseMatrix.build_analytical` is available.\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 8,
+   "id": "23154d5d-bdf8-43ed-8b4a-cfc616d7e659",
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "array([[ 16.94897864,  -7.67022307,   2.50968594, ...,   5.3125638 ,\n",
+       "         -5.57239476,  14.39501938],\n",
+       "       [-10.6549627 ,   3.59085167,  -6.21666755, ...,  -2.46632948,\n",
+       "          8.42841189, -18.50171186],\n",
+       "       [-10.99744814,   3.91741643,  -5.60080281, ...,  -2.73604877,\n",
+       "          7.61448343, -17.01886021],\n",
+       "       ...,\n",
+       "       [-17.0182359 ,   7.61358522,  -2.73824047, ...,  -5.60018031,\n",
+       "          3.92050176, -11.00305834],\n",
+       "       [-18.50166601,   8.42772786,  -2.46888914, ...,  -6.21619686,\n",
+       "          3.59444354, -10.66160249],\n",
+       "       [ 14.38971545,  -5.56943704,   5.31320321, ...,   2.50757509,\n",
+       "         -7.6711941 ,  16.94982252]], shape=(80, 128))"
+      ]
+     },
+     "execution_count": 8,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "resp_dx.build_tracking(use_mp=True)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "6ac09735-b15f-49a4-b44e-4051889e0bae",
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "source": [
+    "### Matrix normalisation\n",
+    "\n",
+    "To be correctly inverted, the response matrix must be correctly normalised: the norms\n",
+    "of its columns must be of the same order of magnitude, and similarly for the rows.\n",
+    "\n",
+    "Normalisation is done by adjusting the weights {math}`w_v` for the variables {math}`\\mathbf{V}`\n",
+    "and {math}`w_o` for the observables {math}`\\mathbf{O}`.\n",
+    "With {math}`\\mathbf{R}` the response matrix:\n",
+    "\n",
+    ":::{math}\n",
+    "\n",
+    "   \\mathbf{O} = \\mathbf{R} . \\mathbf{V}\n",
+    ":::\n",
+    "\n",
+    "The weighted response matrix {math}`\\mathbf{R}_w` is:\n",
+    "\n",
+    ":::{math}\n",
+    "\n",
+    "   \\frac{\\mathbf{O}}{w_o} = \\mathbf{R}_w . \\frac{\\mathbf{V}}{w_v}\n",
+    ":::\n",
+    "The {math}`\\mathbf{R}_w` is dimensionless and should be normalised. This can be checked\n",
+    "using:\n",
+    "\n",
+    "* {py:meth}`~.ResponseMatrix.check_norm` which prints the ratio of the maximum / minimum\n",
+    "  norms for variables and observables. These should be less than 10.\n",
+    "* {py:meth}`~.ResponseMatrix.plot_norm`\n",
+    "\n",
+    "Both natural and weighted response matrices can be retrieved with the\n",
+    "{py:meth}`~.ResponseMatrix.response` and {py:meth}`~.ResponseMatrix.weighted_response`\n",
+    "properties."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 9,
+   "id": "f74e856f-0c6f-4a89-9493-255c1cc3f03f",
+   "metadata": {
+    "editable": true,
+    "jp-MarkdownHeadingCollapsed": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "max/min Observables: 2.8352796928877786\n",
+      "max/min Variables: 4.768846272299542\n"
+     ]
+    },
+    {
+     "data": {
+      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAosAAAHrCAYAAACn9tfQAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjEsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvc2/+5QAAAAlwSFlzAAAPYQAAD2EBqD+naQAAWkZJREFUeJzt3Qm8jHX///GPfc2eNVtSyE5ESkVRp0W6bcmW6C6ylUK2kqh+RCGp6O4ukfuOVFJCWuxb3VqkErJLlgjh+j/eX/+ZZubMdZyD4yzzej4e48xcc8011/Wd65j3+W5XBs/zPAMAAACiyBhtIQAAACCERQAAAPgiLAIAAMAXYREAAAC+CIsAAADwRVgEAACAL8IiAAAAfBEWAQAA4IuwCAAAAF+ERQDn3dy5c6169eqWPXt2y5Ahg+3bty9F9mPo0KHu/c/mtXv27LGU8umnn7p90M9z4bXXXnPbW7lypaVXOr7u3bsnuix++eWX87JfQGpGWARSUOALSaFp69at8Z6/9tprrXLlypae/Pbbb9ayZUvLkSOHjR8/3v79739brly5LD176qmnbNasWSm9GwBwRgiLQCpw9OhRGzlypMWCFStW2MGDB23YsGHWuXNnu/vuuy1Lliwpsi8DBw60P//8M9nfh7AIIC0jLAKpgJpkX375Zdu2bVuyvYfneeclGJ3Orl273M98+fKl9K5Y5syZXa0uUpfUcq4COIWwCKQCAwYMsBMnTiSqdvH48eOuVq5cuXKWLVs2K1OmjHu9aidDafktt9xiH330kdWuXds1+7700kvBfm5vv/22Pf7441aiRAm74IIL7B//+Ift37/fbadXr15WuHBhy507t3Xq1Cnetv3MmDHDatWq5d6rUKFCrtYwtHldzeodOnRw96+44gq3Hx07doy6ra+//to9P3v27OCyVatWuWU1a9YMW/emm26yunXrhi378MMP7eqrr3ZN3Dq+uLg4++abb07bZ1EhpUePHm7/9brbbrvNHYPW0/qR1N9Sx6DwmzdvXldehw8fDj6v1x06dMj+9a9/ufuRx6xt33PPPVakSBH3eV5++eU2efLkeO/z66+/WrNmzdzx6LPp3bt3oj8XWbNmjSunPHnyuM+1UaNGtnTp0qjrav/vu+8+K1iwoFu/ffv29vvvv4eto36NTZo0ceWkz7ts2bLuOEKdPHnSxowZ445JoVzHqO1GbsvvXFUXjOuuuy7e/mm7Om91zgb83//9n9WvX9/ts16v8/A///mPb3m8+eabdtlll7n90rqfffZZosoxMefVjh073Hlw0UUXuc+0WLFidvvtt9P/EWmXByDFTJkyxdOv4YoVK7x77rnHy549u7d169bg8w0bNvQuv/zysNd06NDBveYf//iHN378eK99+/bucbNmzcLWK126tHfJJZd4+fPn9/r16+dNnDjRW7hwobtp/erVq3v16tXznn/+ea9Hjx5ehgwZvNatW3t33XWXd9NNN7ltt2vXzq37+OOPJ/pYrrjiCu+5555z75kjRw6vTJky3u+//+7W+fjjj72uXbu69Z544gnv3//+t7d48eKo2ztx4oSXL18+76GHHgou03YzZszobvv37w+ulydPHu/hhx8Orvf666+742natKn3wgsveE8//bTbD21v48aNwfWGDBni9iVUy5Yt3TIdu8pAj6tVq+aWaf3I19aoUcNr3ry5N2HCBO/ee+91yx555JHgejrGbNmyeVdffbW7H3rMO3bs8C666CKvZMmSrjxefPFF77bbbnPb0LEGHD582Lv00kvd+aFtjxkzxqtVq5ZXtWpVt64+04SsW7fOy5Url1esWDFv2LBh3siRI72yZcu6/Vq6dGm8z7BKlSpuf3VudOvWzZX3Nddc4508edKtt3PnTndeaZ+effZZ7+WXX/Yee+wxr2LFimHvq/LInDmz16VLF3f+Pfroo24/dI4cO3bstOeqykTvvX379rDtLlq0yO3njBkzgstUjg888IA3btw4b/To0V6dOnXcOu+//37Ya7WscuXKXqFChdz2dW7o/XWu/u9//4tXFqHnS2LPq/r163t58+b1Bg4c6L3yyiveU0895V133XVuv4G0iLAIpJKw+NNPP7kvVgU3v7C4du1at76+hEMpKGn5ggULgsv0Bahlc+fODVs3EBb1hRn6hd2mTRv3RaigGEqBUttKiLZTuHBht80///wzuFxf1HqvwYMHRz3m04mLi3Nf+gEKZbplypTJ+/DDD92y1atXu+29++677vHBgwfdl7cCSigFM32Bhy6PDIurVq1yj3v16hX22o4dO/qGRYX8UHfccYdXsGDBsGUKSAr5kTp37uwC3J49e8KWK7RrXxUSReFQ7/X2228H1zl06JALWIkJi/pDImvWrO4cC9i2bZt3wQUXuBAY+dkoiIaeG88880xYGc+cOfO0n+Hnn3/u1nnzzTfDlut8jFzud66uX7/eLVcwC6VQmDt37mD5SOh90f7rfLz++uvDlmt7uq1cuTK4bNOmTS6I67OLLItACEzseaU/jPQ6hWggvaAZGkglLr74YmvXrp1NmjTJtm/fHnWdOXPmuJ99+vQJW/7QQw+5nx988EHYcjUNqqkwGjUthg4sUTOuvksjmxK1fMuWLa7524+aJNUX8YEHHgjrA6gmugoVKsTbr8RSc9/q1atdM6588cUXdvPNN7s+np9//rlbpp9q2m3QoIF7PG/ePNc03KZNGzetTeCWKVMmdywLFy5McEof0XGEevDBB31f889//jPePmvE94EDBxI8NpX1f//7X7v11lvd/dB91WemLgE69sDnrqbM0GbXnDlzWteuXe101L3h448/dk3YOscCtL277rrLlWnkvmq7oefG/fff7/p3Bs6/QH/T999/3/766y/fLglqlr/hhhvCjk1NvmoGj/wcop2rl156qfusp0+fHnY8al5Wuam5OSD0vpq5VX6B8ydSvXr13H4ElCpVyjUTqxlc248mseeV9iNr1qyuu0dkczuQVhEWgVREo3MVyvz6Lm7atMkyZsxol1xySdjyokWLui9wPR/5BexHX5Ch9MUuJUuWjLdcfcT05esn8L7qAxZJYTFyvxJLX/YqjyVLltj69etdINWya665JiwsVqpUyQoUKOAeb9iwwf28/vrr7cILLwy7KTQFBtgkVL6R5RZZ3gmVY/78+d3P0wWF3bt3u/ChPw4i91P93SSwr9ov7UNk/8po5R3tfdQHMdq6FStWdJ+t/hgIVb58+bDHCncKl4E+dw0bNrQ777zT9XlVn0UFrSlTpoT1odTnoHNG/Ssjj++PP/6I9zn4nautWrWyL7/8Mtj3VSFMr9XyUAquV155pftjReeC3ufFF1+Met5GHl8gmKqcVF7RJPa8Uh/Fp59+2vVtVB9NnavPPPOM68cIpFWZU3oHAPxNNT8aFKIA0a9fP9/1EjuRdGhtSyTViCRl+akWvPNLgx305a/BBwplCh76UldgnDBhggsnCot33HFH8DUKP6L5GxWiI6mG7Fw60/IK7Kc+78Cgn0hVq1a11Ejnn2r3NEDmvffeczVyqpEeNWqUW6ZwqePT56WBJNEoZCXmXFUo7N+/v6up1MArDczSHzBNmzYNrqNzQAORFMx0XijYqmZUAXbq1Knn5JiTcl5pP1XzqemSVDaDBg2yESNG2IIFC6xGjRrnZH+A84mwCKTC2sU33njD1U5EKl26tPvSUi2HaoUCdu7c6Wqp9HxKCLyvav9U8xJKy850v9ScV6dOHRcGFBYVEkU/FRQVRHTsCgkBGiUuCiqNGzdO8nGofDdu3BhW+/Tjjz/a2YgW7hWWNJpWzZ6n20/t17p161wADd2WyvZ09D5qso627vfff+9qUiNrk3V+hY5CVk2gukaoC0Ao1eTpNnz4cBfK2rZta9OmTbN7773XfQ6ffPKJXXXVVQn+0XI6qnHUOaCmaF155Z133nFN6qrBC1Bzvv6oUDALXa6wmFAtYagffvjBlVNkiD3T80rrq3uIbno/NacrTOt3G0hraIYGUhl9yai2SVOHRDZdBb6sNR1JqNGjRwf7CKYE1QDqS3TixIlhTZFqivvuu+/Oar8UDJctW+b6hAXCopo+FZYDgTqwXNTvTdO9aCLsaP3p/JoZA68V1U6FeuGFF+xsaJqVyEsaqkZSTbkKOgqCCe2nPnfNwRk6FYyaTFUDfTp6nxtvvNHefffdsKlbFLIV8NTXU+UVStsNLTs156o7gKbeCTSxR9acKgxJ4PPXVXoUhDXNUyRtKymXeFTtomosNaWQ+glGNkHrGBWiQ/sb6lj9JkJXt4bQvoxqhlf5qJz8aooTe17pczly5Ei832n9YZCUqY6A1ISaRSAVeuyxx1xzl2qDNEddQLVq1VyTpb7M9WWrvmPLly93c/iptiXanHTng5r8FNzU1077pEEACiNjx451c+hpTsAzpSComit9oYeGQtUmKlBr+5rPLkBf6Ao3Giyk+Rhbt27taos2b97sBtqopmvcuHFR30uDHhTgFMY1SEW1ZosWLXK1TnKm15HWdlXLplBfvHhxV1umQRHqm6oQrPtdunRxfS/37t3rgozW133Rc9pnDUrSXJNqZtX5oZqwxHjyySfdAA0FQw3eUZOpyk7hRf3pIh07dszNw6jAp3NQ4VmvVVOv6HzTMjX/KwjpijyaVF5lH/iDRueB5lRU8+vatWtdENN5olo2NSnr3AgdsJMQ7cfDDz/sbuqPGFmzpz9GVLZqmtagHfUf1KUk1c9T83VG0vyNCn+aT1M1kYE/DtQH009izyudK4Gy0+epsp45c6b7fdBrgDQppYdjA7EsoWlkAvMpRs6z+Ndff7l5DzVPXpYsWdwcff379/eOHDkStp6mI9HUM5ECU+eEzlGX0L4EpojZvXv3aY9n+vTpbt5Bzd9XoEABr23btt6vv/6a6GOO5sCBA26qHE3zcvz48eDyN954IzgfYjQ6ziZNmrhpTTQtSrly5dwUOKFTpkSbZ1FT0mhuQe2/pmfRtDOBKVw0P+HpyiXa/Hzff/+9m6JGc/npudBpdDRnod5Pn6M+z6JFi3qNGjXyJk2aFLZdTe+iORhz5szp5gjs2bNncBqa002dE5hiSOWhY9I2NO9f5ByXgX3XfICaD1PzHmp9fY6//fZb2LY01VKpUqXcZ61pk2655Zawsg3QcWgqHh27PkPN4ai5IjV1z+nO1VBXXXVV1GmjAl599VWvfPnybn8qVKjgjiXa56vHKm+dP4H1dc5GlmG0zzEx55WmQdL2tQ+aMknr1a1bN2zaIyCtyaB/UjqwAkBqppoxDUxQfzP1ywOAWEKfRQAIEe2axGqW1kCQ0IE0ABAr6LMIACHUh0/9AtX/U/3NNEhHN01UHTlqGABiAc3QABBCA0E00OHbb791U8Zoyh4NatCgo3M9RyMApAWERQAAAPiizyIAAAB8ERYBAADgK2Y64OgSXroCgmbRP9OJdQEAANIL9UTUpPq6WIBmfLBYD4sKioxkBAAACKcrZIVeCStmw6JqFAMFEnkdVAAAgFhz4MABV5EWyEgW62Ex0PSsoEhYBAAAOOV03fMY4AIAAABfhEUAAAD4IiwCAADAF2ERAAAAvgiLAAAA8BUzo6FTizL9Poi6/JeRcZYW9/1s9vtcb+98oRyS51zmdyP5tne+UA7Js+/8biTf9s6XMmn4MxRqFgEAAOCLsAgAAABfhEUAAAD4IiwCAADAF2ERAAAAvgiLAAAA8EVYBAAAgC/CIgAAAHwRFgEAAOCLsAgAAABfhEUAAAD4IiwCAADAV2b/p3A20urFzs/1BdJjqRzE73jT+kXkkyqhzz2Wzgl+N/7G78Yp/G6c2TkRi+WQmhAWU5HUcMKk1n043/uRGvbBbz/Yh9SzD+d7P1LDPvjtB/vAPqTkfqSGfUgtn0dyoBkaAAAAvgiLAAAA8EVYBAAAgC/CIgAAAHwRFgEAAOCLsAgAAABfhEUAAAD4IiwCAADAF2ERAAAAvgiLAAAA8EVYBAAAgC/CIgAAAHwRFgEAAOCLsAgAAABfhEUAAAD4IiwCAADAF2ERAAAAvgiLAAAA8EVYBAAAgC/CIgAAAHwRFgEAAOCLsAgAAABfhEUAAAD4IiwCAADAF2ERAAAA5zYsjh8/3sqUKWPZs2e3unXr2vLlyxNcf8aMGVahQgW3fpUqVWzOnDlhz7/zzjt24403WsGCBS1Dhgy2du3aeNs4cuSIdevWza2TO3duu/POO23nzp1nsvsAAABIrrA4ffp069Onjw0ZMsRWr15t1apVsyZNmtiuXbuirr948WJr06aNde7c2dasWWPNmjVzt3Xr1gXXOXTokDVo0MCefvpp3/ft3bu3vffeey54Llq0yLZt22bNmzdP6u4DAAAgOcPi6NGjrUuXLtapUyerVKmSTZw40XLmzGmTJ0+Ouv7YsWOtadOm1rdvX6tYsaINGzbMatasaePGjQuu065dOxs8eLA1btw46jb2799vr776qnvv66+/3mrVqmVTpkxxQXTp0qVJPQQAAAAkR1g8duyYrVq1KizUZcyY0T1esmRJ1NdoeWQIVE2k3/rR6D3/+uuvsO2oWbtUqVK+2zl69KgdOHAg7AYAAIBkDIt79uyxEydOWJEiRcKW6/GOHTuivkbLk7K+3zayZs1q+fLlS/R2RowYYXnz5g3eSpYsmej3AwAAQDofDd2/f3/XfB24bdmyJaV3CQAAIM3JnJSVCxUqZJkyZYo3ClmPixYtGvU1Wp6U9f22oSbwffv2hdUuJrSdbNmyuRsAAADOU82imoI1uGT+/PnBZSdPnnSP69WrF/U1Wh66vsybN893/Wj0nlmyZAnbzvr1623z5s1J2g4AAACSsWZRNG1Ohw4drHbt2lanTh0bM2aMm/pGo6Olffv2VqJECddnUHr27GkNGza0UaNGWVxcnE2bNs1WrlxpkyZNCm5z7969LvhpOpxAEBTVGuqmPoeaekfvXaBAAcuTJ489+OCDLiheeeWV56osAAAAcLZhsVWrVrZ792431Y0Gl1SvXt3mzp0bHMSi0KcR0gH169e3qVOn2sCBA23AgAFWvnx5mzVrllWuXDm4zuzZs4NhU1q3bu1+ai7HoUOHuvvPPfec264m49ZIZ42onjBhQlJ3HwAAAMkZFqV79+7uFs2nn34ab1mLFi3czU/Hjh3dLSG6+ouuHKMbAAAAzo90OxoaAAAAZ4+wCAAAAF+ERQAAAPgiLAIAAMAXYREAAAC+CIsAAADwRVgEAACAL8IiAAAAfBEWAQAA4IuwCAAAAF+ERQAAAPgiLAIAAMAXYREAAAC+CIsAAADwRVgEAACAL8IiAAAAfBEWAQAA4IuwCAAAAF+ERQAAAPgiLAIAAMAXYREAAAC+CIsAAADwRVgEAACAL8IiAAAAfBEWAQAA4IuwCAAAAF+ERQAAAPgiLAIAAMAXYREAAAC+CIsAAADwRVgEAACAL8IiAAAAfBEWAQAA4Cuz/1NITcr0+yDesl9GxqWa7Z0vlIP/fp/Nvp/r7Z1PnBOnUA6n8LvxN86JUyiHs0fNIgAAAHwRFgEAAOCLsAgAAABfhEUAAAD4IiwCAADAF2ERAAAAvgiLAAAA8EVYBAAAgC/CIgAAAHwRFgEAAOCLsAgAAABfhEUAAACc27A4fvx4K1OmjGXPnt3q1q1ry5cvT3D9GTNmWIUKFdz6VapUsTlz5oQ973meDR482IoVK2Y5cuSwxo0b24YNG8LW0ftlyJAh7DZy5Mgz2X0AAAAkV1icPn269enTx4YMGWKrV6+2atWqWZMmTWzXrl1R11+8eLG1adPGOnfubGvWrLFmzZq527p164LrPPPMM/b888/bxIkTbdmyZZYrVy63zSNHjoRt64knnrDt27cHbw8++GBSdx8AAADJGRZHjx5tXbp0sU6dOlmlSpVcwMuZM6dNnjw56vpjx461pk2bWt++fa1ixYo2bNgwq1mzpo0bNy5YqzhmzBgbOHCg3X777Va1alV7/fXXbdu2bTZr1qywbV1wwQVWtGjR4E2hEgAAAKkkLB47dsxWrVrlmomDG8iY0T1esmRJ1Ndoeej6olrDwPobN260HTt2hK2TN29e17wduU01OxcsWNBq1Khhzz77rB0/ftx3X48ePWoHDhwIuwEAACBpMidl5T179tiJEyesSJEiYcv1+Pvvv4/6GgXBaOtreeD5wDK/daRHjx6uRrJAgQKuabt///6uKVo1ndGMGDHCHn/88aQcHgAAAM4mLKYk9ZMMUFN11qxZ7b777nOhMFu2bPHWV5gMfY1qFkuWLGnpTZl+H0Rd/svIuKjPaXl65Xe8CZVRepTQ5x5L5wS/G3/jd+MUfjdO4Xfjb/xuJEMzdKFChSxTpky2c+fOsOV6rD6E0Wh5QusHfiZlm6JmajVD//LLL1GfV4DMkydP2A0AAADJGBZVm1erVi2bP39+cNnJkyfd43r16kV9jZaHri/z5s0Lrl+2bFkXCkPXUS2gRkX7bVPWrl3r+ksWLlw4KYcAAACA5GyGVtNuhw4drHbt2lanTh03kvnQoUNudLS0b9/eSpQo4ZqHpWfPntawYUMbNWqUxcXF2bRp02zlypU2adIk97zmS+zVq5c9+eSTVr58eRceBw0aZMWLF3dT7IgGuig8XnfddW5EtB737t3b7r77bsufP39SDwEAAADJFRZbtWplu3fvdpNoawBK9erVbe7cucEBKps3b3Y1fgH169e3qVOnuqlxBgwY4AKhpsSpXLlycJ1HHnnEBc6uXbvavn37rEGDBm6bmsQ70KSskDl06FA3ylmBUmExtE8iAAAAUskAl+7du7tbNJ9++mm8ZS1atHA3P6pd1ITbukWjUdBLly49k10FAADAWeDa0AAAAPBFWAQAAIAvwiIAAAB8ERYBAADgi7AIAAAAX4RFAAAA+CIsAgAAwBdhEQAAAL4IiwAAAPBFWAQAAIAvwiIAAAB8ERYBAADgi7AIAAAAX4RFAAAA+CIsAgAAwBdhEQAAAL4IiwAAAPBFWAQAAIAvwiIAAAB8ERYBAADgi7AIAAAAX4RFAAAA+CIsAgAAwBdhEQAAAL4IiwAAAPBFWAQAAIAvwiIAAAB8ERYBAADgi7AIAAAAX4RFAAAA+CIsAgAAwBdhEQAAAL4IiwAAAPBFWAQAAIAvwiIAAAB8ERYBAADgi7AIAAAAX4RFAAAA+CIsAgAAwBdhEQAAAL4IiwAAAPBFWAQAAIAvwiIAAAB8ERYBAADgi7AIAAAAX4RFAAAAnNuwOH78eCtTpoxlz57d6tata8uXL09w/RkzZliFChXc+lWqVLE5c+aEPe95ng0ePNiKFStmOXLksMaNG9uGDRvC1tm7d6+1bdvW8uTJY/ny5bPOnTvbH3/8cSa7DwAAgOQKi9OnT7c+ffrYkCFDbPXq1VatWjVr0qSJ7dq1K+r6ixcvtjZt2rhwt2bNGmvWrJm7rVu3LrjOM888Y88//7xNnDjRli1bZrly5XLbPHLkSHAdBcVvvvnG5s2bZ++//7599tln1rVr16TuPgAAAJIzLI4ePdq6dOlinTp1skqVKrmAlzNnTps8eXLU9ceOHWtNmza1vn37WsWKFW3YsGFWs2ZNGzduXLBWccyYMTZw4EC7/fbbrWrVqvb666/btm3bbNasWW6d7777zubOnWuvvPKKq8ls0KCBvfDCCzZt2jS3HgAAAJJH5qSsfOzYMVu1apX1798/uCxjxoyu2XjJkiVRX6PlqokMpVrDQBDcuHGj7dixw20jIG/evC4U6rWtW7d2P9X0XLt27eA6Wl/vrZrIO+64I977Hj161N0C9u/f734eOHDAzoeTRw/HW6b3jrY8oecC+5uS20sN+3C+t5ca9iG1by817ENq315q2IeU3l5q2IfUvr3UsA+pfXupYR9SYnvJLfA+qrhLkJcEW7du1da8xYsXhy3v27evV6dOnaivyZIlizd16tSwZePHj/cKFy7s7n/55Zdum9u2bQtbp0WLFl7Lli3d/eHDh3uXXnppvG1feOGF3oQJE6K+75AhQ9x2uXHjxo0bN27cuJnvbcuWLQnmvyTVLKYlqv0MrdE8efKkGyRTsGBBy5Ahw3lL7CVLlrQtW7a4gTmxjLI4hXI4hXL4G2VxCuVwCuXwN8oi+ctBNYoHDx604sWLJ7heksJioUKFLFOmTLZz586w5XpctGjRqK/R8oTWD/zUMo2GDl2nevXqwXUiB9AcP37chT+/982WLZu7hVJTdkrQhxvLJ3ooyuIUyuEUyuFvlMUplMMplMPfKIvkLQd1/TunA1yyZs1qtWrVsvnz54fV2OlxvXr1or5Gy0PXF41oDqxftmxZF/hC11GKVl/EwDr6uW/fPtdfMmDBggXuvdW3EQAAAMkjyc3Qatrt0KGDG2xSp04dN5L50KFDbnS0tG/f3kqUKGEjRoxwj3v27GkNGza0UaNGWVxcnBvBvHLlSps0aZJ7Xk3CvXr1sieffNLKly/vwuOgQYNclaim2BGNotaIao3C1ujrv/76y7p37+4Gv5yu6hQAAADnMSy2atXKdu/e7SbR1ihmNRVrWpsiRYq45zdv3uxGKQfUr1/fpk6d6qbGGTBggAuEGglduXLl4DqPPPKIC5yaN1E1iJoaR9vUJN4Bb775pguIjRo1ctu/88473dyMqZmawTUfZWRzeCyiLE6hHE6hHP5GWZxCOZxCOfyNskg95ZBBo1xS7N0BAACQqnFtaAAAAPgiLAIAAMAXYREAAAC+CIsAAADwRVhMRuPHj7cyZcq4Ud2aD3L58uWWnn322Wd26623uumMNCVS4PrfARpLpVH0mnw9R44c7vreGzZssPRG00ZdccUVdsEFF1jhwoXdFFDr168PW+fIkSPWrVs3d0Wh3Llzu9H9kZPXpwcvvviiVa1aNTiZrOZM/fDDD2OuHCKNHDkyOG1YLJXF0KFD3XGH3ipUqBBTZRBq69atdvfdd7vj1f+JVapUcVPLxdL/mfqOjDwndNN5EEvnxIkTJ9y0gZo+UJ91uXLlbNiwYWHXbE7R8yHBiwHijE2bNs3LmjWrN3nyZO+bb77xunTp4uXLl8/buXOnl17NmTPHe+yxx7x33nnHXWty5syZYc+PHDnSy5s3rzdr1izvq6++8m677TavbNmy3p9//umlJ02aNPGmTJnirVu3zlu7dq138803e6VKlfL++OOP4Dr//Oc/vZIlS3rz58/3Vq5c6V155ZVe/fr1vfRm9uzZ3gcffOD98MMP3vr1670BAwa468WrbGKpHEItX77cK1OmjFe1alWvZ8+eweWxUBZDhgzxLr/8cm/79u3B2+7du2OqDAL27t3rlS5d2uvYsaO3bNky7+eff/Y++ugj78cff4yp/zN37doVdj7MmzfPfX8sXLgwps6J4cOHewULFvTef/99b+PGjd6MGTO83Llze2PHjk0V5wNhMZnUqVPH69atW/DxiRMnvOLFi3sjRozwYkFkWDx58qRXtGhR79lnnw0u27dvn5ctWzbvrbfe8tIz/Weo8li0aFHwuBWY9J9BwHfffefWWbJkiZfe5c+f33vllVdishwOHjzolS9f3n0hNmzYMBgWY6UsFBarVasW9blYKYOARx991GvQoIHv87H6f6Z+J8qVK+eOP5bOibi4OO+ee+4JW9a8eXOvbdu2qeJ8oBk6GRw7dsxdmlBVxAGaSFyPlyxZYrFo48aNbhL30DLR9SjVPJ/ey2T//v3uZ4ECBdxPnRu6ClFoWagprlSpUum6LNTMois4aQJ+NUfHYjmoOU1Xsgo9ZomlslCzmbqqXHzxxda2bVt3IYdYKwOZPXu2uxJaixYtXHeVGjVq2MsvvxzT/2fqu/ONN96we+65xzVFx9I5Ub9+fXfZ4x9++ME9/uqrr+yLL76wm266KVWcD0m+ggtOb8+ePe6LMXBVmwA9/v777y0W6SSXaGUSeC490vXL1S/tqquuCl61SMer66zny5cvJsrif//7nwuH6nukPkczZ860SpUq2dq1a2OqHBSUV69ebStWrIj3XKycE/pie+211+yyyy6z7du32+OPP25XX321rVu3LmbKIODnn392fXp1CV1d3UznRY8ePVwZ6JK6sfh/pvq56ypuHTt2dI9j6Zzo16+fHThwwIXhTJkyuQwxfPhw9weVpPT5QFgEkrkmSV+E+gsxVikYKBiqhvU///mP+yJctGiRxZItW7ZYz549bd68eWGXMY01gVoS0cAnhcfSpUvb22+/7TrsxxL9Iamaxaeeeso9Vs2i/q+YOHGi+x2JRa+++qo7R1TzHGvefvttd1ljXR758ssvd/9nqqJBZZEazgeaoZNBoUKF3F8GkSO29Lho0aIWiwLHHUtlomuZv//++7Zw4UK76KKLgst1vGpu0V/QsVAWqhm45JJLrFatWm6keLVq1Wzs2LExVQ5qTtu1a5fVrFnTMmfO7G4KzLq+ve6rdiBWyiKUaowuvfRS+/HHH2PqfBCNaFUNe6iKFSsGm+Vj7f/MTZs22SeffGL33ntvcFksnRN9+/Z1tYutW7d2o+LbtWtnvXv3dv9npobzgbCYTF+O+mJU/4PQvyL1WM1xsUjTAeiEDi0TVbkvW7Ys3ZWJxvcoKKq5dcGCBe7YQ+ncyJIlS1hZaGodfUmkt7KIRr8LR48ejalyaNSokWuOV21B4KZaJTUxBe7HSlmE+uOPP+ynn35ywSmWzgdR15TIKbXUX001rbH2f6ZMmTLF9d1Un96AWDonDh8+7MY2hFKlk/6/TBXnQ7IPoYnhqXM0Sum1117zvv32W69r165u6pwdO3Z46ZVGeq5Zs8bddGqNHj3a3d+0aVNw2L/K4N133/W+/vpr7/bbb09300DI/fff76Y3+PTTT8OmhDh8+HBwHU0Hoel0FixY4KaDqFevnrulN/369XOjwDUVhD5zPc6QIYP38ccfx1Q5RBM6GjpWyuKhhx5yvxc6H7788kuvcePGXqFChdyMAbFSBqFTKGXOnNlNmbJhwwbvzTff9HLmzOm98cYbwXVi5f9MzRaiz10jxCPFyjnRoUMHr0SJEsGpczQFnX43HnnkkVRxPhAWk9ELL7zgTnLNt6ipdJYuXeqlZ5oXSyEx8qZfgsDQ/0GDBnlFihRxQbpRo0Zu7r30JloZ6Ka5FwP0y/3AAw+4aWT0BXHHHXe4QJneaCoIzSWn34ELL7zQfeaBoBhL5ZCYsBgLZdGqVSuvWLFi7nzQF6Meh84rGAtlEOq9997zKleu7P4/rFChgjdp0qSw52Pl/0zNL6n/I6MdW6ycEwcOHHD/HygzZM+e3bv44ovdvMVHjx5NFedDBv2T/PWXAAAASIvoswgAAABfhEUAAAD4IiwCAADAF2ERAAAAvgiLAAAA8EVYBAAAgC/CIgAAAHwRFgEAAOCLsAggJpQpU8bGjBljadG1115rvXr1SrfHByB1IywCSNO2bNli99xzjxUvXtyyZs1qpUuXtp49e9pvv/2W0ruWpj377LN21113uftTp06166+/PqV3CUAKISwCSLN+/vlnq127tm3YsMHeeust+/HHH23ixIk2f/58q1evnu3duzfF9u3EiRN28uRJS6uWLFliV111lbv/+eefB+8DiD2ERQBpVrdu3Vxt4scff2wNGza0UqVK2U033WSffPKJbd261R577LGw9Q8ePGht2rSxXLlyWYkSJWz8+PHB5zzPs6FDh7ptZMuWzdVU9ujRI/j80aNH7eGHH3av0+vr1q1rn376afD51157zfLly2ezZ8+2SpUquW288sorlj17dtu3b1/YfqjmM1BTpxpQ7ZO2mzNnTqtSpYoLvpGOHz9u3bt3t7x581qhQoVs0KBBbp/96D3vvfdeu/DCCy1Pnjzu/b766qszCotffPEFYRGIYYRFAGmSag0/+ugje+CBByxHjhxhzxUtWtTatm1r06dPDwtUalqtVq2arVmzxvr16+dC27x589xz//3vf+25556zl156ydVUzpo1ywW3AAU1Bahp06bZ119/bS1atLCmTZu6dQMOHz5sTz/9tAuJ33zzjdsHBUhtO7TGUful5+TIkSNWq1Yt++CDD2zdunXWtWtXa9eunS1fvjzsmP71r39Z5syZ3fKxY8fa6NGj3fv40f7t2rXLPvzwQ1u1apXVrFnTGjVqlGBt68iRI93+6rZjxw4XwHVf+9WyZUt3X8ERQIzxACANWrp0qVKgN3PmzKjPjx492j2/c+dO97h06dJe06ZNw9Zp1aqVd9NNN7n7o0aN8i699FLv2LFj8ba1adMmL1OmTN7WrVvDljdq1Mjr37+/uz9lyhT3fmvXrg1bp2fPnt71118ffPzRRx952bJl837//XffY4uLi/Meeuih4OOGDRt6FStW9E6ePBlc9uijj7plATq+5557zt3//PPPvTx58nhHjhwJ2265cuW8l156yfd9tU8bN270hgwZ4jVp0sTdHz9+vHfFFVe4+7r9+eefvq8HkD5RswggTUuoKTaS+jFGPv7uu++CNXF//vmnXXzxxdalSxebOXOma/qV//3vf65G8NJLL7XcuXMHb4sWLbKffvopuD01iVetWjXsPVSDqObqbdu2ucdvvvmmxcXFuVo60XaHDRvmajELFCjgtqsa082bN4dt58orr7QMGTKE7btqNfX6SGpu/uOPP6xgwYJh+7tx48aw/Y2kfdKoatVe3nnnne6+amFvu+02d183NasDiC2ZU3oHAOBMXHLJJS48Kezdcccd8Z7X8vz587s+e4lRsmRJW79+vevvqKZpNW+r2VqBUMErU6ZMrjlXP0MphAWoOTw00MkVV1xh5cqVc83X999/vwuh6t8YoPdQs7KmvVFgVH9ITZNz7NgxO1Pa32LFioX1qQwIhNRIGsSi/p6B5nS9tnfv3i5AZ8mSxTVRDxgwwN0AxBbCIoA0SbVmN9xwg02YMMGFmtB+i+pvpxq89u3bh4W3pUuXhm1DjytWrBh8rG3ceuut7qbBMxUqVHC1ijVq1HA1eOoDePXVVyd5X1W7qP256KKLLGPGjK5mMeDLL7+022+/3e6++273WCOof/jhBzdIJtSyZcvi7Xv58uXjhVdR/0SVgfo4qjYwMTSqfO3atS4QP/LII25EuWo3Vau4evVqt9+q+QQQe2iGBpBmjRs3zo1SbtKkiX322WduzsW5c+e6EKnRxcOHDw9bX8HsmWeecWFMI6FnzJjhBrmIavteffVVN5hDU/K88cYbLjxq3kY1PyvwKXy+8847rjlXTbUjRoxwA1NOR69V4NL+/OMf/3AjpQMU+FSTuXjxYlcbet9999nOnTvjbUPBrU+fPq72U6OlX3jhheC+R2rcuLFrpm7WrJkbKf7LL7+47Wt0+MqVK6O+Rseq2lodmyYB1/1ff/3VjYLW8esxYRGITYRFAGmWgpbCj/oZarSumns1mvi6665zI5cjw81DDz3k1ldN4ZNPPulGFCtoBppnX375ZReO1O9QzdHvvfeeq8GUKVOmuLCobVx22WUuiK1YscJNtXM6Clp16tRxo6gDo6ADBg4c6GoCtR8KaRrJrW1H0nurSVjbUa2ngqKONRrVps6ZM8euueYa69Spkwt7rVu3tk2bNlmRIkUS3Fc1P+t1oib4wH0AsSuDRrmk9E4AAAAgdaJmEQAAAL4IiwAAAPBFWAQAAIAvwiIAAAB8ERYBAADgi7AIAAAAX4RFAAAA+CIsAgAAwBdhEQAAAL4IiwAAAPBFWASQJsydO9eqV69u2bNnd9c+3rdvX4rsx9ChQ937n81r9+zZYylF137WPuhnanE2ZVqmTBm75ZZb0uRxA2kFYRFIY1577TX3pafQtHXr1njPX3vttVa5cmVLT3777Tdr2bKl5ciRw8aPH2///ve/LVeuXCm9W8nqqaeeslmzZqX0bgAAYRFIq44ePWojR460WLBixQo7ePCgDRs2zDp37mx33323ZcmSJUX2ZeDAgfbnn38m+/vEUlg8X2UK4MwQFoE0Sk2yL7/8sm3bti3Z3sPzvFTxJb5r1y73M1++fCm9K5Y5c2ZXq4uzd+jQIfeTMgVSN8IikEYNGDDATpw4kajaxePHj7tauXLlylm2bNlcPy+9XrWT0fp/ffTRR1a7dm3X7PvSSy8F+3u9/fbb9vjjj1uJEiXsggsusH/84x+2f/9+t51evXpZ4cKFLXfu3NapU6d42/YzY8YMq1WrlnuvQoUKuVrD0OZ1Nat36NDB3b/iiivcfnTs2DHqtr7++mv3/OzZs4PLVq1a5ZbVrFkzbN2bbrrJ6tatG7bsww8/tKuvvto1cev44uLi7Jtvvjlt/zoF6h49erj91+tuu+02dwxaT+tHUn9LHYPCb968eV15HT58OPi8Xqcg9a9//cvdjzxmbfuee+6xIkWKuM/z8ssvt8mTJ8d7n19//dWaNWvmjkefTe/evRP1ufznP/9x77lo0aJ4z+l80HPr1q0Llrn27eKLL3aBr2jRom7f1HUgWrl9++23dtddd1n+/PmtQYMGvmU6ZcoUu/76691+6xgrVapkL774ou8+f/zxx8E+rVr3nXfescRYtmyZNW3a1H0OOXPmtIYNG9qXX34Zto5qtXV+6/dD+6J9uuGGG2z16tWJeg8grcuc0jsA4MyULVvW2rdv72oX+/XrZ8WLF/dd995773XBQ+HuoYcecl+QI0aMsO+++85mzpwZtu769eutTZs2dt9991mXLl3ssssuCz6n1yjU6f1+/PFHe+GFF1xzcMaMGe333393X/pLly51/Sq1f4MHD07wGLSegpJCoLa9c+dOGzt2rPuyXrNmjQtTjz32mNuHSZMm2RNPPOG2q9Abjfpq6jWfffaZC2zy+eefu/376quv7MCBA5YnTx47efKkLV682Lp27Rp8rfpBKpQ2adLEnn76aRfeFE4UaLQvCgp+FJYUpNu1a2dXXnmlC1kKmn7U/1LHoWNW4HjllVdcANH7BvZFn1mdOnWC+xg4ZpWR3kPhqnv37nbhhRe6kKvmeR2fQk0gwDZq1Mg2b97sgqzOD213wYIFdjrad4V+HZPCU6jp06e7cBroFztv3jz7+eef3eeooKhwrc9KP3UuRIbAFi1aWPny5V0zu2qu/ajs9T76HFXz+N5779kDDzzgPrtu3bqFrbthwwZr1aqV/fOf/3SfoYKm3keDohTq/Kgs9EeD/lgZMmSIO08CIVXnjcpftF0FaJW3gqiC8BdffOF+fyL/CAHSJQ9AmjJlyhR9w3orVqzwfvrpJy9z5sxejx49gs83bNjQu/zyy4OP165d69a/9957w7bz8MMPu+ULFiwILitdurRbNnfu3LB1Fy5c6JZXrlzZO3bsWHB5mzZtvAwZMng33XRT2Pr16tVz20qItlO4cGG3zT///DO4/P3333fvNXjw4KjHfDpxcXFenTp1go+bN2/ubpkyZfI+/PBDt2z16tVue++++657fPDgQS9fvnxely5dwra1Y8cOL2/evGHLhwwZ4l4bsGrVKve4V69eYa/t2LGjW671I197zz33hK17xx13eAULFgxblitXLq9Dhw7xjq9z585esWLFvD179oQtb926tdvXw4cPu8djxoxx7/X2228H1zl06JB3ySWXuOX6TBOiz1afz/Hjx4PLtm/f7mXMmNF74okngssC7xfqrbfecu/x2WefxTt2bTdSZJn6bbdJkybexRdfHLYscM7+97//DS7bv3+/K6MaNWrEO4cDx33y5EmvfPnybpu6H/q+ZcuW9W644YbgMpVrt27dfMsKSO9ohgbSMDX9qTZLNTnbt2+Pus6cOXPczz59+oQtVw2jfPDBB2HLVeOl2rVoVJMZOrBEzbiqHVKzYygt37Jli2v+9rNy5UrXF1G1RaH91VSrVaFChXj7lVhqRlZtXaA/nGqAbr75ZtdEqdoi0U/VeAWaQVU7pqZh1ahqWpvALVOmTO5YFi5c6Pt+qr0SHUeoBx980Pc1qqmK3GfVVqlmMCEq6//+97926623uvuh+6rPTF0CAk2j+tyLFSvmapMD1MwaWpuaENXU6fMJnWpGtWuq2dNzAappDjhy5IjbF9V8SrRm2shj9xO6XR2XtqtaTtVi6nEo1ZrecccdwceqPda5qhrhHTt2RN3+2rVrXY2kmsRV9oFy1HmjGlnVTutYRbXVqo1Pzv7BQGpGWATSwUhShTK/voubNm1yzWuXXHJJ2HI1GepLUM9HhkU/pUqVCnusfl5SsmTJeMv1RRv5pR65XxLazB2gsBi5X4ml4KXyWLJkiWtSV+DRsmuuuSYsLKo5sUCBAu6xQoOo+VHNuqE39YULDLBJqHwjyy2yvBMqR/XfEzXlJ2T37t0u1OqPg8j9VDOwBPZV+6V9iGwGjlbe0QT68anZOUD3FbovvfTS4LK9e/daz549Xf9JBTztS6Ason3+CZ1fodQVoXHjxq6/pc5TbVf9bKNtN9pxBvbxl19+ibr9wGeuZuvIslS3APXtDLzPM8884/po6jxX07S6Wyi0ArGCPotAOqhd1KAQBQj1JfST2EmPQ2t0IqmmLSnLE+qTllw0MEc1laoZUihTX0AFBwXGCRMmuBCgsBhaExWoQVKfPoXoSOozdy6daXkF9lOfd2DQT6SqVauegz00N5BDg2PUp1Xlpr6SCnDqaxjZ/1L9P/v27euCpPo6aj8VNgP7m9jzK+Cnn35ytXv6o2H06NEupGXNmtXVlj733HNRt5tUgW08++yzbr+j0bEEjlHnj8pCfzzoNepfqkE06vMIpHeERSCd1C6+8cYbwQESoUqXLu2+GFWTUrFixeByffmrlkrPp4TA+6r2TzV6obTsTPdLoUK1PwqECov6khf9VFB888033bGrpjEgMHhEwVK1WUk9DpXvxo0b3cCNAA0AOhvRwr1qvTTaWqPgT7ef2i/VhimAhm5LZZtYam7WwKj58+e7wRzaVmgTtGpC9ZxGyIcOZgrU2p0pDWbRZ6VR7aG1sH7dAVTWkcf5ww8/uJ9+A5MCn7marBPzmatJX10NdFPtrQa2DB8+nLCImEAzNJAO6ItPtU2a1iSyj5b668mYMWPClqvGRhIatZvcNYAKZxMnTgybzkUjexVMzma/FAzVx0zhIhAWNa2NwnIgUAeWi/r7KTSo1uyvv/6K2vzrJ9C/U7VvoTRS/Gyo+TXykoaqkbzzzjtdv8XA1DV++6nPXX3s1M8wQCO8VQOdWApRaqpX87NuCuGhzciBGtLIGtHIcy2pom1XTcIaqRyNjjN0VL/6fr7++uuuxjBaTbFoBLR+b/7v//7P/vjjD9+yVDCPbPbWeat+komdHgpI66hZBNIJTTGjZlTVHGnKkYBq1aq5JkuFBIUPDRJYvny5qzFSM+N1112XIvurgTIKbuprp33S4JLA1DmqDdKcgGdKQVC1PhpkExoKVZuoQK3tX3TRRcHlCoqaqkWDhVRj1Lp1a1eLp2lnNNDmqquusnHjxvmGDgU4BSQNlAhMnROo2TrTax5ru5988okL9QomCmkabKO+qQrBuq+pjdT3Uv0GNZhE6+u+6DntswZ6aK5J1Yzp/NAgl6R8Rs2bN7dp06a5gR8KVqFUbipT9elTyNb8m2qmVS3r2bjxxhtdDbEG8mgKJ4U5TRGlkBZtIJe6GWjqIF3pR30nNeekziW/cCnqZ6q+iaoZ1O+LzkPtv+awVPnq2FTDqTkWda5ooJB+l9Q0rXLWe40aNeqsjhNIM1J6ODaApEloGhlNtaLnQqfOkb/++st7/PHH3ZQgWbJk8UqWLOn179/fO3LkSLxpSDT1TKTAtCMzZsxI1L4EpkLZvXv3aY9n+vTpboqTbNmyeQUKFPDatm3r/frrr4k+5mgOHDjgpsq54IILwqZ+eeONN9x22rVrF/V1Ok5NpaKpUrJnz+6VK1fOTYGzcuXKeMcWSlPSaGoV7X/u3Lm9Zs2aeevXr3frjRw58rTlEji+jRs3Bpd9//333jXXXOPlyJHDPRc6jc7OnTvd++lz1OdZtGhRr1GjRt6kSZPCtrtp0ybvtttu83LmzOkVKlTI69mzp5sWKTFT5wTMmzfPra8pkrZs2RLveX1WmvpHUw+p3Fq0aOFt27bNd9qgaOdEtDKdPXu2V7VqVfc5lClTxnv66ae9yZMnxyunwDn70UcfufV1HlWoUCHeuRo5dU7AmjVr3NRKmrpIr9X2WrZs6c2fP989f/ToUa9v375etWrV3PmkKY10f8KECYkqPyA9yKB/UjqwAkB6o6lZatSo4fqStm3bNqV3BwDOGH0WAeAsRbt+tpql1dQZOpAGANIi+iwCwFlSnz31C1T/T02zo0E6umkC7Mg5KAEgraEZGgDOkq4Ao+ljvv32WzcYQ9O9aLCMBh2d6zkaAeB8IywCAADAF30WAQAAcG7D4vjx4908Zbqklub60pxtCZkxY4a7bJPWr1KlirtkUyhdMknzahUsWNDNSaZRhJGuvfZa91zoLbEXpAcAAMCZSXJnGs3i36dPH3fVBQVFjfjTFQw0EbAmTI2ka4Zqst0RI0bYLbfcYlOnTnUTAWsC2cqVK7t1NNlrgwYN3PU3NZGsHz33xBNPBB8nZXJZXY5Ls/zrUllnOkkuAABAeqGeiJp4XhP/a/aGhFZMkjp16rjJYANOnDjhFS9e3BsxYkTU9TW5aeQkv3Xr1vXuu+++eOtqolXtkiZJjdSwYUM3oeyZ0mSy2jY3bty4cePGjRs3C96iTbgfKkk1i8eOHXPTQ/Tv3z+4TElU1w9dsmRJ1NdouWoiQ6kmctasWZZUb775ppvgVtf61GWgBg0a5Fu7qGt2hl63MzCOR5f/0mWcAAAAYtmBAwfc9F5qdU1IksLinj173EXVde3NUHr8/fffR33Njh07oq6v5Ulx1113WenSpV1V6ddff22PPvqoa/pWf8do1OytqSwiKSgSFgEAAE45Xfe8NDMBmCa3DdAgmWLFilmjRo3sp59+snLlysVbX7WfoTWagfQMAACAxEtSWCxUqJBlypTJdu7cGbZcj9U0HI2WJ2X9xNLgGvnxxx+jhsVs2bK5GwAAAM7T1DlZs2a1WrVq2fz588NGGetxvXr1or5Gy0PXD1ztwG/9xApMr6MaRgAAACSPJDdDq2m3Q4cOVrt2batTp46bOkdT33Tq1Mk93759eytRooTrMyg9e/a0hg0b2qhRoywuLs6mTZtmK1eutEmTJgW3uXfvXtu8ebOb2kbUF1FU+6ibmpo15c7NN9/s5mJUn8XevXvbNddcY1WrVj1XZQEAAICzDYutWrWy3bt32+DBg90glerVq9vcuXODg1gU+kLn6qlfv74LegMHDrQBAwZY+fLl3UjowByLMnv27GDYlNatW7ufQ4YMsaFDh7oazU8++SQYTNX38M4773TbBAAAQPKJmWtDa4BL3rx5bf/+/el2NHSZfh8E7/8yMi7e41hBOZwSetzRxGpZcE6cQjlEL4fAslgQ7bg5J2KrHA4kMhtxbWgAAAD4IiwCAADAF2ERAAAAvgiLAAAA8EVYBAAAgC/CIgAAAHwRFgEAAOCLsAgAAABfhEUAAAD4IiwCAADAF2ERAAAAvgiLAAAA8EVYBAAAgC/CIgAAAHwRFgEAAOCLsAgAAABfhEUAAAD4IiwCAADAF2ERAAAAvgiLAAAA8EVYBAAAgC/CIgAAAHwRFgEAAOCLsAgAAABfhEUAAAD4IiwCAADAF2ERAAAAvgiLAAAA8EVYBAAAgC/CIgAAAHwRFgEAAOCLsAgAAABfhEUAAAD4yuz/FFK7Mv0+CN7/ZWTcOV8/Fsohsa9JC87k802P58SZfL7psRyE/yNO4f+I6MeV1Nekl3IQfjeShppFAAAA+CIsAgAAwBdhEQAAAL4IiwAAAPBFWAQAAIAvwiIAAAB8ERYBAADgi7AIAAAAX4RFAAAAnNuwOH78eCtTpoxlz57d6tata8uXL09w/RkzZliFChXc+lWqVLE5c+aEPf/OO+/YjTfeaAULFrQMGTLY2rVr423jyJEj1q1bN7dO7ty57c4777SdO3eeye4DAAAgucLi9OnTrU+fPjZkyBBbvXq1VatWzZo0aWK7du2Kuv7ixYutTZs21rlzZ1uzZo01a9bM3datWxdc59ChQ9agQQN7+umnfd+3d+/e9t5777nguWjRItu2bZs1b948qbsPAACA5AyLo0ePti5dulinTp2sUqVKNnHiRMuZM6dNnjw56vpjx461pk2bWt++fa1ixYo2bNgwq1mzpo0bNy64Trt27Wzw4MHWuHHjqNvYv3+/vfrqq+69r7/+eqtVq5ZNmTLFBdGlS5cm9RAAAACQHGHx2LFjtmrVqrBQlzFjRvd4yZIlUV+j5ZEhUDWRfutHo/f866+/wrajZu1SpUr5bufo0aN24MCBsBsAAACSMSzu2bPHTpw4YUWKFAlbrsc7duyI+hotT8r6ftvImjWr5cuXL9HbGTFihOXNmzd4K1myZKLfDwAAAOl8NHT//v1d83XgtmXLlpTeJQAAgDQnc1JWLlSokGXKlCneKGQ9Llq0aNTXaHlS1vfbhprA9+3bF1a7mNB2smXL5m4AAAA4TzWLagrW4JL58+cHl508edI9rlevXtTXaHno+jJv3jzf9aPRe2bJkiVsO+vXr7fNmzcnaTsAAABIxppF0bQ5HTp0sNq1a1udOnVszJgxbuobjY6W9u3bW4kSJVyfQenZs6c1bNjQRo0aZXFxcTZt2jRbuXKlTZo0KbjNvXv3uuCn6XACQVBUa6ib+hxq6h29d4ECBSxPnjz24IMPuqB45ZVXJvUQAAAAkFxhsVWrVrZ792431Y0Gl1SvXt3mzp0bHMSi0KcR0gH169e3qVOn2sCBA23AgAFWvnx5mzVrllWuXDm4zuzZs4NhU1q3bu1+ai7HoUOHuvvPPfec264m49ZIZ42onjBhQlJ3HwAAAMkZFqV79+7uFs2nn34ab1mLFi3czU/Hjh3dLSG6+ouuHKMbAAAAzo90OxoaAAAAZ4+wCAAAAF+ERQAAAPgiLAIAAODcDnDB6ZXp90Hw/i8j4yxWj5tyiF8OgWWxINpxc05QDgGUg79YLQvOidR53NQsAgAAwBdhEQAAAL4IiwAAAPBFWAQAAIAvwiIAAAB8ERYBAADgi7AIAAAAX4RFAAAA+CIsAgAAwBdhEQAAAL4IiwAAAPBFWAQAAICvzP5PIdYvIB55ofvzsV8p8Z6p9fNJ7ecE5fA3fjdO4Zw4hXI4hXKwVLdfZ4qwmEqlhi+j1IL/fE6hHE6hHFLP/xGpsSw4J06hHP4Wy9+f5wrN0AAAAPBFWAQAAIAvwiIAAAB8ERYBAADgi7AIAAAAX4RFAAAA+CIsAgAAwBdhEQAAAL4IiwAAAPBFWAQAAIAvwiIAAAB8ERYBAADgi7AIAAAAX4RFAAAA+CIsAgAAwBdhEQAAAL4IiwAAAPBFWAQAAIAvwiIAAAB8ERYBAADgi7AIAAAAX4RFAAAA+CIsAgAAwBdhEQAAAL4IiwAAADi3YXH8+PFWpkwZy549u9WtW9eWL1+e4PozZsywChUquPWrVKlic+bMCXve8zwbPHiwFStWzHLkyGGNGze2DRs2hK2j98uQIUPYbeTIkWey+wAAAEiusDh9+nTr06ePDRkyxFavXm3VqlWzJk2a2K5du6Kuv3jxYmvTpo117tzZ1qxZY82aNXO3devWBdd55pln7Pnnn7eJEyfasmXLLFeuXG6bR44cCdvWE088Ydu3bw/eHnzwwaTuPgAAAJIzLI4ePdq6dOlinTp1skqVKrmAlzNnTps8eXLU9ceOHWtNmza1vn37WsWKFW3YsGFWs2ZNGzduXLBWccyYMTZw4EC7/fbbrWrVqvb666/btm3bbNasWWHbuuCCC6xo0aLBm0IlAAAAUklYPHbsmK1atco1Ewc3kDGje7xkyZKor9Hy0PVFtYaB9Tdu3Gg7duwIWydv3ryueTtym2p2LliwoNWoUcOeffZZO378uO++Hj161A4cOBB2AwAAQNJkTsrKe/bssRMnTliRIkXCluvx999/H/U1CoLR1tfywPOBZX7rSI8ePVyNZIECBVzTdv/+/V1TtGo6oxkxYoQ9/vjjSTk8AAAAnE1YTEnqJxmgpuqsWbPafffd50JhtmzZ4q2vMBn6GtUslixZ8rztLwAAQMw1QxcqVMgyZcpkO3fuDFuux+pDGI2WJ7R+4GdStilqplYz9C+//BL1eQXIPHnyhN0AAACQjGFRtXm1atWy+fPnB5edPHnSPa5Xr17U12h56Poyb9684Pply5Z1oTB0HdUCalS03zZl7dq1rr9k4cKFk3IIAAAASM5maDXtdujQwWrXrm116tRxI5kPHTrkRkdL+/btrUSJEq55WHr27GkNGza0UaNGWVxcnE2bNs1WrlxpkyZNcs9rvsRevXrZk08+aeXLl3fhcdCgQVa8eHE3xY5ooIvC43XXXedGROtx79697e6777b8+fMn9RAAAACQXGGxVatWtnv3bjeJtgagVK9e3ebOnRscoLJ582ZX4xdQv359mzp1qpsaZ8CAAS4QakqcypUrB9d55JFHXODs2rWr7du3zxo0aOC2qUm8A03KCplDhw51o5wVKBUWQ/skAgAAIJUMcOnevbu7RfPpp5/GW9aiRQt386PaRU24rVs0GgW9dOnSM9lVAAAAnAWuDQ0AAABfhEUAAACk/XkWY12Zfh+EPf5lZNw53ea52N75cq73OznK9nxIjs8vLZ4TkZ9frJbD+fjdSCv43TiFcjiF78+zR1g8T9LriZXU46Ic4q+f2NekBWcSLjgnzmz9tILfjTP/fNPjOXEmn296LIe0dlw0QwMAAMAXYREAAAC+CIsAAADwRVgEAACAL8IiAAAAfBEWAQAA4IuwCAAAAF+ERQAAAPgiLAIAAMAXYREAAAC+CIsAAADwxbWhU0hauiZkcu73mVxLOL1+fmnxnKAcku96xmmxHM7H/xFppSz43Ui+/+PTYjmk5f0WahYBAADgi7AIAAAAX4RFAAAA+CIsAgAAwBdhEQAAAL4IiwAAAPBFWAQAAIAvwiIAAAB8ERYBAADgi7AIAAAAX4RFAAAA+CIsAgAAwBdhEQAAAL4IiwAAAPBFWAQAAIAvwiIAAAB8ERYBAADgi7AIAAAAX4RFAAAA+CIsAgAAwBdhEQAAAL4IiwAAAPBFWAQAAIAvwiIAAAB8ERYBAADgi7AIAAAAX4RFAAAA+CIsAgAA4NyGxfHjx1uZMmUse/bsVrduXVu+fHmC68+YMcMqVKjg1q9SpYrNmTMn7HnP82zw4MFWrFgxy5EjhzVu3Ng2bNgQts7evXutbdu2lidPHsuXL5917tzZ/vjjjzPZfQAAACRXWJw+fbr16dPHhgwZYqtXr7Zq1apZkyZNbNeuXVHXX7x4sbVp08aFuzVr1lizZs3cbd26dcF1nnnmGXv++edt4sSJtmzZMsuVK5fb5pEjR4LrKCh+8803Nm/ePHv//ffts88+s65duyZ19wEAAJCcYXH06NHWpUsX69Spk1WqVMkFvJw5c9rkyZOjrj927Fhr2rSp9e3b1ypWrGjDhg2zmjVr2rhx44K1imPGjLGBAwfa7bffblWrVrXXX3/dtm3bZrNmzXLrfPfddzZ37lx75ZVXXE1mgwYN7IUXXrBp06a59QAAAJA8Midl5WPHjtmqVausf//+wWUZM2Z0zcZLliyJ+hotV01kKNUaBoLgxo0bbceOHW4bAXnz5nWhUK9t3bq1+6mm59q1awfX0fp6b9VE3nHHHfHe9+jRo+4WsH//fvfzwIEDdj6cPHo4eF/veTaPoznbbfKevCfvyXvynrwn75l63/N8CLyPKu4S5CXB1q1btTVv8eLFYcv79u3r1alTJ+prsmTJ4k2dOjVs2fjx473ChQu7+19++aXb5rZt28LWadGihdeyZUt3f/jw4d6ll14ab9sXXnihN2HChKjvO2TIELddbty4cePGjRs3buZ727JlS4L5L0k1i2mJaj9DazRPnjzpBskULFjQMmTIcF7SesmSJW3Lli1uUA7OHGV57lCW5w5lee5QlucOZXnuxEJZep5nBw8etOLFiye4XpLCYqFChSxTpky2c+fOsOV6XLRo0aiv0fKE1g/81DKNhg5dp3r16sF1IgfQHD9+3IU/v/fNli2bu4VSU/b5phMsvZ5k5xtlee5QlucOZXnuUJbnDmV57uRJ52Wprn/ndIBL1qxZrVatWjZ//vywGjs9rlevXtTXaHno+qIRzYH1y5Yt6wJf6DpK8+qLGFhHP/ft2+f6SwYsWLDAvbf6NgIAACB5JLkZWk27HTp0cINN6tSp40YyHzp0yI2Olvbt21uJEiVsxIgR7nHPnj2tYcOGNmrUKIuLi3MjmFeuXGmTJk1yz6tJuFevXvbkk09a+fLlXXgcNGiQqxLVFDuiUdQaUa1R2Bp9/ddff1n37t3d4JfTVZ0CAADgPIbFVq1a2e7du90k2hrFrKZiTWtTpEgR9/zmzZvdKOWA+vXr29SpU93UOAMGDHCBUCOhK1euHFznkUcecYFT8yaqBlFT42ibmsQ74M0333QBsVGjRm77d955p5ubMbVSE7jmooxsCkfSUZbnDmV57lCW5w5lee5QlucOZfm3DBrlEvIYAAAACOLa0AAAAPBFWAQAAIAvwiIAAAB8ERYBAADgi7CYTMaPH29lypRxI7o1F+Ty5ctTepdSNU21dMUVV9gFF1xghQsXdtMmrV+/PmydI0eOWLdu3dxVeHLnzu1GxEdO+I74Ro4cGZyiKoCyTLytW7fa3Xff7coqR44cVqVKFTf9V4DGCGp2CF1UQM/ruvUbNmxI0X1OjU6cOOGmRdP0aCqncuXK2bBhw8KuSUtZRvfZZ5/Zrbfe6qaK0++yZhQJlZhy00Us2rZt6yaX1gUqOnfubH/88YfFmoTKUtPyPfroo+53PFeuXG4dTQe4bds2i/WyJCwmg+nTp7v5KDXkfvXq1VatWjVr0qRJvKvQ4G+LFi1y4WXp0qVu0nb90t54441uSqWA3r1723vvvWczZsxw6+sXuHnz5im636ndihUr7KWXXrKqVauGLacsE+f333+3q666yrJkyWIffvihffvtt27O2Pz58wfXeeaZZ9w0XpoDVhcT0JeMft8VyPG3p59+2l588UUbN26cfffdd+6xyu6FF14IrkNZRqf/B/U9okqIaBJTbgo333zzjfv/9f3333ehSdPVxZqEyvLw4cPuO1t/1OjnO++84yotbrvttrD1YrIsE7xyNM5InTp1vG7dugUfnzhxwitevLg3YsSIFN2vtGTXrl3u4uaLFi1yj/ft2+dlyZLFmzFjRnCd7777zq2zZMmSFNzT1OvgwYNe+fLlvXnz5nkNGzb0evbs6ZZTlon36KOPeg0aNPB9/uTJk17RokW9Z599NrhM5ZstWzbvrbfeOk97mTbExcV599xzT9iy5s2be23btnX3KcvE0e/pzJkzg48TU27ffvute92KFSuC63z44YdehgwZvK1bt3qxKrIso1m+fLlbb9OmTTFdltQsnmPHjh1zlyVUM0CAJhHX4yVLlqTovqUl+/fvdz8LFCjgfqpMVdsYWq4VKlSwUqVKUa4+VFOrqyaFlplQlok3e/Zsd7WqFi1auO4RNWrUsJdffjn4/MaNG93FCULLUtdZVdcTyjKcLtCgy7r+8MMP7vFXX31lX3zxhd10003uMWV5ZhJTbvqp5lKdywFaX99NqolEwt9FGTJkcOUXy2WZ5Cu4IGF79uxxfXMCV7QJ0OPvv/8+xfYrLdE1v9W/Ts1/gSv96D9DXZs88AsbWq56DuF0WU01o6gZOhJlmXg///yzazpVtxJdgUrl2aNHD1d+uuxpoLyi/b5TluH69etnBw4ccH+YZMqUyf0/OXz4cNekJ5TlmUlMuemn/tgJlTlzZvfHOGXrT8346sPYpk0b1z8xlsuSsIhUWSO2bt06V+uApNuyZYu7Jrv604ReMhNn9oeLahCeeuop91g1izo31TdMYRGJ9/bbb7vLturyr5dffrmtXbvW/VGoQQSUJVIbtb60bNnSDR568cUXLdbRDH2OFSpUyP3VHDmyVI+LFi2aYvuVVuj63+owvHDhQrvooouCy1V2auLXtcNDUa7xqZlZg6lq1qzp/uLVTYNY1AFe91XjQFkmjkaXVqpUKWxZxYoVbfPmze5+oLz4fT+9vn37utrF1q1bu9Gm7dq1cwOtNBOCUJZnJjHlpp+RAyyPHz/uRvVStv5BcdOmTe6P7jz/v1YxlsuSsHiOqXmqVq1arm9OaO2EHterVy9F9y01019vCoozZ860BQsWuOk1QqlMNSI1tFw1Sk1f2pRruEaNGtn//vc/V3MTuKl2TM19gfuUZeKoK0TkFE7qc1e6dGl3X+epviBCy1JNreq7RFnGH2mqfl2h9Ie1/n8UyvLMJKbc9FN/HOoPyQD9P6uyV99GxA+Kmnrok08+cVNmhYrZskzpETbp0bRp09xItNdee82NnOratauXL18+b8eOHSm9a6nW/fff7+XNm9f79NNPve3btwdvhw8fDq7zz3/+0ytVqpS3YMECb+XKlV69evXcDacXOhpaKMvE0UjIzJkze8OHD/c2bNjgvfnmm17OnDm9N954I7jOyJEj3e/3u+++63399dfe7bff7pUtW9b7888/U3TfU5sOHTp4JUqU8N5//31v48aN3jvvvOMVKlTIe+SRR4LrUJb+MxusWbPG3fS1PXr0aHc/MEI3MeXWtGlTr0aNGt6yZcu8L774ws2U0KZNGy/WJFSWx44d82677Tbvoosu8tauXRv2XXT06NGYLkvCYjJ54YUX3Jdx1qxZ3VQ6S5cuTeldStX0SxvtNmXKlOA6+o/vgQce8PLnz+++sO+44w73S4ykh0XKMvHee+89r3Llyu4PwAoVKniTJk0Ke15TlwwaNMgrUqSIW6dRo0be+vXrU2x/U6sDBw64c1D/L2bPnt27+OKLvcceeyzsS5iyjG7hwoVR/39UAE9suf32228u0OTOndvLkyeP16lTJxecYk1CZak/Yvy+ixYuXBjTZZlB/6R07SYAAABSJ/osAgAAwBdhEQAAAL4IiwAAAPBFWAQAAIAvwiIAAAB8ERYBAADgi7AIAAAAX4RFAAAA+CIsAsBZypAhg82aNSvR6w8dOtSqV6+e4DodO3a0Zs2anYO9A4CzQ1gEkO7deuut1rRp06jPff755y7sff3112e8/e3bt9tNN91kqdmKFSusePHi7v62bdssR44cduzYsZTeLQBpAGERQLrXuXNnmzdvnv3666/xnpsyZYrVrl3bqlatmuTtBsJW0aJFLVu2bJaaLVmyxK666qpgQNYxZ82aNaV3C0AaQFgEkO7dcsstduGFF9prr70WtvyPP/6wGTNmuDD522+/WZs2baxEiRKWM2dOq1Klir311lth61977bXWvXt369WrlxUqVMiaNGkStRn60UcftUsvvdRt5+KLL7ZBgwbZX3/9FW+/XnrpJStZsqRbr2XLlrZ//37fYzh58qSNGDHCypYt62oFq1WrZv/5z38SXQaLFy8OhsUvvvgieB8AToewCCDdy5w5s7Vv396FRc/zgssVFE+cOOFC4pEjR6xWrVr2wQcf2Lp166xr167Wrl07W758edi2/vWvf7kauS+//NImTpwY9f0uuOAC917ffvutjR071l5++WV77rnnwtb58ccf7e2337b33nvP5s6da2vWrLEHHnjA9xgUFF9//XX3nt9884317t3b7r77blu0aJHvaxQK8+XL524Klo899pi7r208//zz7v7IkSOTUJIAYlEGL/R/TgBIp77//nurWLGiLVy40NUQyjXXXGOlS5e2f//73741khUqVLD/+7//c4/1ugMHDtjq1avD1lPN4syZM30HpOj106ZNs5UrVwYHuDz55JO2adMmV5MpCoxxcXG2detW16ytAS779u1zNZZHjx61AgUK2CeffGL16tULbvfee++1w4cP29SpU6O+rwLwjh073LHfddddtmrVKtu7d6/Vr1/fvvrqK8uePXswTAKAn8y+zwBAOqLQp5A0efJkF/pUs6e+e0888YR7XjWMTz31lKvtU2BTf0SFNDURh1Lt4+lMnz7d1dz99NNPrqn7+PHjlidPnrB1SpUqFQyKohCopub169e7sBhK+6pQeMMNN4Qt1z7WqFHDdz8UBsuUKeOOSQNw1ISt5uirr77alQcAJAZhEUDMUN/EBx980MaPH+8GtpQrV84aNmzonnv22Wddk/GYMWNcf8VcuXK5vomRI4a1/HQDSdq2bWuPP/6469OYN29eV6s4atSoM95vBU5RE3lowJSEBtbkzp3b/VTozZgxo7377rvueNSgpOcUGj/88MMz3i8AsYGwCCBmaBBJz549XbOt+v/df//9rglZ1Afx9ttvd/0ARbV8P/zwg1WqVClJ76GaOzVtq39ggJqbI23evNlNYROYzmbp0qUu0F122WXx1tU+KBTqNYFwmxhr1651tZqa01FN2KqxVECcMGGCC8QaKAMAp0NYBBAzVJvWqlUr69+/v+t7qH6BAeXLl3eDQBT28ufPb6NHj7adO3cmOSxqOwp1qk284oorXG2g+jNGayLu0KGD68+ofenRo4cLs5FN0IEBMw8//LAb1KIQ26BBAzdyWgFXzdvaTjSXXHKJC6FFihRxr9F+HTx40M07qUE/AJAYjIYGEHNN0b///rtrIg7U6snAgQOtZs2abrn6NCq0nckVVG677TYX6jTFjmr0FD41dU60INe8eXO7+eab7cYbb3TzPKrGz8+wYcPcdjQqWgN1NMm4gqj6ISbk008/dQN5RCOn1TeSoAggKRgNDQAAAF/ULAIAAMAXYREAAAC+CIsAAADwRVgEAACAL8IiAAAAfBEWAQAA4IuwCAAAAF+ERQAAAPgiLAIAAMAXYREAAAC+CIsAAAAwP/8PmS03sGeKIIIAAAAASUVORK5CYII=",
+      "text/plain": [
+       "<Figure size 640x480 with 2 Axes>"
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    }
+   ],
+   "source": [
+    "resp_dx.plot_norm()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "c348df35-e03e-478e-8d03-8b7656015b45",
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "source": [
+    "### Matrix pseudo-inversion\n",
+    "\n",
+    "The {py:meth}`~.ResponseMatrix.solve` method computes the singular values of the\n",
+    "weighted response matrix and its pseudo-inverse."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 10,
+   "id": "64e02b6c-2b89-4f5a-b7e1-b9e4e8af0a8b",
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "outputs": [],
+   "source": [
+    "resp_dx.solve()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "45546624-1884-4792-a9c1-d144a62488c9",
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "source": [
+    "We can plot the singular values:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 11,
+   "id": "bee36f70-ccab-4f5e-a746-c4d86cdbb3fc",
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "outputs": [
+    {
+     "data": {
+      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAi8AAAGzCAYAAADnmPfhAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjEsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvc2/+5QAAAAlwSFlzAAAPYQAAD2EBqD+naQAAJQxJREFUeJzt3Q2UV2WBP/AHIcDXIaBA4s2X1JDEggFR20DZZZWlTGtp13UBNzytg2tRtlBH0LOruFvL0XXZOJXInrXypU10JVEji3JRBgjTJlEKjA15W4PhJcGG+z/P/Z+ZZnidwRl/89zf53POBe6dO/feZ+7w+31/z9vtkGVZFgAAEnFcqS8AAKAlhBcAICnCCwCQFOEFAEiK8AIAJEV4AQCSIrwAAEkRXgCApAgvAEBShBcouIEDB4ZJkyaF9mLUqFH5kqL29rOEciW8QKJeeOGF8PGPfzwMGDAgdO3aNbznPe8Jf/zHfxzuvvvuUl8aQJvq1LaHB9rC//zP/4TRo0eH/v37hylTpoTevXuHDRs2hGeffTbcdddd4YYbbmjYd82aNeG443xOAYpDeIEE3XbbbaGioiJUV1eHbt26Nfnali1bmqx36dIlFNXvf//7sH///tC5c+dSXwrwNvJxDBL0y1/+Mpx77rkHBZfo3e9+9xH7aSxYsCB06NAhPPPMM2HatGnhXe96VzjxxBPDxz72sbB169Ym3xuDwS233BL69OkTTjjhhLy2p6am5qBjxn3iMQ9Uf67169cftiz79u0LM2fODEOHDs0DWbyWD33oQ+Hpp59usl88RjzWV77ylXDnnXeGM844Iw9m8XoOZfDgwfn1HiiWKTaxxSa3evGYF154YejRo0c4/vjj82v5zne+c9hrPtZyP/7443nZYhlPPvnkMG7cuPDzn/+8yT6bNm0KkydPDn379s3Ld+qpp4aPfvSjR/wZQrlR8wIJiv1cli1bFl588cX8TfpYxKald77znWHWrFn5G2MMBFOnTg0PPPBAwz4zZswI//zP/xzGjx8fxo4dG55//vn87zfeeKPVylJbWxu+8Y1vhL/4i7/Im8B27twZ7rnnnvw8y5cvD+eff36T/e+99978/Nddd13+5t69e/dDHnfChAl5uIhhIDar1fvJT34SNm7cGD75yU82bItNbR/5yEfC1VdfnYep+++/P3ziE58Ijz32WB4wWsN//ud/hokTJ+bl+qd/+qewZ8+e8NWvfjVcfPHF4ac//WkeCKOrrroqDzTx/sRtsSbtqaeeCr/+9a8b9oGylwHJefLJJ7OOHTvmy8iRI7MvfOEL2RNPPJHt27fvoH0HDBiQTZw4sWH93nvvzeJ//TFjxmT79+9v2P7Zz342P9727dvz9U2bNmWdOnXKrrjiiibHu+WWW/Lvb3zMWbNm5dsOVH+udevWNWz78Ic/nC/1fv/732d79+5t8n2//e1vs169emXXXnttw7Z4jHisU045JduyZctRf0Zr1qzJ97/77rubbL/++uuzk046KduzZ0/Dtsb/juLPcfDgwdkll1xyxJ9lc8u9c+fOrFu3btmUKVOa7Bd/xhUVFQ3bY7nj9335y18+avmgnGk2ggTFUUWx5iXWFsTakFg7Ej/Rx+aQRx99tFnHiDUXjZs8YnNGXV1dePXVV/P1JUuW5H1Krr/++ibf17gzcGvo2LFjQ5+V2KTz+uuv5+cdNmxYWLVq1UH7x5qJ2NR1NGeddVZea9O4JimWLzYHxZqk2DxUr/G/f/vb34YdO3bkP49Dnf9YxJqT7du357VL27Zta1hi2UeMGNHQRBavI/4sfvjDH+bXARya8AKJqqysDN/97nfzN7nYvBKbeGKTS+zLcbh+II3FkUqNxSakqP5Nsz7EnHnmmU32i8009fu2lv/4j/8I5513Xj7kO/Y7ieFk0aJFeYg40Gmnndbs48amo9i35ze/+U2+HkNBbIaJ2xuLzUMXXHBBfv5Yvnj+2KRzqPMfi1deeSX/+5JLLsmP3Xh58sknGzpZx2aw2KQU+8b06tUr/NEf/VEeTGPTF/AHwgskLn5Sj0Hm9ttvz99w33zzzfDQQw8d9fvip/5DybLYctEyh+q0Wl/TcTT33Xdf3vk3dsCNfV0WL16c11TEN/pYE3OgxrUkRxNDSixP/c/jwQcfzDsF/+mf/mnDPj/+8Y/zGqwYXP793/89fO9738vP/5d/+ZdH/Vk0t9z15Yj9XuKxD1weeeSRhn0/85nPhJdffjnMnj07v6abb745vO9978v7xQD/nw67UCCxqSV67bXXWqVTcLR27domtR3/93//d1CTRn1NTGwaaTwCqr725khiM87pp5+e1yI1DgOxI/FbFa97+PDhedNR7Iwcz3HFFVc0GT7+X//1X3lIeOKJJ5psjx2Dj6a55Y7BrH4k2JgxY4563Lj/5z73uXyJtTax+etf/uVf8qAHqHmBJMU+EoeqFYi1BtHZZ5/9ls9x6aWXhk6dOuW1OY3927/920H71r85L126tGHb7t278+ag5tYANS7Pc889l/fpaQ2x9iVO3jd//vy8n8mBTUbx/DE0Na4tiaOvFi5ceNRjN7fcsT/SKaeckteOxZqxA9UPUY8jkA4cyRXPEYdV7927t9llhqJT8wIJip1m4xtdnJvlnHPOyYf3xll3Yw1DHE4b5wl5q2KfixtvvDH/xB+bVWJTS+wcHPtj9OzZs0ktyZ/8yZ/kfWj+5m/+Jtx00015IIhhIfbpiEN8j+TP/uzP8hqRWJY4LHndunVh3rx5YdCgQWHXrl1vuRx//ud/Hj7/+c/nS+zPcmDNRzznnDlz8vLFpqLY/2Tu3Ll5X5+f/exnRzx2c8sdg0sMgddcc0344Ac/mA/Trt8n9u256KKL8lAYm4tiaIzXHMsfw+PDDz8cNm/e3GRoN5S9Ug93Alru8ccfz4cRn3POOfmw386dO2dnnnlmdsMNN2SbN29u1lDp6urqJvs9/fTT+fb4d+NhzDfffHPWu3fv7Pjjj8+HDv/iF7/IevTokX36059u8v0rV67MRowYkV9L//79szlz5jRrqHQcrn377bfn19mlS5fsAx/4QPbYY4/l1xy3HThU+liGEV900UX5937qU5865Nfvueee7L3vfW9+/vgzjdd9qGHQB/4sW1Lu+p/x2LFj8+HRXbt2zc4444xs0qRJ2YoVK/Kvb9u2Lauqqsqv4cQTT8z3i8d+8MEHW1xmKLIO8Y9SByggHbF/R+zr8Y//+I/hS1/6UqkvByhD+rwAh/W73/3uoG1xJt5o1KhRJbgiAH1egCOIfWjic3ouv/zycNJJJ+VT63/729/O+3rEfhoApdAua15ix71YLd34wWnA2y9OHBc7jcaJ0uL8I3FOlNiJNw4vBiiVdtnnJc6CGWcKjcMNm/NkVwCgfLTLmpfYlh7nNQAAeMvhJU7GFB9q1qdPn3yeh0NN5BTnSIhzTcRZK+NDx+JzVwAAStJhN84eOWTIkHDttdeGK6+88pAd/KZNm5ZPMhWDSxyZEGeXXLNmTT41dhSnuo5PjT1QfEBZDEVvRXyGyMaNG/Oam8M9dwQAaF9iL5bYZSTmgOOOO0rdyluZJCZ++8MPP9xk2/Dhw/NJlurV1dVlffr0yWbPnt2iY8fJnK666qqj7vfGG29kO3bsaFhqamry67JYLBaLxRKSWzZs2HDU9/5WHSodpyhfuXJlmDFjRsO2mJ7idNyt9ZySA8Unr956660Hbd+wYUM+JTcA0P7V1taGfv36NavPa6uGl/jQs/hws/hMlMbi+ksvvdTs48SwE5+hEpuo+vbtmz/OfuTIkYfcNwal2Ex1YOFjcBFeACAtzeny0S4nqfv+97/f7H3jI+wbP8YeACi2Vh0qHZ80G5+qGp+A2lhc7927d2ueCgAoU60aXjp37hyGDh0alixZ0mT0T1w/XLNPa4nDs+Mj5CsrK9v0PABAabW42WjXrl1h7dq1Devr1q0Lq1evDt27dw/9+/fP+59MnDgxDBs2LAwfPjwfKh37rkyePDm0paqqqnyJfV4qKira9FwAQELhZcWKFWH06NEN6/WdZWNgiQ9wmzBhQti6dWuYOXNm2LRpUz6ny+LFiw/qxAsAUJhnG70V9TUvO3bsMNoIAAr4/t0un20EAFD48KLDLgCUB81GAEDJaTYCAApLeAEAkiK8AABJKUx40WEXAMqDDrstNHD6oibr6+8Yd9C2+u0AQPPosAsAFJbwAgAkRXgBAJIivAAASSlMeDHaCADKQ2HCS1VVVaipqQnV1dWlvhQAoA0VJrwAAOVBeAEAkiK8AABJEV4AgKQILwBAUoQXACAphQkv5nkBgPJQmPBinhcAKA+FCS8AQHkQXgCApAgvAEBShBcAICnCCwCQFOEFAEiK8AIAJKUw4cUkdQBQHgoTXkxSBwDloTDhBQAoD8ILAJAU4QUASIrwAgAkRXgBAJIivAAASRFeAICkCC8AQFKEFwAgKcILAJCUwoQXzzYCgPJQmPDi2UYAUB4KE14AgPIgvAAASRFeAICkCC8AQFKEFwAgKcILAJAU4QUASIrwAgAkRXgBAJIivAAASRFeAICkCC8AQFKEFwAgKcILAJAU4QUASEphwsvcuXPDoEGDQmVlZakvBQBoQ4UJL1VVVaGmpiZUV1eX+lIAgDZUmPACAJQH4QUASIrwAgAkRXgBAJIivAAASRFeAICkCC8AQFKEFwAgKcILAJAU4QUASIrwAgAkRXgBAJIivAAASRFeAICkCC8AQFKEFwAgKcILAJAU4QUASIrwAgAkRXgBAJLS7sLLhg0bwqhRo8KgQYPCeeedFx566KFSXxIA0I50Cu1Mp06dwp133hnOP//8sGnTpjB06NBw+eWXhxNPPLHUlwYAtAPtLryceuqp+RL17t079OzZM7z++uvCCwBwbM1GS5cuDePHjw99+vQJHTp0CAsXLjxon7lz54aBAweGrl27hhEjRoTly5eHY7Fy5cpQV1cX+vXrd0zfDwAUT4vDy+7du8OQIUPygHIoDzzwQJg2bVqYNWtWWLVqVb7v2LFjw5YtWxr2iU1CgwcPPmjZuHFjwz6xtuWv//qvw9e+9rVjLRsAUEAtbja67LLL8uVw5syZE6ZMmRImT56cr8+bNy8sWrQozJ8/P0yfPj3ftnr16iOeY+/eveGKK67I97/wwguPum9c6tXW1rawRABA2Y422rdvX97UM2bMmD+c4Ljj8vVly5Y16xhZloVJkyaFSy65JFxzzTVH3X/27NmhoqKiYdHEBADF1qrhZdu2bXkflV69ejXZHtfjyKHmeOaZZ/Kmp9iXJjYvxeWFF1447P4zZswIO3bsaFjiUGsAoLja3Wijiy++OOzfv7/Z+3fp0iVfAIDy0Ko1L3FYc8eOHcPmzZubbI/rcdgzAEC7Ci+dO3fOJ5VbsmRJw7ZYixLXR44cGdpSHP0UZ+WtrKxs0/MAAIk1G+3atSusXbu2YX3dunX56KHu3buH/v3758OkJ06cGIYNGxaGDx+ez5Ybh1fXjz5qK1VVVfkSRxvFjrsAQDG1OLysWLEijB49umE9hpUoBpYFCxaECRMmhK1bt4aZM2fmnXRjh9vFixcf1IkXAOBtCS/xoYlxOPORTJ06NV8AAAr/VGkAgLIILzrsAkB5KEx4iZ11a2pqQnV1dakvBQBoQ4UJLwBAeRBeAICkCC8AQFKEFwAgKe3uwYxvZbRRXOJTrduLgdMXHbRt/R3jSnItAFAUhal5MdoIAMpDYcILAFAehBcAICnCCwCQFOEFAEhKYcKLZxsBQHkoTHgx2ggAykNhwgsAUB6EFwAgKcILAJAU4QUASIrwAgAkpTDhxVBpACgPhQkvhkoDQHkoTHgBAMqD8AIAJEV4AQCS0qnUF1COBk5f1GR9/R3jDtpWvx0AaErNCwCQFOEFAEiK8AIAJEV4AQCSUpjwYoZdACgPhQkvZtgFgPJQmPACAJQH4QUASIrwAgAkRXgBAJIivAAASRFeAICkCC8AQFI8Vbqd87RpAGhKzQsAkJTChBePBwCA8lCY8OLxAABQHgoTXgCA8iC8AABJEV4AgKQILwBAUoQXACApwgsAkBThBQBIivACACRFeAEAkiK8AABJEV4AgKQILwBAUoQXACApwgsAkBThBQBISmHCy9y5c8OgQYNCZWVlqS8FAGhDnUJBVFVV5UttbW2oqKgIRTdw+qIm6+vvGFeyawGAt1Nhal4AgPJQmJoXDq6Nqa+RUUsDQJGoeQEAkiK8AABJEV4AgKQILwBAUoQXACApwgsAkBThBQBIinleytih5n853FwxANBeqHkBAJIivAAASRFeAICkCC8AQFKEFwAgKcILAJAU4QUASIrwAgAkRXgBAJIivAAASRFeAICktLvwsn379jBs2LBw/vnnh8GDB4evf/3rpb4kAKAdaXcPZjz55JPD0qVLwwknnBB2796dB5grr7wy9OjRo9SXVtYO9RDHI20HgLIJLx07dsyDS7R3796QZVm+kI7DPZla0AGgJM1GsVZk/PjxoU+fPqFDhw5h4cKFB+0zd+7cMHDgwNC1a9cwYsSIsHz58hY3HQ0ZMiT07ds33HTTTaFnz54tvUwSEkNN4wUAWjW8xKacGCxiQDmUBx54IEybNi3MmjUrrFq1Kt937NixYcuWLQ371PdnOXDZuHFj/vVu3bqF559/Pqxbty5861vfCps3bz7s9cTamdra2iYLAFBcLW42uuyyy/LlcObMmROmTJkSJk+enK/PmzcvLFq0KMyfPz9Mnz4937Z69epmnatXr155+Pnxj38cPv7xjx9yn9mzZ4dbb721pcUAABLVqqON9u3bF1auXBnGjBnzhxMcd1y+vmzZsmYdI9ay7Ny5M//3jh078maqs88++7D7z5gxI9+vftmwYUMrlAQAKIsOu9u2bQt1dXV5jUljcf2ll15q1jFeffXVcN111zV01L3hhhvC+9///sPu36VLl3wBAMpDuxttNHz48GY3KwEA5adVm43iqKA41PnADrZxvXfv3q15KgCgTLVqeOncuXMYOnRoWLJkScO2/fv35+sjR44MbSmOfho0aFCorKxs0/MAAIk1G+3atSusXbu2YT0OZ47NPN27dw/9+/fPh0lPnDgxn+I/NgHdeeed+fDq+tFHbaWqqipf4lDpioqKNj0Xbc9EdwC0WnhZsWJFGD16dMN6DCtRDCwLFiwIEyZMCFu3bg0zZ84MmzZtyud0Wbx48UGdeAEA3pbwMmrUqKNO1z916tR8gVLW0gBQTO3uqdLHSp8XACgPhQkvsb9LTU1NqK6uLvWlAABtqDDhBQAoD8ILAJAU4QUASIrwAgAkpTDhxWgjACgPhQkvRhsBQHkoTHgBAMqD8AIAJEV4AQCSIrwAAMV+MGN7Hm0Ul7q6ulJfCu38gY0HbvcQR4C0FKbmxWgjACgPhQkvAEB5EF4AgKQUps8LtEX/GADaHzUvAEBShBcAICmajeAIDLcGaH8KU/PiqdIAUB46FWmel7jU1taGioqKUl8OZUinX4C3R2FqXgCA8iC8AABJEV4AgKQILwBAUoQXACApwgsAkBThBQBISqciTVIXl7q6ulJfChzzLL0tndHXTL9AOSpMeDFJHfx/AhBQdJqNAICkCC8AQFKEFwAgKcILAJCUwnTYBUrX6Rfg7aTmBQBIivACACRFeAEAkqLPC/CWtcZswQDNpeYFAEhKYcJLfK7RoEGDQmVlZakvBQBoQ4UJL/G5RjU1NaG6urrUlwIAtCF9XoCSa27/mPrtQHkrTM0LAFAehBcAICnCCwCQFOEFAEiKDrtAoSfFA4pHzQsAkBQ1L0ChteQxBS0dsq0GCEpDzQsAkBThBQBIivACACRFeAEAkiK8AABJEV4AgKQUZqj03Llz86Wurq7UlwJwkJYMtwbKpOalqqoq1NTUhOrq6lJfCgDQhgoTXgCA8lCYZiOAotDEBEcmvAAk7q086qA+AHnUASnRbAQAJEXNCwAt0lYPtYTmUvMCACRFeAEAkiK8AABJEV4AgKQILwBAUoQXACApwgsAkBThBQBIivACACRFeAEAkiK8AABJEV4AgKQILwBAUoQXACApwgsAkBThBQBISrsNL3v27AkDBgwIn//850t9KQBAO9Juw8ttt90WLrjgglJfBgDQzrTL8PLKK6+El156KVx22WWlvhQAIPXwsnTp0jB+/PjQp0+f0KFDh7Bw4cKD9pk7d24YOHBg6Nq1axgxYkRYvnx5i84Rm4pmz57d0ksDAMpAi8PL7t27w5AhQ/KAcigPPPBAmDZtWpg1a1ZYtWpVvu/YsWPDli1bGvY5//zzw+DBgw9aNm7cGB555JFw1lln5QsAwIE6hRaKTTlHas6ZM2dOmDJlSpg8eXK+Pm/evLBo0aIwf/78MH369Hzb6tWrD/v9zz77bLj//vvDQw89FHbt2hXefPPNcMopp4SZM2cecv+9e/fmS73a2tqWFgkAKNc+L/v27QsrV64MY8aM+cMJjjsuX1+2bFmzjhGbizZs2BDWr18fvvKVr+RB6HDBpX7/ioqKhqVfv36tUhYAoAzCy7Zt20JdXV3o1atXk+1xfdOmTaEtzJgxI+zYsaNhicEHACiuFjcbvZ0mTZp01H26dOmSLwBAeWjVmpeePXuGjh07hs2bNzfZHtd79+7dmqcCAMpUq4aXzp07h6FDh4YlS5Y0bNu/f3++PnLkyNCW4uinQYMGhcrKyjY9DwCQWLNRHAG0du3ahvV169blo4e6d+8e+vfvnw+TnjhxYhg2bFgYPnx4uPPOO/Ph1fWjj9pKVVVVvsTRRrHjLgBQTC0OLytWrAijR49uWI9hJYqBZcGCBWHChAlh69at+Qih2Ek3zumyePHigzrxAgC8LeFl1KhRIcuyI+4zderUfAEAKKvRRgCUj4HTFx20bf0d4w67nfLVLh/MeCx02AWA8lCY8BI769bU1ITq6upSXwoA0IYKE14AgPIgvAAASRFeAICkFCa86LALAOWhMOFFh10AKA+FCS8AQHkQXgCApAgvAEBShBcAICnCCwCQlMKEF0OlAaA8FCa8GCoNAOWhMOEFACgPwgsAkBThBQBIivACACRFeAEAklKY8GKoNACUh8KEF0OlAaA8FCa8AADlQXgBAJIivAAASRFeAICkCC8AQFKEFwAgKcILAJCUwoQXk9QBQHkoTHgxSR0AlIfChBcAoDwILwBAUoQXACApwgsAkJROpb4AADgWA6cvarK+/o5xJbsW3l5qXgCApAgvAEBShBcAICnCCwCQFOEFAEhKYcKLZxsBQHkoTHjxbCMAKA+FCS8AQHkQXgCApAgvAEBShBcAICnCCwCQFOEFAEiK8AIAJEV4AQCSIrwAAEkRXgCApAgvAEBShBcAICnCCwCQFOEFAEiK8AIAJKUw4WXu3Llh0KBBobKystSXAgC0oU6hIKqqqvKltrY2VFRUlPpyACiBgdMXHbRt/R3jmr29Jfse6Ri0rcKEFwBoL1oadN5KiFpfomOUUmGajQCA8iC8AABJEV4AgKQILwBAUoQXACApwgsAkBThBQBIivACACRFeAEAkiK8AABJEV4AgKQILwBAUoQXACApwgsAkBThBQBISqdQMFmW5X/X1ta2yfH3793TZD2e58BtLd3+dh+j/mdTlGOk9vNvjWO0p59/axwjtZ9/ezmGe5j+MdrTPWzpMVpb/THr38ePpEPWnL0S8r//+7+hX79+pb4MAOAYbNiwIfTt27e8wsv+/fvDxo0bw8knnxw6dOjQJueI6TAGpPgDPuWUU0IRFb2MRS9fOZSx6OUrhzIWvXzlUMbaVixfjCM7d+4Mffr0Cccdd1x5NRvFAh8tsbWWeKOK+MtYTmUsevnKoYxFL185lLHo5SuHMp7SSuWrqKho1n467AIASRFeAICkCC/HoEuXLmHWrFn530VV9DIWvXzlUMail68cylj08pVDGbuUqHyF67ALABSbmhcAICnCCwCQFOEFAEiK8AIAJEV4AQCSIrwcg7lz54aBAweGrl27hhEjRoTly5eHFC1dujSMHz8+n4o5Pkph4cKFTb4eB6LNnDkznHrqqeH4448PY8aMCa+88kpIxezZs0NlZWX+qIh3v/vd4Yorrghr1qxpss8bb7wRqqqqQo8ePcJJJ50UrrrqqrB58+aQiq9+9avhvPPOa5jdcuTIkeHxxx8vTPkOdMcdd+S/q5/5zGcKU8ZbbrklL1Pj5ZxzzilM+aLf/OY34a/+6q/yMsTXkve///1hxYoVhXmtie8HB97DuMT7VoR7WFdXF26++eZw2mmn5ffnjDPOCP/wD//Q5AGKb/s9jEOlab77778/69y5czZ//vzs5z//eTZlypSsW7du2ebNm7PUfO9738u+9KUvZd/97nfjb2D28MMPN/n6HXfckVVUVGQLFy7Mnn/++ewjH/lIdtppp2W/+93vshSMHTs2u/fee7MXX3wxW716dXb55Zdn/fv3z3bt2tWwz6c//emsX79+2ZIlS7IVK1ZkF1xwQXbhhRdmqXj00UezRYsWZS+//HK2Zs2a7Itf/GL2jne8Iy9zEcrX2PLly7OBAwdm5513XnbjjTc2bE+9jLNmzcrOPffc7LXXXmtYtm7dWpjyvf7669mAAQOySZMmZc8991z2q1/9KnviiSeytWvXFua1ZsuWLU3u31NPPZW/pj799NOFuIe33XZb1qNHj+yxxx7L1q1blz300EPZSSedlN11110lu4fCSwsNHz48q6qqalivq6vL+vTpk82ePTtL2YHhZf/+/Vnv3r2zL3/5yw3btm/fnnXp0iX79re/naUovsDEcv7oRz9qKE98o4//Eev94he/yPdZtmxZlqp3vvOd2Te+8Y1ClW/nzp3Ze9/73vxN4cMf/nBDeClCGWN4GTJkyCG/VoTy/f3f/3128cUXH/brRXytib+fZ5xxRl62ItzDcePGZddee22TbVdeeWV29dVXl+weajZqgX379oWVK1fm1WGNHwQZ15ctWxaKZN26dWHTpk1NyhofmBWbyVIt644dO/K/u3fvnv8d7+Wbb77ZpIyxur5///5JljFW7d5///1h9+7defNRkcoXq9zHjRvXpCxRUcoYq9dj8+3pp58err766vDrX/+6MOV79NFHw7Bhw8InPvGJvPn2Ax/4QPj6179e2Nea+D5x3333hWuvvTZvOirCPbzwwgvDkiVLwssvv5yvP//88+EnP/lJuOyyy0p2Dwv3VOm2tG3btvwNolevXk22x/WXXnopFEn8RYwOVdb6r6Vk//79eT+Jiy66KAwePDjfFsvRuXPn0K1bt6TL+MILL+RhJbarx/b0hx9+OAwaNCisXr26EOWLgWzVqlWhurr6oK8V4R7GF/gFCxaEs88+O7z22mvh1ltvDR/60IfCiy++WIjy/epXv8r7Zk2bNi188YtfzO/j3/3d3+XlmjhxYuFea2Lfwe3bt4dJkybl60W4h9OnTw+1tbV56OrYsWP+PnjbbbflQTsqxT0UXigL8ZN7fDOInxaKJr7pxaASa5a+853v5G8IP/rRj0IRbNiwIdx4443hqaeeyjvIF1H9p9codr6OYWbAgAHhwQcfzDs+pi5+cIg1L7fffnu+Hmte4v/FefPm5b+rRXPPPffk9zTWpBXFgw8+GL75zW+Gb33rW+Hcc8/NX2/ih8FYxlLdQ81GLdCzZ888dR7YSzyu9+7dOxRJfXmKUNapU6eGxx57LDz99NOhb9++DdtjOWIVb/yUlHIZ46e6M888MwwdOjQfYTVkyJBw1113FaJ8scp9y5Yt4YMf/GDo1KlTvsRg9q//+q/5v+Mnu9TLeKD4Cf2ss84Ka9euLcQ9jKNPYk1gY+973/samsaK9Frz6quvhu9///vhU5/6VMO2ItzDm266Ka99+eQnP5mPFLvmmmvCZz/72fz1plT3UHhp4ZtEfIOIbX+NP1XE9VhtXyRxSFz8pWtc1lht+NxzzyVT1tgPOQaX2Izygx/8IC9TY/FevuMd72hSxjiUOr6oplLGQ4m/k3v37i1E+S699NK8WSx+0qtf4qf4WF1d/+/Uy3igXbt2hV/+8pf5m34R7mFsqj1wioLYdyLWLhXltabevffem/frif2z6hXhHu7Zsyfv39lY/CAfX2tKdg/bpBtwwYdKxx7UCxYsyGpqarLrrrsuHyq9adOmLDVxBMdPf/rTfIm/CnPmzMn//eqrrzYMfYtle+SRR7Kf/exn2Uc/+tGkhi/+7d/+bT5074c//GGTYYx79uxp2CcOYYzDp3/wgx/kQxhHjhyZL6mYPn16PnoqDl+M9yiud+jQIXvyyScLUb5DaTzaqAhl/NznPpf/jsZ7+Mwzz2RjxozJevbsmY+OK0L54hD3Tp065cNtX3nlleyb3/xmdsIJJ2T33Xdfwz6pv9bUjzyN9ymOrjpQ6vdw4sSJ2Xve856GodJxeo34O/qFL3yhZPdQeDkGd999d/6LGOd7iUOnn3322SxFcQ6CGFoOXOIvav3wt5tvvjnr1atXHtguvfTSfC6RVByqbHGJc7/Ui/+xrr/++nx4cXxB/djHPpYHnFTE4YtxDo34u/iud70rv0f1waUI5WtOeEm9jBMmTMhOPfXU/B7GN4i43ngOlNTLF/33f/93Nnjw4Px15Jxzzsm+9rWvNfl66q81UZy7Jr6+HOq6U7+HtbW1+f+5+L7XtWvX7PTTT8/nCNu7d2/J7mGH+Efb1OkAALQ+fV4AgKQILwBAUoQXACApwgsAkBThBQBIivACACRFeAEAkiK8AABJEV4AgKQILwBAUoQXACCk5P8BJ1iW2NTo1MkAAAAASUVORK5CYII=",
+      "text/plain": [
+       "<Figure size 640x480 with 1 Axes>"
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    }
+   ],
+   "source": [
+    "resp_dx.plot_singular_values()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "5ceff349-0253-40c8-ae08-4a03aa935913",
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "source": [
+    "After solving, correction is available, for instance with\n",
+    "\n",
+    "* {py:meth}`~.ResponseMatrix.correction_matrix` which returns the correction matrix (pseudo-inverse of\n",
+    "  the response matrix),\n",
+    "* {py:meth}`~.ResponseMatrix.get_correction` which returns a correction vector when given observed values,\n",
+    "* {py:meth}`~.ResponseMatrix.correct` which computes and optionally applies a correction\n",
+    "  for the provided {py:class}`.Lattice`.\n",
+    "\n",
+    "### Exclusion of variables and observables\n",
+    "\n",
+    "Variables may be added to a set of excluded values, and similarly for observables.\n",
+    "Excluding an item does not change the response matrix. The values are excluded from the\n",
+    "pseudo-inversion of the response, possibly reducing the number of singular values.\n",
+    "After inversion, the correction matrix is expanded to its original size by inserting\n",
+    "zero lines and columns at the location of excluded items. This way:\n",
+    "\n",
+    "- error and correction vectors keep the same size independently of excluded values,\n",
+    "- excluded error values are ignored,\n",
+    "- excluded corrections are set to zero.\n",
+    "\n",
+    "#### Exclusion of variables\n",
+    "\n",
+    "Excluded variables are selected by their name or their index in the variable list:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 12,
+   "id": "317e2a0f-0e81-479e-89fb-71c49fe03d70",
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "outputs": [],
+   "source": [
+    "resp_dx.exclude_vars(0, \"dx_9\", \"dx_47\", -1)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "d2a5eddb-56ec-4009-a7e8-d748214ee90c",
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "source": [
+    "Where *-1* refers to the last variable.\n",
+    "\n",
+    "#### Exclusion of observables\n",
+    "\n",
+    "Observables are selected by their name or their index in the observable list. In addition, for\n",
+    "{py:class}`.ElementObservable` observables, we need to specify a _refpts_ to identify which item\n",
+    "in the array will be excluded.\n",
+    "\n",
+    "Let's exclude all Monitors with name _BPM\\_07_:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 13,
+   "id": "4223c1bf-0108-4068-b938-4d9313c48684",
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "outputs": [],
+   "source": [
+    "resp_dx.exclude_obs(obsid=\"orbit[x]\", refpts=\"BPM_07\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "d09b2106-ee9b-4171-8c4c-36e855fb4958",
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "source": [
+    "Or by using the observable index:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 14,
+   "id": "6ab3809a-86e3-486f-a688-544c875c9550",
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": [
+     "remove-output"
+    ]
+   },
+   "outputs": [
+    {
+     "name": "stderr",
+     "output_type": "stream",
+     "text": [
+      "/var/folders/9d/tcctx2j125xd3nzkr5bp1zq4000103/T/ipykernel_63123/1399543609.py:1: AtWarning: No new excluded value\n",
+      "  resp_dx.exclude_obs(obsid=0, refpts=\"BPM_07\")\n"
+     ]
+    }
+   ],
+   "source": [
+    "resp_dx.exclude_obs(obsid=0, refpts=\"BPM_07\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "9402e021-90b2-442a-9d7e-a1c02a1ec235",
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "source": [
+    "Or even, since there is a single observable and *obsid* defaults to *0*:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 15,
+   "id": "7acb0cf6-118b-4157-8bc1-6e464ffdaa71",
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": [
+     "remove-output"
+    ]
+   },
+   "outputs": [
+    {
+     "name": "stderr",
+     "output_type": "stream",
+     "text": [
+      "/var/folders/9d/tcctx2j125xd3nzkr5bp1zq4000103/T/ipykernel_63123/3319419291.py:1: AtWarning: No new excluded value\n",
+      "  resp_dx.exclude_obs(refpts=\"BPM_07\")\n"
+     ]
+    }
+   ],
+   "source": [
+    "resp_dx.exclude_obs(refpts=\"BPM_07\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "2e806ea3-3edc-4265-9a6f-558720f0befb",
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "source": [
+    "After excluding items, the pseudo-inverse is discarded so one must recompute it again:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 16,
+   "id": "208f3e79-7ab5-4647-95bd-bb4d8787bea0",
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "outputs": [
+    {
+     "data": {
+      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAi8AAAGzCAYAAADnmPfhAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjEsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvc2/+5QAAAAlwSFlzAAAPYQAAD2EBqD+naQAAJl9JREFUeJzt3QtUVWX+//EvSmBpokiCpEh3I1KLi1E2qTEx5jiZ1TAzTaFOtipsLKsZXLPSWmsKpwvLpmFyuqitqZVmM9KFSS2y7KIJmN0okwaNSVGZErwUJOz/+j7/3znDTQXleM6z9/u11g73ZnP2fs6hcz48z/fZO8xxHEcAAAAs0SPYJwAAANAVhBcAAGAVwgsAALAK4QUAAFiF8AIAAKxCeAEAAFYhvAAAAKsQXgAAgFUILwAAwCqEF8DlEhMTZcqUKRIqxowZYxYbhdpzCXgV4QWw1McffyxXX321DB06VHr16iUnn3yy/PjHP5ZHH3002KcGAAEVHtiHBxAI7733nowdO1YSEhJk+vTpEhcXJ9XV1bJu3Tp55JFH5NZbb/Xvu2nTJunRg79TALgH4QWw0H333SdRUVFSWloq/fr1a/W9nTt3tlqPjIwUtzpw4IA0NzdLREREsE8FwDHEn2OAhb788ks555xz2gUXNXDgwEPWaSxevFjCwsLk3XfflVmzZslJJ50kvXv3liuvvFJ27drV6mc1GNxzzz0SHx8vJ5xwguntqaioaPeYuo8+Zlu+Y23ZsuWgbWlsbJQ5c+ZISkqKCWR6LhdffLGsXr261X76GPpYDz30kMyfP19OO+00E8z0fDqSnJxszrctbZMOsemQm48+5oUXXigDBgyQ448/3pzLCy+8cNBzPtJ2v/rqq6Zt2sYTTzxRJkyYIJ9++mmrfWpqamTq1KkyePBg075BgwbJFVdcccjnEPAael4AC2mdy9q1a+WTTz4xH9JHQoeW+vfvL3PnzjUfjBoIZsyYIUuXLvXvM3v2bHnggQdk4sSJkpWVJR9++KH5+v3333dbW+rr6+XJJ5+UX/7yl2YIbM+ePfLUU0+Z46xfv15GjhzZav9FixaZ4994443mwz06OrrDx83OzjbhQsOADqv5vPPOO7Jt2zb5xS9+4d+mQ20/+9nP5NprrzVhasmSJXLNNdfIK6+8YgJGd/j73/8uOTk5pl1/+tOfZP/+/fLYY4/J6NGj5YMPPjCBUF111VUm0Ojro9u0J+21116Tr776yr8P4HkOAOusWrXK6dmzp1kyMjKc3/3ud87KlSudxsbGdvsOHTrUycnJ8a8vWrTI0f/1MzMznebmZv/222+/3Tze7t27zXpNTY0THh7uTJo0qdXj3XPPPebnWz7m3Llzzba2fMeqqqryb7vkkkvM4nPgwAGnoaGh1c99++23TmxsrDNt2jT/Nn0Mfay+ffs6O3fuPOxztGnTJrP/o48+2mr7Lbfc4vTp08fZv3+/f1vLfyt9HpOTk51x48Yd8rnsbLv37Nnj9OvXz5k+fXqr/fQ5joqK8m/XduvPPfjgg4dtH+BlDBsBFtJZRdrzor0F2huivSP6F70Oh7z00kudegztuWg55KHDGU1NTbJ161azXlJSYmpKbrnlllY/17IYuDv07NnTX7OiQzrffPONOW5qaqps2LCh3f7aM6FDXYdz5plnml6blj1J2j4dDtKeJB0e8mn572+//Vbq6urM89HR8Y+E9pzs3r3b9C7V1tb6F237qFGj/ENkeh76XLz55pvmPAB0jPACWCotLU3++c9/mg85HV7RIR4dctFajoPVgbSkM5Va0iEk5fvQ9IWY008/vdV+Okzj27e7PP300zJ8+HAz5VvrTjScFBcXmxDR1imnnNLpx9WhI63t+frrr826hgIdhtHtLenw0AUXXGCOr+3T4+uQTkfHPxKbN282X8eNG2ceu+WyatUqf5G1DoPpkJLWxsTGxsqPfvQjE0x16AvA/xBeAMvpX+oaZO6//37zgfvDDz/IsmXLDvtz+ld/RxxHRy66pqOiVV9Px+E888wzpvhXC3C11mXFihWmp0I/6LUnpq2WvSSHoyFF2+N7Pp5//nlTFPyTn/zEv8/bb79terA0uPz1r3+Vf/3rX+b4v/rVrw77XHS23b52aN2LPnbb5cUXX/Tve9ttt8kXX3wh+fn55pzuvvtuOfvss01dDID/j4JdwEV0qEVt3769W4qCVWVlZavejv/+97/thjR8PTE6NNJyBpSv9+ZQdBjn1FNPNb1ILcOAFhIfLT3v9PR0M3Skxch6jEmTJrWaPv6Pf/zDhISVK1e22q6FwYfT2XZrMPPNBMvMzDzs4+r+d9xxh1m010aHvx5++GET9ADQ8wJYSWskOuoV0F4DddZZZx31MS699FIJDw83vTkt/eUvf2m3r+/Dec2aNf5t+/btM8NBne0Batme999/39T0dAftfdGL9y1cuNDUmbQdMtLja2hq2Vuis6+KiooO+9idbbfWI/Xt29f0jmnPWFu+Keo6A6ntTC49hk6rbmho6HSbAbej5wWwkBbN6gedXptl2LBhZnqvXnVXexh0Oq1eJ+Roac3FzJkzzV/8OqyiQy1aHKz1GDExMa16SS677DJTQ/Ob3/xG7rrrLhMINCxoTYdO8T2Un/70p6ZHRNui05KrqqpkwYIFkpSUJHv37j3qdvz85z+XO++80yxaz9K250OPWVBQYNqnQ0Vaf1JYWGhqfT766KNDPnZn263BRUPgddddJ+eff76Zpu3bR2t7LrroIhMKdbhIQ6Oes7Zfw+Py5ctlx44draZ2A54X7OlOALru1VdfNdOIhw0bZqb9RkREOKeffrpz6623Ojt27OjUVOnS0tJW+61evdps168tpzHffffdTlxcnHP88cebqcOfffaZM2DAAOemm25q9fPl5eXOqFGjzLkkJCQ4BQUFnZoqrdO177//fnOekZGRznnnnee88sor5px1W9up0kcyjfiiiy4yP3vDDTd0+P2nnnrKOeOMM8zx9TnV8+5oGnTb57Ir7fY9x1lZWWZ6dK9evZzTTjvNmTJlilNWVma+X1tb6+Tm5ppz6N27t9lPH/v555/vcpsBNwvT/wQ7QAGwh9Z3aK3HH//4R/nDH/4Q7NMB4EHUvAA4qO+++67dNr0SrxozZkwQzggAqHkBcAhaQ6P36bn88sulT58+5tL6zz33nKn10DoNAAiGkOx50cI97ZZueeM0AMeeXjhOi0b1Qml6/RG9JooW8er0YgAIlpCsedGrYOqVQnW6YWfu7AoAALwjJHtedCxdr2sAAABw1OFFL8akNzWLj48313no6EJOeo0EvdaEXrVSbzqm910BAAAISsGuXj1yxIgRMm3aNJk8eXKHBX6zZs0yF5nS4KIzE/Tqkps2bTKXxlZ6qWu9a2xbeoMyDUVHQ+8hsm3bNtNzc7D7jgAAgNCiVSxaMqI5oEePw/StHM1FYvTHly9f3mpbenq6uciST1NTkxMfH+/k5+d36bH1Yk5XXXXVYff7/vvvnbq6Ov9SUVFhzouFhYWFhYVFrFuqq6sP+9nfrVOl9RLl5eXlMnv2bP82TU96Oe7uuk9JW3rn1Xvvvbfd9urqanNJbgAAEPrq6+tlyJAhnap57dbwojc905ub6T1RWtL1zz//vNOPo2FH76GiQ1SDBw82t7PPyMjocF8NSjpM1bbxGlwILwAA2KUzJR8heZG6119/vdP76i3sW97GHgAAuFu3TpXWO83qXVX1Dqgt6XpcXFx3HgoAAHhUt4aXiIgISUlJkZKSklazf3T9YMM+3UWnZ+st5NPS0gJ6HAAAEFxdHjbau3evVFZW+terqqpk48aNEh0dLQkJCab+JCcnR1JTUyU9Pd1MldbalalTp0og5ebmmkVrXqKiogJ6LAAAYFF4KSsrk7Fjx/rXfcWyGlj0Bm7Z2dmya9cumTNnjtTU1JhruqxYsaJdES8AAIBr7m10NHw9L3V1dcw2AgDAhZ/fIXlvoyNBzQsAAN5AzwsAAAg6T/a8AAAAbyC8AAAAqxBeAACAVVwTXijYBQDAGyjYBQAAVn1+h+SNGUNZYl5xu21b5k0IyrkAAOBFhJduQqgBAODYcE3NCwAA8AbCCwAAsIprwguzjQAA8AbXhJfc3FypqKiQ0tLSYJ8KAAAIINeEFwAA4A2EFwAAYBXCCwAAsArhBQAAWIXwAgAArOKa8MJUaQAAvME14YWp0gAAeINrwgsAAPAGwgsAALAK4QUAAFiF8AIAAKxCeAEAAFYhvAAAAKu4JrxwnRcAALzBNeGF67wAAOANrgkvAADAGwgvAADAKoQXAABgFcILAACwCuEFAABYhfACAACsQngBAABWIbwAAACrEF4AAIBVXBNeuD0AAADe4Jrwwu0BAADwBteEFwAA4A2EFwAAYBXCCwAAsArhBQAAWIXwAgAArEJ4AQAAViG8AAAAqxBeAACAVQgvAADAKoQXAABgFcILAACwCuEFAABYhfACAACs4prwUlhYKElJSZKWlhbsUwEAAAHkmvCSm5srFRUVUlpaGuxTAQAAAeSa8AIAALyB8AIAAKxCeAEAAFYhvAAAAKsQXgAAgFUILwAAwCqEFwAAYBXCCwAAsArhBQAAWIXwAgAArEJ4AQAAViG8AAAAqxBeAACAVQgvAADAKoQXAABgFcILAACwCuEFAABYhfACAACsQngBAABWCbnwUl1dLWPGjJGkpCQZPny4LFu2LNinBAAAQki4hJjw8HCZP3++jBw5UmpqaiQlJUUuv/xy6d27d7BPDQAAhICQCy+DBg0yi4qLi5OYmBj55ptvCC8AAODIho3WrFkjEydOlPj4eAkLC5OioqJ2+xQWFkpiYqL06tVLRo0aJevXr5cjUV5eLk1NTTJkyJAj+nkAAOA+XQ4v+/btkxEjRpiA0pGlS5fKrFmzZO7cubJhwwazb1ZWluzcudO/jw4JJScnt1u2bdvm30d7W66//np5/PHHj7RtAADAhbo8bDR+/HizHExBQYFMnz5dpk6datYXLFggxcXFsnDhQsnLyzPbNm7ceMhjNDQ0yKRJk8z+F1544WH31cWnvr6+iy0CAACenW3U2NhohnoyMzP/d4AePcz62rVrO/UYjuPIlClTZNy4cXLdddcddv/8/HyJioryLwwxAQDgbt0aXmpra02NSmxsbKvtuq4zhzrj3XffNUNPWkujw0u6fPzxxwfdf/bs2VJXV+dfdKo1AABwr5CbbTR69Ghpbm7u9P6RkZFmAQAA3tCtPS86rblnz56yY8eOVtt1Xac9AwAAhFR4iYiIMBeVKykp8W/TXhRdz8jIkEDS2U96Vd60tLSAHgcAAFg2bLR3716prKz0r1dVVZnZQ9HR0ZKQkGCmSefk5Ehqaqqkp6ebq+Xq9Grf7KNAyc3NNYvONtLCXQAA4E5dDi9lZWUyduxY/7qGFaWBZfHixZKdnS27du2SOXPmmCJdLbhdsWJFuyJeAACAYxJe9KaJOp35UGbMmGEWAAAA199V+khR8wIAgDe4JrxovUtFRYWUlpYG+1QAAEAAuSa8AAAAbyC8AAAAqxBeAACAVVwTXijYBQDAG0Lu3kZuu0hdYl5xu21b5k0IyrkAAOAGrul5AQAA3kB4AQAAViG8AAAAqxBeAACAVVwTXphtBACAN7gmvHB7AAAAvME14QUAAHgD4QUAAFiF8AIAAKxCeAEAAFYhvAAAAKu4JrwwVRoAAG9wTXhhqjQAAN7gmvACAAC8gfACAACsQngBAABWIbwAAACrEF4AAIBVwoN9Al6WmFfcbtuWeRMOuh0AALio54XrvAAA4A2uCS9c5wUAAG9wTXgBAADeQHgBAABWIbwAAACrEF4AAIBVCC8AAMAqhBcAAGAVwgsAALAK4QUAAFiF8AIAAKzimvDC7QEAAPAG14QXbg8AAIA3cFdpi3C3aQAAXNTzAgAAvIHwAgAArEJ4AQAAViG8AAAAqxBeAACAVQgvAADAKoQXAABgFcILAACwCuEFAABYhfACAACsQngBAABWIbwAAACruCa8FBYWSlJSkqSlpQX7VAAAQAC5Jrzk5uZKRUWFlJaWBvtUAABAALkmvAAAAG8gvAAAAKuEB/sEcPQS84rbbdsyb0JQzgUAgEAjvLgcwQYA4DaEF48i1AAAbEXNCwAAsArhBQAAWIXwAgAArEJ4AQAAViG8AAAAqxBeAACAVZgqjVaYQg0ACHX0vAAAAKvQ84JOo1cGABAK6HkBAABWIbwAAACrEF4AAIBVCC8AAMAqhBcAAGCVkAsvu3fvltTUVBk5cqQkJyfLE088EexTAgAAISTkpkqfeOKJsmbNGjnhhBNk3759JsBMnjxZBgwYEOxTAwAAISDkel569uxpgotqaGgQx3HMAgAAcEThRXtFJk6cKPHx8RIWFiZFRUXt9iksLJTExETp1auXjBo1StavX9/loaMRI0bI4MGD5a677pKYmBherRC/eF3bBQCAkBk20qEcDRbTpk0zwzltLV26VGbNmiULFiwwwWX+/PmSlZUlmzZtkoEDB5p9tJ7lwIED7X521apVJhT169dPPvzwQ9mxY4c5xtVXXy2xsbFH2kYECVfkBQCERHgZP368WQ6moKBApk+fLlOnTjXrGmKKi4tl4cKFkpeXZ7Zt3LixU8fSwKJB6e233zYBpiM6tKSLT319fRdbhGONUAMACJmC3cbGRikvL5fZs2f7t/Xo0UMyMzNl7dq1nXoM7W3Rmhct3K2rqzPDVDfffPNB98/Pz5d77723W84foRtsCDwAgIAU7NbW1kpTU1O7IR5dr6mp6dRjbN26VS6++GLT46Jfb731Vjn33HMPur8GJQ05vqW6uvqo2wEAAEJXyE2VTk9P7/SwkoqMjDQLAADwhm7tedFZQTrVWYd+WtL1uLi47jwUAADwqG4NLxEREZKSkiIlJSX+bc3NzWY9IyNDAkmnZyclJUlaWlpAjwMAACwbNtq7d69UVlb616uqqswwT3R0tCQkJJhp0jk5OeYS/zoEpFOldXq1b/ZRoOTm5ppFZxtFRUUF9FgAAMCi8FJWViZjx471r2tYURpYFi9eLNnZ2bJr1y6ZM2eOKdLVa7qsWLGC67QAAIDghJcxY8Yc9nL9M2bMMAsAAIDrZxsdTc2LLjpVG97B9V8AwHtcE16oeUFbBBsAcCfXhBegswg1AGC3bp0qDQAAEGiEFwAAYBXXhBcuUgcAgDe4JrxosW5FRYWUlpYG+1QAAEAAuSa8AAAAbyC8AAAAqxBeAACAVQgvAADAKq4JL8w2AgDAG1wTXphtBACAN7gmvAAAAG8gvAAAAKsQXgAAgFW4qzTwf7jbNADYgfACHAahBgBCS7ibpkrr0tTUFOxTgYcQbADg2At301RpXerr6yUqKirYpwOPO1ioIewAwNGjYBcAAFiF8AIAAKzimmEjwGYdDScphpQAoD16XgAAgFUILwAAwCqEFwAAYBXXhBe9xktSUpKkpaUF+1QAAEAAuSa86DVeKioqpLS0NNinAgAAAojZRkCI44J3ANAa4QVwmUOFGoIQADcgvAA4KEINgFBEeAHQZYQaAMHkmoJdAADgDfS8AOhW9MoACDR6XgAAgFXoeQFwTNAjA6C7EF4AWDe1G4C3uWbYiNsDAADgDeFuuj2ALvX19RIVFRXs0wEQQPTIAN7mmp4XAADgDa7peQEAemQAb6DnBQAAWIXwAgAArEJ4AQAAViG8AAAAqxBeAACAVQgvAADAKoQXAABgFcILAACwChepA+B6R3Lzx65uB3DsEF4AoJsQbIBjg/ACAAFGLw7Qvah5AQAAVnFNz0thYaFZmpqagn0qAHBU6JEBPNLzkpubKxUVFVJaWhrsUwEAAAHkmp4XAPACemUAwgsAuEJ3TPn2fQ8IdYQXAEAr9O4g1Lmm5gUAAHgD4QUAAFiF8AIAAKxCeAEAAFahYBcA0CkU8iJUEF4AAEeFUINjjWEjAABgFcILAACwCuEFAABYhfACAACsQsEuACBgKOZFINDzAgAArEJ4AQAAViG8AAAAqxBeAACAVUI2vOzfv1+GDh0qd955Z7BPBQAAhJCQDS/33XefXHDBBcE+DQAAEGJCMrxs3rxZPv/8cxk/fnywTwUAANgeXtasWSMTJ06U+Ph4CQsLk6Kionb7FBYWSmJiovTq1UtGjRol69ev79IxdKgoPz+/q6cGAAA8oMvhZd++fTJixAgTUDqydOlSmTVrlsydO1c2bNhg9s3KypKdO3f69xk5cqQkJye3W7Zt2yYvvviinHnmmWYBAAA46ivs6lDOoYZzCgoKZPr06TJ16lSzvmDBAikuLpaFCxdKXl6e2bZx48aD/vy6detkyZIlsmzZMtm7d6/88MMP0rdvX5kzZ06H+zc0NJjFp76+vqtNAgAAXq15aWxslPLycsnMzPzfAXr0MOtr167t1GPocFF1dbVs2bJFHnroIROEDhZcfPtHRUX5lyFDhnRLWwAAgAfCS21trTQ1NUlsbGyr7bpeU1MjgTB79mypq6vzLxp8AACAe4X0jRmnTJly2H0iIyPNAgAAvKFbe15iYmKkZ8+esmPHjlbbdT0uLq47DwUAADyqW8NLRESEpKSkSElJiX9bc3OzWc/IyJBA0tlPSUlJkpaWFtDjAAAAy4aNdAZQZWWlf72qqsrMHoqOjpaEhAQzTTonJ0dSU1MlPT1d5s+fb6ZX+2YfBUpubq5ZdLaRFu4CAAB36nJ4KSsrk7Fjx/rXNawoDSyLFy+W7Oxs2bVrl5khpEW6ek2XFStWtCviBQAAOCbhZcyYMeI4ziH3mTFjhlkAAAA8cW+jI0HNCwAA3uCa8KL1LhUVFVJaWhrsUwEAAAHkmvACAAC8gfACAACsEtJX2AUAuFNiXnG7bVvmTQjKucA+rul5oWAXAABvcE14oWAXAABvcE14AQAA3kB4AQAAViG8AAAAqxBeAACAVVwTXphtBACAN7gmvDDbCAAAb3BNeAEAAN5AeAEAAFYhvAAAAKsQXgAAgFUILwAAwCquCS9MlQYAwBtcE16YKg0AgDe4JrwAAABvILwAAACrEF4AAIBVCC8AAMAqhBcAAGAVwgsAALCKa8IL13kBAMAbXBNeuM4LAADe4JrwAgAAvIHwAgAArEJ4AQAAViG8AAAAqxBeAACAVQgvAADAKoQXAABglfBgnwAAAD6JecXttm2ZNyEo54LQRXgBAIQ8Qg1cGV709gC6NDU1BftUAADHEMHGe1xT88LtAQAA8AbX9LwAANCZHhl6auxHeAEA4P8QbOzgmmEjAADgDYQXAABgFcILAACwCuEFAABYhfACAACsQngBAABWIbwAAACrEF4AAIBVCC8AAMAqXGEXAIDD4FYDoYWeFwAAYBXCCwAAsIprwkthYaEkJSVJWlpasE8FAAAEkGtqXnJzc81SX18vUVFRwT4dAICHHaoWpqv1M0fyWG7nmvACAACk2wNSKHLNsBEAAPAGwgsAALAK4QUAAFiF8AIAAKxCeAEAAFYhvAAAAKsQXgAAgFUILwAAwCqEFwAAYBXCCwAAsArhBQAAWIXwAgAArEJ4AQAAViG8AAAAqxBeAACAVQgvAADAKoQXAABglXAJQYmJidK3b1/p0aOH9O/fX1avXh3sUwIAACEiJMOLeu+996RPnz7BPg0AABBiGDYCAADuDi9r1qyRiRMnSnx8vISFhUlRUVG7fQoLC83QT69evWTUqFGyfv36Lh1DH/eSSy6RtLQ0efbZZ7t6igAAwMW6PGy0b98+GTFihEybNk0mT57c7vtLly6VWbNmyYIFC0xwmT9/vmRlZcmmTZtk4MCBZp+RI0fKgQMH2v3sqlWrTCh655135OSTT5bt27dLZmamnHvuuTJ8+PAjbSMAAPByeBk/frxZDqagoECmT58uU6dONesaYoqLi2XhwoWSl5dntm3cuPGQx9DgogYNGiSXX365bNiw4aDhpaGhwSw+9fX1XW0SAADwas1LY2OjlJeXm94S/wF69DDra9eu7XTPzp49e8y/9+7dK2+88Yacc845B90/Pz9foqKi/MuQIUO6oSUAAMAT4aW2tlaampokNja21XZdr6mp6dRj7NixQ0aPHm2Gpi644AK5/vrrTe3LwcyePVvq6ur8S3V19VG3AwAAhK6Qmyp96qmnyocfftjp/SMjI80CAAC8oVvDS0xMjPTs2dP0nrSk63Fxcd15KAAAcAwk5hW327Zl3gRxzbBRRESEpKSkSElJiX9bc3OzWc/IyJBA0unZSUlJhxxiAgAAHux50SLayspK/3pVVZWZPRQdHS0JCQlmmnROTo6kpqZKenq6mSqtRbi+2UeBkpubaxadbaSFuwAAwJ26HF7Kyspk7Nix/nUNK0oDy+LFiyU7O1t27dolc+bMMUW6ek2XFStWtCviBQAAOCbhZcyYMeI4ziH3mTFjhlkAAAC6m2vubUTNCwAA3uCa8KL1LhUVFVJaWhrsUwEAAAHkmvACAAC8gfACAACsQngBAABWcU14oWAXAABvcE14oWAXAABvcE14AQAA3kB4AQAAViG8AAAAqxBeAACAu+9tFMqzjXQ5cOCAWde7SwdCc8P+dtv0WF3d3p2PdSyOwfke+2PYdr5efw1tO19+5zjfoz1Gd/M95uHun6jCnM7sZZH//Oc/MmTIkGCfBgAAOALV1dUyePBgb4WX5uZm2bZtm5x44okSFhYWkGNoOtSApE9w3759xUtoO233Utu92m5F22l732Pcdo0je/bskfj4eOnRo4c3ho18tMGHS2zdRV9Yr/1i+9B22u4lXm23ou20/ViKiorq1H4U7AIAAKsQXgAAgFUIL0cgMjJS5s6da756DW2n7V7i1XYr2k7bQ5nrCnYBAIC70fMCAACsQngBAABWIbwAAACrEF4AAIBVCC8AAMAqhJcjoDeATExMlF69esmoUaNk/fr14jZr1qyRiRMnmss0620WioqKWn1fJ6nNmTNHBg0aJMcff7xkZmbK5s2bxXb5+fmSlpZmbi8xcOBAmTRpkmzatKnVPt9//73k5ubKgAEDpE+fPnLVVVfJjh07xHaPPfaYDB8+3H9lzYyMDHn11Vdd3+625s2bZ37nb7vtNk+0/Z577jHtbbkMGzbME23/+uuv5de//rVpm76PnXvuuVJWVub697nExMR2r7ku+jrb8poTXrpo6dKlMmvWLDMPfsOGDTJixAjJysqSnTt3ipvs27fPtE2DWkceeOAB+fOf/ywLFiyQ999/X3r37m2eB/2lt9lbb71l/qddt26dvPbaa/LDDz/IZZddZp4Pn9tvv11efvllWbZsmdlf76U1efJksZ3eVkM/uMvLy80b+Lhx4+SKK66QTz/91NXtbqm0tFT+9re/mRDXktvbfs4558j27dv9yzvvvOP6tn/77bdy0UUXyXHHHWdCekVFhTz88MPSv39/17/PlZaWtnq99b1OXXPNNfa85nqdF3Reenq6k5ub619vampy4uPjnfz8fMet9Ndk+fLl/vXm5mYnLi7OefDBB/3bdu/e7URGRjrPPfec4yY7d+407X/rrbf87TzuuOOcZcuW+ff57LPPzD5r16513KZ///7Ok08+6Yl279mzxznjjDOc1157zbnkkkucmTNnmu1ub/vcuXOdESNGdPg9N7f997//vTN69OiDft9L73MzZ850TjvtNNNmW15zel66oLGx0fxVql2HLW8Eqetr164Vr6iqqpKamppWz4PeTEuH0Nz2PNTV1Zmv0dHR5qu+/tob07Lt2sWekJDgqrY3NTXJkiVLTI+TDh95od3a4zZhwoRWbVReaLsOhegQ8amnnirXXnutfPXVV65v+0svvSSpqammt0GHiM877zx54oknPPc+19jYKM8884xMmzbNDB3Z8poTXrqgtrbWvKnHxsa22q7r+kvuFb62uv15aG5uNnUP2rWcnJxstmn7IiIipF+/fq5s+8cff2zGuPXS4DfddJMsX75ckpKSXN9uDWo6DKw1T225ve36Ybx48WJZsWKFqXvSD+2LL75Y9uzZ4+q2//vf/zbtPeOMM2TlypVy8803y29/+1t5+umnPfU+V1RUJLt375YpU6aYdVte8/BgnwAQyn+Jf/LJJ63G/93urLPOko0bN5oepxdeeEFycnLMmLebVVdXy8yZM824vxbhe8348eP9/9ZaHw0zQ4cOleeff94UqbqV/nGiPS/333+/WdeeF/3/Xetb9PfeK5566inzO6A9bzah56ULYmJipGfPnu2qrnU9Li5OvMLXVjc/DzNmzJBXXnlFVq9ebQpZfbR92s2qf6m4se36F9fpp58uKSkpphdCi7YfeeQRV7dbu8m14P7888+X8PBws2hg00JN/bf+xenWtndE/+I+88wzpbKy0tWvu84g0l7Fls4++2z/kJkX3ue2bt0qr7/+utxwww3+bba85oSXLr6x65t6SUlJq/Su61oX4BWnnHKK+SVu+TzU19ebanzbnwetT9bgosMlb7zxhmlrS/r66+yElm3XqdT6hmd72zuiv98NDQ2ubvell15qhsu0x8m36F/kWvvh+7db296RvXv3ypdffmk+3N38uutwcNvLIHzxxRem18nt73M+ixYtMvU+WuvlY81rHuyKYdssWbLEVJsvXrzYqaiocG688UanX79+Tk1NjeMmOvPigw8+MIv+mhQUFJh/b9261Xx/3rx5pt0vvvii89FHHzlXXHGFc8oppzjfffedY7Obb77ZiYqKct58801n+/bt/mX//v3+fW666SYnISHBeeONN5yysjInIyPDLLbLy8szs6qqqqrMa6rrYWFhzqpVq1zd7o60nG3k9rbfcccd5vddX/d3333XyczMdGJiYsxMOze3ff369U54eLhz3333OZs3b3aeffZZ54QTTnCeeeYZ/z5ufZ/zzZTV11VnXbVlw2tOeDkCjz76qHlhIyIizNTpdevWOW6zevVqE1raLjk5Oeb7OqXu7rvvdmJjY02Yu/TSS51NmzY5tuuozbosWrTIv4++cd1yyy1mGrG+2V155ZUm4Nhu2rRpztChQ83v9UknnWReU19wcXO7OxNe3Nz27OxsZ9CgQeZ1P/nkk816ZWWlJ9r+8ssvO8nJyeY9bNiwYc7jjz/e6vtufZ9TK1euNO9tHbXHhtc8TP8T7N4fAACAzqLmBQAAWIXwAgAArEJ4AQAAViG8AAAAqxBeAACAVQgvAADAKoQXAABgFcILAACwCuEFAABYhfACAACsQngBAABik/8Hvdh6uI73u3MAAAAASUVORK5CYII=",
+      "text/plain": [
+       "<Figure size 640x480 with 1 Axes>"
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    }
+   ],
+   "source": [
+    "resp_dx.solve()\n",
+    "resp_dx.plot_singular_values()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "47f6baaa-1d3d-45c4-829d-0a129d849c5b",
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "source": [
+    "There are now only 72 singular values instead of 80 (number of active monitors).\n",
+    "\n",
+    "The excluded items can be retrieved with the {py:attr}`~.ResponseMatrix.excluded_obs` and\n",
+    "{py:attr}`~.ResponseMatrix.excluded_vars` properties:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 17,
+   "id": "1cb0b249-8cca-4409-9277-217ba1442a78",
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "{'orbit[x]': array([ 79, 200, 321, 442, 563, 684, 805, 926], dtype=uint32)}\n",
+      "['dx_5', 'dx_9', 'dx_47', 'dx_964']\n"
+     ]
+    }
+   ],
+   "source": [
+    "print(resp_dx.excluded_obs)\n",
+    "print(resp_dx.excluded_vars)"
+   ]
+  },
+  {
+   "attachments": {},
+   "cell_type": "markdown",
+   "id": "29c9fde1-cc6d-4f62-8dba-cff7645a473f",
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "source": [
+    "The exclusion masks can be reset using {py:meth}`~.ResponseMatrix.reset_vars` and\n",
+    "{py:meth}`~.ResponseMatrix.reset_obs`.\n",
+    "\n",
+    "## Orbit response matrix\n",
+    "\n",
+    "An {py:class}`.OrbitResponseMatrix` defines its observables as instances of\n",
+    "{py:class}`.OrbitObservable` and its variables as _KickAngle_ attributes of elements.\n",
+    "\n",
+    "### Instantiation\n",
+    "\n",
+    "By default, the observables are all the {py:class}`.Monitor` elements, and the\n",
+    "variables are all the elements having a *KickAngle* attribute. This is equivalent to:\n",
+    "```python\n",
+    "resp_v = at.OrbitResponseMatrix(ring, \"v\", bpmrefs = at.Monitor,\n",
+    "                                steerrefs = at.checkattr(\"KickAngle\"))\n",
+    "```\n",
+    "The variable elements must have the *KickAngle* attribute used for correction.\n",
+    "It's available for all magnets, though not present by default except in\n",
+    "{py:class}`.Corrector` magnets. For other magnets, the attribute should be\n",
+    "explicitly created.\n",
+    "\n",
+    "There are options in {py:class}`.OrbitResponseMatrix` to include the RF frequency in the\n",
+    "variable list, and the sum of correction angles in the list of observables:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 18,
+   "id": "8e703ee8-0aba-452c-808d-0621ef4e5f1e",
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "(81, 49)\n"
+     ]
+    }
+   ],
+   "source": [
+    "resp_h = at.OrbitResponseMatrix(ring, \"h\", cavrefs=at.RFCavity, steersum=True)\n",
+    "print(resp_h.shape)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "ca0dca48-4b0d-4aa4-96c7-b5ed24280ec2",
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "source": [
+    "### Matrix building\n",
+    "\n",
+    "{py:class}`.OrbitResponseMatrix` has a {py:meth}`~.OrbitResponseMatrix.build_analytical` build method,\n",
+    "using formulas from [^Franchi].\n",
+    "\n",
+    "[^Franchi]: A. Franchi, S.M. Liuzzo, Z. Marti, _\"Analytic formulas for the rapid evaluation of the orbit\n",
+    "response matrix and chromatic functions from lattice parameters in circular accelerators\"_,\n",
+    "arXiv:1711.06589 [physics.acc-ph]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 19,
+   "id": "65a21d36-6d0d-46f0-ba72-5092fec3ab17",
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "array([[-8.54870374e+00, -1.81104969e+01, -1.21497454e+01, ...,\n",
+       "        -2.85830326e+01, -1.60772960e+01, -5.75838151e-08],\n",
+       "       [ 1.85290756e+01,  3.16574224e+01,  1.82799745e+01, ...,\n",
+       "         1.90656557e+01,  8.85865108e+00, -2.68128309e-06],\n",
+       "       [ 1.68213604e+01,  2.87634523e+01,  1.65301788e+01, ...,\n",
+       "         1.94726063e+01,  9.44097149e+00, -2.44436213e-06],\n",
+       "       ...,\n",
+       "       [ 8.84897619e+00,  1.90656922e+01,  1.29411619e+01, ...,\n",
+       "         3.16578009e+01,  1.85390265e+01, -2.68138514e-06],\n",
+       "       [-1.60775574e+01, -2.85833742e+01, -1.65903209e+01, ...,\n",
+       "        -1.81113028e+01, -8.54900471e+00, -5.76356367e-08],\n",
+       "       [ 1.00000000e+00,  1.00000000e+00,  1.00000000e+00, ...,\n",
+       "         1.00000000e+00,  1.00000000e+00,  0.00000000e+00]],\n",
+       "      shape=(81, 49))"
+      ]
+     },
+     "execution_count": 19,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "resp_h.build_analytical()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "538b8409-99f1-4f9c-9f00-861a0bdb3ac4",
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "source": [
+    "### Matrix normalisation\n",
+    "\n",
+    "This is critical when including the RF frequency response which is not commensurate\n",
+    "with steerer responses. Similarly for rows, the sum of steerers is not commensurate with\n",
+    "monitor readings.\n",
+    "\n",
+    "By default, the normalisation is done automatically by adjusting the RF frequency step\n",
+    "and the weight of the steerer sum based on an approximate analytical response matrix.\n",
+    "Explicitly specifying the *cavdelta* and *stsumweight* prevents this automatic normalisation.\n",
+    "\n",
+    "After building the response matrix, and before solving, normalisation may be applied\n",
+    "with the {py:meth}`~.OrbitResponseMatrix.normalise` method. The default normalisation\n",
+    "gives a higher priority to RF response and steerer sum."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "eb624b1d-1624-4196-909e-35c055f3aa34",
+   "metadata": {
+    "editable": true,
+    "raw_mimetype": "",
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "source": [
+    "## Trajectory response matrix\n",
+    "\n",
+    "A {py:class}`.TrajectoryResponseMatrix` defines its observables as instances of\n",
+    "{py:class}`.TrajectoryObservable` and its variables as _KickAngle_ attributes of elements.\n",
+    "\n",
+    "### Instantiation\n",
+    "\n",
+    "By default, the observables are all the {py:class}`.Monitor` elements, and the\n",
+    "variables are all the elements having a *KickAngle* attribute. This is equivalent to:\n",
+    "```python\n",
+    "resp_v = at.TrajectoryResponseMatrix(lattice, \"v\", bpmrefs = at.Monitor,\n",
+    "                                     steerrefs = at.checkattr(\"KickAngle\"))\n",
+    "```\n",
+    "The variable elements must have the *KickAngle* attribute used for correction.\n",
+    "It's available for all magnets, though not present by default except in\n",
+    "{py:class}`.Corrector` magnets. For other magnets, the attribute should be\n",
+    "explicitly created."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 20,
+   "id": "8c7de11a-d8e5-4fd9-9784-25a9244a0305",
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "(80, 48)\n"
+     ]
+    }
+   ],
+   "source": [
+    "resp_h = at.TrajectoryResponseMatrix(ring, \"h\")\n",
+    "print(resp_h.shape)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "896f42fa-6c67-44b4-8e54-3ff39589949a",
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "source": [
+    "## References"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3 (ipykernel)",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.11.11"
+  },
+  "toc": {
+   "base_numbering": 2
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
diff --git a/pyat/at/lattice/axisdef.py b/pyat/at/lattice/axisdef.py
index 6f51608cf..5fa4222be 100644
--- a/pyat/at/lattice/axisdef.py
+++ b/pyat/at/lattice/axisdef.py
@@ -1,23 +1,23 @@
 """Helper functions for axis and plane descriptions"""
 
 from __future__ import annotations
-from typing import Optional, Union
 
-# For sys.version_info.minor < 9:
+# Necessary for type aliases in python <= 3.8 :
 from typing import Tuple
+from typing import Union
 
 AxisCode = Union[str, int, slice, None, type(Ellipsis)]
 AxisDef = Union[AxisCode, Tuple[AxisCode, AxisCode]]
 
-_axis_def = dict(
-    x=dict(index=0, label="x", unit=" [m]"),
-    px=dict(index=1, label=r"$p_x$", unit=" [rad]"),
-    y=dict(index=2, label="y", unit=" [m]"),
-    py=dict(index=3, label=r"$p_y$", unit=" [rad]"),
-    dp=dict(index=4, label=r"$\delta$", unit=""),
-    ct=dict(index=5, label=r"$\beta c \tau$", unit=" [m]"),
-)
-for xk, xv in [it for it in _axis_def.items()]:
+_axis_def = {
+    "x": {"index": 0, "label": "x", "unit": " [m]"},
+    "px": {"index": 1, "label": r"$p_x$", "unit": " [rad]"},
+    "y": {"index": 2, "label": "y", "unit": " [m]"},
+    "py": {"index": 3, "label": r"$p_y$", "unit": " [rad]"},
+    "dp": {"index": 4, "label": r"$\delta$", "unit": ""},
+    "ct": {"index": 5, "label": r"$\beta c \tau$", "unit": " [m]"},
+}
+for xk, xv in list(_axis_def.items()):
     xv["code"] = xk
     _axis_def[xv["index"]] = xv
     _axis_def[xk.upper()] = xv
@@ -26,15 +26,15 @@
 _axis_def["yp"] = _axis_def["py"]  # For backward compatibility
 _axis_def["s"] = _axis_def["ct"]
 _axis_def["S"] = _axis_def["ct"]
-_axis_def[None] = dict(index=None, label="", unit="", code=":")
-_axis_def[Ellipsis] = dict(index=Ellipsis, label="", unit="", code="...")
-
-_plane_def = dict(
-    x=dict(index=0, label="x", unit=" [m]"),
-    y=dict(index=1, label="y", unit=" [m]"),
-    z=dict(index=2, label="z", unit=""),
-)
-for xk, xv in [it for it in _plane_def.items()]:
+_axis_def[None] = {"index": None, "label": "", "unit": "", "code": ":"}
+_axis_def[Ellipsis] = {"index": Ellipsis, "label": "", "unit": "", "code": "..."}
+
+_plane_def = {
+    "x": {"index": 0, "label": "x", "unit": " [m]"},
+    "y": {"index": 1, "label": "y", "unit": " [m]"},
+    "z": {"index": 2, "label": "z", "unit": ""},
+}
+for xk, xv in list(_plane_def.items()):
     xv["code"] = xk
     _plane_def[xv["index"]] = xv
     _plane_def[xk.upper()] = xv
@@ -42,25 +42,27 @@
 _plane_def["v"] = _plane_def["y"]
 _plane_def["H"] = _plane_def["x"]
 _plane_def["V"] = _plane_def["y"]
-_plane_def[None] = dict(index=None, label="", unit="", code=":")
-_plane_def[Ellipsis] = dict(index=Ellipsis, label="", unit="", code="...")
+_plane_def[None] = {"index": None, "label": "", "unit": "", "code": ":"}
+_plane_def[Ellipsis] = {"index": Ellipsis, "label": "", "unit": "", "code": "..."}
 
 
-def _descr(dd: dict, arg: AxisDef, key: Optional[str] = None):
-    if isinstance(arg, tuple):
-        return tuple(_descr(dd, a, key=key) for a in arg)
-    else:
-        try:
-            descr = dd[arg]
-        except (TypeError, KeyError):
-            descr = dict(index=arg, code=arg, label="", unit="")
-        if key is None:
-            return descr
+def _descr(dd: dict, *args: AxisDef, key: str | None = None):
+    for arg in args:
+        if isinstance(arg, tuple):
+            for a in arg:
+                yield from _descr(dd, a, key=key)
         else:
-            return descr[key]
+            if isinstance(arg, slice):
+                descr = {"index": arg, "code": arg, "label": "", "unit": ""}
+            else:
+                descr = dd[arg]
+            if key is None:
+                yield descr
+            else:
+                yield descr[key]
 
 
-def axis_(axis: AxisDef, key: Optional[str] = None):
+def axis_(*axis: AxisDef, key: str | None = None):
     r"""Return axis descriptions
 
     Parameters:
@@ -100,28 +102,32 @@ def axis_(axis: AxisDef, key: Optional[str] = None):
 
     Examples:
 
-        >>> axis_(('x','dp'), key='index')
+        >>> axis_("x", "dp", key="index")
         (0, 4)
 
         returns the indices in the standard coordinate vector
 
-        >>> dplabel = axis_('dp', key='label')
+        >>> dplabel = axis_("dp", key="label")
         >>> print(dplabel)
         $\delta$
 
         returns the coordinate label for plot annotation
 
-        >>> axis_((0,'dp'))
+        >>> axis_(0, "dp")
         ({'plane': 0, 'label': 'x', 'unit': ' [m]', 'code': 'x'},
          {'plane': 4, 'label': '$\\delta$', 'unit': '', 'code': 'dp'})
 
         returns the entire description directories
 
     """
-    return _descr(_axis_def, axis, key=key)
+    ret = tuple(_descr(_axis_def, *axis, key=key))
+    if len(ret) > 1:
+        return ret
+    else:
+        return ret[0]
 
 
-def plane_(plane: AxisDef, key: Optional[str] = None):
+def plane_(*plane: AxisDef, key: str | None = None):
     r"""Return plane descriptions
 
     Parameters:
@@ -154,16 +160,20 @@ def plane_(plane: AxisDef, key: Optional[str] = None):
 
     Examples:
 
-        >>> plane_('v', key='index')
+        >>> plane_("v", key="index")
         1
 
         returns the indices in the standard coordinate vector
 
-        >>> plane_(('x','y'))
-        ({'plane': 0, 'label': 'h', 'unit': ' [m]', 'code': 'h'},
-         {'plane': 1, 'label': 'v', 'unit': ' [m]', 'code': 'v'})
+        >>> plane_("x", "y")
+        ({'plane': 0, 'label': 'x', 'unit': ' [m]', 'code': 'h'},
+         {'plane': 1, 'label': 'y', 'unit': ' [m]', 'code': 'v'})
 
         returns the entire description directories
 
     """
-    return _descr(_plane_def, plane, key=key)
+    ret = tuple(_descr(_plane_def, *plane, key=key))
+    if len(ret) > 1:
+        return ret
+    else:
+        return ret[0]
diff --git a/pyat/at/lattice/utils.py b/pyat/at/lattice/utils.py
index ff87986a8..179cd97b2 100644
--- a/pyat/at/lattice/utils.py
+++ b/pyat/at/lattice/utils.py
@@ -30,47 +30,63 @@
    is :py:obj:`True` for selected elements.
 
 """
-import numpy
+
+from __future__ import annotations
+
 import functools
 import re
-from typing import Callable, Optional, Sequence, Iterator
-from typing import Union, Tuple, List, Type
 from enum import Enum
+from fnmatch import fnmatch
 from itertools import compress
 from operator import attrgetter
-from fnmatch import fnmatch
+from typing import Optional, Union
+# Necessary for type aliases in python <= 3.8 :
+# from collections.abc import Callable, Sequence, Iterator
+from typing import Callable, Sequence, Iterator, Type
+
+import numpy
+import numpy.typing as npt
+
 from .elements import Element, Dipole
 
 _GEOMETRY_EPSIL = 1.0e-3
 
-ElementFilter = Callable[[Element], bool]
-BoolRefpts = numpy.ndarray
-Uint32Refpts = numpy.ndarray
-
-
-__all__ = ['All', 'End', 'AtError', 'AtWarning', 'axis_descr',
-           'check_radiation', 'check_6d',
-           'set_radiation', 'set_6d',
-           'make_copy', 'uint32_refpts', 'bool_refpts',
-           'get_uint32_index', 'get_bool_index',
-           'checkattr', 'checktype', 'checkname',
-           'get_elements', 'get_s_pos',
-           'refpts_count', 'refpts_iterator',
-           'set_shift', 'set_tilt', 'set_rotation',
-           'tilt_elem', 'shift_elem', 'rotate_elem',
-           'get_value_refpts', 'set_value_refpts', 'Refpts',
-           'get_geometry', 'setval', 'getval']
-
-_axis_def = dict(
-    x=dict(index=0, label="x", unit=" [m]"),
-    xp=dict(index=1, label="x'", unit=" [rad]"),
-    y=dict(index=2, label="y", unit=" [m]"),
-    yp=dict(index=3, label="y'", unit=" [rad]"),
-    dp=dict(index=4, label=r"$\delta$", unit=""),
-    ct=dict(index=5, label=r"$\beta c \tau$", unit=" [m]"),
-)
-for vvv in [vv for vv in _axis_def.values()]:
-    _axis_def[vvv['index']] = vvv
+__all__ = [
+    "All",
+    "End",
+    "AtError",
+    "AtWarning",
+    "BoolRefpts",
+    "Uint32Refpts",
+    "check_radiation",
+    "check_6d",
+    "set_radiation",
+    "set_6d",
+    "make_copy",
+    "uint32_refpts",
+    "bool_refpts",
+    "get_uint32_index",
+    "get_bool_index",
+    "checkattr",
+    "checktype",
+    "checkname",
+    "get_elements",
+    "get_s_pos",
+    "refpts_count",
+    "refpts_iterator",
+    "set_shift",
+    "set_tilt",
+    "set_rotation",
+    "tilt_elem",
+    "shift_elem",
+    "rotate_elem",
+    "get_value_refpts",
+    "set_value_refpts",
+    "Refpts",
+    "get_geometry",
+    "setval",
+    "getval",
+]
 
 
 class AtError(Exception):
@@ -83,14 +99,17 @@ class AtWarning(UserWarning):
 
 _typ1 = "None, All, End, int, bool"
 
-_typ2 = "None, All, End, int, bool, str, Type[Element], ElementFilter"
+_typ2 = "None, All, End, int, bool, str, type[Element], ElementFilter"
 
 
 class RefptsCode(Enum):
-    All = 'All'
-    End = 'End'
+    All = "All"
+    End = "End"
 
 
+ElementFilter = Callable[[Element], bool]
+BoolRefpts = npt.NDArray[bool]
+Uint32Refpts = npt.NDArray[numpy.uint32]
 RefIndex = Union[None, int, Sequence[int], bool, Sequence[bool], RefptsCode]
 Refpts = Union[Type[Element], Element, ElementFilter, str, RefIndex]
 
@@ -105,19 +124,44 @@ class RefptsCode(Enum):
 End = RefptsCode.End
 
 
+def _chkattr(attrname: str, el):
+    return hasattr(el, attrname)
+
+
+def _chkattrval(attrname: str, attrvalue, el):
+    try:
+        v = getattr(el, attrname)
+    except AttributeError:
+        return False
+    else:
+        return v == attrvalue
+
+
+def _chkpattern(pattern: str, el):
+    return fnmatch(el.FamName, pattern)
+
+
+def _chkregex(pattern: str, el):
+    rgx = re.compile(pattern)
+    return rgx.fullmatch(el.FamName)
+
+
+def _chktype(eltype: type, el):
+    return isinstance(el, eltype)
+
+
 def _type_error(refpts, types):
     if isinstance(refpts, numpy.ndarray):
         tp = refpts.dtype.type
     else:
         tp = type(refpts)
-    return TypeError(
-        "Invalid refpts type {0}. Allowed types: {1}".format(tp, types))
+    return TypeError(f"Invalid refpts type {tp}. Allowed types: {types}")
 
 
 # setval and getval return pickleable functions: no inner, nested function
 # are allowed. So nested functions are replaced be module-level callable
 # class instances
-class _AttrItemGetter(object):
+class _AttrItemGetter:
     __slots__ = ["attrname", "index"]
 
     def __init__(self, attrname: str, index: int):
@@ -133,7 +177,7 @@ def getval(attrname: str, index: Optional[int] = None) -> Callable:
     attribute *attrname* of its operand. Examples:
 
     - After ``f = getval('Length')``, ``f(elem)`` returns ``elem.Length``
-    - After ``f = getval('PolynomB, index=1)``, ``f(elem)`` returns
+    - After ``f = getval('PolynomB', index=1)``, ``f(elem)`` returns
       ``elem.PolynomB[1]``
 
     """
@@ -143,7 +187,7 @@ def getval(attrname: str, index: Optional[int] = None) -> Callable:
         return _AttrItemGetter(attrname, index)
 
 
-class _AttrSetter(object):
+class _AttrSetter:
     __slots__ = ["attrname"]
 
     def __init__(self, attrname: str):
@@ -153,7 +197,7 @@ def __call__(self, elem, value):
         setattr(elem, self.attrname, value)
 
 
-class _AttrItemSetter(object):
+class _AttrItemSetter:
     __slots__ = ["attrname", "index"]
 
     def __init__(self, attrname: str, index: int):
@@ -170,7 +214,7 @@ def setval(attrname: str, index: Optional[int] = None) -> Callable:
 
     - After ``f = setval('Length')``, ``f(elem, value)`` is equivalent to
       ``elem.Length = value``
-    - After ``f = setval('PolynomB, index=1)``, ``f(elem, value)`` is
+    - After ``f = setval('PolynomB', index=1)``, ``f(elem, value)`` is
       equivalent to ``elem.PolynomB[1] = value``
 
     """
@@ -180,56 +224,6 @@ def setval(attrname: str, index: Optional[int] = None) -> Callable:
         return _AttrItemSetter(attrname, index)
 
 
-# noinspection PyIncorrectDocstring
-def axis_descr(*args, key=None) -> Tuple:
-    r"""axis_descr(axis [ ,axis], key=None)
-
-    Return a tuple containing for each input argument the requested information
-
-    Parameters:
-        axis (Union[int, str]):    either an index in 0:6 or a string in
-          ['x', 'xp', 'y', 'yp', 'dp', 'ct']
-        key:                        key in the coordinate description
-          dictionary, selecting the desired information. One of :
-
-          'index'
-            index in the standard AT coordinate vector
-          'label'
-            label for plot annotation
-          'unit'
-            coordinate unit
-          :py:obj:`None`
-            entire description dictionary
-
-    Returns:
-        descr (Tuple): requested information for each input argument.
-
-    Examples:
-
-        >>> axis_descr('x','dp', key='index')
-        (0, 4)
-
-        returns the indices in the standard coordinate vector
-
-        >>> dplabel, = axis_descr('dp', key='label')
-        >>> print(dplabel)
-        $\delta$
-
-        returns the coordinate label for plot annotation
-
-        >>> axis_descr('x','dp')
-        ({'index': 0, 'label': 'x', 'unit': ' [m]'},
-         {'index': 4, 'label': '$\\delta$', 'unit': ''})
-
-        returns the entire description directories
-
-    """
-    if key is None:
-        return tuple(_axis_def[k] for k in args)
-    else:
-        return tuple(_axis_def[k][key] for k in args)
-
-
 def check_radiation(rad: bool) -> Callable:
     r"""Deprecated decorator for optics functions (see :py:func:`check_6d`).
 
@@ -271,15 +265,17 @@ def check_6d(is_6d: bool) -> Callable:
     See Also:
         :py:func:`set_6d`
     """
+
     def radiation_decorator(func):
         @functools.wraps(func)
         def wrapper(ring, *args, **kwargs):
-            ringrad = getattr(ring, 'is_6d', is_6d)
+            ringrad = getattr(ring, "is_6d", is_6d)
             if ringrad != is_6d:
-                raise AtError('{0} needs "ring.is_6d" {1}'.format(
-                    func.__name__, is_6d))
+                raise AtError(f'{func.__name__} needs "ring.is_6d" {is_6d}')
             return func(ring, *args, **kwargs)
+
         return wrapper
+
     return radiation_decorator
 
 
@@ -311,21 +307,26 @@ def set_6d(is_6d: bool) -> Callable:
         See Also:
             :py:func:`check_6d`, :py:meth:`.Lattice.enable_6d`,
             :py:meth:`.Lattice.disable_6d`
-        """
+    """
     if is_6d:
+
         def setrad_decorator(func):
             @functools.wraps(func)
             def wrapper(ring, *args, **kwargs):
                 rg = ring if ring.is_6d else ring.enable_6d(copy=True)
                 return func(rg, *args, **kwargs)
+
             return wrapper
     else:
+
         def setrad_decorator(func):
             @functools.wraps(func)
             def wrapper(ring, *args, **kwargs):
                 rg = ring.disable_6d(copy=True) if ring.is_6d else ring
                 return func(rg, *args, **kwargs)
+
             return wrapper
+
     return setrad_decorator
 
 
@@ -346,6 +347,7 @@ def make_copy(copy: bool) -> Callable:
           :pycode:`ring`
     """
     if copy:
+
         def copy_decorator(func):
             @functools.wraps(func)
             def wrapper(ring, refpts, *args, **kwargs):
@@ -353,20 +355,22 @@ def wrapper(ring, refpts, *args, **kwargs):
                     ring = ring.replace(refpts)
                 except AttributeError:
                     check = get_bool_index(ring, refpts)
-                    ring = [el.deepcopy() if ok else el 
-                            for el, ok in zip(ring, check)]
+                    ring = [el.deepcopy() if ok else el for el, ok in zip(ring, check)]
                 func(ring, refpts, *args, **kwargs)
                 return ring
+
             return wrapper
     else:
+
         def copy_decorator(func):
             return func
+
     return copy_decorator
 
 
-def uint32_refpts(refpts: RefIndex, n_elements: int,
-                  endpoint: bool = True,
-                  types: str = _typ1) -> Uint32Refpts:
+def uint32_refpts(
+    refpts: RefIndex, n_elements: int, endpoint: bool = True, types: str = _typ1
+) -> Uint32Refpts:
     r"""Return a :py:obj:`~numpy.uint32` array of element indices selecting
     ring elements.
 
@@ -390,7 +394,7 @@ def uint32_refpts(refpts: RefIndex, n_elements: int,
     """
     refs = numpy.ravel(refpts)
     if refpts is RefptsCode.All:
-        stop = n_elements+1 if endpoint else n_elements
+        stop = n_elements + 1 if endpoint else n_elements
         return numpy.arange(stop, dtype=numpy.uint32)
     elif refpts is RefptsCode.End:
         if not endpoint:
@@ -401,23 +405,22 @@ def uint32_refpts(refpts: RefIndex, n_elements: int,
     elif numpy.issubdtype(refs.dtype, numpy.bool_):
         return numpy.flatnonzero(refs).astype(numpy.uint32)
     elif numpy.issubdtype(refs.dtype, numpy.integer):
-
         # Handle negative indices
         if endpoint:
-            refs = numpy.array([i if (i == n_elements) else i % n_elements
-                                for i in refs], dtype=numpy.uint32)
+            refs = numpy.array(
+                [i if (i == n_elements) else i % n_elements for i in refs],
+                dtype=numpy.uint32,
+            )
         else:
-            refs = numpy.array([i % n_elements
-                                for i in refs], dtype=numpy.uint32)
+            refs = numpy.array([i % n_elements for i in refs], dtype=numpy.uint32)
         # Check ascending
         if refs.size > 1:
             prev = refs[0]
             for nxt in refs[1:]:
                 if nxt < prev:
-                    raise IndexError('Index out of range or not in ascending'
-                                     ' order')
+                    raise IndexError("Index out of range or not in ascending order")
                 elif nxt == prev:
-                    raise IndexError('Duplicated index')
+                    raise IndexError("Duplicated index")
                 prev = nxt
 
         return refs
@@ -426,9 +429,9 @@ def uint32_refpts(refpts: RefIndex, n_elements: int,
 
 
 # noinspection PyIncorrectDocstring
-def get_uint32_index(ring: Sequence[Element], refpts: Refpts,
-                     endpoint: bool = True,
-                     regex: bool = False) -> Uint32Refpts:
+def get_uint32_index(
+    ring: Sequence[Element], refpts: Refpts, endpoint: bool = True, regex: bool = False
+) -> Uint32Refpts:
     # noinspection PyUnresolvedReferences, PyShadowingNames
     r"""Returns an integer array of element indices, selecting ring elements.
 
@@ -456,7 +459,7 @@ def get_uint32_index(ring: Sequence[Element], refpts: Refpts,
 
         numpy array([:pycode:`len(ring)+1`])
 
-        >>> get_uint32_index(ring, at.checkattr('Frequency'))
+        >>> get_uint32_index(ring, at.checkattr("Frequency"))
         array([0], dtype=uint32)
 
         numpy array of indices of all elements having a 'Frequency'
@@ -473,13 +476,14 @@ def get_uint32_index(ring: Sequence[Element], refpts: Refpts,
     else:
         return uint32_refpts(refpts, len(ring), endpoint=endpoint, types=_typ2)
 
-    return numpy.fromiter((i for i, el in enumerate(ring) if checkfun(el)),
-                          dtype=numpy.uint32)
+    return numpy.fromiter(
+        (i for i, el in enumerate(ring) if checkfun(el)), dtype=numpy.uint32
+    )
 
 
-def bool_refpts(refpts: RefIndex, n_elements: int,
-                endpoint: bool = True,
-                types: str = _typ1) -> BoolRefpts:
+def bool_refpts(
+    refpts: RefIndex, n_elements: int, endpoint: bool = True, types: str = _typ1
+) -> BoolRefpts:
     r"""Returns a :py:class:`bool` array of element indices, selecting ring
     elements.
 
@@ -502,7 +506,7 @@ def bool_refpts(refpts: RefIndex, n_elements: int,
           :py:class:`.Element`\ s in a lattice.
     """
     refs = numpy.ravel(refpts)
-    stop = n_elements+1 if endpoint else n_elements
+    stop = n_elements + 1 if endpoint else n_elements
     if refpts is RefptsCode.All:
         return numpy.ones(stop, dtype=bool)
     elif refpts is RefptsCode.End:
@@ -528,8 +532,9 @@ def bool_refpts(refpts: RefIndex, n_elements: int,
 
 
 # noinspection PyIncorrectDocstring
-def get_bool_index(ring: Sequence[Element], refpts: Refpts,
-                   endpoint: bool = True, regex: bool = False) -> BoolRefpts:
+def get_bool_index(
+    ring: Sequence[Element], refpts: Refpts, endpoint: bool = True, regex: bool = False
+) -> BoolRefpts:
     # noinspection PyUnresolvedReferences, PyShadowingNames
     r"""Returns a bool array of element indices, selecting ring elements.
 
@@ -557,7 +562,7 @@ def get_bool_index(ring: Sequence[Element], refpts: Refpts,
         Returns a numpy array of booleans where all elements whose *FamName*
         matches "Q[FD]*" are :py:obj:`True`
 
-        >>> refpts = get_bool_index(ring, at.checkattr('K', 0.0))
+        >>> refpts = get_bool_index(ring, at.checkattr("K", 0.0))
 
         Returns a numpy array of booleans where all elements whose *K*
         attribute is 0.0 are :py:obj:`True`
@@ -583,8 +588,7 @@ def get_bool_index(ring: Sequence[Element], refpts: Refpts,
     return boolrefs
 
 
-def checkattr(attrname: str, attrvalue: Optional = None) \
-        -> ElementFilter:
+def checkattr(attrname: str, attrvalue: Optional = None) -> ElementFilter:
     # noinspection PyUnresolvedReferences
     r"""Checks the presence or the value of an attribute
 
@@ -605,27 +609,23 @@ def checkattr(attrname: str, attrvalue: Optional = None) \
 
     Examples:
 
-        >>> cavs = filter(checkattr('Frequency'), ring)
+        >>> cavs = filter(checkattr("Frequency"), ring)
 
         Returns an iterator over all elements in *ring* that have a
         :pycode:`Frequency` attribute
 
-        >>> elts = filter(checkattr('K', 0.0), ring)
+        >>> elts = filter(checkattr("K", 0.0), ring)
 
         Returns an iterator over all elements in ring that have a
         :pycode:`K` attribute equal to 0.0
     """
-    def testf(el):
-        try:
-            v = getattr(el, attrname)
-            return (attrvalue is None) or (v == attrvalue)
-        except AttributeError:
-            return False
-
-    return testf
+    if attrvalue is None:
+        return functools.partial(_chkattr, attrname)
+    else:
+        return functools.partial(_chkattrval, attrname, attrvalue)
 
 
-def checktype(eltype: Union[type, Tuple[type, ...]]) -> ElementFilter:
+def checktype(eltype: Union[type, tuple[type, ...]]) -> ElementFilter:
     # noinspection PyUnresolvedReferences
     r"""Checks the type of an element
 
@@ -646,7 +646,7 @@ def checktype(eltype: Union[type, Tuple[type, ...]]) -> ElementFilter:
 
         Returns an iterator over all quadrupoles in ring
     """
-    return lambda el: isinstance(el, eltype)
+    return functools.partial(_chktype, eltype)
 
 
 def checkname(pattern: str, regex: bool = False) -> ElementFilter:
@@ -669,19 +669,19 @@ def checkname(pattern: str, regex: bool = False) -> ElementFilter:
 
     Examples:
 
-        >>> qps = filter(checkname('QF*'), ring)
+        >>> qps = filter(checkname("QF*"), ring)
 
         Returns an iterator over all with name starting with ``QF``.
     """
     if regex:
-        rgx = re.compile(pattern)
-        return lambda el: rgx.fullmatch(el.FamName)
+        return functools.partial(_chkregex, pattern)
     else:
-        return lambda el: fnmatch(el.FamName, pattern)
+        return functools.partial(_chkpattern, pattern)
 
 
-def refpts_iterator(ring: Sequence[Element], refpts: Refpts,
-                    regex: bool = False) -> Iterator[Element]:
+def refpts_iterator(
+    ring: Sequence[Element], refpts: Refpts, regex: bool = False
+) -> Iterator[Element]:
     r"""Return an iterator over selected elements in a lattice
 
     Parameters:
@@ -722,9 +722,9 @@ def refpts_iterator(ring: Sequence[Element], refpts: Refpts,
 
 
 # noinspection PyUnusedLocal,PyIncorrectDocstring
-def refpts_count(refpts: RefIndex, n_elements: int,
-                 endpoint: bool = True,
-                 types: str = _typ1) -> int:
+def refpts_count(
+    refpts: RefIndex, n_elements: int, endpoint: bool = True, types: str = _typ1
+) -> int:
     r"""Returns the number of reference points
 
     Parameters:
@@ -745,7 +745,7 @@ def refpts_count(refpts: RefIndex, n_elements: int,
     """
     refs = numpy.ravel(refpts)
     if refpts is RefptsCode.All:
-        return n_elements+1 if endpoint else n_elements
+        return n_elements + 1 if endpoint else n_elements
     elif refpts is RefptsCode.End:
         if not endpoint:
             raise IndexError('"End" index out of range')
@@ -760,8 +760,9 @@ def refpts_count(refpts: RefIndex, n_elements: int,
         raise _type_error(refpts, types)
 
 
-def _refcount(ring: Sequence[Element], refpts: Refpts,
-              endpoint: bool = True, regex: bool = False) -> int:
+def _refcount(
+    ring: Sequence[Element], refpts: Refpts, endpoint: bool = True, regex: bool = False
+) -> int:
     # noinspection PyUnresolvedReferences, PyShadowingNames
     r"""Returns the number of reference points
 
@@ -792,7 +793,7 @@ def _refcount(ring: Sequence[Element], refpts: Refpts,
         121
 
         Returns *len(ring)*
-        """
+    """
     if isinstance(refpts, type):
         checkfun = checktype(refpts)
     elif callable(refpts):
@@ -808,8 +809,7 @@ def _refcount(ring: Sequence[Element], refpts: Refpts,
 
 
 # noinspection PyUnusedLocal,PyIncorrectDocstring
-def get_elements(ring: Sequence[Element], refpts: Refpts,
-                 regex: bool = False) -> list:
+def get_elements(ring: Sequence[Element], refpts: Refpts, regex: bool = False) -> list:
     r"""Returns a list of elements selected by *key*.
 
     Deprecated: :pycode:`get_elements(ring, refpts)` is :pycode:`ring[refpts]`
@@ -827,9 +827,13 @@ def get_elements(ring: Sequence[Element], refpts: Refpts,
     return list(refpts_iterator(ring, refpts, regex=regex))
 
 
-def get_value_refpts(ring: Sequence[Element], refpts: Refpts,
-                     attrname: str, index: Optional[int] = None,
-                     regex: bool = False):
+def get_value_refpts(
+    ring: Sequence[Element],
+    refpts: Refpts,
+    attrname: str,
+    index: Optional[int] = None,
+    regex: bool = False,
+):
     r"""Extracts attribute values from selected
         lattice :py:class:`.Element`\ s.
 
@@ -847,14 +851,21 @@ def get_value_refpts(ring: Sequence[Element], refpts: Refpts,
         attrvalues: numpy Array of attribute values.
     """
     getf = getval(attrname, index=index)
-    return numpy.array([getf(elem) for elem in refpts_iterator(ring, refpts,
-                                                               regex=regex)])
-
-
-def set_value_refpts(ring: Sequence[Element], refpts: Refpts,
-                     attrname: str, attrvalues, index: Optional[int] = None,
-                     increment: bool = False,
-                     copy: bool = False, regex: bool = False):
+    return numpy.array(
+        [getf(elem) for elem in refpts_iterator(ring, refpts, regex=regex)]
+    )
+
+
+def set_value_refpts(
+    ring: Sequence[Element],
+    refpts: Refpts,
+    attrname: str,
+    attrvalues,
+    index: Optional[int] = None,
+    increment: bool = False,
+    copy: bool = False,
+    regex: bool = False,
+):
     r"""Set the values of an attribute of an array of elements based on
     their refpts
 
@@ -885,13 +896,11 @@ def set_value_refpts(ring: Sequence[Element], refpts: Refpts,
     """
     setf = setval(attrname, index=index)
     if increment:
-        attrvalues += get_value_refpts(ring, refpts,
-                                       attrname, index=index,
-                                       regex=regex)
+        attrvalues += get_value_refpts(ring, refpts, attrname, index=index, regex=regex)
     else:
-        attrvalues = numpy.broadcast_to(attrvalues,
-                                        (_refcount(ring, refpts,
-                                                   regex=regex),))
+        attrvalues = numpy.broadcast_to(
+            attrvalues, (_refcount(ring, refpts, regex=regex),)
+        )
 
     # noinspection PyShadowingNames
     @make_copy(copy)
@@ -902,8 +911,9 @@ def apply(ring, refpts, values, regex):
     return apply(ring, refpts, attrvalues, regex)
 
 
-def get_s_pos(ring: Sequence[Element], refpts: Refpts = All,
-              regex: bool = False) -> Sequence[float]:
+def get_s_pos(
+    ring: Sequence[Element], refpts: Refpts = All, regex: bool = False
+) -> Sequence[float]:
     # noinspection PyUnresolvedReferences
     r"""Returns the locations of selected elements
 
@@ -923,16 +933,21 @@ def get_s_pos(ring: Sequence[Element], refpts: Refpts = All,
         array([26.37428795])
 
         Position at the end of the last element: length of the lattice
-        """
+    """
     # Positions at the end of each element.
-    s_pos = numpy.cumsum([getattr(el, 'Length', 0.0) for el in ring])
+    s_pos = numpy.cumsum([getattr(el, "Length", 0.0) for el in ring])
     # Prepend position at the start of the first element.
     s_pos = numpy.concatenate(([0.0], s_pos))
     return s_pos[get_bool_index(ring, refpts, regex=regex)]
 
 
-def rotate_elem(elem: Element, tilt: float = 0.0, pitch: float = 0.0,
-                yaw: float = 0.0, relative: bool = False) -> None:
+def rotate_elem(
+    elem: Element,
+    tilt: float = 0.0,
+    pitch: float = 0.0,
+    yaw: float = 0.0,
+    relative: bool = False,
+) -> None:
     r"""Set the tilt, pitch and yaw angle of an :py:class:`.Element`.
     The tilt is a rotation around the *s*-axis, the pitch is a
     rotation around the *x*-axis and the yaw is a rotation around
@@ -964,13 +979,14 @@ def rotate_elem(elem: Element, tilt: float = 0.0, pitch: float = 0.0,
         relative:       If :py:obj:`True`, the rotation is added to the
           previous one
     """
+
     # noinspection PyShadowingNames
     def _get_rm_tv(le, tilt, pitch, yaw):
         tilt = numpy.around(tilt, decimals=15)
         pitch = numpy.around(pitch, decimals=15)
         yaw = numpy.around(yaw, decimals=15)
         ct, st = numpy.cos(tilt), numpy.sin(tilt)
-        ap, ay = 0.5*le*numpy.tan(pitch), 0.5*le*numpy.tan(yaw)
+        ap, ay = 0.5 * le * numpy.tan(pitch), 0.5 * le * numpy.tan(yaw)
         rr1 = numpy.asfortranarray(numpy.diag([ct, ct, ct, ct, 1.0, 1.0]))
         rr1[0, 2] = st
         rr1[1, 3] = st
@@ -979,10 +995,10 @@ def _get_rm_tv(le, tilt, pitch, yaw):
         rr2 = rr1.T
         t1 = numpy.array([ay, numpy.sin(-yaw), -ap, numpy.sin(pitch), 0, 0])
         t2 = numpy.array([ay, numpy.sin(yaw), -ap, numpy.sin(-pitch), 0, 0])
-        rt1 = numpy.eye(6, order='F')
+        rt1 = numpy.eye(6, order="F")
         rt1[1, 4] = t1[1]
         rt1[3, 4] = t1[3]
-        rt2 = numpy.eye(6, order='F')
+        rt2 = numpy.eye(6, order="F")
         rt2[1, 4] = t2[1]
         rt2[3, 4] = t2[3]
         return rr1 @ rt1, rt2 @ rr2, t1, t2
@@ -992,17 +1008,17 @@ def _get_rm_tv(le, tilt, pitch, yaw):
     yaw0 = 0.0
     t10 = numpy.zeros(6)
     t20 = numpy.zeros(6)
-    if hasattr(elem, 'R1') and hasattr(elem, 'R2'):
-        rr10 = numpy.eye(6, order='F')
+    if hasattr(elem, "R1") and hasattr(elem, "R2"):
+        rr10 = numpy.eye(6, order="F")
         rr10[:4, :4] = elem.R1[:4, :4]
         rt10 = rr10.T @ elem.R1
         tilt0 = numpy.arctan2(rr10[0, 2], rr10[0, 0])
         yaw0 = numpy.arcsin(-rt10[1, 4])
         pitch0 = numpy.arcsin(rt10[3, 4])
         _, _, t10, t20 = _get_rm_tv(elem.Length, tilt0, pitch0, yaw0)
-    if hasattr(elem, 'T1') and hasattr(elem, 'T2'):
-        t10 = elem.T1-t10
-        t20 = elem.T2-t20
+    if hasattr(elem, "T1") and hasattr(elem, "T2"):
+        t10 = elem.T1 - t10
+        t20 = elem.T2 - t20
     if relative:
         tilt += tilt0
         pitch += pitch0
@@ -1011,8 +1027,8 @@ def _get_rm_tv(le, tilt, pitch, yaw):
     r1, r2, t1, t2 = _get_rm_tv(elem.Length, tilt, pitch, yaw)
     elem.R1 = r1
     elem.R2 = r2
-    elem.T1 = t1+t10
-    elem.T2 = t2+t20
+    elem.T1 = t1 + t10
+    elem.T2 = t2 + t20
 
 
 def tilt_elem(elem: Element, rots: float, relative: bool = False) -> None:
@@ -1037,8 +1053,9 @@ def tilt_elem(elem: Element, rots: float, relative: bool = False) -> None:
     rotate_elem(elem, tilt=rots, relative=relative)
 
 
-def shift_elem(elem: Element, deltax: float = 0.0, deltaz: float = 0.0,
-               relative: bool = False) -> None:
+def shift_elem(
+    elem: Element, deltax: float = 0.0, deltaz: float = 0.0, relative: bool = False
+) -> None:
     r"""Sets the transverse displacement of an :py:class:`.Element`
 
     The translation vectors are stored in the :pycode:`T1` and :pycode:`T2`
@@ -1052,7 +1069,7 @@ def shift_elem(elem: Element, deltax: float = 0.0, deltaz: float = 0.0,
           existing one
     """
     tr = numpy.array([deltax, 0.0, deltaz, 0.0, 0.0, 0.0])
-    if relative and hasattr(elem, 'T1') and hasattr(elem, 'T2'):
+    if relative and hasattr(elem, "T1") and hasattr(elem, "T2"):
         elem.T1 -= tr
         elem.T2 += tr
     else:
@@ -1060,8 +1077,9 @@ def shift_elem(elem: Element, deltax: float = 0.0, deltaz: float = 0.0,
         elem.T2 = tr
 
 
-def set_rotation(ring: Sequence[Element], tilts=0.0,
-                 pitches=0.0, yaws=0.0, relative=False) -> None:
+def set_rotation(
+    ring: Sequence[Element], tilts=0.0, pitches=0.0, yaws=0.0, relative=False
+) -> None:
     r"""Sets the tilts of a list of elements.
 
     Parameters:
@@ -1115,12 +1133,13 @@ def set_shift(ring: Sequence[Element], dxs, dzs, relative=False) -> None:
         shift_elem(el, dx, dy, relative=relative)
 
 
-def get_geometry(ring: List[Element],
-                 refpts: Refpts = All,
-                 start_coordinates: Tuple[float, float, float] = (0, 0, 0),
-                 centered: bool = False,
-                 regex: bool = False
-                 ):
+def get_geometry(
+    ring: list[Element],
+    refpts: Refpts = All,
+    start_coordinates: tuple[float, float, float] = (0, 0, 0),
+    centered: bool = False,
+    regex: bool = False,
+):
     # noinspection PyShadowingNames
     r"""Compute the 2D ring geometry in cartesian coordinates
 
@@ -1148,15 +1167,13 @@ def get_geometry(ring: List[Element],
        >>> geomdata, radius = get_geometry(ring)
     """
 
-    geom_dtype = [("x", numpy.float64),
-                  ("y", numpy.float64),
-                  ("angle", numpy.float64)]
+    geom_dtype = [("x", numpy.float64), ("y", numpy.float64), ("angle", numpy.float64)]
     boolrefs = get_bool_index(ring, refpts, endpoint=True, regex=regex)
     nrefs = refpts_count(boolrefs, len(ring))
-    geomdata = numpy.recarray((nrefs, ), dtype=geom_dtype)
-    xx = numpy.zeros(len(ring)+1)
-    yy = numpy.zeros(len(ring)+1)
-    angle = numpy.zeros(len(ring)+1)
+    geomdata = numpy.recarray((nrefs,), dtype=geom_dtype)
+    xx = numpy.zeros(len(ring) + 1)
+    yy = numpy.zeros(len(ring) + 1)
+    angle = numpy.zeros(len(ring) + 1)
     x0, y0, t0 = start_coordinates
     x, y = 0.0, 0.0
     t = t0
@@ -1168,30 +1185,30 @@ def get_geometry(ring: List[Element],
         ll = el.Length
         if isinstance(el, Dipole) and el.BendingAngle != 0:
             ang = 0.5 * el.BendingAngle
-            ll *= numpy.sin(ang)/ang
+            ll *= numpy.sin(ang) / ang
         else:
             ang = 0.0
         t -= ang
         x += ll * numpy.cos(t)
         y += ll * numpy.sin(t)
         t -= ang
-        xx[ind+1] = x
-        yy[ind+1] = y
-        angle[ind+1] = t
+        xx[ind + 1] = x
+        yy[ind + 1] = y
+        angle[ind + 1] = t
 
     dff = (t + _GEOMETRY_EPSIL) % (2.0 * numpy.pi) - _GEOMETRY_EPSIL
     if abs(dff) < _GEOMETRY_EPSIL:
         xcenter = numpy.mean(xx)
         ycenter = numpy.mean(yy)
-    elif abs(dff-numpy.pi) < _GEOMETRY_EPSIL:
-        xcenter = 0.5*x
-        ycenter = 0.5*y
+    elif abs(dff - numpy.pi) < _GEOMETRY_EPSIL:
+        xcenter = 0.5 * x
+        ycenter = 0.5 * y
     else:
-        num = numpy.cos(t)*x + numpy.sin(t)*y
-        den = numpy.sin(t-t0)
-        xcenter = -num*numpy.sin(t0)/den
-        ycenter = num*numpy.cos(t0)/den
-    radius = numpy.sqrt(xcenter*xcenter + ycenter*ycenter)
+        num = numpy.cos(t) * x + numpy.sin(t) * y
+        den = numpy.sin(t - t0)
+        xcenter = -num * numpy.sin(t0) / den
+        ycenter = num * numpy.cos(t0) / den
+    radius = numpy.sqrt(xcenter * xcenter + ycenter * ycenter)
     if centered:
         xx -= xcenter
         yy -= ycenter
diff --git a/pyat/at/lattice/variables.py b/pyat/at/lattice/variables.py
index 5b90e5f5a..08c7a0897 100644
--- a/pyat/at/lattice/variables.py
+++ b/pyat/at/lattice/variables.py
@@ -94,6 +94,7 @@ def _getfun(self, **kwargs):
 from typing import Union
 
 import numpy as np
+import numpy.typing as npt
 
 Number = Union[int, float]
 
@@ -403,6 +404,12 @@ class VariableList(list):
     appending, insertion or concatenation with the "+" operator.
     """
 
+    def __getitem__(self, index):
+        if isinstance(index, slice):
+            return VariableList(super().__getitem__(index))
+        else:
+            return super().__getitem__(index)
+
     def get(self, ring=None, **kwargs) -> Sequence[float]:
         r"""Get the current values of Variables
 
@@ -453,6 +460,12 @@ def __str__(self) -> str:
         return self.status()
 
     @property
-    def deltas(self) -> Sequence[Number]:
+    def deltas(self) -> npt.NDArray[Number]:
         """delta values of the variables"""
         return np.array([var.delta for var in self])
+
+    @deltas.setter
+    def deltas(self, value: Number | Sequence[Number]) -> None:
+        deltas = np.broadcast_to(value, len(self))
+        for var, delta in zip(self, deltas):
+            var.delta = delta
diff --git a/pyat/at/latticetools/__init__.py b/pyat/at/latticetools/__init__.py
index b1b34e2fb..a634efcb2 100644
--- a/pyat/at/latticetools/__init__.py
+++ b/pyat/at/latticetools/__init__.py
@@ -1,5 +1,6 @@
-"""Defines classes for modifying a lattice and observing its parameters"""
+"""Defines classes for modifying a lattice and observing its parameters."""
 
 from .observables import *
 from .observablelist import *
-from .matching import *
+# from .matching import *
+from .response_matrix import *
diff --git a/pyat/at/latticetools/observables.py b/pyat/at/latticetools/observables.py
index e0ba3c726..8f3012a3c 100644
--- a/pyat/at/latticetools/observables.py
+++ b/pyat/at/latticetools/observables.py
@@ -312,6 +312,9 @@ def evaluate(self, *data, initial: bool = False):
               sent to the evaluation function
             initial:    It :py:obj:`None`, store the result as the initial
               value
+
+        Returns:
+            value:      The value of the observable.
         """
         for d in data:
             if isinstance(d, Exception):
@@ -339,6 +342,10 @@ def weight(self):
         """Observable weight."""
         return np.broadcast_to(self.w, np.asarray(self._value).shape)
 
+    @weight.setter
+    def weight(self, w):
+        self.w = w
+
     @property
     def weighted_value(self):
         """Weighted value of the Observable, computed as
@@ -616,8 +623,8 @@ def __init__(
 
             Observe the horizontal closed orbit at monitor locations
         """
-        name = self._set_name(name, "orbit", axis_(axis, "code"))
-        fun = _ArrayAccess(axis_(axis, "index"))
+        name = self._set_name(name, "orbit", axis_(axis, key="code"))
+        fun = _ArrayAccess(axis_(axis, key="index"))
         needs = {Need.ORBIT}
         super().__init__(fun, refpts, needs=needs, name=name, **kwargs)
 
@@ -666,8 +673,8 @@ def __init__(
             Observe the transfer matrix from origin to monitor locations and
             extract T[0,1]
         """
-        name = self._set_name(name, "matrix", axis_(axis, "code"))
-        fun = _ArrayAccess(axis_(axis, "index"))
+        name = self._set_name(name, "matrix", axis_(axis, key="code"))
+        fun = _ArrayAccess(axis_(axis, key="index"))
         needs = {Need.MATRIX}
         super().__init__(fun, refpts, needs=needs, name=name, **kwargs)
 
@@ -700,12 +707,12 @@ def __init__(
         shape of *value*.
         """
         needs = {Need.GLOBALOPTICS}
-        name = self._set_name(name, param, plane_(plane, "code"))
+        name = self._set_name(name, param, plane_(plane, key="code"))
         if callable(param):
             fun = param
             needs.add(Need.CHROMATICITY)
         else:
-            fun = _RecordAccess(param, plane_(plane, "index"))
+            fun = _RecordAccess(param, plane_(plane, key="index"))
             if param == "chromaticity":
                 needs.add(Need.CHROMATICITY)
         super().__init__(fun, needs=needs, name=name, **kwargs)
@@ -803,11 +810,11 @@ def __init__(
             ax_ = plane_
 
         needs = {Need.LOCALOPTICS}
-        name = self._set_name(name, param, ax_(plane, "code"))
+        name = self._set_name(name, param, ax_(plane, key="code"))
         if callable(param):
             fun = param
         else:
-            fun = _RecordAccess(param, _all_rows(ax_(plane, "index")))
+            fun = _RecordAccess(param, _all_rows(ax_(plane, key="index")))
             if param == "mu" or all_points:
                 needs.add(Need.ALL_POINTS)
             if param in {"W", "Wp", "dalpha", "dbeta", "dmu", "ddispersion", "dR"}:
@@ -843,7 +850,9 @@ def __init__(
 
         Example:
 
-            >>> obs = LatticeObservable(at.Sextupole, "KickAngle", index=0, statfun=np.sum)
+            >>> obs = LatticeObservable(
+            ...     at.Sextupole, "KickAngle", index=0, statfun=np.sum
+            ... )
 
             Observe the sum of horizontal kicks in Sextupoles
         """
@@ -888,8 +897,8 @@ def __init__(
         The *target*, *weight* and *bounds* inputs must be broadcastable to the
         shape of *value*.
         """
-        name = self._set_name(name, "trajectory", axis_(axis, "code"))
-        fun = _ArrayAccess(axis_(axis, "index"))
+        name = self._set_name(name, "trajectory", axis_(axis, key="code"))
+        fun = _ArrayAccess(axis_(axis, key="index"))
         needs = {Need.TRAJECTORY}
         super().__init__(fun, refpts, needs=needs, name=name, **kwargs)
 
@@ -941,11 +950,11 @@ def __init__(
 
             Observe the horizontal emittance
         """
-        name = self._set_name(name, param, plane_(plane, "code"))
+        name = self._set_name(name, param, plane_(plane, key="code"))
         if callable(param):
             fun = param
         else:
-            fun = _RecordAccess(param, plane_(plane, "index"))
+            fun = _RecordAccess(param, plane_(plane, key="index"))
         needs = {Need.EMITTANCE}
         super().__init__(fun, needs=needs, name=name, **kwargs)
 
@@ -1008,10 +1017,10 @@ def GlobalOpticsObservable(
     """
     if param == "tune" and use_integer:
         # noinspection PyProtectedMember
-        name = ElementObservable._set_name(name, param, plane_(plane, "code"))
+        name = ElementObservable._set_name(name, param, plane_(plane, key="code"))
         return LocalOpticsObservable(
             End,
-            _Tune(plane_(plane, "index")),
+            _Tune(plane_(plane, key="index")),
             name=name,
             summary=True,
             all_points=True,
diff --git a/pyat/at/latticetools/response_matrix.py b/pyat/at/latticetools/response_matrix.py
new file mode 100644
index 000000000..ddfad19b9
--- /dev/null
+++ b/pyat/at/latticetools/response_matrix.py
@@ -0,0 +1,1195 @@
+# noinspection PyUnresolvedReferences
+r"""Definition of :py:class:`.ResponseMatrix` objects.
+
+A :py:class:`ResponseMatrix` object defines a general-purpose response matrix, based
+on a :py:class:`.VariableList` of attributes which will be independently varied, and an
+:py:class:`.ObservableList` of attributes which will be recorded for each
+variable step.
+
+:py:class:`ResponseMatrix` objects can be combined with the "+" operator to define
+combined responses. This concatenates the variables and the observables.
+
+This module also defines two commonly used response matrices:
+:py:class:`OrbitResponseMatrix` for circular machines and
+:py:class:`TrajectoryResponseMatrix` for beam lines. Other matrices can be easily
+defined by providing the desired Observables and Variables to the
+:py:class:`ResponseMatrix` base class.
+
+Generic response matrix
+-----------------------
+
+The :py:class:`ResponseMatrix` class defines a general-purpose response matrix, based
+on a :py:class:`.VariableList` of quantities which will be independently varied, and an
+:py:class:`.ObservableList` of quantities which will be recorded for each step.
+
+For instance let's take the horizontal displacements of all quadrupoles as variables:
+
+>>> variables = VariableList(
+...     RefptsVariable(ik, "dx", name=f"dx_{ik}", delta=0.0001)
+...     for ik in ring.get_uint32_index(at.Quadrupole)
+... )
+
+The variables are the horizontal displacement ``dx`` of all quadrupoles. The variable
+name is set to *dx_nnnn* where *nnnn* is the index of the quadrupole in the lattice.
+The step is set to 0.0001 m.
+
+Let's take the horizontal positions at all beam position monitors as observables:
+
+>>> observables = at.ObservableList([at.OrbitObservable(at.Monitor, axis="x")])
+
+This is a single observable named *orbit[x]* by default, with multiple values.
+
+Instantiation
+^^^^^^^^^^^^^
+
+>>> resp_dx = at.ResponseMatrix(ring, variables, observables)
+
+At that point, the response matrix is empty.
+
+Matrix Building
+^^^^^^^^^^^^^^^
+
+The response matrix may be filled by several means:
+
+#. Direct assignment of an array to the :py:attr:`~.ResponseMatrix.response` property.
+   The shape of the array is checked.
+#. :py:meth:`~ResponseMatrix.load` loads data from a file containing previously
+   saved values or experimentally measured values,
+#. :py:meth:`~ResponseMatrix.build_tracking` computes the matrix using tracking,
+#. For some specialized response matrices a
+   :py:meth:`~OrbitResponseMatrix.build_analytical` method is available.
+
+Matrix normalisation
+^^^^^^^^^^^^^^^^^^^^
+
+To be correctly inverted, the response matrix must be correctly normalised: the norms
+of its columns must be of the same order of magnitude, and similarly for the rows.
+
+Normalisation is done by adjusting the weights :math:`w_v` for the variables
+:math:`\mathbf{V}` and :math:`w_o` for the observables :math:`\mathbf{O}`.
+With :math:`\mathbf{R}` the response matrix:
+
+.. math::
+
+   \mathbf{O} = \mathbf{R} . \mathbf{V}
+
+The weighted response matrix :math:`\mathbf{R}_w` is:
+
+.. math::
+
+   \frac{\mathbf{O}}{w_o} = \mathbf{R}_w . \frac{\mathbf{V}}{w_v}
+
+The :math:`\mathbf{R}_w` is dimensionless and should be normalised. This can be checked
+using:
+
+* :py:meth:`~ResponseMatrix.check_norm` which prints the ratio of the maximum / minimum
+  norms for variables and observables. These should be less than 10.
+* :py:meth:`~.ResponseMatrix.plot_norm`
+
+Both natural and weighted response matrices can be retrieved with the
+:py:attr:`~ResponseMatrix.response` and :py:attr:`~ResponseMatrix.weighted_response`
+properties.
+
+Matrix pseudo-inversion
+^^^^^^^^^^^^^^^^^^^^^^^
+
+The :py:meth:`~ResponseMatrix.solve` method computes the singular values of the
+weighted response matrix.
+
+After solving, correction is available, for instance with
+
+* :py:meth:`~ResponseMatrix.correction_matrix` which returns the correction matrix
+  (pseudo-inverse of the response matrix),
+* :py:meth:`~ResponseMatrix.get_correction` which returns a correction vector when
+  given error values,
+* :py:meth:`~ResponseMatrix.correct` which computes and optionally applies a correction
+  for the provided :py:class:`.Lattice`.
+
+Exclusion of variables and observables
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+Variables may be added to a set of excluded values, and similarly for observables.
+Excluding an item does not change the response matrix. The values are excluded from the
+pseudo-inversion of the response, possibly reducing the number of singular values.
+After inversion the correction matrix is expanded to its original size by inserting
+zero lines and columns at the location of excluded items. This way:
+
+- error and correction vectors keep the same size independently of excluded values,
+- excluded error values are ignored,
+- excluded corrections are set to zero.
+
+Variables can be added to the set of excluded variables using
+:py:meth:`~.ResponseMatrix.exclude_vars` and observables using
+:py:meth:`~.ResponseMatrix.exclude_obs`.
+
+After excluding items, the pseudo-inverse is discarded so one must recompute it again
+by calling :py:meth:`~ResponseMatrix.solve`.
+
+The exclusion masks can be reset with  :py:meth:`~.ResponseMatrix.reset_vars` and
+:py:meth:`~.ResponseMatrix.reset_obs`.
+"""
+
+from __future__ import annotations
+
+__all__ = [
+    "sequence_split",
+    "ResponseMatrix",
+    "OrbitResponseMatrix",
+    "TrajectoryResponseMatrix",
+]
+
+import os
+import multiprocessing
+import concurrent.futures
+import abc
+import warnings
+from collections.abc import Sequence, Generator, Callable
+from typing import Any, ClassVar
+from itertools import chain
+from functools import partial
+import math
+
+import numpy as np
+import numpy.typing as npt
+
+from .observables import ElementObservable
+from .observables import TrajectoryObservable, OrbitObservable, LatticeObservable
+from .observables import LocalOpticsObservable, GlobalOpticsObservable
+from .observablelist import ObservableList
+from ..lattice import AtError, AtWarning, Refpts, Uint32Refpts, All
+from ..lattice import AxisDef, plane_, Lattice, Monitor, checkattr
+from ..lattice.lattice_variables import RefptsVariable
+from ..lattice.variables import VariableList
+
+FloatArray = npt.NDArray[np.float64]
+
+_orbit_correctors = checkattr("KickAngle")
+
+_globring: Lattice | None = None
+_globobs: ObservableList | None = None
+
+warnings.filterwarnings("always", category=AtWarning, module=__name__)
+
+
+def sequence_split(seq: Sequence, nslices: int) -> Generator[Sequence, None, None]:
+    """Split a sequence into multiple sub-sequences.
+
+    The length of *seq* does not have to be a multiple of *nslices*.
+
+    Args:
+        seq: sequence to split
+        nslices: number of sub-sequences
+
+    Returns:
+        subseqs: Iterator over sub-sequences
+    """
+
+    def _split(seqsizes):
+        beg = 0
+        for size in seqsizes:
+            end = beg + size
+            yield seq[beg:end]
+            beg = end
+
+    lna = len(seq)
+    sz, rem = divmod(lna, nslices)
+    lsubseqs = [sz] * nslices
+    for k in range(rem):
+        lsubseqs[k] += 1
+    return _split(lsubseqs)
+
+
+def _resp(
+    ring: Lattice, observables: ObservableList, variables: VariableList, **kwargs
+):
+    def _resp_one(variable: RefptsVariable):
+        """Single response"""
+        variable.step_up(ring=ring)
+        observables.evaluate(ring, **kwargs)
+        op = observables.flat_values
+        variable.step_down(ring=ring)
+        observables.evaluate(ring, **kwargs)
+        om = observables.flat_values
+        variable.reset(ring=ring)
+        return (op - om) / (2.0 * variable.delta)
+
+    return [_resp_one(v) for v in variables]
+
+
+def _resp_fork(variables: VariableList, **kwargs):
+    """Response for fork parallel method."""
+    return _resp(_globring, _globobs, variables, **kwargs)
+
+
+class _SvdSolver(abc.ABC):
+    """SVD solver for response matrices."""
+
+    _shape: tuple[int, int]
+    _obsmask: npt.NDArray[bool]
+    _varmask: npt.NDArray[bool]
+    _response: FloatArray | None = None
+    _v: FloatArray | None = None
+    _uh: FloatArray | None = None
+    #: Singular values of the response matrix
+    singular_values: FloatArray | None = None
+
+    def __init__(self, nobs: int, nvar: int):
+        self._shape = (nobs, nvar)
+        self._obsmask = np.ones(nobs, dtype=bool)
+        self._varmask = np.ones(nvar, dtype=bool)
+
+    def reset_vars(self):
+        """Reset the variable exclusion mask: enable all variables"""
+        self._varmask = np.ones(self.shape[1], dtype=bool)
+        self._v = None
+        self._uh = None
+        self.singular_values = None
+
+    def reset_obs(self):
+        """Reset the observable exclusion mask: enable all observables"""
+        self._obsmask = np.ones(self.shape[0], dtype=bool)
+        self._v = None
+        self._uh = None
+        self.singular_values = None
+
+    @property
+    @abc.abstractmethod
+    def varweights(self) -> np.ndarray: ...
+
+    @property
+    @abc.abstractmethod
+    def obsweights(self) -> np.ndarray: ...
+
+    @property
+    def shape(self) -> tuple[int, int]:
+        """Shape of the response matrix."""
+        return self._shape
+
+    def solve(self) -> None:
+        """Compute the singular values of the response matrix."""
+        resp = self.weighted_response
+        selected = np.ix_(self._obsmask, self._varmask)
+        u, s, vh = np.linalg.svd(resp[selected], full_matrices=False)
+        self._v = vh.T * (1.0 / s) * self.varweights[self._varmask].reshape(-1, 1)
+        self._uh = u.T / self.obsweights[self._obsmask]
+        self.singular_values = s
+
+    def check_norm(self) -> tuple[FloatArray, FloatArray]:
+        """Display the norm of the rows and columns of the weighted response matrix.
+
+        Adjusting the variables and observable weights to equalize the norms
+        of rows and columns is important.
+
+        Returns:
+            obs_norms:      Norms of observables (rows)
+            var_norms:      Norms of Variables (columns)
+        """
+        resp = self.weighted_response
+        obs = np.linalg.norm(resp, axis=1)
+        var = np.linalg.norm(resp, axis=0)
+        print(f"max/min Observables: {np.amax(obs) / np.amin(obs)}")
+        print(f"max/min Variables: {np.amax(var) / np.amin(var)}")
+        return obs, var
+
+    @property
+    def response(self) -> FloatArray:
+        """Response matrix."""
+        resp = self._response
+        if resp is None:
+            raise AtError("No matrix yet: run build() or load() first")
+        return resp
+
+    @response.setter
+    def response(self, response: FloatArray) -> None:
+        l1, c1 = self._shape
+        l2, c2 = response.shape
+        if l1 != l1 or c1 != c2:
+            raise ValueError(
+                f"Input matrix has incompatible shape. Expected: {self.shape}"
+            )
+        self._response = response
+
+    @property
+    def weighted_response(self) -> FloatArray:
+        """Weighted response matrix."""
+        return self.response * (self.varweights / self.obsweights.reshape(-1, 1))
+
+    def correction_matrix(self, nvals: int | None = None) -> FloatArray:
+        """Return the correction matrix (pseudo-inverse of the response matrix).
+
+        Args:
+            nvals:  Desired number of singular values. If :py:obj:`None`, use
+              all singular values
+
+        Returns:
+            cormat: Correction matrix
+        """
+        if self.singular_values is None:
+            self.solve()
+        if nvals is None:
+            nvals = len(self.singular_values)
+        cormat = np.zeros(self._shape[::-1])
+        selected = np.ix_(self._varmask, self._obsmask)
+        cormat[selected] = self._v[:, :nvals] @ self._uh[:nvals, :]
+        return cormat
+
+    def get_correction(
+        self, observed: FloatArray, nvals: int | None = None
+    ) -> FloatArray:
+        """Compute the correction of the given observation.
+
+        Args:
+            observed:   Vector of observed deviations,
+            nvals:      Desired number of singular values. If :py:obj:`None`, use
+              all singular values
+
+        Returns:
+            corr:       Correction vector
+        """
+        return -self.correction_matrix(nvals=nvals) @ observed
+
+    def save(self, file) -> None:
+        """Save a response matrix in the NumPy .npy format.
+
+        Args:
+            file:   file-like object, string, or :py:class:`pathlib.Path`: File to
+              which the data is saved. If file is a file-object, it must be opened in
+              binary mode. If file is a string or Path, a .npy extension will
+              be appended to the filename if it does not already have one.
+        """
+        if self._response is None:
+            raise AtError("No response matrix: run build_tracking() or load() first")
+        np.save(file, self._response)
+
+    def load(self, file) -> None:
+        """Load a response matrix saved in the NumPy .npy format.
+
+        Args:
+            file:   file-like object, string, or :py:class:`pathlib.Path`: the file to
+              read. A file object must always be opened in binary mode.
+        """
+        self.response = np.load(file)
+
+
+class ResponseMatrix(_SvdSolver):
+    r"""Base class for response matrices.
+
+    It is defined by any arbitrary set of :py:class:`~.variables.VariableBase`\ s and
+    :py:class:`.Observable`\s
+
+    Addition is defined on :py:class:`ResponseMatrix` objects as the addition
+    of their :py:class:`~.variables.VariableBase`\ s and :py:class:`.Observable`\s to
+    produce combined responses.
+    """
+
+    ring: Lattice
+    variables: VariableList  #: List of matrix :py:class:`Variable <.VariableBase>`\ s
+    observables: ObservableList  #: List of matrix :py:class:`.Observable`\s
+    _eval_args: dict[str, Any] = {}
+
+    def __init__(
+        self,
+        ring: Lattice,
+        variables: VariableList,
+        observables: ObservableList,
+    ):
+        r"""
+        Args:
+            ring:           Design lattice, used to compute the response
+            variables:      List of :py:class:`Variable <.VariableBase>`\ s
+            observables:    List of :py:class:`.Observable`\s
+        """
+
+        def limits(obslist):
+            beg = 0
+            for obs in obslist:
+                end = beg + obs.value.size
+                yield beg, end
+                beg = end
+
+        # for efficiency of parallel computation, the variable's refpts must be integer
+        for var in variables:
+            var.refpts = ring.get_uint32_index(var.refpts)
+        self.ring = ring
+        self.variables = variables
+        self.observables = observables
+        variables.get(ring=ring, initial=True)
+        observables.evaluate(ring=ring, initial=True)
+        super().__init__(len(observables.flat_values), len(variables))
+        self._ob = [self._obsmask[beg:end] for beg, end in limits(self.observables)]
+
+    def __add__(self, other: ResponseMatrix):
+        if not isinstance(other, ResponseMatrix):
+            raise TypeError(
+                f"Cannot add {type(other).__name__} and {type(self).__name__}"
+            )
+        return ResponseMatrix(
+            self.ring,
+            VariableList(self.variables + other.variables),
+            self.observables + other.observables,
+        )
+
+    def __str__(self):
+        no, nv = self.shape
+        return f"{type(self).__name__}({no} observables, {nv} variables)"
+
+    @property
+    def varweights(self) -> np.ndarray:
+        """Variable weights."""
+        return self.variables.deltas
+
+    @property
+    def obsweights(self) -> np.ndarray:
+        """Observable weights."""
+        return self.observables.flat_weights
+
+    def correct(
+        self, ring: Lattice, nvals: int = None, niter: int = 1, apply: bool = False
+    ) -> FloatArray:
+        """Compute and optionally apply the correction.
+
+        Args:
+            ring:       Lattice description. The response matrix observables
+              will be evaluated for *ring* and the deviation from target will
+              be corrected
+            apply:      If :py:obj:`True`, apply the correction to *ring*
+            niter:      Number of iterations. For more than one iteration,
+              *apply* must be :py:obj:`True`
+            nvals:      Desired number of singular values. If :py:obj:`None`,
+              use all singular values. *nvals* may be a scalar or an iterable with
+              *niter* values.
+
+        Returns:
+            correction: Vector of correction values
+        """
+        if niter > 1 and not apply:
+            raise ValueError("If niter > 1, 'apply' must be True")
+        obs = self.observables
+        if apply:
+            self.variables.get(ring=ring, initial=True)
+        sumcorr = np.array([0.0])
+        for it, nv in zip(range(niter), np.broadcast_to(nvals, (niter,))):
+            print(f'step {it+1}, nvals = {nv}')
+            obs.evaluate(ring, **self._eval_args)
+            err = obs.flat_deviations
+            if np.any(np.isnan(err)):
+                raise AtError(
+                    f"Step {it + 1}: Invalid observables, cannot compute correction"
+                )
+            corr = self.get_correction(obs.flat_deviations, nvals=nv)
+            sumcorr = sumcorr + corr  # non-broadcastable sumcorr
+            if apply:
+                self.variables.increment(corr, ring=ring)
+        return sumcorr
+
+    def build_tracking(
+        self,
+        use_mp: bool = False,
+        pool_size: int | None = None,
+        start_method: str | None = None,
+        **kwargs,
+    ) -> FloatArray:
+        """Build the response matrix.
+
+        Args:
+            use_mp:             Use multiprocessing
+            pool_size:          number of processes. If None,
+              :pycode:`min(len(self.variables, nproc)` is used
+            start_method:       python multiprocessing start method.
+              :py:obj:`None` uses the python default that is considered safe.
+              Available values: ``'fork'``, ``'spawn'``, ``'forkserver'``.
+              Default for linux is ``'fork'``, default for macOS and  Windows
+              is ``'spawn'``. ``'fork'`` may be used on macOS to speed up the
+              calculation, however it is considered unsafe.
+
+        Keyword Args:
+            dp (float):     Momentum deviation. Defaults to :py:obj:`None`
+            dct (float):    Path lengthening. Defaults to :py:obj:`None`
+            df (float):     Deviation from the nominal RF frequency.
+              Defaults to :py:obj:`None`
+            r_in (Orbit):   Initial trajectory, used for
+              :py:class:`TrajectoryResponseMatrix`, Default: zeros(6)
+
+        Returns:
+            response:       Response matrix
+        """
+        self._eval_args = kwargs
+        self.observables.evaluate(self.ring)
+        ring = self.ring.deepcopy()
+
+        if use_mp:
+            global _globring
+            global _globobs
+            ctx = multiprocessing.get_context(start_method)
+            if pool_size is None:
+                pool_size = min(len(self.variables), os.cpu_count())
+            obschunks = sequence_split(self.variables, pool_size)
+            if ctx.get_start_method() == "fork":
+                _globring = ring
+                _globobs = self.observables
+                _single_resp = partial(_resp_fork, **kwargs)
+            else:
+                _single_resp = partial(_resp, ring, self.observables, **kwargs)
+            with concurrent.futures.ProcessPoolExecutor(
+                max_workers=pool_size,
+                mp_context=ctx,
+            ) as pool:
+                results = list(chain(*pool.map(_single_resp, obschunks)))
+            _globring = None
+            _globobs = None
+        else:
+            results = _resp(ring, self.observables, self.variables, **kwargs)
+
+        resp = np.stack(results, axis=-1)
+        self.response = resp
+        return resp
+
+    def build_analytical(self) -> FloatArray:
+        """Build the response matrix."""
+        raise NotImplementedError(
+            f"build_analytical not implemented for {self.__class__.__name__}"
+        )
+
+    def _on_obs(self, fun: Callable, *args, obsid: int | str = 0):
+        """Apply a function to the selected observable"""
+        if not isinstance(obsid, str):
+            return fun(self.observables[obsid], *args)
+        else:
+            for obs in self.observables:
+                if obs.name == obsid:
+                    return fun(obs, *args)
+            else:
+                raise ValueError(f"Observable {obsid} not found")
+
+    def get_target(self, *, obsid: int | str = 0) -> FloatArray:
+        r"""Return the target of the specified observable
+
+        Args:
+            obsid:  :py:class:`.Observable` name or index in the observable list.
+
+        Returns:
+            target: observable target
+        """
+        def _get(obs):
+            return obs.target
+
+        return self._on_obs(_get, obsid=obsid)
+
+    def set_target(self, target: npt.ArrayLike, *, obsid: int | str = 0) -> None:
+        r"""Set the target of the specified observable
+
+        Args:
+            target: observable target. Must be broadcastable to the shape of the
+              observable value.
+            obsid:  :py:class:`.Observable` name or index in the observable list.
+        """
+
+        def _set(obs, targ):
+            obs.target = targ
+
+        return self._on_obs(_set, target, obsid=obsid)
+
+    def exclude_obs(self, *, obsid: int | str = 0, refpts: Refpts = None) -> None:
+        # noinspection PyUnresolvedReferences
+        r"""Add an observable item to the set of excluded values
+
+        After excluding observation points, the matrix must be inverted again using
+        :py:meth:`solve`.
+
+        Args:
+            obsid:      :py:class:`.Observable` name or index in the observable list.
+            refpts:     location of elements to exclude for
+              :py:class:`.ElementObservable` objects, otherwise ignored.
+
+        Raises:
+            ValueError: No observable with the given name.
+            IndexError: Observableindex out of range.
+
+        Example:
+            >>> resp = OrbitResponseMatrix(ring, "h", Monitor, Corrector)
+            >>> resp.exclude_obs(obsid="x_orbit", refpts="BPM_02")
+
+            Create an horizontal :py:class:`OrbitResponseMatrix` from
+            :py:class:`.Corrector` elements to :py:class:`.Monitor` elements,
+            and exclude the monitor with name "BPM_02"
+        """
+
+        def exclude(ob, msk):
+            inimask = msk.copy()
+            if isinstance(ob, ElementObservable) and not ob.summary:
+                boolref = self.ring.get_bool_index(refpts)
+                # noinspection PyProtectedMember
+                msk &= np.logical_not(boolref[ob._boolrefs])
+            else:
+                msk[:] = False
+            if np.all(msk == inimask):
+                warnings.warn(AtWarning("No new excluded value"), stacklevel=3)
+            # Force a new computation
+            self.singular_values = None
+
+        if not isinstance(obsid, str):
+            exclude(self.observables[obsid], self._ob[obsid])
+        else:
+            for obs, mask in zip(self.observables, self._ob):
+                if obs.name == obsid:
+                    exclude(obs, mask)
+                    break
+            else:
+                raise ValueError(f"Observable {obsid} not found")
+
+    @property
+    def excluded_obs(self) -> dict:
+        """Directory of excluded observables.
+
+        The dictionary keys are the observable names, the values are the integer
+        indices of excluded items (empty list if no exclusion).
+        """
+
+        def ex(obs, mask):
+            if isinstance(obs, ElementObservable) and not obs.summary:
+                refpts = self.ring.get_bool_index(None)
+                # noinspection PyProtectedMember
+                refpts[obs._boolrefs] = np.logical_not(mask)
+                refpts = self.ring.get_uint32_index(refpts)
+            else:
+                refpts = np.arange(0 if np.all(mask) else mask.size, dtype=np.uint32)
+            return refpts
+
+        return {ob.name: ex(ob, mask) for ob, mask in zip(self.observables, self._ob)}
+
+    def exclude_vars(self, *varid: int | str) -> None:
+        # noinspection PyUnresolvedReferences
+        """Add variables to the set of excluded variables.
+
+        Args:
+            *varid:  :py:class:`Variable <.VariableBase>` names or variable indices
+              in the variable list
+
+        After excluding variables, the matrix must be inverted again using
+        :py:meth:`solve`.
+
+        Examples:
+            >>> resp.exclude_vars(0, "var1", -1)
+
+            Exclude the 1st variable, the variable named "var1" and the last variable.
+        """
+        nameset = set(nm for nm in varid if isinstance(nm, str))
+        varidx = [nm for nm in varid if not isinstance(nm, str)]
+        mask = np.array([var.name in nameset for var in self.variables])
+        mask[varidx] = True
+        miss = nameset - {var.name for var, ok in zip(self.variables, mask) if ok}
+        if miss:
+            raise ValueError(f"Unknown variables: {miss}")
+        self._varmask &= np.logical_not(mask)
+
+    @property
+    def excluded_vars(self) -> list:
+        """List of excluded variables"""
+        return [var.name for var, ok in zip(self.variables, self._varmask) if not ok]
+
+
+class OrbitResponseMatrix(ResponseMatrix):
+    # noinspection PyUnresolvedReferences
+    r"""Orbit response matrix.
+
+    An :py:class:`OrbitResponseMatrix` applies to a single plane, horizontal or
+    vertical. A combined response matrix is obtained by adding horizontal and
+    vertical matrices. However, the resulting matrix has the :py:class:`ResponseMatrix`
+    class, which implies that the :py:class:`OrbitResponseMatrix` specific methods are
+    not available.
+
+    Variables are a set of steerers and optionally the RF frequency. Steerer
+    variables are named ``xnnnn`` or ``ynnnn`` where nnnn is the index in the
+    lattice. The RF frequency variable is named ``RF frequency``.
+
+    Observables are the closed orbit position at selected points, named
+    ``orbit[x]`` for the horizontal plane or ``orbit[y]`` for the vertical plane,
+    and optionally the sum of steerer angles named ``sum(h_kicks)`` or
+    ``sum(v_kicks)``
+
+    The variable elements must have the *KickAngle* attribute used for correction.
+    It's available for all magnets, though not present by default
+    except in :py:class:`.Corrector` magnets. For other magnets, the attribute
+    should be explicitly created.
+
+    By default, the observables are all the :py:class:`.Monitor` elements, and the
+    variables are all the elements having a *KickAngle* attribute.
+    This is equivalent to:
+
+    >>> resp_v = OrbitResponseMatrix(
+    ...     ring, "v", bpmrefs=at.Monitor, steerrefs=at.checkattr("KickAngle")
+    ... )
+    """
+
+    bpmrefs: Uint32Refpts  #: location of position monitors
+    steerrefs: Uint32Refpts  #: location of steerers
+
+    def __init__(
+        self,
+        ring: Lattice,
+        plane: AxisDef,
+        bpmrefs: Refpts = Monitor,
+        steerrefs: Refpts = _orbit_correctors,
+        *,
+        cavrefs: Refpts = None,
+        bpmweight: float | Sequence[float] = 1.0,
+        bpmtarget: float | Sequence[float] = 0.0,
+        steerdelta: float | Sequence[float] = 0.0001,
+        cavdelta: float | None = None,
+        steersum: bool = False,
+        stsumweight: float | None = None,
+    ):
+        """
+        Args:
+            ring:       Design lattice, used to compute the response.
+            plane:      One out of {0, 'x', 'h', 'H'} for horizontal orbit, or
+              one of {1, 'y', 'v', 'V'} for vertical orbit.
+            bpmrefs:    Location of closed orbit observation points.
+              See ":ref:`Selecting elements in a lattice <refpts>`".
+              Default: all :py:class:`.Monitor` elements.
+            steerrefs:  Location of orbit steerers. Their *KickAngle* attribute
+              is used and must be present in the selected elements.
+              Default: All Elements having a *KickAngle* attribute.
+            cavrefs:    Location of RF cavities. Their *Frequency* attribute
+              is used. If :py:obj:`None`, no cavity is included in the response.
+              Cavities must be active. Cavity variables are appended to the steerer
+              variables.
+            bpmweight:  Weight of position readings. Must be broadcastable to the
+              number of BPMs.
+            bpmtarget:  Target orbit position. Must be broadcastable to the number of
+              observation points.
+            cavdelta:   Step on RF frequency for matrix computation [Hz]. This
+              is also the cavity weight. Default: automatically computed.
+            steerdelta: Step on steerers for matrix computation [rad]. This is
+              also the steerer weight. Must be broadcastable to the number of steerers.
+            steersum:   If :py:obj:`True`, the sum of steerers is appended to the
+              Observables.
+            stsumweight: Weight on steerer summation. Default: automatically computed.
+
+        :ivar VariableList variables: matrix variables
+        :ivar ObservableList observables: matrix observables
+
+        By default, the weights of cavities and steerers summation are set to give
+        a factor 2 more efficiency than steerers and BPMs
+
+        """
+
+        def steerer(ik, delta):
+            name = f"{plcode}{ik:04}"
+            return RefptsVariable(ik, "KickAngle", index=pl, name=name, delta=delta)
+
+        def set_norm():
+            bpm = LocalOpticsObservable(bpmrefs, "beta", plane=pl)
+            sts = LocalOpticsObservable(steerrefs, "beta", plane=pl)
+            dsp = LocalOpticsObservable(bpmrefs, "dispersion", plane=2 * pl)
+            tun = GlobalOpticsObservable("tune", plane=pl)
+            obs = ObservableList([bpm, sts, dsp, tun])
+            result = obs.evaluate(ring=ring)
+            alpha = ring.disable_6d(copy=True).get_mcf(0)
+            freq = ring.get_rf_frequency(cavpts=cavrefs)
+            nr = np.outer(
+                np.sqrt(result[0]) / bpmweight, np.sqrt(result[1]) * steerdelta
+            )
+            vv = np.mean(np.linalg.norm(nr, axis=0))
+            vo = np.mean(np.linalg.norm(nr, axis=1))
+            korb = 0.25 * math.sqrt(2.0) / math.sin(math.pi * result[3])
+            cd = vv * korb * alpha * freq / np.linalg.norm(result[2] / bpmweight)
+            sw = np.linalg.norm(deltas) / vo / korb
+            return cd, sw
+
+        pl = plane_(plane, key="index")
+        plcode = plane_(plane, key="code")
+        ids = ring.get_uint32_index(steerrefs)
+        nbsteers = len(ids)
+        deltas = np.broadcast_to(steerdelta, nbsteers)
+        if steersum and stsumweight is None or cavrefs and cavdelta is None:
+            cavd, stsw = set_norm()
+
+        # Observables
+        bpms = OrbitObservable(bpmrefs, axis=2 * pl, target=bpmtarget, weight=bpmweight)
+        observables = ObservableList([bpms])
+        if steersum:
+            # noinspection PyUnboundLocalVariable
+            sumobs = LatticeObservable(
+                steerrefs,
+                "KickAngle",
+                name=f"{plcode}_kicks",
+                target=0.0,
+                index=pl,
+                weight=stsumweight if stsumweight else stsw / 2.0,
+                statfun=np.sum,
+            )
+            observables.append(sumobs)
+
+        # Variables
+        variables = VariableList(steerer(ik, delta) for ik, delta in zip(ids, deltas))
+        if cavrefs is not None:
+            active = (el.longt_motion for el in ring.select(cavrefs))
+            if not all(active):
+                raise ValueError("Cavities are not active")
+            # noinspection PyUnboundLocalVariable
+            cavvar = RefptsVariable(
+                cavrefs,
+                "Frequency",
+                name="RF frequency",
+                delta=cavdelta if cavdelta else 2.0 * cavd,
+            )
+            variables.append(cavvar)
+
+        super().__init__(ring, variables, observables)
+        self.plane = pl
+        self.steerrefs = ids
+        self.nbsteers = nbsteers
+        self.bpmrefs = ring.get_uint32_index(bpmrefs)
+
+    def exclude_obs(self, *, obsid: int | str = 0, refpts: Refpts = None) -> None:
+        # noinspection PyUnresolvedReferences
+        r"""Add an observable item to the set of excluded values.
+
+        After excluding observation points, the matrix must be inverted again using
+        :py:meth:`solve`.
+
+        Args:
+            obsid:      If 0 (default), act on Monitors. Otherwise,
+              it must be 1 or "sum(x_kicks)"  or "sum(y_kicks)"
+            refpts:    location of Monitors to exclude
+
+        Raises:
+            ValueError: No observable with the given name.
+            IndexError: Observableindex out of range.
+
+        Example:
+            >>> resp = OrbitResponseMatrix(ring, "h")
+            >>> resp.exclude_obs("BPM_02")
+
+            Create an horizontal :py:class:`OrbitResponseMatrix` from
+            :py:class:`.Corrector` elements to :py:class:`.Monitor` elements,
+            and exclude all monitors with name "BPM_02"
+        """
+        super().exclude_obs(obsid=obsid, refpts=refpts)
+
+    def exclude_vars(self, *varid: int | str, refpts: Refpts = None) -> None:
+        # noinspection PyUnresolvedReferences
+        """Add correctors to the set of excluded variables.
+
+        Args:
+            *varid:  :py:class:`Variable <.VariableBase>` names or variable indices
+              in the variable list
+            refpts:  location of correctors to exclude
+
+        After excluding correctors, the matrix must be inverted again using
+        :py:meth:`solve`.
+
+        Examples:
+            >>> resp.exclude_vars(0, "x0097", -1)
+
+            Exclude the 1st variable, the variable named "x0097" and the last variable.
+
+            >>> resp.exclude_vars(refpts="SD1E")
+
+            Exclude all variables associated with the element named "SD1E".
+        """
+        plcode = plane_(self.plane, key="code")
+        names = [f"{plcode}{ik:04}" for ik in self.ring.get_uint32_index(refpts)]
+        super().exclude_vars(*varid, *names)
+
+    def normalise(
+        self, cav_ampl: float | None = 2.0, stsum_ampl: float | None = 2.0
+    ) -> None:
+        """Normalise the response matrix
+
+        Adjust the RF cavity delta and/or the weight of steerer summation so that the
+        weighted response matrix is normalised.
+
+        Args:
+            cav_ampl: Desired ratio between the cavity response and the average of
+              steerer responses. If :py:obj:`None`, do not normalise.
+            stsum_ampl: Desired inverse ratio between the weight of the steerer
+              summation and the average of Monitor responses. If :py:obj:`None`,
+              do not normalise.
+
+        By default, the normalisation gives to the RF frequency and steerer summation
+        a factor 2 more efficiency than steerers and BPMs
+        """
+        resp = self.weighted_response
+        normvar = np.linalg.norm(resp, axis=0)
+        normobs = np.linalg.norm(resp, axis=1)
+        if len(self.variables) > self.nbsteers and cav_ampl is not None:
+            self.cavdelta *= np.mean(normvar[:-1]) / normvar[-1] * cav_ampl
+        if len(self.observables) > 1 and stsum_ampl is not None:
+            self.stsumweight = (
+                self.stsumweight * normobs[-1] / np.mean(normobs[:-1]) / stsum_ampl
+            )
+
+    def build_analytical(self, **kwargs) -> FloatArray:
+        """Build analytically the response matrix.
+
+        Keyword Args:
+            dp (float):     Momentum deviation. Defaults to :py:obj:`None`
+            dct (float):    Path lengthening. Defaults to :py:obj:`None`
+            df (float):     Deviation from the nominal RF frequency.
+              Defaults to :py:obj:`None`
+
+        Returns:
+            response:       Response matrix
+
+        References:
+            .. [#Franchi] A. Franchi, S.M. Liuzzo, Z. Marti, *"Analytic formulas for
+               the rapid evaluation of the orbit response matrix and chromatic functions
+               from lattice parameters in circular accelerators"*,
+               arXiv:1711.06589 [physics.acc-ph]
+        """
+
+        def tauwj(muj, muw):
+            tau = muj - muw
+            if tau < 0.0:
+                tau += 2.0 * pi_tune
+            return tau - pi_tune
+
+        ring = self.ring
+        pl = self.plane
+        _, ringdata, elemdata = ring.linopt6(All, **kwargs)
+        pi_tune = math.pi * ringdata.tune[pl]
+        dataw = elemdata[self.steerrefs]
+        dataj = elemdata[self.bpmrefs]
+        dispj = dataj.dispersion[:, 2 * pl]
+        dispw = dataw.dispersion[:, 2 * pl]
+        lw = np.array([elem.Length for elem in ring.select(self.steerrefs)])
+        taufunc = np.frompyfunc(tauwj, 2, 1)
+
+        sqbetaw = np.sqrt(dataw.beta[:, pl])
+        ts = lw / sqbetaw / 2.0
+        tc = sqbetaw - dataw.alpha[:, pl] * ts
+        twj = np.astype(taufunc.outer(dataj.mu[:, pl], dataw.mu[:, pl]), np.float64)
+        jcwj = tc * np.cos(twj) + ts * np.sin(twj)
+        coefj = np.sqrt(dataj.beta[:, pl]) / (2.0 * np.sin(pi_tune))
+        resp = coefj[:, np.newaxis] * jcwj
+        if ring.is_6d:
+            alpha_c = ring.disable_6d(copy=True).get_mcf()
+            resp += np.outer(dispj, dispw) / (alpha_c * ring.circumference)
+            if len(self.variables) > self.nbsteers:
+                rfrsp = -dispj / (alpha_c * ring.rf_frequency)
+                resp = np.concatenate((resp, rfrsp[:, np.newaxis]), axis=1)
+        if len(self.observables) > 1:
+            sumst = np.ones(resp.shape[1], np.float64)
+            if len(self.variables) > self.nbsteers:
+                sumst[-1] = 0.0
+            resp = np.concatenate((resp, sumst[np.newaxis]), axis=0)
+        self.response = resp
+        return resp
+
+    @property
+    def bpmweight(self) -> FloatArray:
+        """Weight of position readings."""
+        return self.observables[0].weight
+
+    @bpmweight.setter
+    def bpmweight(self, value: npt.ArrayLike):
+        self.observables[0].weight = value
+
+    @property
+    def stsumweight(self) -> FloatArray:
+        """Weight of steerer summation."""
+        return self.observables[1].weight
+
+    @stsumweight.setter
+    def stsumweight(self, value: float):
+        self.observables[1].weight = value
+
+    @property
+    def steerdelta(self) -> FloatArray:
+        """Step and weight of steerers."""
+        return self.variables[: self.nbsteers].deltas
+
+    @steerdelta.setter
+    def steerdelta(self, value: npt.ArrayLike):
+        self.variables[: self.nbsteers].deltas = value
+
+    @property
+    def cavdelta(self) -> FloatArray:
+        """Step and weight of RF frequency deviation."""
+        return self.variables[self.nbsteers].delta
+
+    @cavdelta.setter
+    def cavdelta(self, value: float):
+        self.variables[self.nbsteers].delta = value
+
+
+class TrajectoryResponseMatrix(ResponseMatrix):
+    """Trajectory response matrix.
+
+    A :py:class:`TrajectoryResponseMatrix` applies to a single plane, horizontal or
+    vertical. A combined response matrix is obtained by adding horizontal and vertical
+    matrices. However, the resulting matrix has the :py:class:`ResponseMatrix`
+    class, which implies that the :py:class:`OrbitResponseMatrix` specific methods are
+    not available.
+
+    Variables are a set of steerers. Steerer variables are named ``xnnnn`` or
+    ``ynnnn`` where *nnnn* is the index in the lattice.
+
+    Observables are the trajectory position at selected points, named ``trajectory[x]``
+    for the horizontal plane or ``trajectory[y]`` for the vertical plane.
+
+    The variable elements must have the *KickAngle* attribute used for correction.
+    It's available for all magnets, though not present by default
+    except in :py:class:`.Corrector` magnets. For other magnets, the attribute
+    should be explicitly created.
+
+    By default, the observables are all the :py:class:`.Monitor` elements, and the
+    variables are all the elements having a *KickAngle* attribute.
+
+    """
+
+    bpmrefs: Uint32Refpts
+    steerrefs: Uint32Refpts
+    _default_twiss_in: ClassVar[dict] = {"beta": np.ones(2), "alpha": np.zeros(2)}
+
+    def __init__(
+        self,
+        ring: Lattice,
+        plane: AxisDef,
+        bpmrefs: Refpts = Monitor,
+        steerrefs: Refpts = _orbit_correctors,
+        *,
+        bpmweight: float = 1.0,
+        bpmtarget: float = 0.0,
+        steerdelta: float = 0.0001,
+    ):
+        """
+        Args:
+            ring:       Design lattice, used to compute the response
+            plane:      One out of {0, 'x', 'h', 'H'} for horizontal orbit, or
+              one of {1, 'y', 'v', 'V'} for vertical orbit
+            bpmrefs:    Location of closed orbit observation points.
+              See ":ref:`Selecting elements in a lattice <refpts>`".
+              Default: all :py:class:`.Monitor` elements.
+            steerrefs:  Location of orbit steerers. Their *KickAngle* attribute
+              is used and must be present in the selected elements.
+              Default: All Elements having a *KickAngle* attribute.
+            bpmweight:  Weight on position readings. Must be broadcastable to the
+              number of BPMs
+            bpmtarget:  Target position
+            steerdelta: Step on steerers for matrix computation [rad]. This is
+              also the steerer weight. Must be broadcastable to the number of steerers.
+        """
+
+        def steerer(ik, delta):
+            name = f"{plcode}{ik:04}"
+            return RefptsVariable(ik, "KickAngle", index=pl, name=name, delta=delta)
+
+        pl = plane_(plane, key="index")
+        plcode = plane_(plane, key="code")
+        ids = ring.get_uint32_index(steerrefs)
+        nbsteers = len(ids)
+        deltas = np.broadcast_to(steerdelta, nbsteers)
+        # Observables
+        bpms = TrajectoryObservable(
+            bpmrefs, axis=2 * pl, target=bpmtarget, weight=bpmweight
+        )
+        observables = ObservableList([bpms])
+        # Variables
+        variables = VariableList(steerer(ik, delta) for ik, delta in zip(ids, deltas))
+
+        super().__init__(ring, variables, observables)
+        self.plane = pl
+        self.steerrefs = ids
+        self.nbsteers = nbsteers
+        self.bpmrefs = ring.get_uint32_index(bpmrefs)
+
+    def build_analytical(self, **kwargs) -> FloatArray:
+        """Build analytically the response matrix.
+
+        Keyword Args:
+            dp (float):     Momentum deviation. Defaults to :py:obj:`None`
+            dct (float):    Path lengthening. Defaults to :py:obj:`None`
+            df (float):     Deviation from the nominal RF frequency.
+              Defaults to :py:obj:`None`
+
+        Returns:
+            response:       Response matrix
+        """
+        ring = self.ring
+        pl = self.plane
+        twiss_in = self._eval_args.get("twiss_in", self._default_twiss_in)
+        _, _, elemdata = ring.linopt6(All, twiss_in=twiss_in, **kwargs)
+        dataj = elemdata[self.bpmrefs]
+        dataw = elemdata[self.steerrefs]
+        lw = np.array([elem.Length for elem in ring.select(self.steerrefs)])
+
+        sqbetaw = np.sqrt(dataw.beta[:, pl])
+        ts = lw / sqbetaw / 2.0
+        tc = sqbetaw - dataw.alpha[:, pl] * ts
+        twj = dataj.mu[:, pl].reshape(-1, 1) - dataw.mu[:, pl]
+        jswj = tc * np.sin(twj) - ts * np.cos(twj)
+        coefj = np.sqrt(dataj.beta[:, pl])
+        resp = coefj[:, np.newaxis] * jswj
+        resp[twj < 0.0] = 0.0
+        self.response = resp
+        return resp
+
+    def exclude_obs(self, *, obsid: int | str = 0, refpts: Refpts = None) -> None:
+        # noinspection PyUnresolvedReferences
+        r"""Add a monitor to the set of excluded values.
+
+        After excluding observation points, the matrix must be inverted again using
+        :py:meth:`solve`.
+
+        Args:
+            refpts:    location of Monitors to exclude
+
+        Raises:
+            ValueError: No observable with the given name.
+            IndexError: Observableindex out of range.
+
+        Example:
+            >>> resp = TrajectoryResponseMatrix(ring, "v")
+            >>> resp.exclude_obs("BPM_02")
+
+            Create a vertical :py:class:`TrajectoryResponseMatrix` from
+            :py:class:`.Corrector` elements to :py:class:`.Monitor` elements,
+            and exclude all monitors with name "BPM_02"
+        """
+        super().exclude_obs(obsid=0, refpts=refpts)
+
+    def exclude_vars(self, *varid: int | str, refpts: Refpts = None) -> None:
+        # noinspection PyUnresolvedReferences
+        """Add correctors to the set of excluded variables.
+
+        Args:
+            *varid:  :py:class:`Variable <.VariableBase>` names or variable indices
+              in the variable list
+            refpts:  location of correctors to exclude
+
+        After excluding correctors, the matrix must be inverted again using
+        :py:meth:`solve`.
+
+        Examples:
+            >>> resp.exclude_vars(0, "x0103", -1)
+
+            Exclude the 1st variable, the variable named "x0103" and the last variable.
+
+            >>> resp.exclude_vars(refpts="SD1E")
+
+            Exclude all variables associated with the element named "SD1E".
+        """
+        plcode = plane_(self.plane, key="code")
+        names = [f"{plcode}{ik:04}" for ik in self.ring.get_uint32_index(refpts)]
+        super().exclude_vars(*varid, *names)
+
+    @property
+    def bpmweight(self) -> FloatArray:
+        """Weight of position readings."""
+        return self.observables[0].weight
+
+    @bpmweight.setter
+    def bpmweight(self, value: npt.ArrayLike):
+        self.observables[0].weight = value
+
+    @property
+    def steerdelta(self) -> np.ndarray:
+        """Step and weight on steerers."""
+        return self.variables.deltas
+
+    @steerdelta.setter
+    def steerdelta(self, value):
+        self.variables.deltas = value
diff --git a/pyat/at/plot/__init__.py b/pyat/at/plot/__init__.py
index d810d35cb..3975fbd87 100644
--- a/pyat/at/plot/__init__.py
+++ b/pyat/at/plot/__init__.py
@@ -12,3 +12,4 @@
     from .specific import *
     from .standalone import *
     from .resonances import *
+    from .response_matrix import *
diff --git a/pyat/at/plot/response_matrix.py b/pyat/at/plot/response_matrix.py
new file mode 100644
index 000000000..8533a3ee7
--- /dev/null
+++ b/pyat/at/plot/response_matrix.py
@@ -0,0 +1,113 @@
+from __future__ import annotations
+from ..lattice import Lattice
+from ..latticetools import ResponseMatrix
+from typing import Optional
+import matplotlib.pyplot as plt
+from matplotlib.axes import Axes
+
+
+def plot_norm(resp: ResponseMatrix, ax: Optional[tuple[Axes, Axes]] = None) -> None:
+    r"""Plot the norm of the lines and columns of the weighted response matrix
+
+    For a stable solution, the norms should have the same order of magnitude.
+    If not, the weights of observables and variables should be adjusted.
+
+    Args:
+        resp:           Response matrix object
+        ax:             tuple of :py:class:`~.matplotlib.axes.Axes`. If given,
+          plots will be drawn in these axes.
+    """
+    obs, var = resp.check_norm()
+    if ax is None:
+        fig, (ax1, ax2) = plt.subplots(nrows=2, layout="constrained")
+    else:
+        ax1, ax2 = ax[:2]
+    ax1.bar(range(len(obs)), obs)
+    ax1.set_title("Norm of weighted observables")
+    ax1.set_xlabel("Observable #")
+    ax2.bar(range(len(var)), var)
+    ax2.set_title("Norm of weighted variables")
+    ax2.set_xlabel("Variable #")
+
+
+def plot_singular_values(
+    resp: ResponseMatrix, ax: Axes = None, logscale: bool = True
+) -> None:
+    r"""Plot the singular values of a response matrix
+
+    Args:
+        resp:           Response matrix object
+        logscale:       If :py:obj:`True`, use log scale
+        ax:             If given, plots will be drawn in these axes.
+    """
+    if resp.singular_values is None:
+        resp.solve()
+    singvals = resp.singular_values
+    if ax is None:
+        fig, ax = plt.subplots()
+    ax.bar(range(len(singvals)), singvals)
+    if logscale:
+        ax.set_yscale("log")
+    ax.set_title("Singular values")
+
+
+def plot_obs_analysis(
+    resp: ResponseMatrix, lattice: Lattice, ax: Axes = None, logscale: bool = True
+) -> None:
+    """Plot the decomposition of an error vector on the basis of singular
+    vectors
+
+    Args:
+        resp:           Response matrix object
+        lattice:        Lattice description. The response matrix observables
+          will be evaluated for this :py:class:`.Lattice` and the deviation
+          from   target will be decomposed on the basis of singular vectors,
+        logscale:       If :py:obj:`True`, use log scale
+        ax:             If given, plots will be drawn in these axes.
+    """
+    if resp.singular_values is None:
+        resp.solve()
+    obs = resp.observables
+    # noinspection PyProtectedMember
+    obs.evaluate(lattice, **resp._eval_args)
+    corr = resp._uh @ obs.flat_deviations
+    if ax is None:
+        fig, ax = plt.subplots()
+    ax.bar(range(len(corr)), corr)
+    if logscale:
+        ax.set_yscale("log")
+    ax.set_title("SVD decomposition")
+    ax.set_xlabel("Singular vector #")
+
+
+def plot_var_analysis(
+    resp: ResponseMatrix, lattice: Lattice, ax: Axes = None, logscale: bool = False
+) -> None:
+    """Plot the decomposition of a correction vector on the basis of singular
+    vectors
+
+    Args:
+        resp:           Response matrix object
+        lattice:        Lattice description. The variables will be evaluated
+          for this :py:class:`.Lattice` and will be decomposed on the basis
+          of singular vectors,
+        logscale:       If :py:obj:`True`, use log scale
+        ax:             If given, plots will be drawn in these axes.
+    """
+    if resp.singular_values is None:
+        resp.solve()
+    var = resp.variables
+    if ax is None:
+        fig, ax = plt.subplots()
+    corr = (resp._v * resp.singular_values).T @ var.get(lattice)
+    ax.bar(range(len(corr)), corr)
+    if logscale:
+        ax.set_yscale("log")
+    ax.set_title("SVD decomposition")
+    ax.set_xlabel("Singular vector #")
+
+
+ResponseMatrix.plot_norm = plot_norm
+ResponseMatrix.plot_singular_values = plot_singular_values
+ResponseMatrix.plot_obs_analysis = plot_obs_analysis
+ResponseMatrix.plot_var_analysis = plot_var_analysis