{ "cells": [ { "cell_type": "markdown", "id": "7057dc7c", "metadata": {}, "source": [ "# 23 · Fast Complementary Dynamics — secondary motion that never fights the rig\n", "\n", "An animator poses a character with a **rig**: a few handles driving a smooth Linear Blend\n", "Skinning (LBS) deformation. *Complementary dynamics* (Zhang et al. 2020) layers physical\n", "secondary motion — jiggle, sag, follow-through — **on top** of that rig, under one\n", "rule: the simulation may never reproduce a motion the rig can already make itself.\n", "\n", "Split every vertex's displacement into a rig part and a complementary part,\n", "\n", "$$\n", "u \\;=\\; \\underbrace{u_r}_{\\text{rig (the animator)}} \\;+\\; \\underbrace{u_c}_{\\text{complementary (the physics)}},\n", "$$\n", "\n", "and constrain the complementary part to be **orthogonal to the rig**:\n", "\n", "$$\n", "\\boxed{\\,D\\,M\\,J^{\\top}\\,u_c \\;=\\; 0\\,}.\n", "$$\n", "\n", "* $J$ — the LBS Jacobian of the rig. Its columns span every displacement the handles can\n", " produce, i.e. the **control space**.\n", "* $M$ — the mass matrix, so \"orthogonal\" is measured in the physical *momentum* inner product:\n", " $u_c$ carries no momentum that the rig already accounts for.\n", "* $D$ — a diagonal **momentum-leaking** field. Where $D=1$ momentum is conserved against the\n", " rig; where $D\\to 0$ momentum is allowed to *leak*, letting that tissue ride along with the\n", " handle instead of resisting it.\n", "\n", "We build this twice. First a **full-order** solver, to watch the constraint do its job. Then the\n", "**fast** version of *Fast Complementary Dynamics via Skinning Eigenmodes* (Benchekroun et al.,\n", "SIGGRAPH 2023), which makes it real-time with two `simkit` ingredients: **rig-orthogonal skinning\n", "eigenmodes** and **rotation clustering**." ] }, { "cell_type": "code", "execution_count": 1, "id": "3fa923db", "metadata": { "execution": { "iopub.execute_input": "2026-06-09T05:43:18.352668Z", "iopub.status.busy": "2026-06-09T05:43:18.352480Z", "iopub.status.idle": "2026-06-09T05:43:18.899384Z", "shell.execute_reply": "2026-06-09T05:43:18.899126Z" } }, "outputs": [], "source": [ "%matplotlib inline\n", "import time, warnings\n", "import numpy as np\n", "import scipy as sp\n", "import scipy.sparse as sps\n", "import igl\n", "import matplotlib.pyplot as plt\n", "from matplotlib.animation import FuncAnimation\n", "import utils\n", "from simkit import (massmatrix, deformation_jacobian, lbs_jacobian, volume,\n", " ympr_to_lame, polar_svd, orthonormalize,\n", " project_into_subspace, cluster_grouping_matrices)\n", "from simkit.diffuse_field import diffuse_field\n", "from simkit.normalize_and_center import normalize_and_center\n", "from simkit.lbs_weight_space_constraint import lbs_weight_space_constraint\n", "from simkit.skinning_eigenmodes import skinning_eigenmodes\n", "from simkit.spectral_cubature import spectral_cubature\n", "from simkit.solvers import block_coord\n", "warnings.filterwarnings(\"ignore\"); np.seterr(all=\"ignore\")\n", "\n", "YM, RHO, DT = 1e5, 1e3, 1e-2 # Young's modulus, density, timestep\n", "RIG_C, FULL_C, RED_C = \"#d62728\", \"#1f77b4\", \"#2ca02c\" # rig / full / reduced colors" ] }, { "cell_type": "markdown", "id": "fc292850", "metadata": {}, "source": [ "## Setup: the beam, the shared FEM operators, and a side-by-side renderer\n", "\n", "Not the focus of the tutorial. These are the standard ARAP/FEM operators of tutorials 1–3\n", "(mass matrix, deformation Jacobian, ARAP stiffness), the metronome **beam** mesh, the rig-pose\n", "solve, and a small helper that animates several deforming meshes side by side. Run it, then skip\n", "to the next cell." ] }, { "cell_type": "code", "execution_count": 2, "id": "ad5c7640", "metadata": { "execution": { "iopub.execute_input": "2026-06-09T05:43:18.900696Z", "iopub.status.busy": "2026-06-09T05:43:18.900599Z", "iopub.status.idle": "2026-06-09T05:43:18.903933Z", "shell.execute_reply": "2026-06-09T05:43:18.903754Z" } }, "outputs": [], "source": [ "def beam(width=40, height=5, thickness=0.1):\n", " \"\"\"A thin vertical beam (a metronome arm), centered and normalized.\"\"\"\n", " X, T = igl.triangulated_grid(width, height)\n", " X[:, 1] *= thickness\n", " Y = X.copy(); X[:, 0] = -Y[:, 1]; X[:, 1] = Y[:, 0]\n", " return normalize_and_center(X), T\n", "\n", "def elastic_operators(X, T):\n", " \"\"\"Standard ARAP/FEM operators (covered in tutorials 1-3).\"\"\"\n", " dim = X.shape[1]\n", " Mv = sps.kron(massmatrix(X, T, rho=RHO), sps.identity(dim)).tocsc() # mass\n", " G = deformation_jacobian(X, T) # x -> per-elem F\n", " mu, _ = ympr_to_lame(YM, 0); vol = volume(X, T)\n", " AMu = sps.diags((mu * vol).flatten()) # elastic weights\n", " AMue = sps.kron(AMu, sps.identity(dim * dim))\n", " L = (G.T @ AMue @ G).tocsc() # ARAP stiffness\n", " return dict(Mv=Mv, G=G, AMue=AMue, L=L, AMu=AMu, dim=dim)\n", "\n", "def rig_rest(J, Mv, X):\n", " \"\"\"Rig parameters p whose pose J p best matches the rest shape.\"\"\"\n", " JMJ, JMy = J.T @ Mv @ J, J.T @ Mv @ X.reshape(-1, 1)\n", " return project_into_subspace(X.reshape(-1, 1), J, M=Mv, BMB=JMJ, BMy=JMy)\n", "\n", "# left/right sweep of the handle's x-translation (the only thing the animator does)\n", "PERIOD, AMP = 80, 3.0\n", "def handle(p0, i):\n", " p = p0.copy(); p[-2] = AMP * np.sin(2 * np.pi * i / PERIOD); return p\n", "\n", "def thin(states, target=110):\n", " k = max(1, len(states) // target)\n", " return states[::k]\n", "\n", "def animate_panels(panels, *, fps=30, suptitle=None, lims=None):\n", " \"\"\"Deforming meshes side by side, one color each. panels: list of dicts\n", " {states, T, title, color}. Returns (fig, anim).\"\"\"\n", " if lims is None:\n", " lims = utils.auto_limits([s for p in panels for s in p[\"states\"]], pad=0.6)\n", " npan = len(panels)\n", " fig, axes = plt.subplots(1, npan, figsize=(5.2 * npan, 4.4), squeeze=False)\n", " axes = axes[0]; arts = []\n", " nframes = max(len(p[\"states\"]) for p in panels)\n", " for ax, p in zip(axes, panels):\n", " utils.setup_axes(ax, lims[0], lims[1], title=p[\"title\"])\n", " arts.append(utils.PolyMeshArtist(ax, p[\"states\"][0], p[\"T\"],\n", " facecolor=p[\"color\"], edgecolor=\"0.3\", lw=0.5, alpha=0.95))\n", " def update(i):\n", " for art, p in zip(arts, panels):\n", " art.update(p[\"states\"][min(i, len(p[\"states\"]) - 1)])\n", " return ()\n", " fig.tight_layout(rect=[0, 0, 1, 0.92])\n", " if suptitle: fig.suptitle(suptitle, fontsize=13)\n", " return fig, FuncAnimation(fig, update, frames=nframes, interval=1000 / fps, blit=False)" ] }, { "cell_type": "markdown", "id": "bd6cbb30", "metadata": {}, "source": [ "## 1 · The rig: one global handle\n", "\n", "The simplest possible rig is a **single affine handle** — one LBS bone with weight $1$\n", "everywhere. `lbs_jacobian(X, ones)` builds its Jacobian $J$, whose $6$ columns (in 2D) span\n", "every affine map: translation, rotation, scale, shear. The animator's pose is a vector of rig\n", "parameters $p$, and the rigged shape is simply $x = J\\,p$.\n", "\n", "Here the animation is a metronome: sweep the handle's $x$-translation left and right. On its own\n", "the beam stays perfectly rigid — there is no physics yet." ] }, { "cell_type": "code", "execution_count": 3, "id": "4d1cbf05", "metadata": { "execution": { "iopub.execute_input": "2026-06-09T05:43:18.904806Z", "iopub.status.busy": "2026-06-09T05:43:18.904757Z", "iopub.status.idle": "2026-06-09T05:43:21.910392Z", "shell.execute_reply": "2026-06-09T05:43:21.910161Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "mesh: 200 verts, 312 triangles | rig: J is (400, 6) -> 6 DOFs\n" ] }, { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "X, T = beam(); n, dim = X.shape\n", "op = elastic_operators(X, T)\n", "J = lbs_jacobian(X, np.ones((n, 1))) # control space: one global affine handle\n", "p0 = rig_rest(J, op[\"Mv\"], X)\n", "print(f\"mesh: {n} verts, {T.shape[0]} triangles | rig: J is {J.shape} -> {J.shape[1]} DOFs\")\n", "\n", "NSTEPS = 200\n", "rig_states = [(J @ handle(p0, i)).reshape(n, dim) for i in range(NSTEPS)]\n", "\n", "fig, anim = animate_panels(\n", " [dict(states=thin(rig_states), T=T, title=\"rig pose $x = J\\\\,p$\", color=RIG_C)],\n", " suptitle=\"The animator's input: a rigid beam swept left and right\")\n", "utils.show_video(fig, anim, \"media/23_rig.mp4\", fps=30)" ] }, { "cell_type": "markdown", "id": "59763a34", "metadata": {}, "source": [ "## 2 · Full-order complementary dynamics\n", "\n", "Now add physics. Each timestep is a **backward-Euler** implicit step that minimizes the usual\n", "inertia-plus-elasticity objective over the *total* position $x = J\\,p_{\\text{next}} + u_c$,\n", "\n", "$$\n", "\\min_{u_c}\\;\\; \\tfrac{1}{2h^2}\\,\\lVert x - y\\rVert_M^2 \\;+\\; E_{\\text{ARAP}}(x)\n", " \\qquad\\text{subject to}\\qquad J^{\\top} M D\\, u_c = 0,\n", "$$\n", "\n", "where $y = x_{\\text{curr}} + h\\,\\dot x$ is the inertial prediction. The only thing that makes this\n", "*complementary* dynamics rather than ordinary elastodynamics is that one linear constraint.\n", "\n", "We solve the ARAP minimization with **block coordinate descent** (tutorials 21 & 22):\n", "\n", "* **local step** — for each triangle, the closest rotation to its current deformation gradient\n", " (one `polar_svd` per element);\n", "* **global step** — a single linear solve for $u_c$. With the constraint, that solve becomes a\n", " **KKT system** $\\begin{bmatrix} H & A^{\\top}\\\\ A & 0\\end{bmatrix}$ with $A = J^{\\top} M D$. Both\n", " $H = L + M/h^2$ and $A$ are constant, so we **prefactor once** and reuse the factorization every\n", " step and every inner iteration.\n", "\n", "The momentum-leaking field $D$ is built by diffusing a value of $1$ inward from the boundary\n", "(`diffuse_field`) and using $1-\\text{(that)}$: near the handle/boundary momentum leaks (tissue\n", "follows the rig), deep inside it is conserved (tissue resists, producing secondary motion).\n", "\n", "`FullCD` does all of this precompute in `__init__` and exposes a `step(...)` method." ] }, { "cell_type": "code", "execution_count": 4, "id": "8104ca3b", "metadata": { "execution": { "iopub.execute_input": "2026-06-09T05:43:21.911569Z", "iopub.status.busy": "2026-06-09T05:43:21.911487Z", "iopub.status.idle": "2026-06-09T05:43:21.914262Z", "shell.execute_reply": "2026-06-09T05:43:21.914071Z" } }, "outputs": [], "source": [ "class FullCD:\n", " \"\"\"Full-order complementary dynamics: per-element rotations + a prefactored\n", " KKT global step enforcing the rig-orthogonality constraint J^T M D u_c = 0.\"\"\"\n", "\n", " def __init__(self, X, T):\n", " n, dim = X.shape\n", " op = elastic_operators(X, T)\n", " self.Mv, self.G, self.AMue, self.L = op[\"Mv\"], op[\"G\"], op[\"AMue\"], op[\"L\"]\n", " self.dim = dim\n", "\n", " # rig J: one global affine handle, and the rest pose that matches X\n", " self.J = lbs_jacobian(X, np.ones((n, 1)))\n", " self.p0 = rig_rest(self.J, self.Mv, X)\n", "\n", " # momentum-leaking field D: ~0 near the boundary (leaks), ->1 in the interior (conserved)\n", " bnd = np.unique(igl.boundary_facets(T)[0])\n", " leak = diffuse_field(X, T, bnd, np.ones((bnd.shape[0], 1)), dt=1, normalize=True)\n", " D = sps.kron(sps.diags((1 - leak).flatten()), sps.identity(dim))\n", "\n", " # rig-orthogonality operator A = J^T M D, and the constant backward-Euler matrix H\n", " self.A = sps.csc_matrix(self.J.T @ self.Mv @ D)\n", " self.H = (self.L + self.Mv / DT**2).tocsc()\n", "\n", " # prefactor the KKT system [[H, A^T], [A, 0]] once, reuse every step / inner iter\n", " KKT = sps.bmat([[self.H, self.A.T],\n", " [self.A, sps.csc_matrix((self.A.shape[0],) * 2)]]).tocsc()\n", " self.kkt = sps.linalg.factorized(KKT)\n", " self.nc = self.A.shape[0]\n", "\n", " self.state_dim = n * dim\n", " self.n_rot = T.shape[0]\n", " self.B = None # full order: u_c lives in full coordinates\n", "\n", " def step(self, uc, uc_dot, p, p_dot, p_next, iters=10):\n", " Jp = self.J @ p_next\n", " y = (uc + DT * uc_dot) + self.J @ (p + DT * p_dot) # inertial prediction\n", "\n", " def local(u): # closest rotation, per triangle\n", " return polar_svd((self.G @ (Jp + u)).reshape(-1, self.dim, self.dim))[0].reshape(-1, 1)\n", "\n", " def glob(u, R): # constrained (rig-orthogonal) solve\n", " E_inertia = self.Mv @ y / DT**2 # backward-Euler inertia term\n", " E_elastic = self.G.T @ (self.AMue @ R) # ARAP (local-rotation) elastic term\n", " E_rigpose = -self.H @ Jp # pin the rig pose into the rhs\n", " rhs = E_inertia + E_elastic + E_rigpose\n", " return self.kkt(np.vstack([rhs, np.zeros((self.nc, 1))]))[:uc.shape[0]]\n", "\n", " return block_coord(uc.copy(), glob, local, 0.0, iters)" ] }, { "cell_type": "markdown", "id": "8580a6a9", "metadata": {}, "source": [ "### Rolling out a solver under the metronome rig\n", "\n", "The roll-out is the same for both solvers: at each step, advance the handle pose, take one\n", "implicit step for the complementary state, finite-difference the velocities, and record both the\n", "deformed mesh and the rig-orthogonality residual $\\lVert J^{\\top} M D\\,u_c\\rVert$. For the reduced\n", "solver the state lives in subspace coordinates $z$, so we lift $u_c = B\\,z$ before measuring." ] }, { "cell_type": "code", "execution_count": 5, "id": "b5eeed11", "metadata": { "execution": { "iopub.execute_input": "2026-06-09T05:43:21.915223Z", "iopub.status.busy": "2026-06-09T05:43:21.915157Z", "iopub.status.idle": "2026-06-09T05:43:22.941040Z", "shell.execute_reply": "2026-06-09T05:43:22.940820Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "full-order: 5.1 ms/step, 312 rotations/step\n", "max rig-orthogonality residual ||J^T M D u_c|| = 1.1e-13 (constraint holds)\n" ] } ], "source": [ "def run(solver, nsteps=NSTEPS):\n", " \"\"\"Roll out a solver under the left/right handle sweep. Returns (states, violation).\"\"\"\n", " J, p0, B = solver.J, solver.p0, solver.B\n", " s = np.zeros((solver.state_dim, 1)); s_dot = np.zeros_like(s)\n", " p = p0.copy(); p_dot = np.zeros_like(p)\n", " states, viol = [], []\n", " for i in range(nsteps):\n", " p_next = handle(p0, i)\n", " s_new = solver.step(s, s_dot, p, p_dot, p_next)\n", " s_dot = (s_new - s) / DT; p_dot = (p_next - p) / DT\n", " s, p = s_new, p_next.copy()\n", " uc = (B @ s) if B is not None else s # u_c in full coordinates\n", " states.append((J @ p + uc).reshape(-1, dim))\n", " viol.append(float(np.linalg.norm(solver.A @ uc)))\n", " return states, viol\n", "\n", "full = FullCD(X, T)\n", "t0 = time.perf_counter()\n", "full_states, full_viol = run(full)\n", "t_full = (time.perf_counter() - t0) / NSTEPS * 1000\n", "print(f\"full-order: {t_full:.1f} ms/step, {full.n_rot} rotations/step\")\n", "print(f\"max rig-orthogonality residual ||J^T M D u_c|| = {max(full_viol):.1e} (constraint holds)\")" ] }, { "cell_type": "markdown", "id": "e16c92ab", "metadata": {}, "source": [ "### Rig vs. rig + complementary\n", "\n", "Left: the rig alone (rigid). Right: the same rig **plus** the complementary physics. The beam now\n", "whips and trails as the handle reverses — secondary motion the animator never authored — yet\n", "the constraint guarantees that motion adds nothing the handle could have done itself." ] }, { "cell_type": "code", "execution_count": 6, "id": "061d48f7", "metadata": { "execution": { "iopub.execute_input": "2026-06-09T05:43:22.942113Z", "iopub.status.busy": "2026-06-09T05:43:22.942057Z", "iopub.status.idle": "2026-06-09T05:43:28.036427Z", "shell.execute_reply": "2026-06-09T05:43:28.036147Z" } }, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "fig, anim = animate_panels([\n", " dict(states=thin(rig_states), T=T, title=\"rig only $J\\\\,p$\", color=RIG_C),\n", " dict(states=thin(full_states), T=T, title=\"rig + complementary $J\\\\,p+u_c$\", color=FULL_C),\n", "], suptitle=\"Full-order complementary dynamics\")\n", "utils.show_video(fig, anim, \"media/23_full.mp4\", fps=30)" ] }, { "cell_type": "markdown", "id": "e8b5c1c5", "metadata": {}, "source": [ "## 3 · Why the full solve is slow\n", "\n", "Two costs scale with the mesh:\n", "\n", "1. **The local step** computes one `polar_svd` **per triangle** — thousands of tiny rotations\n", " every inner iteration.\n", "2. **The global step** solves a sparse KKT system whose size is the **full vertex count** plus the\n", " constraint rows.\n", "\n", "Refine the mesh and both grow. *Fast* complementary dynamics removes both dependencies with two\n", "ideas, each one line of `simkit`." ] }, { "cell_type": "markdown", "id": "e05fc85a", "metadata": {}, "source": [ "## 4 · Fast complementary dynamics\n", "\n", "**Idea 1 — rig-orthogonal skinning eigenmodes.** Instead of solving for the full field $u_c$, restrict\n", "it to a small subspace $u_c = B\\,z$. *Skinning eigenmodes* are the natural rotation-aware modes of\n", "the shape (tutorial 17). The trick: pass the rig-orthogonality requirement as an equality constraint\n", "when building them (`lbs_weight_space_constraint` → `skinning_eigenmodes(..., Aeq=...)`), so that\n", "**every column of $B$ already satisfies $J^{\\top} M D\\,B = 0$.** The constraint is now *structural* —\n", "true for any $z$ — so the KKT block vanishes and the global solve shrinks to a tiny dense system.\n", "\n", "**Idea 2 — rotation clustering.** The local step doesn't need a distinct rotation per triangle.\n", "`spectral_cubature` groups elements that move together into a handful of **clusters**; the local step\n", "then fits **one rotation per cluster**. `cluster_grouping_matrices` assembles the (mass-weighted)\n", "grouping operator used in the local/global steps.\n", "\n", "Together: the global solve is independent of mesh resolution (it is the size of the subspace), and\n", "the local step is independent of mesh resolution (it is the number of clusters).\n", "\n", "`ReducedCD` does both precomputes in `__init__` (subspace $B$, clusters, the tiny Cholesky factor)\n", "and exposes the same `step(...)` interface as `FullCD`." ] }, { "cell_type": "code", "execution_count": 7, "id": "a8e23617", "metadata": { "execution": { "iopub.execute_input": "2026-06-09T05:43:28.037941Z", "iopub.status.busy": "2026-06-09T05:43:28.037871Z", "iopub.status.idle": "2026-06-09T05:43:28.041420Z", "shell.execute_reply": "2026-06-09T05:43:28.041193Z" } }, "outputs": [], "source": [ "class ReducedCD:\n", " \"\"\"Fast complementary dynamics: a rig-orthogonal skinning-eigenmode subspace\n", " + clustered rotations. Same step() interface as FullCD, but tiny and dense.\"\"\"\n", "\n", " def __init__(self, X, T, k_modes=10, n_clusters=20):\n", " n, dim = X.shape\n", " op = elastic_operators(X, T)\n", " Mv, G, AMu, L = op[\"Mv\"], op[\"G\"], op[\"AMu\"], op[\"L\"]\n", " self.dim = dim\n", "\n", " # rig J and its rest pose\n", " self.J = lbs_jacobian(X, np.ones((n, 1)))\n", " self.p0 = rig_rest(self.J, Mv, X)\n", "\n", " # momentum-leaking field D (same construction as the full-order solver)\n", " bnd = np.unique(igl.boundary_facets(T)[0])\n", " leak = diffuse_field(X, T, bnd, np.ones((bnd.shape[0], 1)), dt=1, normalize=True)\n", " D = sps.kron(sps.diags((1 - leak).flatten()), sps.identity(dim))\n", "\n", " # Idea 1: subspace B with rig-orthogonality (J^T M D) B = 0 baked in\n", " Aeq = lbs_weight_space_constraint(X, T, (D @ Mv @ self.J).T)\n", " W, _E, B = skinning_eigenmodes(X, T, k_modes, Aeq=Aeq)\n", " self.B = orthonormalize(B, M=Mv)\n", "\n", " # Idea 2: cluster elements that rotate together; one rotation per cluster\n", " _cI, _cW, labels = spectral_cubature(X, T, W, n_clusters, return_labels=True)\n", " self.labels = np.asarray(labels).flatten()\n", " P, _Pm = cluster_grouping_matrices(self.labels, X, T)\n", " PAMue = sps.kron(P @ AMu, sps.identity(dim * dim))\n", " self.GB = (PAMue @ G) @ self.B # clustered deformation operators\n", " self.GJ = (PAMue @ G) @ self.J\n", "\n", " # reduced stiffness / mass and rig-coupling blocks\n", " self.BLB, self.BMB = self.B.T @ L @ self.B, self.B.T @ Mv @ self.B\n", " self.BMJ, self.BLJ = self.B.T @ Mv @ self.J, self.B.T @ L @ self.J\n", " self.chol = sp.linalg.cho_factor(self.BLB + self.BMB / DT**2) # tiny dense factor, reused\n", "\n", " # rig-orthogonality operator (used only to report the residual)\n", " self.A = sps.csc_matrix(self.J.T @ Mv @ D)\n", " self.state_dim = self.B.shape[1]\n", " self.n_rot = P.shape[0]\n", "\n", " def step(self, z, z_dot, p, p_dot, p_next, iters=10):\n", " y = z + DT * z_dot\n", " q = p_next - (p + DT * p_dot)\n", "\n", " def local(zz): # one rotation per cluster\n", " return polar_svd((self.GB @ zz + self.GJ @ p_next).reshape(-1, self.dim, self.dim))[0].reshape(-1, 1)\n", "\n", " def glob(zz, R): # dense solve in subspace coordinates\n", " E_inertia = self.BMB @ y / DT**2 # reduced inertia term\n", " E_rigvel = -self.BMJ @ q / DT**2 # rig momentum coupling\n", " E_elastic = self.GB.T @ R # clustered ARAP elastic term\n", " E_rigpose = -self.BLJ @ p_next # rig pose coupling\n", " rhs = E_inertia + E_rigvel + E_elastic + E_rigpose\n", " return sp.linalg.cho_solve(self.chol, rhs)\n", "\n", " return block_coord(z.copy(), glob, local, 0.0, iters)" ] }, { "cell_type": "code", "execution_count": 8, "id": "d4d881d6", "metadata": { "execution": { "iopub.execute_input": "2026-06-09T05:43:28.042343Z", "iopub.status.busy": "2026-06-09T05:43:28.042283Z", "iopub.status.idle": "2026-06-09T05:43:28.194654Z", "shell.execute_reply": "2026-06-09T05:43:28.194445Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "subspace B: (400, 60) -> 60 DOFs\n", "structural rig-orthogonality ||(J^T M D) B|| = 2.4e-08\n", "rotation clusters: 20 (vs 312 triangles)\n", "reduced: 0.68 ms/step -> 7.5x faster than full-order\n" ] } ], "source": [ "red = ReducedCD(X, T)\n", "print(f\"subspace B: {red.B.shape} -> {red.state_dim} DOFs\")\n", "print(f\"structural rig-orthogonality ||(J^T M D) B|| = {np.linalg.norm((red.A @ red.B)):.1e}\")\n", "print(f\"rotation clusters: {red.n_rot} (vs {T.shape[0]} triangles)\")\n", "\n", "t0 = time.perf_counter()\n", "red_states, _ = run(red)\n", "t_red = (time.perf_counter() - t0) / NSTEPS * 1000\n", "print(f\"reduced: {t_red:.2f} ms/step -> {t_full / t_red:.1f}x faster than full-order\")" ] }, { "cell_type": "markdown", "id": "ff11b5ef", "metadata": {}, "source": [ "### The rotation clusters\n", "\n", "Each color is one cluster of triangles that share a single rotation in the local step. A few dozen\n", "clusters capture the bending of the whole beam, so the local step cost no longer depends on how\n", "finely the beam is meshed." ] }, { "cell_type": "code", "execution_count": 9, "id": "c258b173", "metadata": { "execution": { "iopub.execute_input": "2026-06-09T05:43:28.195713Z", "iopub.status.busy": "2026-06-09T05:43:28.195656Z", "iopub.status.idle": "2026-06-09T05:43:28.225435Z", "shell.execute_reply": "2026-06-09T05:43:28.225211Z" } }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEiCAYAAAD+lrlwAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjksIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvJkbTWQAAAAlwSFlzAAAPYQAAD2EBqD+naQAAT7hJREFUeJztnQeYXGXZ/u/pM9t7b9leso1NDyQBAkFRQZEmShFQUSyAfyTKB4L65VMsKOKHHVD5BCyoSA9ICYHU7T2bbLb3ma3Tzpz/9byzc3ZmsykbdrPJyfPjOmTnnjMzp7znPu95nuecVyPLsgyGYRjmrEG71AvAMAzDnFrY+BmGYc4y2PgZhmHOMtj4GYZhzjLY+BmGYc4y2PgZhmHOMtj4GYZhzjLY+BmGYc4y2PgZhmHOMlRl/DfeeCMyMjJwuqLRaPDtb38bZyNf/OIXcdFFF+FM5vHHHxf7cM+ePTgTUEN7823zQ4cO4Wzg29/+tljfwcHBY87ncrmQmpqKX/ziF2eG8dOC0s48Wbq7u8XGqaysXNDlOtM5nbfLwYMH8Zvf/Abf/OY3l3pRmEXghRde+MAnmP/+7//Gc889t2DLpHYMBgPuvPNOfO9734Pdbj87jP+BBx6Y0+B+/etfo6mpCWcjx9ouS81Pf/pTLFu2DOeff/5SLwqzSMZPbW8xjP8zn/kMpqamkJ6e/oG+X43cdNNN4srgqaeeOvXGPzExgdPpLGgymZZ6MVTFB92/dEn6pz/9CVddddWCLZPamZycXNLfpx6kx+PB6YBOp4PZbBbhDyaQiIgIXHzxxSfVkdaeTPypvr4en/rUpxAZGYlzzz1XvOd2u/Gd73wHWVlZwnwp1k6X9g6HQ/k8aXV1dXjzzTfF99C0adMm8d7w8DC+/vWvo7i4GCEhIQgLC8OHPvQhVFVVKZ//z3/+g5UrVypnO993+FZ8rhg/Gdddd90l4mG0XHl5efjhD3+I2Q8lpe+5/fbbRa9j+fLlYt6ioiK89NJLJ3yw0PbJzc0VDTUxMRGf+MQncODAgXnnJHzb2Z9XX31VbGva2bR9aD18oZPjbRfi/fffxyWXXILw8HAEBQVh48aN2LFjxwnv397eXvHdKSkpYtvQ+l122WXHjb2+8847oleyefPmI97r7+/HzTffjPj4eLHNSktL8cQTTwTMQ99Py0T77Fe/+pXSvmh9d+/efcR3NjY24pOf/CSioqLEd65YsQL//Oc/cSL8+c9/RkVFBUJDQ0X7o7ZIVyuzoTZNl9mxsbEIDg7Gxz/+cQwMDATM849//AOXXnopkpKSxPLSctPxIUlSwHzU/qm97d27Fxs2bBD7xrdf6Xfuv/9+ZGdni++gNnz33XcHHFO++e644w6xPLTsH/vYx9DZ2XlC60xth7Yvrfu9996L5ORksQyjo6Pi/WeffVZsE4vFgpiYGHz6059GV1dXQBt+9NFHxd++duffdmm/rVu3DtHR0eI76Lv+8pe/BCwDzU/HKe173+fpe48V46fIAR2ftF1oG3/pS1+C1Wqdc9tSe6arTVovWr8f/OAHOFH++Mc/KutPbeqaa65BR0dHwDxvv/02rrzySqSlpSn7ifYHXanM1T6pE0T7ir6TjuNvfetbR8xH60LbgI53Ombp2JurQ0B5MzrGyD/ngx4nAa1kTk6OuDzzGegtt9widhwddGS0ZDTbtm1DQ0MD/v73v4t5Hn74YXz5y18WxuVbWTroiba2NmG69N0UFujr68Mvf/lLYVC042jnFhQU4MEHH8R9992Hz33uczjvvPPEZ6lhzQUtGx0Eb7zxhjCYsrIyvPzyy/h//+//icb7k5/8JGB+2oB/+9vfRCKSDqCf/exnuOKKK3D48GHRcI8GHcwf+chHsH37dtEwvvrVr2JsbEyYdW1trTjoPwh0sqTvLykpEetPjau1tVUx7uNtl9dff12cRKkBk5FotVr8/ve/xwUXXCAa7apVq467f2k70HLQ/qOTFZk2rR9tm2Ml1N99911x4JaXlwfodFDQgUnrQSdc2udkMtTYqdHTNvSHLmdpm37+858X30cHL51Yqd3QlZ5vO61fv14c3Pfcc48w5WeeeQaXX345/vrXvwqDPhq0Ltdeey0uvPBCfP/73xcatV3axrOXhbYBnRRpW5IhUbumdXj66aeVeciwqJ3TCYL+pX1A+4cM9aGHHgr4vqGhIbF/qO2QsdIxQT1uarvUJmmf0j6uqakRbba5uTkgLELHHhkUnaxpn9Nv0UlnPtBJyWg0is4XnUjob1oHMhw6ydKxTMcknQhpm+zfv1+YEu0PCjPS9vvDH/5wxPfS/LQe1113HZxOpzjBUPt6/vnnlWWkz9E6UDukdSWOdcxQB4VCS9SZuO2220R493//939FR4CWzdceiJGREdHhobZChksnnW984xvipE7b/FhQ/Py//uu/xOdo+ejk/sgjj4gTtG/9CWq3ZMq0LOQTu3btEvPRyZfe81FdXS2OTVo+Wk86bqhj+K9//Uv8lj/0m3RM0Hbft2+fyJHFxcUpbdMHHdN0jNJxRh5xwsjz4P777ycXkK+99toAvbKyUui33HJLgP71r39d6K+//rqiFRUVyRs3bjziu+12uyxJUoB28OBB2WQyyQ8++KCi7d69W3zn73//+yO+44YbbpDT09OV188995yY97vf/W7AfJ/85CdljUYjt7a2KhrNZzQaA7SqqiqhP/LII8fcLr/73e/EfD/+8Y+PeM/j8QT8Bm3Doy3v7O3s4yc/+Yl4PTAwcNRlONp2od/PycmRt2zZErAsk5OT8rJly+SLLrrouPt3ZGRE6A899JA8Xz796U/L0dHRR+gPP/yw+M4//vGPiuZ0OuW1a9fKISEh8ujoqNIGaD76juHhYWXef/zjH0L/17/+pWgXXnihXFxcLNqS//qvW7dObINj8dWvflUOCwuT3W73UeehbUu/uXnz5oBteccdd8g6nU62Wq0B23c2n//85+WgoKCA5aNjgb7zscceC5j3D3/4g6zVauW33347QKf5aP4dO3YEHHtf/OIXA+b71Kc+dUR7m4s33nhDzJeZmRmwzLQv4uLi5OXLl8tTU1OK/vzzz4v577vvPkX70pe+FNBe/Zm9Heh76TsvuOCCAD04OFgcD0fb5tQOiP7+fnGcXnzxxQF+8fOf/1zMR8fi7G375JNPKprD4ZATEhLkK6644pjb5dChQ2Kffu973wvQa2pqZL1eH6DPta+3bdsmPKa9vV3RNmzYIIeGhgZohH9b8h2Dn/3sZwPm+fjHPz7ncdTd3S3m//73vy/Ph5OK8X/hC184IrlDUO/GH+r5E//+97+P+53Ui6WeqK8HTb0gX0iDzngnAy0XxQi/8pWvHLFc5MMvvvhigE49CP+eBvWw6ZKfepXHgnqTdBlMPcHZLERs0tezoPDBfGOvlOxtaWkRvUHaphR2oYkural3+9Zbbx3xnbP3L12SUg+QwgLUg5oP9JvUO55r3yQkJIhetg/qCdG+Gh8fF+FAf66++uqA7/Fd1fj2DV3qUk+Xekp0ZeBbT/r9LVu2iG3gH6KYaxvTNqGe6/Gg3pr/fqVloTbb3t4esM18+JaH5qOeIV3uz2771LP2h3qK1MvPz89X1oUmukoj6CrWtx2J2W38a1/7GubDDTfcELDMVLJKV3V09UshMx/US6dlOpFjmvD/Tmo7NptNbIeTPaZfe+01ceVA6+fzC+LWW28Vx+rs5SIPoasoH9SOV61addxjmq786big9uS//anN0tWwb/vPXkdqQzQfXXmRx9CVAUFXC3SsffaznxUhoeN5xOxjkLYZtWVfCM6H75g4XvnngoR66BLEH2rwtBMoFukPbSQ6oPwPiKNBG5kuCyl2R+V//rHQY4VZjgX9LoWIKGzjDx1Qvvf9mb1DfBv2eGZHl2t0gtLrT2pzHhcyPbrUo8tNCmGQYdOlK4XV/Bv/XJDh+Q7so0EHo7+pzt6/ZEx0iUknTApDrFmzRlxWXn/99WIfH4+5BnmjbU8H0OzlP9F941te376hkBH9Dl2a0zQXZGQUBpoLMjgKC9HlP81DSTM66ClMMJvjLYsv7EQxczoZzT5YaXv7Q79HhjR7v1GoiWLBR1sX/2NvdmiE2uMHPaaP9j1k/BSCOhEopPPd735XdED8cxMn2yE62nLR9svMzDyi3VBOavZvRUZGirDLsaDtT+2J2uhc+IeTKNxJYTzKJc32Ct++9p1oKOdwIhyrjdEJbvaxNd/teVJO5X+GW6jeLcWT6YClMyLFGymRQg2azuynqsKArg7mYrFGpzza9pqdAKTtTb0F6mVQj4YSzhRPpt7fK6+8ctTlJnzbjuLKlOOYC+oVzf692dB++OhHPypiy5QnoX1F8Ucyttnxe3/opD3fq4ST2Te+9aQYNfXw52J2x8Qfip+SOdG60ZUgTZQHoZPb7ITz8ZaFchSUm6IDlHIvZMrUa6ZeLsWXZ7fnubY3zUNx6B//+Mdz/hYlEBeSox3THwTKH1F8n2Li1KGjggAyTNquJ1OCeCqPaY/HI45PagdzfYfvmKFjlRKsdMVJ+5ZOipRboqtLyledrHed6HL7ji2KOMyHBemiUo0trSCdJX09NoKSQXQQ+NfgHs3sKOlCmfff/va3ATp93n+l5nNyod+lS0O61Pbv9fsutReqNpgObEpmU+mif0/geNBZfHYlAjHXFRKdBKmnTxOZAZ0oKUFOJwMKUR1tu/h6gmRCc1XWzAf6Lur100T7mk4kP/rRj0Ri8WjQgUDlnNTzoeoEH7TtqddF7ca/13+y+4Z6ewRt/5NdT+o10smNJlouugqgAgM6yR3rpDEbConRZTmFC8j0fNCV7Hy2NVW00f4+Vpv3HXu+q04fH/R+Ft/2p+/xhZf8v/tEjmkKgdIJj06m/mXWZPyzOdHj2n+5fPucoPAPbd8P2sb9tz+ZLF0JUaXe0aCEOyXbqXNAnQQfs0OGvmWlYo+FxNem/H33lN3A9eEPf1j8S9UN/vh6K/4VBnQ2nMvs6Aw3+2xGcc7ZcVn6PDHXd8y1XHRG/vnPfx6gU2UENbTjZfVPFKp4oRjb7N85Xs+CGhcZov9lZ09Pj1IF5WOuUi1f7913+Xy07UJZf/odKquj2PlsZpchzgXFpWffHUjfSSfT2aWFs1m7dq3YBlSuOHvfUImofyUMlQRTNQT1pqjHPB+ox05VQmTUtA3nu55k1P7QyYhyPMTx1vFovTX/fU/GNJ/b6ynMRG2fbkqcDVVE+e6v8LVhqkDzZ/axOF+oDJa26WOPPRaw/tQDphDU7GN6rrZH24GOM/8rWKqCmutGraP5wmzI2OkETevrv32pw0jH0nyrmY4GhVJp+al6aPYxTK997WWufU1/zy4DppAddQJ+97vfidDQ7O87Wei4om1Mx9kp7/FT/TXFkKnO2neZSyVNdBakUjr/OzbJiKj0iuJ+1IuixkU9CooZ02UxJbkoMUJnUuop+p/VfYZDeQNqkGQ81GBWr159RIySoJ4b/Tb1jKnB0XJSaISSpBS6+KBllj7oTP/kk0+K5DatNyVi6MCkqw3qNVK9+1xQ+R5dHlKZISXnyGBp21APwz/5RduFQj3UqKnHQ/FdMhGKX/rq7I+1XSg/QAZBdc+0fSmmTKZCVwt0JUDlZMeCejTU8yQzKiwsFLkMOjnRFR2tw7Gg5aNwD20L/54jJUjJpOlymBovlbbRVR+V45Fpzc7LnAhUT06/RyESSvZR26Fl3Llzpyit878nZDaUP6ETLC0jbVe66qKTEJ1g59ubovZLV3N0TNB+pQOTShbnc4DTHauUc6AkH+0nKlMlA6UrItKpF03mTMtHCXJqD2R89NtUVkw5jw8CXTlRXofaCx3P9Bu+ck7aV1Sn7n9ME7SuFGYjM6R2Qe2VOn+UJ6HiAmq3tI/ouJ8dY6fvoDZC81Nejtottd/ZkIFu3bpVGDJ9L4WSqPdP609lp/6J3A9CVlaW8Cj6LfIO8jFqk9TDprZP7ZfCinRFS/PS33RM0fFEVzpzhTfpZEXt85xzzhGfp3Wk76bw7cnecU9XFtQ25p0HnU8JkK/UaK6yQpfLJT/wwAOiRNBgMMipqany1q1bA0rXiN7eXvnSSy8VZU30Xb7STprvrrvukhMTE2WLxSKvX79e3rlzp3h/dvknlfIVFhaKsir/Esa5yiPHxsZEuV1SUpJYLirro7JE/xIqgr6HytJmQ983V5nZbKik61vf+pay/lQyRmWjBw4cCPiN2eV1r7zyiihvoxK1vLw8Ud44u5xz+/bt8mWXXSbWgeajf6nksrm5+YS2C7F//375E5/4hCgJoxJZWq+rrrpKfPfx9u/g4KDYNvn5+aLsLjw8XF69erX8zDPPyCfCV77yFTk7O/sIva+vT77pppvkmJgYsV5Uijm7HNVXzjlXKelc25O29/XXXy+2P+2H5ORk+SMf+Yj8l7/85ZjLSO9TiSCVMNKypKWlifLLnp6eI0oLqXR2rpJI+tcHlVuuWbNGtGXaX3fffbf88ssvHzEftW0qcZ4LKn2kMj16n/ZZZGSkXFFRIY4zm82mzEfllrSNad/S/vnoRz8qd3R0zKuc89lnn53z/aefflouLy8Xvx8VFSVfd911cmdnZ8A8VAL75S9/WY6NjRUljP5t97e//a045ujz1H5oG85u30RjY6Mod6TtRe/5jrnZ5Zz+5Zv0fbSP4+Pj5dtuu02UHftztG17w1HKqOfir3/9q3zuueeK7UoT/SYdC01NTco89fX1osSXypCpLd96661KKfjs9lxbWytKMyMiImSz2SyO+f/6r/867jE413ag8mFqq7/5zW/k+aKh/53UqYZhThCqaKCeEYUJ6MqBYZgPDl0Z042MlN+Zb3KejZ85JdBdjRR+OJE6eYZhjg0VklCIicq7KZw8X9j4GYZhzjJUNRALwzAMc3zY+BmGYc4y2PgZhmHOMtj4GYZhzjIW56liDLME0KML6NnwdKPNQo7YRBXP9NgPurHoeA/FY5gzATZ+RjWQ6S/0w8v8oZGX6K5ehjnTYeNnVIPvMQ9k0P6Prv2g0GOV6YRyMo+RYJjTETZ+RjX4wjtk+gtp/LO/n2HOdDhgyTAMc5bBxs8wDHOWwcbPMAxzlsHGzywKNH4AjYdAJZAUG59r8I25Rq6iZ5XTaE30zPbHH3+c9w7DLAJs/MyiQAPR0MA3NPDGiUADXNDAHTRwDg1KQQPl0OAoNOAIwzALC1f1MIsCjfg1n6EtaeQwGpGIxvAlaNSrd955RwyTebTB0xmGOTm4x8+cFtDwiLMHyibDJ51hmIWFe/zMaQENvB4fHx+g0Wu6eYoGF59rhCEaBNx/IHCal6Dxcv1vttqxfz9KsrMDtN319ciMjw8Yq/RQSwtMZjMS/e7+7e/upmc2LOCaMszSw8bPnLFs27ZNDLo9mxve3I3JBO+jFfRuF/Lra/CLrmEMJaQqWl7DPnxadiA6IlJoDo+Eyv4BJBmNSIyKUrT3BoaQtXLlKV0vhllsONTDnBYkJCSgr68vQKPXdAfu0cYT3bp1K2w2mzLRoxqIsbBIWKPjxL+pXQfRl1OE8fBoRUvvbENDQQXiTRZkBAchwWzCLtsYrk5KQrTRqGjvjY7hs6nJ0POD2RiVwT1+5rRg7dq1eOGFFwI0Gp+X9KNBZZ80zYXO5UJp9XuoKlmDjMFuQJKFVl69E/tL1iJ0dETMZ5ckPNPbh6sS4zEpSYr2594+XB0fx6bPqBLu8TOLwvj4uCjLpMlXrkl/Hz58WOmtX3/99cr8X/jCF9DW1oa7774bjY2N+MUvfoFnnnkGd9xxx7x/W+92K6YvGQxC00kunFPlNX2f5pIlPNPTh6sS4mHW6YTmlmXF9C167hcx6oRbNrMo7NmzR9Tk+7jzzjvFvzfccIO4Maunp0c5CRBUyvnvf/9bGP1Pf/pT8fjj3/zmNydVylnYsA/63EKc29kiXmsGeqEdGYSUVaBoknUYrw8MYXVqCt6yexPEtskpNAyPoCIlBW87nIDDCbfHA8fU1AfeHgxzOsHGzywKmzZtEgOYHI257sqlz+zfv/8D//ah1ExMRMYridzCsWZ0Ji3DoL92sAXFoSEonI7fUyL3eZsNOSFBKNLNaC8ODiEnY9kHXiaGOZ3gUA+jOkYjYpREblZHK3ryijEWGatoOYdbUFu0AkmWICWR++7oGK5JSkScyaxoO0bHcWNyEsf5GdXBxs+oEkrkrqjagT2l6+DR6RVtdeUO7CJNb1ASuf9HyV0R09cp2lN9/bgyLgbBHOdnVAgbP6M6KLm7snIH9pSsC0jukum/X7oOktGoJHefmjZ9n8FTcpdM/5Mx0Qid/izDqA2O8TOqo6BhHzpKVyPeOiBehw71I63zEFqLVypakG0Y/xkcwkeSk3B4OpFrczqxy2bFlvgEdLncYpJkD+x+dwczjBpg42dUR29sIiS3G3C7oZPciOo6hL6ouAAtq60J5sh0NIZkis+4PS4099UgJCwZLWE5itbWXY/iNO8VAsOoBTZ+RnUMJKRiLHWZiOmvrXwHnQVlsGkN6EvJENqW/TtwSdY62B0TWJaUC5fbiT3Nb2FlwQYMWvsUbV/L21hVcD4kg/cZQAyjFjjGz6gSMvh1le/gPXHDlrfHLky/cgc+lbEeZqNZaGTwe5vfwvKMVQgyhyjavtZ3kJ9ahrCgiCVdD4ZZDNj4GVUmd32m7zaZleQumf61aesQZA4Wmtvjxt6Wt1GYvhKhQWFCkzxu7D/wDvKSSxAeMvPkToZRE2z8jOooaq7E7sIVkHR6aEVc34XSllpclbwKJqMZbskNySOhtbMeecml4kQgNElCZ/9BLIsvQEhQhDIf6QyjJjjGz6iOnLZWrG7zDvkoabXoiopEQnQ6DnY3eTVZQt9IJ4x6Cw73HlC0fms3tFoduvsPi4m0wdFe5BVkL+n6MMxCw8bPqI71rS3ImpyCU6vFWytXYH17OzRllyA1PkuJ3+emFEOWPchOLRJaZdu7KMgow8TkGAqXlQutpv09FGdVwBDMPX5GXXCoh1ElwvRXrcSKqmqEOp1CIzP3xu9LER4SqWhVB3diWXw+osPiFK2m/X0kRWUgLjx5SdeDYRYDNn5Gdbi0Wry9ehUqqmsQYbd7NYmStjuQk1iMiNBoJblbfWgn0uNyEBOeoCR3aw/vQmJUOhIi05Z0PRhmsWDjZ1TH/rw8lFfXIHJyUryWPR60dNYgK6EQkWGxynwHe5uQHJ2J2PAkResePozYsGQk+pm+x+M5xWvAMIsLx/gZ1eGemsC7Gd4xd+nB0JPmIDhsw6ipelmZx6k3QoYWvYNdYiKs40OQXBIGR/rERIxOjmBZLvf8GXXBxs+oDp3ODd8wuRPRiaIc0zLSAZPb+8ydqdhkuK1TyC3/qEjuEo2d+xGti0VcRJJI7hJNXZUICkqCxTz3mL8Mc6bCoR5GtUyFRgH2KVgkl6LZw6MhO+ywTI0pWnt/M43ThWWJ+Yp2eKAFbklCVoL3xMAwaoKNn1ElTnMwnB4PgiesgZqsQdDokKIN2LphHRtGfkrZjDbagyFbP4rSKk75cjPMqYCNn1Edbq0BY3oTwv1M363TY9xoQdjooKJNOiZxqKcFpVlrFM3utKOtqxHl2etP+XIzzKmCY/yM6hgLDkdmQjw0Gm+1zqRtFOM6N9JiIqGJ95ZtOmw2HOppQlFeCdoGaoRmtVnR2X8IRXmliuZyuRAcwQOyMOqCjZ9RHUETNkwdtIm/PRoNJiwRsDgnMHVwQNHGLOEwjh7E4e4GRRsJioDeacfh7kZFs4VEYdV5G5dwbRhm4eFQD6M6TB4XzJILRsmNSUsYokGv3TNaUDiiJm0w+WvBkYgZG4bR46eFRCHW2gedzjsWL8OoBTZ+RpV4oMFwaDSixoaURi608FhEjQ5CKyr8p7WIOETa+gO0kcgERFr7FI1h1AQbP6M6KEQzHOYzfXlGI9O3DSiarJk2fWs/dLNMP2KkV9EYRm1wjJ9RHSOWcCTpDXBFexO54x5g3GJEnE6raFNuCTZoEa/VwB2dADclfCUJ4xIQq4GiyZAxZZ9a4jVimIWFjZ9RHcH2CciuCUjTPf3JoEgEOcYhTzgVbTwoEmYnzWdXtLGQKBhdk8Dk1Mx8obHInh66kWHUAod6GNVh9rgQJHkTuvaQSMTITlgkt6I5KPY/MRKgOcNiED06BLNnRnNFxCHK2gs9J3cZlcHGz6gSb9I2XsTv/ZO7s+P3pFmjkxA23Aud7FE0W2wKQoa6oZd5EBZGfbDxM6pDJHKF6ffNGDzV6UclBpg+JXeF6Q/1QI8Z0x+N85q+YfpEwDBqg42fUR3WoAgEjw5B0urg1Orh8MiwBkcgyDqgaC6tDragcJitA/BotV5Np8OkORT64X7I0EzPp4fL5R3Bi2HUAid3GdWRb8hFUniM+NslS6h1tSNLDkNyXIWi1bnbkeYJwbLYcxSt3tOBJJcZWTEz8zVpuqGDeQnXhmEWHjZ+RnUUmTORa0mDU3bhJWkf1lqKYXHpkWtM9WryPmx0l8Ehu1BqyhLaK5pKXOiqwKDeitXmAqG9qq/Gxe6VGLHwnbuMuuBQD6NKhMF79uF8bTEsWpOivYz9OA9FCNMGKdor2kqs8uQiRhemaK/pa1AuLUMivIOyM4yaYONnVIdLdgvT36RZjtBpg6ewDZn+erkAkdoQZb5XtFVY6clBnDZcaG7Zg+2GGpRI6UhC1JKuB8MsFmz8jOrYKTViI5YjTBssXkseCbulZqyV8xGlDVXmq5IPYYUnE/GaCEVrRTcK3SlIQbSi0dCNDKMmOMbPqI4wyYJKZ6v4W5I96PAMwCLpUKc5qGjd8jCMbi2a0CEm0no0w9A4PDiEHjGR1q+zId2Ru8RrxDALCxs/ozrKTFlKcpdCOTlSKmLcIUpy9zVDDfImU2HU65Tk7uumOhROLYNT41CSu/8JakCJMwvuIB5snVEXHOphVImoytFVYYUnGzG6cEXbbqhFsTsNydooRXvDXI9cKRGZmnhFezO4EamuaOS6vaN4MYyaYONnVAclbV/VVeMcKQsJ0/F7Su6+bqxFkZSCVPhq/N34j6UeWe44ZErxSnL3zZAmJLuiUOBOXtL1YJjFgo2fUR3voxllUgYSNd5STI/Hgz2aA8h1JSFNjlXmq9EeRqojGtlSoqK16noRbw9BoWvG9Dm5y6gNjvEzi8ajjz6Khx56CL29vSgtLcUjjzyCVatWzTnv448/jptuuilAM5lMsNvt8/5do0OLNl032tAtXnd5hqBxSeiS+tGFfqH1yVZ4nG4MSiN4GyNCG8Qo3A4XbNIY3sb0AOzaCcQ60ua9DAxzOsPGzywKTz/9NO6880489thjWL16NR5++GFs2bIFTU1NiIuLm/MzYWFh4n0fGo3mpH57lTlfJHeJHdpGhDgtyDVkieQu8a6+GSFTFuTokkRyl3jP2IIQu4QkTZRI7hLvm1sRIRkQxMldRmVwqIdZFH784x/j1ltvFb34wsJCcQIICgrC7373u6N+how+ISFBmeLjvXH3k6VW0w6zpEeuIXVG03ZA65ZRoslQtDpdB1ySC6uQo2j1hi5Muu1Y65zRGEYtsPEzC47T6cTevXuxefPmmYam1YrXO3fuPOrnxsfHkZ6ejtTUVFx22WWoq6s75u84HA6Mjo4GTD46MIAR9zgqNNl+2iAG3VasQV6A1isN4zypUNE6NUPocPdjk6vopNafYU532PiZBWdwcFAkRGf32Ok1xfvnIi8vT1wN/OMf/8Af//hHkZBdt24dOjs7j/o727ZtQ3h4uDLRCYOwyZNodHbgPN2McY/Kk6hzH8ImbbGijclTqJEO4UKUKtqk7MB+6QAu9pSddKiJYU53OMbPnBasXbtWTD7I9AsKCvDLX/4S3/nOd+b8zNatW0UewQf1+Mn897pbUB5TiF2aNqF3jfei3zmCsth8RRucGMYhVw8qImbmG50awwHXYZRGFmL39F2+bo8bWrv3eT8MoxbY+JkFJyYmBjqdDn19fQE6vabY/YlgMBhQXl6O1lbvoxfmgqp+aJpNbOHb8CS+Kv52u3Xoq8tBXGofPIkvK1pHfTaiE4bhSZ2Zr60hB6HLRiBnvCYGWyetqSkfJcEXzmv9GeZ0h0M9zIJjNBpRUVGB7du3KxqFbui1f6/+WFCoqKamBomJMzX2J0pE+ASiokcRGjaBjs405Ob1IzLST+vIQEFhC8LDJhXtcEcWysprERoy5aflYMXKvdDp+Xn8jLpg42cWBQrB/PrXv8YTTzyBhoYG3HbbbZiYmFBq9a+//noRqvHx4IMP4pVXXkFbWxv27duHT3/602hvb8ctt9xyUr/vculQXVWMkpJa6HSyolVVlaKkpBoGvWdGqyxHScl+6PWS33wVKC7eD6PRtQBbg2FOLzjUwywKV199NQYGBnDfffeJhG5ZWRleeuklJeF7+PBhUenjY2RkRJR/0ryRkZHiiuHdd98VpaDzhUI01VUlKCmthsHgNXO3W4uqqjKUllTCYPRqkkTaOSgpJYOX4HLpIXm0qK5ageXL98FkckGSOMHLqA+NLMve7hDDnOFQcpeqe77+9StQUjIAnd7btAcGItHdFYLlyzsVzWYz4WBbAoqWd8Bg8Pb+7VMU009CQUE3jCbvycHj0cBo+DSuvvpzsNls4iYzhjnT4VAPozpiYwfg8TjgcjoxNSmhpycEUdFDAVp7exwiIoYB2a5ora0pCA2zQqOZUrTmpkwu62RUBxs/ozri40eQkjqA+IRhDA6moCC/D0mJftpQOpYXNSMuzqZoQ8PLUFpWi+iocUUbGclGWfk+GE3GpV4lhllQ2PgZVeJN0FL8vhL66VCON+Hri9+7Fa2mejXy8/fBYnEoWm3NOmTn7EdIyMSSrgfDLAZs/IzqoOTuTFWOW0nuVletRGHRfpjNLiW5S6afl78PwcFe06fkbl3tOmRlVyI0lE2fUSds/IzqaGnJQX5+FbRaWfTeXU6N0HJyq0TJJml0cjh4aBmWZVaLkk2vpkdPdwqSU+pgsUx6P0uVPm7vyYNh1AKXczKqwxI0gsPt3juEPR4t+vriYLaMoLMjTtH6++NgMtrQ0x2taAMDCdBqxzHQHyom0oaHk5Cfz8bPqAs2fkZ1ZKT3ISV1QonfL8vsh9lkRVLyoBK/z8puhUcClmX2CK2+7jzk5jVhYlyH3LxOoTU1bkJ+fhOCgmeeMsowaoBDPYwq8Zr+GuTm7UNQkEvR6mrXIyurCqHTSVvSGuo3IDWtzlveOa01NZ2P+PhmRMd0Lel6MMxiwMbPqA6K39fUrEVOLlXl2JXkbl3demQsq0ZY+KiS3G1s2IiU1HpERQ14NY8WzU0XIC6uBbFxHUu6HgyzWLDxM6qjtTUf2dn7EBo6KV57PDLa2vKRnl6DiAirMl9HRw4SE8n0Z54i2t+XgahoqvFvVzSP5L2Ll2HUAsf4GdWhNwyj4/DMuL7DQ9HQaIbR3RUuJmLEGg2PbEVfnxl9fd5hGG22aHikCQwPaTA85NXGx2ORm8VPNWHUBRs/ozqys3pEcpdoaToHhRGpcMS8I5K7RGvLCpSGJ2EweLdI7hIHWldimS4FvdoqZOYdElrbgdXItyyDbA5dwrVhmIWHQz2Mauk8nI3lphTkxs48cqGrMwc5uhSUJc70ebo685CGVKxPnTkcurvyEC+l4vw0DvMw6oONn1ElQ4NxiJrKxrlp5gAtdCwbF2TMnAiGhuJhGc3GJZkzJ4Lh4QToRnLwsWw+PBh1wqEeRnVMTpmhHy7CVcUzj1CenLLAM5iDq4pnwjZ2uxna0Sx8tiRY0RxOI0bas/C5Uh5nl1EvbPyM6jjcloMNBdl4qdM7iEp77xAOW5fh3LwsResbHEFnfyo2FixTtBHbKHq7E7A+P0PRnC4XPGFTS7g2DLPwsPEzqmNDcjgKgofE3w63jLpxO9YmhQZo9TY7VsSHBGgNQw4URQcFaC8dciE7OnsJ14ZhFh4OYjKqIyVcj4woPRLCdNjTJWFDTghSIwyKtqtTwkcKzEgKn9He6/Dg6lIz4kJ0irazQ8ZnV5qh1/Fg64y6YONnVIndLeMvNQ5cUWyCQadRtGerHbiyxASjfkZ7usqJq0oMMBtmtP+rcuGqEj3M0/MxjJpg42dUh5NMv9qBTxabFOMm7ZlqpzB9n+Zyy/jztOlbDN5DwS1BmP7VJXoETWsMozY4xs+ojn83OrAmLww7u70Gf7B/Cp0jdqzJjlS0fqsL9d0OrM+NxHs9Xm10QkJVhwOrsqLw/rTmlmRMhHJyl1EXbPyM6jgnWYesUIeSoN016kBpYqC2e8iBgjgtssJmtL1tTmRGaZETPqO93CIhp8SwhGvDMAsPX8syqiM90oCsGAOSI/TY1SXj/LwQZEQZFe39LuBjhWak+WnvdWlwVakZiWEzn32vW4vrzzHCoOf+EaMu2PgZVWKfjt9fXWKAUTcraVusD0jkPlXlxpXLdQgx+WnVHnw8X4MwMx8ijPrgVs2oDkrk/t+spC1pT02bfrBRqyR3n6qScGWRDqGmmeQumf7l+UBkEJdxMuqEjZ9RHS82O7ExXYuRSQ+6bG4MjLnwSrMDG9I0sE5Na+MS3mxzYl2KjFGHV+sd82BvtwsrEiVMuWSh9YxKcLqcS71KDLOgcPCSUR3p4Rp02bxP1XRJMt5rdyE1DMLEfdrODjfiQ2RxAhgYn56vQ0KEyQPbpHcibXcPkFPkWeI1YpiFhY2fUR1lyUYUp5iV+P2FeWYE610oTTIq8ftLC3SYdMpYl2Hyxv5rgY8vd6PbpsN5Wd7PPl2nxVXFHjSbZp7wyTBqgEM9jCrxmr6ETxbpEGQMTNpelgdETCdtRRK4VoNLsmTEh85oT9frsDHdgxTvgF0MoyrY+BnVIRK51RKuKNQqVTki4Vsj42N5GkQH62bu3K3V4KJMDxJDZ5K7zzTosSHFg4wIflwDo07Y+BnV8XKrG58o0CDc4m3eHo8H29vc+HAOEBM80+TfOezBhRkeJIfNaJV9MtYmurEsaub7JA+PwsWoC47xM6ojyt6JN9/zJmQlaNDqisdURBD+NOgBBj2QPTKmul2wmz141uoBrF7N3uPGhMmDf41LwDiE5uzzICPPtdSrxDALChs/ozrW66pRbLDC7tHjz9IWlOrbULliA4JzQyA5PRh+YwqaEgkWPRBcEiY025tO6M6RYLY7EbwyUmhjb7thXqWBSeLRuBh1waEeRpV4Tf8SXKTZiXj9uNDIzEfesCOk2AhLgnfcXZ/pW3L0CE6b0cbfkWBK1cKcxX0jRn2w8TOqwynr8WfPh3Ch9j0k64eF5nF7YP2PA8FFBliSvA9dk90ejL7lEuZuyfAavOyWMfGOB8ZkLSy5bPqMOmHjZ1THq64V2IT3kaobFK89HhnW/XZY8nSwpMw8aXOs1gljugZBmTMGP9nmhjbeI+b14ZE4ucuoC+7SMKoj2DWI+qkY1CMGMj2PX4qD2dOHoD3dwB4IzeaKhMOjwWj/OEb7qasPOIfccLjtGLOOY2yHV3NbZYQlJC31KjHMgsLGz6iOTcGtKA61ir//LW1AuEOD86TdKB3rFtoL2ITIyRG0nluBmJIwoY3tcsEeNgVzkgHxKyOFNr5bgidKgiXIsoRrwzALD4d6GNXynlSMEMmKYnPXjCaXwuCexGp9g6JNNLjgliSEnzNj8JONbrimXAhZxYOwMOqDjZ9RJS1SMoZcJmw0VitaqycFfc5gXKTfpWj2Tjcm+xyIXDvzPB57lxsTnXZEnGc65cvNMKcCNn5m0Xj00UeRkZEBs9mM1atXY9euGcOdi2effRb5+fli/uLiYrzwwgsn9btDnmBU2lNxqXmPoo3IQdhtz8BlpncVzT0hYaRuHLGbQhRNmvRgpGocMRcGQ6PhRzYw6oSNn1kUnn76adx55524//77sW/fPpSWlmLLli3o76dM6pG8++67uPbaa3HzzTdj//79uPzyy8VUW1s7799+y16AsIhIvGL8kJh2OzLx9lQeosJDFe0duRzDNeMIjw7D1F5ZTKP7nLA1TyA8Lgz2fRDa5F4P7FM82DqjLjSyLFORA8MsKNTDX7lyJX7+858rz8tJTU3Fl7/8Zdxzzz1HzH/11VdjYmICzz//vKKtWbMGZWVleOyxx07oN0dHRxEeHo4XvlyMgpAxoTlkPf45UYoCfSeWm/sU7fmJEoyUZiGqINS7fC4P+nbYYIzRIXp5mKIN7BzDysJVuPKyq2Gz2RAW5n2PYc5kuMfPLDhOpxN79+7F5s2bZxqaVite79y5c87PkO4/P0FXCEebn3A4HMLs/ScixWBDhnEYCfpR7HLlYlNYB1L9tPdc+fhIcBUssSYEJ5lhjjFitNaOpIuiYIo0BGgpH4mGVs9DMDLqgo2fWXAGBwchSRLi4+MDdHrd29s752dIn8/8xLZt20QP3zfRFYX/IxuemVyDK4Peg0HrmXmMw+R6XGV5FyaNpDyeofdVK+IvDId+eixer2ZD3Pnh0JvZ9Bn1wcbPnLFs3bpVhF98U0dHh/LIhmcm1+KqoPdg1rqnNR3+PLUeV1t2wKJzK49sUEx/2uDpkQ19r5Dph8LAg60zKoVv4GIWnJiYGOh0OvT1eWPqPuh1QkLCnJ8hfT7zEyaTSUyzeWGyGJfE9KARywAP0OKIQdNkGC6OOIQmjVfrlSwY2DOKlHNj4eiU4IAEx7gT1voJJK6KhrPbAyec8EgyHFP2k94WDHM6wsbPLDhGoxEVFRXYvn27qMxRBkPZvh233377nJ9Zu3ateP9rX/uaor366qtCny+F+k4YXBNK77/ZHo4cXReMrlFF2zuRg8z2g8h8c5/Q3FotKpctQ+rICHJe36toNdnZMJ6beBJbgWFOXzjUwywKVMr561//Gk888QQaGhpw2223iaqdm266Sbx//fXXi1CNj69+9at46aWX8KMf/QiNjY349re/jT179hz1RHEsck0DKLN0I9/Ujzo5C5tDDyHfPKho1c58VOw7hOSxMWQPDSNtxIqOjAxsrG9AwsSEonVmZeGCqmqYDHz3LqMuuMfPLApUnjkwMID77rtPJGipLJOM3ZfAPXz4sKj08bFu3To89dRTuPfee/HNb34TOTk5eO6557B8+fKT+n1K5P6ffQOuMr+NFgrvuL3aH8cuQPlrjRi1mDFuNMKp1eLNFRVYs78S8vTykPb2qlVYUVWFELsd3qf+MIx64Dp+RjX46vj33JGJKn0prja9hWCdC5WeXLjsU6hyF6LslQYEu93oDQ1BX0gIBtLSsHp/JUKdToyaTGiKjcFoahrKq6sRNTEBt0aD9q334KM33MB1/Ixq4FAPozpenirBpYb3IGs0GPcYMenW4pWJEuS90Soa/JReD4dej7r0dBTX1UPv8ShaW1IS8hsaYHE4hGY3GERpKsOoCQ71MKojqHMcr3Xki789Wi0OxicgqG8UB6OjhCZptWiPjUWQzYbDEeEAwoV2OD4e+rExdIcEi4m0rqQk5Lh4sHVGXXCPn1Ed2V1DOKe1B8vb+jCUmI2o865F/uAgVnZ2obS7BxPJySjt6EDO8LCiTaaloaK9HZk2m6I5li3DykPtCDbPPLmTYdQAGz+jSihB+5/zzkPsJ78Gc1DITNJ2zWqU1dQiYvrBa6TtWL8ORQ0NiLVNl3tqtdi54TxkNrcgada9BQyjBtj4GdXhEqa/AdFXfBVBod7RtNw6Hd5eswaltXWIHh/3alot3j13PQoamxA/NKxo723YgMyWVqT29CzpejDMYsHGz6iOfUXFiPrE7QgOi1JuHqvOzkFJXS1iph/kRjRmZSGnsREJAwOK1p6ahrSWFqR2zYza5ZG8z/phGLXAyV1GdbizStDbViP+liHDNtABfUgkqjZ+SNHGhnuhsQ6hw2xGR1oqZGgwEBUFl92OXp0WvdPacEw0UnX8oDZGXbDxM6ojLa8CCel54u/2mv8gIjoRYfmrEZfmrfTpqHsLJmc+8v78GIr6vAPD7FlRAVx8KbKf/T+sPOx92Nu+lSugu/RyGEJ5sHVGXXCoh1EtvW2VMAVHIjI+TdH6DlZBZwxCTIr3xEA05eRg4KJLEV5UpmjNuTno27QZYaUrT/lyM8xiw8bPqBJrXzsc9kkkZJYGaFMTNiTlrFC07rg4tF54MSLXbVS0nvg4tGy8AJEbLjrly80wpwIO9TCqY2rcirHBMeRUzIzoZZ+0wdrfjtwVWxRtwmTGwbXrEXPpJ2bmM5mwd8P5iPnYVad8uRnmVMHGz6iO7ta9yCsqw9AB77CNvd0dGLUOIaegVNEGB/phzchAWtFy2CvfF5p1ZAQjcfHIKC5RNBpGMsjEhwmjLrhFM6ojNm05DOHp4m+324HJ8QbEpOQHaKPWRphWroWU4B2uUXI5YW1phqmkPEDr27UDyYXepDDDqAWO8TOqIywqHuGxyQiOiEF/ez3SC1YiPDpR0foO1SKr9HwExcQhJDEZluhYWBtrkHbBJbBERCraSH01ln34E9DquX/EqAs2fkaVuF0OtOx5BTkVF0GnNyha8+6XRZzfp0lOJzrfeAnJG7dAZzDOaK+/iJTzL4HO6NUYRk2w8TOqw+12oplMf8XF0BtMM9rul5C7couiedwudJDpb9oC/fTYvR63W5i+v8YwaoOvYRnV0Vb9JjLzyzFwsFK8Hu7rgHWoD+l5ZYo2OjKIvt4DSC0px3hDldCmxscw1NqM5KJijDdWC80jSbDwjbuMymDjZ1RHfHoxQqJnErnjLTWITSsI0KytNQjNyIIhKW0mudtUj+DkNBiTMxRtcNcOpHFyl1EZHOphVEd4dAIi4lIQEhmL/oM1WFa0GhExSX5aNbLLL0BIXAJCk1JEktfWUI208y9BUFS0oo021CDt4o9Cx4OtMyqDjZ9RJZTIbdr1InJWXhKQ3G3a9QJyVlwCvd7kl8h9AUkbLobebFG0rtdfRPz6C2AM9j7Ln2HUBBs/ozpEInfXC8hdeQkMRnOARqZvMFmU5C4ZfNJ5F8EQFDStudH9xouIX3c+TKGhS7oeDLNYsPEzquNQzdtIyDoHU2MjGB3qwZh1EO217yI+swz2cavQJkYHMVC1B2F5y+EcH8VEXw8mB/sx0lyPkOx8uKcmhDYx0AuX07nUq8QwCwondxnVUaPZiaqmHeJvjaxFUn8SRizD2Nf8tqIlDybDahhBdeVORUsZSEFETCZcExNiktwu2JpqEVtUuKTrwzALDff4GdUxEjWCvoQ+DMQOIMERj8HUAVijrYqW5ExEW/IBjEaOoi++DwMxA0hxJqMltRmhKWmIKzkH0fnL4RnqR+oFH4JpOvbPMGqBjZ9RJVpJi5KOEtQm1cKj9yhaWWcZahJq4DK4FK2iqwL1sfWYskwpyd2+N19BVOkKWKLjlnQ9GGYxYONnVIdW1qK0s1SYvtvo9moeLco7y1EXXweXadr0ZS1WdK9AQ2wDJi2TQpMlCf1vvYKIknNgiUtc0vVgmMWCjZ9RHZm9mahJrFFMnwZbX9a/DLVxtXCYHcp86UPpqI+ux7hlXNFG25oRVliKoPgkRZN5sHVGZXByl1EdNr0NiT3e3rpG1iByPBJWrRVx/XEz2kQkbLINkQORoP9Ii56MRniIE7q212BvA8jvB/olBGdlL/EaMczCwsbPqA5K4Nrj7N74fXcFhmOHMawfxkj0iNBW9a5Ce2w73LIbA/EDQts4eC4uLIuH3e5ERYUFTqcH777rxNpyIySJH9bGqAsO9TCqhMxcxO+jG+CwOBRtde9qNEc2Y8wypmgbh87D5StSkZrqNXiv6buQkqJDRgb3jRj1wcbPqA5K2q7sWYnGmEaMB40ryd01fWvQHNGM4aBhZb5NQxvwsfJUpKV5n7vvdsvYuZNMX4vsbO+jHhhGbbDxM6ojtz8XDZENGLWMKsndnKEcNIc2Yyh4SJkv25qFD5ckISNjxuAPHnQiPh4Bpi9xcpdRGXwdy6iOEYwgdDhUTETERARG3CMwWU1ItiZ7talwxBtN6O8fE5MsA0NDbvHI5pGRMezc6dWsVg+SkhKWeI0YZmFh42dUB92NS8ldYvnAckAG+oP7RXKXKB8qw/q4VMTF6lFcHCG0PXscCA+3IynJiIqKKKHt2+dETIwHFgvfucuoCw71MKolYyQD45pxWCOsipZly8SFiZkoLAhWtMZGBzweN8rKvE/oJJqaXLDb3TjnHB5zl1EfbPyMKokej4bOrsOhmEOKFjsRi01BuVi/MlzRurtdGBhwYNWq4ACtu3sK69bNnAgYRk1wqIdRHWanGVH2KLSktQRo63R5+MjFM49hmJiQ0Ns7hosu8oZ2iMlJCbW147joIm8IiGHUCBs/ozqSB5IRlhWGmLEY8do17IJhXIvClWnYv9/7wLbBQQc6OkZRVjaj2Wx2HDo4gZLSVFRWyvSwBrjdEgwG78PbGEYtsPEzquO8sEykeCLF326PB1Xj/YgPCYKphxRZaL3dE8jMC0ZioneQFZdLRlPDKJITLQHa++/bUFi4bClXh2EWHI7xM6ojKtiCmJAghFtMOGQdQ35iHKJnaSVJcYiONiIx0YLoaBNqKydRHpmEiAiDotXVTeJDH4qFXq9b6lVimAWFjZ9ZcIaHh3HdddchLCwMERERuPnmmzE+PvMEzLnYtGkTNBpNwPSFL3zhpJfBJUnY29WPiuQ46LVaRdvT1Y9zkmIVjR7P8NpLwygOjodeN6Nt3z6MTZsiYTbzIcKoDw71MAsOmX5PTw9effVVuFwu3HTTTfjc5z6Hp5566pifu/XWW/Hggw8qr4OmB0CfLxTK8Zm+Qeftrbs9sjD9iqRYGPXU7J3i8QyvvTSC5cHxQnNKTkgSsH37CDZujIDFooPb7Y3/M4yaYONnFpSGhga89NJL2L17N1asWCG0Rx55BB/+8Ifxwx/+EElJM8+5nw0ZfULCB79LtrZnAMUZqeid9D6crX90DN0jNixPSULfFMXvnbBOTOLwPhsKYpLQP61NORxo7RjDilUx6OqS6GEN4nEPdrv3ZjCGUQt8HcssKDt37hThHZ/pE5s3b4ZWq8X7779/zM/+6U9/QkxMDJYvX46tW7dictI7KtbRcDgcGB0dDZiIhJAQGDweaNxuMYxir20M8cHBMMgzWseIDRGaIBhlWdEODI4gDCY4D+gxWgcM13iw+y0bZDf3jxh1wS2aWVB6e3sRFxc4Tq1er0dUVJR472h86lOfQnp6urgiqK6uxje+8Q00NTXhb3/721E/s23bNjzwwANH6AnhwUiJDFNi+kVJCdDJ0ozWPYDSpDhMOV0zWs8gViYnoH9iQtH29g5hVWo83AY+TBh1wT1+5oS45557jki+zp4aGxtPemtSDmDLli0oLi4WOYInn3wSf//733HgwIGjfoauCmw2mzJ1dHQo7/kncg3TSVuf6ZcnRMM0HfsXBt8ziJLYSJiNekXb1zuMopgIhJj4kQ2M+uCuDHNC3HXXXbjxxhuPOU9mZqaI0ff39wfobrdbVPrMJ36/evVq8W9rayuysrLmnMdkMolpruTunq4BnJMYC5MvkUvJ3e4BlMVHw2wwYNzhgkRJ4J4hLI+NRLDJiHGHEx6PLEy/IDocoSajmIdh1AYbP3NCxMbGiul4rF27FlarFXv37kVFRYXQXn/9dZEk9Zn5iVBZWSn+TUycecTCiVLfN4TzstLhgYwplwt2pwtNfUNYm54EaCA0h1tC66AVq9OToNNqpjU3OmxjKE9OgEHn1SSPDDeV+jCMimDjZxaUgoICXHLJJaI087HHHhPlnLfffjuuueYapaKnq6sLF154oQjnrFq1SoRzqNSTKn+io6NFjP+OO+7Ahg0bUFJSMu9lqHOEoeaghoZggVb2IFEawZAnBFWHtDOaZxgj7iBUt+sULRmDGHSYUdmhV7Q07QjKguO5lTCqgo2fWXCoOofMnsydqnmuuOIK/OxnP1Pep5MBJW59VTtGoxGvvfYaHn74YUxMTCA1NVV85t577z2p3z+gT8eEMQk62Y01nga0GbIw4tZj0JgotPVoRLUrmwbcRY8xTWgbdM3Y78yDwTiODmOW0M43tGKfKx8rg2ae3MkwaoCNn1lwqILnWDdrZWRkQKbhraYho3/zzTcXdBnIuNfKDdiDLCTqxgG3JLRzNU3Y68mAVutCsDQmtI26FlRKyZjUGhCLcaFdYDyAOnc8ehAmnu/DMGqCq3oY1aGXJaxDI/YiE3ad9+5fvezBudom7PekYVwbosy3Sd+Kak8ihjXeZ/TrZA8uNLWhUYpFt+x90BvDqA02fkZ1FHgOil79lNYboqEEbb6nHfulVIxqvePwEjnoRrU7HoOYefZ+umYQDa4odHhmntHv4eQuozI41MOojmG3HgnOTuV1jGcEI24NYtAtJiJatsHqBELlPoSizzufPIpRpwSDPIhMDAotVjMBlzN3idaEYRYHNn5GdXSZMjBh8VYQleMA7O4odBvCRHKXqNAcxKgjEqM6M3pMaUJbqTuEUacRE1oZHWbvfQOrDYdhl0Jg4sHWGZXBoR5GtWTL3RhwWzBmmInV56AHVrcePdqZexJyNb2YkHQ4pJm5wSxf1wen24Nqaf73ETDM6Q4bP6NK4jzDkF12dOpnngYaLw9D557EAV2qoiVgBCb3OOrhr1kR7B7FfjnjlC83w5wKONTDqA6LPIUo1ygaTfl+mh3hTitqTQV+mgOhrglUGWY0k+xEvLsLe/RFp3y5GeZUwcbPqI4sqQNxafnI1XiTtmMDvXC7bYhIzUXmtGYbGoTGPoDQpHyk++azWuEZ70dYch4+Nq05HU647VzWyagLNn5GdWTZUxF5gG68ogcvuNGvmUKQnIwwP82gGYdZTkSEnzagbYfBE4soP21I34GwzMDHTDPMmQ7H+BnVYUE4gjTRMCEMY7ouxBlyYUHkjKbvRJxcADPC/LQOJHrKYELIjGboQJJ7FbQ67h8x6oKNn1ElkuxCv64GcVIxtBrdjKavRpy7BNrpi12hGaoQ5y4N1IxViHWVQMcXxYwKYeNnVIcI7+iqhenrNAahyaRNm75P80ASpi8MXpnPg35jJWKdJdBreBAWRp3wNSyjOga0DUgNzoVD471712YfwBhGkBSUpWjjdiusjh4kWXLhDOoCDbc+5ZzAyFQnksxZcFq8mkf2wD7ljfkzjFpg42dUR5gnDdrRSKX3P6UdRxiSA7Rx7SGEIClAG9W1Ixix0I1GKdqwsRmZxuglXBuGWXg41MOoO7lr6ECcflZy13gYcZ7CgOTuuKkDCVJgcnfC3IEEZwUndxnVwcbPqBJf0pbi91rNrKStMzCRO2CqQrSjGDoYFG3QUoVIexH0OHJMX4Y502HjZ9SZ3J02fV+ClpK7A8L0ZzRK7g6YvaZv0JiU5O6gpRqRU0UwaixLuh4Ms1iw8TOqY0DbhFBXKpwYw6Q8hCnJigE0IdiZomh2jGIYBxBkT4IL4975YIUNHTBPxcONyWltGC4npXkZRj1wcpdRHen5sUhIcANwQ5IkNDV1IDkqGgkJEgVyhNbcfAAJoZFISaFhFe1Ca2lpQWxQKFLTNIrW1tYGvWX+A74zzOkMGz+jOtLS0pCeni4Gdd+1axeWL18Og8EgdJ9WVlYGh8OBnJwcoe3evRsrV67E0NCQou3duxdr1qyB200nEYZRDxzqYVSJz+DPOeccmEze+L3P4EtLS2GZHlyFtD179qCoqAghISGKtm/fPuTl5SE83DsWL8OoCTZ+RnVQD91n+j6DJ41Mv6SkRDF4CuVQr76goEAxeNLI9KnXHxnJT+Vk1AkbP6M66urqUF5erpi+x+NBQ0MDiouLFdMnWltbkZ+fj4iImcHWOzs7kZWVhagov8HWPZ5TvAYMs7hwjJ9RHWT4ZPQ+0+7v74fRaERzc7OiDQwMQK/X48CBA4o2ODgIrVaL9vZ2MZE2MjIiQj4MoybY+BnVkZubqyR3KX6fmZmJ4OBgJblLoZzs7GxoNBolkVtZWSkMfmpqSsT7SaupqUFhYaE4aTCMmuBQD6NKfFU5FL8PDQ1VNJ/p++L3pFVVVYmTQlycd8AVn+nHx8cjMZEHW2fUBxs/ozookUum7x+/J23//v0ifh8dHa0kcsn0U1NThcn7tNraWvE6OTl5SdeDYRYLNn5GddTX14twj8/0KVZP8f2MjAzExMQo89HNWdSj95k+0d3dLa4G/E2fk7uM2uAYP6M6wrRmjLR1iYnuy+0Y7lNMnSZieHhYSfLSRFitVoyNjYkTAk0+li1btkRrwjCLAxs/ozo2p5egKDlT/P1y2z6kWCJhSo8VcXyivq4emcYouFJTRXJXaPX1GB0dFQlf5/SzeajCh0I/VP3DMGqCQz2MatnV3YwwvQV5cV7DJw4dOoQElxnlSdkB2uHDh6HTecfm9cGPamDUChs/o0oODPfANjmOtWkFikb1/NrBcVyQXhygNTU1BXxWlmVR2UM9foZRI3wNy6iOkakxDNl6cVnhOkWbnJzEeOcAPlWwQdGoZp9i+WT0PijuTw9vM5vNp3y5GeZUwcbPqI4d7XVYkV+C17vrxOvDfd3oGBvEuqxiResfHkRLf7t4NAPF9Ynx8XER36fKH19vn8I9dCJgGDXBxs+ojvLEbGQHxYq/nW4XaieasSI2K0CrszYLU6cqHoJ6/XRVQHF+OgH4NLvdLm74Yhg1wUFMRnUkh8cgPTJeDLSyt7sV67OKkRIRq2i7u1twSf4q8Yx+ehwDVe3QSYAe4DZbo0c9zE76MsyZDhs/o0ocbif+XvcuPl60DgadXtH+WvcOPlG4DsZpjWL6FOunmL4vvEMa9fRJY9Nn1AgbP6M6nJIbf6vbIUzfpPc+YM0lufDXWjL99TAbpgdWl2XF9H0G7wvv0OAtbPqMWuEYP6M6Xmneg00FFagb7BCvDwx24eBANzbklqF+qFNofbZhEdOnxzNQT58Mn0I7ZPr+mu9EwDBqgo2fUR1FcekI03vLMZ1uNw4P9SIvJhXheouiNfQdEj16MnuafAY/W6MqH75zl1EbHOphFpzvfe97WLduHYKCggJGtzoWZLL33XefeGgaDaSyefNmtLS0nNTvL4tOREF8OjKjE1E70I6NOeXIjElStJqBQ7g0f40I8dDke94+JXLpb59GvX5aBzZ+Rm2w8TMLDvWSr7zyStx2220n/Jkf/OAH+NnPfobHHnsM77//vjDhLVu2nHSYhRK5f6GYftF6GPUGRXu29m18vIDi/MaAG7aomscX0yfNd+cumz6jRjjUwyw4DzzwgPj38ccfP+He/sMPP4x7770Xl112mdCefPJJ8bjk5557Dtdcc828k7t/IYMvPBeW6UQuJXefrXsblxesQ7DJjMFJmxLKIdMng6cHsvli/XQSIM3/rl6GUQvc42eWnIMHD6K3t1eEd3yEh4dj9erV2Llz51E/Rz11eqKm/0RsP7AfFywrg1NywTo1jgnHFN48VINN6cVweyShjTkmlV493blLpu+bCH+NYdQG9/iZJYdMn/AfEMX32vfeXGzbtk25uvBn/cUXwBoeCuv0a2N0Ii4oSsWY0QjvfbqAOy4al2VcFvBMHjoR0ERxfX94+EVGbbDxMyfEPffcg+9///vHnKehoUEMd3iq2Lp1K+68807lNfX4aRjFilUrEBYWtmC/47uSYBi1wMbPnBB33XUXbrzxxmPOk5npHfxkviQkJIh/+/r6AnrX9LqsrOyon6ObrGhiGGZ+sPEzJ0RsbKyYFgMa2pDMf/v27YrRUy+bqnvmUxnEMMyJwcldZsGh0awqKyvFv5Qcpb9p8j31kqCQ0N///nclkfq1r30N3/3ud/HPf/4TNTU1uP7665GUlITLL7+c9xDDLDDc42cWHLoR64knnlBel5eXi3/feOMNbNq0SfxNo17ZbDZlnrvvvhsTExP43Oc+JwY9P/fcc/HSSy/xgCgMswhoZC5UZlQChYeoDJROKAud3F2M72WYpYJDPQzDMGcZbPwMwzBnGRzjZ1SDL2q50HX3vu/jqCijFtj4GdXgGz+XbuJarO+nWD/DnOlwcpdRDfRUTaoWKiwsREdHxzETsb67fE90vvr6euTl5SnDMzLMmQz3+BnVQKacnJws/iYzP5EKnBOdj76XTZ9RC9x9YRiGOctg42cYhjnLYONnVAU9tO3+++8/7sPbFno+hjmT4OQuwzDMWQb3+BmGYc4y2PgZhmHOMtj4GYZhzjLY+BmGYc4y2PiZM5rh4WFcd9114iasiIgI3HzzzcqAL48++igyMjLEM/1Xr16NXbt2Cf1Xv/qVGBeAPkODwDz00ENiPqPRCIvFIv4tLi7GCy+8cMTv1dXV4YorrhDz02cffvjhU77ODPNBYeNnzmjI9MmMX331VTz//PN46623xGAuTz/9tBiInUox9+3bh9LSUmzZsgX9/f2YnJzEJZdcgm9+85viO+69917xPTRaWElJiTD/iy66SIz+VVtbG/B79FkaW/h//ud/lLGCGeaMgwZiYZgzkfr6enocp7x7925Fe/HFF2WNRiOXl5fLX/rSlxRdkiQ5KSlJ3rZtm6K98cYb4vO33HKLfNVVV8mXXnppwHyrV6+WP//5zx/199PT0+Wf/OQni7iGDLM4cI+fOWPZuXOnCO+sWLFC0TZv3ixCMFVVVeJvH/ScHXpNn/HhcrnEvxT2IZ3e95+PrhD852cYtcDGz5yx9Pb2Ii4uLkDT6/XiZEBP6oyPjw94j17TZ3z4xvyNjY0Vum9+33yz52cYtcDGz5x23HPPPaLXfqypsbFxqReTYc5Y+LHMzGnHXXfdhRtvvPGY81CClZKrlKz1x+12w2q1ipBNX19fwHv02j8h6xtUZWBgQOi++X3zzZ6fYdQC9/iZ0w4KveTn5x9zopLLtWvXCpPfu3ev8tnXX39dDJFIVTzbt29XdAr90Gv6jA+DwSD+ffPNN4VO7/vPR5VC/vMzjFrgHj9zxlJQUCDKMm+99VY89thjIll7++2345prrsFll12GG264AVlZWfjpT3+KsrIyTExM4KabbhJx+1tuuUWYPPHEE0/gM5/5DP76179i/fr1YojFwcFB7NmzBzExMdi6dSu2bdsm5nU6nWI0Lt/fXV1dqKysREhICLKzs5d0ezDMCbNI1UIMc0oYGhqSr732WjkkJEQOCwuTb7rpJnlsbEy898gjj4jSTGrm+fn58nvvvSf0+++/X2izJ61WK5vNZlmv18tFRUXyv//9b3njxo3yDTfcoPzewYMH5/wszccwZwr8WGaGYZizDI7xMwzDnGWw8TMMw5xlsPEzDMOcZbDxMwzDnGWw8TMMw5xlsPEzDMOcZbDxMwzDnGWw8TMMw5xlsPEzDMOcZbDxMwzDnGWw8TMMw5xlsPEzDMPg7OL/A7tPKx+x/Pd2AAAAAElFTkSuQmCC", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "fig, ax = plt.subplots(figsize=(8, 3))\n", "utils.setup_axes(ax, title=\"rotation clusters (one shared rotation each)\", grid=False)\n", "ax.tripcolor(X[:, 0], X[:, 1], T, facecolors=red.labels.astype(float),\n", " cmap=\"tab20\", edgecolors=\"0.4\", linewidth=0.3)\n", "ax.set_aspect(\"equal\"); fig.tight_layout(); plt.show()" ] }, { "cell_type": "markdown", "id": "bdc5b514", "metadata": {}, "source": [ "## 5 · Rig vs. full vs. reduced\n", "\n", "The reduced solve reproduces the full-order secondary motion closely — the same whip and trail —\n", "while solving a system that is orders of magnitude smaller." ] }, { "cell_type": "code", "execution_count": 10, "id": "f9cb5684", "metadata": { "execution": { "iopub.execute_input": "2026-06-09T05:43:28.226520Z", "iopub.status.busy": "2026-06-09T05:43:28.226468Z", "iopub.status.idle": "2026-06-09T05:43:35.826593Z", "shell.execute_reply": "2026-06-09T05:43:35.826337Z" } }, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "fig, anim = animate_panels([\n", " dict(states=thin(rig_states), T=T, title=\"rig\", color=RIG_C),\n", " dict(states=thin(full_states), T=T, title=f\"full ({full.n_rot} rot, {full.state_dim} DOF)\", color=FULL_C),\n", " dict(states=thin(red_states), T=T, title=f\"reduced ({red.n_rot} rot, {red.state_dim} DOF)\", color=RED_C),\n", "], suptitle=\"Complementary dynamics: full-order vs. fast (skinning eigenmodes)\")\n", "utils.show_video(fig, anim, \"media/23_compare.mp4\", fps=30)" ] }, { "cell_type": "markdown", "id": "9ad803eb", "metadata": {}, "source": [ "## 6 · Timing: how the gap grows with resolution\n", "\n", "The reduced solver's cost is set by the subspace size and the cluster count — **both fixed**, both\n", "independent of the mesh. So as we refine the beam, the full-order cost climbs while the reduced cost\n", "stays nearly flat, and the speedup widens.\n", "\n", "We sweep a handful of mesh resolutions explicitly, building a fresh `FullCD` and `ReducedCD` at each\n", "one and timing the median per-step cost." ] }, { "cell_type": "code", "execution_count": 11, "id": "07c55cd2", "metadata": { "execution": { "iopub.execute_input": "2026-06-09T05:43:35.828432Z", "iopub.status.busy": "2026-06-09T05:43:35.828358Z", "iopub.status.idle": "2026-06-09T05:43:37.745486Z", "shell.execute_reply": "2026-06-09T05:43:37.745268Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " 200 verts / 312 tris: full 5.0 ms reduced 0.67 ms ( 7.5x)\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ " 490 verts / 828 tris: full 12.7 ms reduced 0.67 ms ( 19.0x)\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ " 800 verts / 1386 tris: full 20.3 ms reduced 0.68 ms ( 30.0x)\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ " 1260 verts / 2224 tris: full 32.2 ms reduced 0.67 ms ( 47.9x)\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ " 1800 verts / 3222 tris: full 46.3 ms reduced 0.67 ms ( 69.0x)\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAABQkAAAGZCAYAAAAq17elAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjksIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvJkbTWQAAAAlwSFlzAAAPYQAAD2EBqD+naQAAoMZJREFUeJzt3Qd809UWwPFDS3eh7LL3BtmgKIpsUNmKgoPhYooMB8+JW5QlMhVxIIKgoKgMQQRBRJAtU/Yss7u0pc37nAsJSQe0pW2a5vf15ZXc/Jv8c5v8c3P+956Tx2KxWAQAAAAAAACA2/Jw9g4AAAAAAAAAcC6ChAAAAAAAAICbI0gIAAAAAAAAuDmChAAAAAAAAICbI0gIAAAAAAAAuDmChAAAAAAAAICbI0gIAAAAAAAAuDmChAAAAAAAAICbI0gIAAAAAAAAuDmChAAAAACAG7r77rvNBc73+uuvS548edK17blz5yQ30Oeizykzff755+Z+Dx8+nKn3C7gagoQAkIlOnjxpBi1bt251m36dM2eOTJgwwdm7AQAA4NbeeecdWbRokbN3I0ejj4DrI0gIAJkcJBw9ejRBQgAAAGSZl19+WWJiYhzaCIDdWGp99Oijj5r+LFeuXCb+lQDXQ5AQAJKIjo6mTwAAAHKJqKgoyW3y5s0rvr6+Tnv83DZe9vT0NP2Z1iXcQG5FkBCAy7DmU9mzZ4/06NFD8ufPL4ULF5ahQ4fKpUuXkm0/e/Zsadiwofj5+UmhQoXkoYcekmPHjjlso3l1ateuLf/884/cdddd4u/vL//73/9S3Ydff/1VmjVrJgUKFJDAwECpVq2abfvff/9dGjdubP7dt29fs6960RwnVhs2bJD27dtLUFCQeazmzZvLunXrbup5pkQf55577pGCBQtKQECA1KlTRyZOnOiwzW+//SZ33nmnuV2fT+fOnWX37t0O20RERMizzz4r5cuXFx8fHylWrJi0adNGNm/ebOu/n3/+WY4cOWJ7vrotAABI3Y0+X5OOUW6//XYznqlQoYJMmzYt2f3FxsbKa6+9JpUrVzb3V6ZMGXn++edNe0bGR2rGjBlSqVIls12TJk3kjz/+SHMeNx0Tabv+zMjzSYnO8nrmmWekSJEiki9fPunUqZOcOHEiWX466zhq165d0qtXLzMW0rGbunz5srz55pvmeWk/af/rOM6+n4YPH27GXRaLxdY2ZMgQc58fffSRrS0kJMS0TZ061dY2adIkqVWrlhnj6eM2atTIpGVJjT6GPh99TKvExEQzLtOgVWhoqK39/fffN4HByMhIh+dppf/WYOgXX3xhG5P16dPH4fH0/rRN71/HojpeTUuw73rj5bS+9q43hrY6c+aMPP744xIcHGwCdnXr1jXP50b0OaU0/kxPH6X2Wp4yZYr5m+pzK1mypAwaNMjh72LfP/qaa9GihemfUqVKyZgxY26470BOk9fZOwAA6aWBMx0IvPvuu/LXX3+ZAdvFixflyy+/tG3z9ttvyyuvvGK2feKJJ+Ts2bNm4KYDmy1btpgBitX58+elQ4cOZpD8yCOPmIFJSv7991+57777TMDtjTfeMIOF//77zxbkq1Gjhml/9dVX5amnnjIBOKUDYWtQTh9HB+Y6mPLw8JBZs2ZJy5YtzcBbB+DpfZ4p0UGY7meJEiVMYLF48eIm+PfTTz+Z62rFihVmXypWrGgGUDrw1v654447zBcU60Crf//+smDBAhk8eLDUrFnT9NXatWvN/TVo0EBeeuklCQsLk+PHj8v48ePN7+jADwAApO5Gn69W+rmvJ/10TNCzZ0/59ttvZcCAAeLt7S39+vWzBZU0YKa/r+MPHY/s2LHDfC7v27fPYWllWsdHM2fOlKefftqMYTSYefDgQfMYGlTUIFBGpeX5pEaDObq9Lgu97bbbZPXq1XLvvfemuv0DDzwgVapUMctLrQE/fc4aILr//vtlxIgR5qSqjrO03xcuXGi20fGb9p2O+zTwo3ScpuM2/amBSmub0r5Tn3zyiblN79t6Ynf79u3mMTRYmRINSunYa82aNbY2/R0dW+nj6RjT+hz18erXr5/qOOurr74yz0/Hk/o6UBoMtaf9roFZfc463vv0009NgFoDkDeS0ng5ra+9G42hlY5FNdim7fq+0P2cP3+++btrUM46hr0ZaekjezpG1jRCrVu3Nq/TvXv3mqDwxo0bzb57eXk5vLZ1IkC3bt1MP+v7+4UXXpBbbrnF9BvgMiwA4CJee+01HeFZOnXq5NA+cOBA075t2zZz/fDhwxZPT0/L22+/7bDdjh07LHnz5nVob968ufndadOm3fDxx48fb7Y9e/Zsqtts3LjRbDNr1iyH9sTEREuVKlUs7dq1M/+2io6OtlSoUMHSpk2bdD/PlFy+fNncX7ly5SwXL15Mtg9W9erVsxQrVsxy/vx5W5ver4eHh+Wxxx6ztQUFBVkGDRp0nV6xWO69917zeAAAIG3S8vlqHaOMHTvW1hYbG2v7DI+LizNtX331lfn8/uOPPxx+X8c2+vvr1q1L1/hI71fvXx9HH89qxowZ5v50v6x0vKNthw4dcrjPVatWmXb9md7nk5J//vnH/O6zzz7r0N6nTx/TrmOnpOOonj17Omy7detW0/7EE084tI8cOdK0//bbb+b6mTNnzPUpU6aY66GhoaZ/H3jgAUtwcLDt95555hlLoUKFbOOrzp07W2rVqmVJrw8++MD8XcLDw831jz76yIyrmjRpYnnhhRdMW0JCgqVAgQKWYcOGJXue9gICAiy9e/dO9hjWbfv16+fQ3rVrV0vhwoVvuI+pjZfT+tpLyxh6woQJZpvZs2fb2vQ10bRpU0tgYKCtf1TSv7k+55TGounpo6SvZX0deHt7W9q2bWv63+rjjz8223322WfJ+ufLL790eG0XL17c0r1791SfM5ATsdwYgMvRaf72dAmI+uWXX8zP77//3pzZ1LN4586ds110Rp2eUV61apXD7+vZTF1ucSPWs+s//PCDuf/00GrH+/fvN2eS9UysdZ90yUOrVq3MGeSk93mj55kSnQVw6NAhc9bffraksi63OHXqlNkfPTOrMwKs9OyuLnWyv3+9Dz0DrgVZAABA5kjr56suL9UZfVY6406v67JMXfqpdLaVzuCqXr26w7hHVyoo67gnreOjTZs2mfvX2Y76eFY6btAlqjcjLc8nJUuXLjU/Bw4cmOLYKCW6//as4xv7pb1KZxQqTZ+iihYtavrSOrtPZ4zp0t/nnnvOLDHW8Zx1Zp8un7WOr/RvqisrdJZZeujMxYSEBPnzzz9t96tterHOVty5c6eZTWddpZJRSftE70/HpeHh4Tf83ZTGy2l97aVlDK1/H30t6gxTK52pp7MzdYm1zhzNTrrqJi4uzoypdVan1ZNPPmlSAVlfL1Y6w1NnWNq/tnXGos7CBVwJQUIALkcHsvZ0mYB+eFtziOjgTU8y6nY60LO/6HISHYja05wh9oNgXeJx+vRp2+XChQum/cEHHzRLQnSZgi6x0OUWuuwlLQFD64Cyd+/eyfZJl3po3hZ93PQ8z5QcOHDA/LQuj0mJ5g9UmgsmKR3oWYOXSnOp6MBUlxbpQEeXXTDYAQDg5qT181VzoGnuYHtVq1Y1P+3HPbqcM+n4wrqdddyT1vGRdZyQdByiARtNU3Iz0vJ8UqL7pGMgXYJqT/PgpSbpttb7SPo7GpjSIJb1eSv7AJ3+1NyCetGTq3pdg2rbtm1zCNrp0lINFOnfU/tOT/YmzTudEl1erjns7B9P71eXMWvAVpctW2+z5lbMqLJlyzpc17yJ1qWyN5J0vJye115axtDa/9pv9gE569jUent2Sm28rH2g74Ok+1O6dOlkRU+0f9PSt0BOQk5CAC4v6QeyDji0bcmSJebMb1JJc7lo4mx7mvPEPkmyFhfRxNu6nZ5V1rOievZQz2rPmzfPnC1dvnx5io9lv0/qgw8+kHr16qW4zY1y+Tmj2prONtCBqubp0eeo+695a3Q2AvlVAABw/uerjjE079m4ceNSvN2aQzC946O0SG1sojPjnC3p+C494ykNxmmOQQ3cWoN2+nvartc12Kn9aR8k1GCW5qzTHNA6Rvzuu+9M0QvNVa157VKjwddbb73VjDE1H5+eoNb71WBafHy8mXGqj6mz9TQAdzNSG6vaF2lJT3+m9bV3M2NoV3kd3kzfAjkJQUIALkfPWtqfHdYBlQ5SrMU2dMadfiDrNtYzmemhFdnslwtYz7IqPbupy4P1ogMiTYatxTt00KNJjVMbpFiTIuvyBN0uM57n9R5HZyek9jjlypUzP3Ugm5RWVNYqe/Zn+bUAii7v0YueEdYz3pr43PolxhnBSwAAXN2NPl+VLkfW2f32n8taEELZj3t0VpuOTa73mZzW8ZF1nKDjEOuyUaUBK01pohVnk46RklZ7TW3WV1qeT2r7pGMgfXz7GY46Nkor633o87LOTlO6hFj33/q8lTX4p8XgdPnwiy++aK7r7D4tXGGdEanF6Oxpm86a04suVdUiFvo3HTVqlKnWmxp9PA0S6xJXHYdpQFD/llpVVwOEetHCHzeS3WOytL720jKG1v7Xoi36N7KfTahjU2X/90lKX4dJX4OpvQ7T2kf242X7GbT6d9XXYVrH84CrYbkxAJczefJkh+talU9ZB9U6INOzeXrWNunZO72uuVeuR6sM6ge/9WIdAFqXHduzzgrU5cLKOuhNOlDR+9CB1IcffmjyqiSl1QXT+zxTol8wdPA/YcKEZPtg7Qv9UqL7rbMl7bfRwKKezdWqg9azr0mXQGsFPB0YW5+v9Tkn3Q4AAKQsrZ+v6vLlyzJ9+nSHAIVe1xll1vGJzko8ceKEmfmWlFaMtaYQSev4SJfV6v1PmzbNPJ7V559/nmxsYT05aV+dV5/fjBkzUnzuaXk+KWnXrp35qTPzUhobpYV1fKNjJHvWWXD2lZJ1LKXLa7VKrwZHdamsNZinqV20cq1WWNYci1ZJx5e6LFXHlNq3eh/Xo/erf3vdN/s8h9quFXk1uJqWfIQ6JkspWJZV0vraS8sYWv8+OotSZxjav170b6yzXHVlT2r0dajvKQ0yWmkObmvF6oz0kX4H0L/hRx995PB+0crf+ljXq6wNuDJmEgJwOXr2rlOnTtK+fXtZv369zJ492xQEsZ7Z1oHCW2+9Zc7aan6bLl26SL58+czv6WDhqaeekpEjR6b7cd944w0zCNZBgZ5d1LP+OljVHCTWHDH62JrXRgfW+pg6ENElJDrY1NyDGuDTs8Ka+FkHnzqw0jOoOsNw8eLF6XqeKdEzr3qGu2PHjmbwpY+jQUE9C6s5Y5YtW2a202VNui9NmzaVxx9/3AzkdBCmCck1L5KKiIgwz+3+++83j6kDND3DrWfUx44da3tMHdTrgE4TgTdu3Nhsp48PAACSS+vnq9LAoc4w0/GMzv7Tz1stPqZBOF2mqh599FGT302LUuiYQgNaGqjTz35t189+DfyldXyk96vbaUERnUmos+J0m1mzZiXLSahjGg2W6X1qIEhz9s2dO9cEd1KSlueTEh1rdO/e3QTRNBinj6mFLKyzENMyO0z7WnND62NpkEiDTn///bc5aap90aJFC4ftNSinz0WX01pnTOrJWB3b6ePqmMxe27ZtTX5D7X9dKqx5Hj/++GMzbtR+vh4dj2nAUWet6d/Byjpz0bo/N6L9pK8lDXxqX+v4U8ehWSWtr720jKH1eWvAWAvkaBEbnVmqwVjN66h/9+v1oeY41JyQXbt2NYVOoqOjTb/pa2zz5s0Z6iMNXOvrWoPqOhbXMbn+fXS/dbxrv+oIyFWcXV4ZANLqtdde09N4ll27dlnuv/9+S758+SwFCxa0DB482BITE5Ns+++++87SrFkzS0BAgLlUr17dMmjQIMvevXtt2zRv3txSq1atND3+ypUrLZ07d7aULFnS4u3tbX727NnTsm/fPoftfvjhB0vNmjUtefPmNfs7a9Ys221btmyxdOvWzVK4cGGLj4+PpVy5cpYePXqY+87o80zJ2rVrLW3atDG/q8+9Tp06lkmTJjlss2LFCssdd9xh8fPzs+TPn9/SsWNH85hWsbGxlueee85St25d2/3ov6dMmeJwP5GRkZZevXpZChQoYPZbnxMAAEhZWj9frWOUTZs2WZo2bWrx9fU1n7Eff/xxsvuMi4uzvP/++2Z7HV/ouKFhw4aW0aNHW8LCwtI9PlK6PxUqVDD316hRI8uaNWvMPunF3oEDByytW7c22wUHB1v+97//WX799VczJli1alWGnk9KoqKizH4WKlTIEhgYaOnSpYvZZ32c9957L9k46uzZs8nuIz4+3vSJPi8vLy9LmTJlLKNGjbJcunQp2baTJ0829zNgwACHdn2u2m4/dlPTp0+33HXXXbYxXqVKlczfOWn/p6Zx48bmfjds2GBrO378uGnT/UzK+jzt7dmzx+yDju30tt69e1+3T3SMqu2HDh267r5db7ycltdeWsfQISEhlr59+1qKFClitrvlllscxtFWus/6nOwtX77cUrt2bfN71apVs8yePTtdfZRaX+jrU98j+nrR17e+Hi5evJim/tH7ZlwMV5NH/8/ZgUoASAud4aZn83RpruZrya3c5XkCAIDU3X333XLu3DmTDiQ3yIrno7MQ69evb1ZbPPzww5l2vwDgrshJCAAAAADI0TQ1SlK6DFVTreiyXADAzSMnIQAAAAAgRxszZozJVae5AzV/35IlS8xFc9mVKVPG2bsHALkCQUIAAAAAQI52++23y6+//ipvvvmmREZGStmyZU2KlpdeesnZuwYAuQY5CQEAAAAAAAA3R05CAAAAAAAAwM25/XLjxMREOXnypOTLl0/y5Mnj7L8HAABAlrNYLBIRESElS5Y0Sf/dBeM+AADgjixpHPu5fZBQA4QkugUAAO7o2LFjUrp0aXEXjPsAAIA7O3aDsZ/bBwl1BqG1o/Lnz5+lZ67PnDkjxYoVc6sz9vboA/qB1wLvCY4LHB/5nMgZn5Xh4eHmJKl1HOQuGPdlL8Z+9AOvBd4THBc4PvIZ4VpjP7cPElqXGGuAMKuDhDExMeYx3DlI6O59oOgH+oDXAu8HjgscH3PKZ4S7pVph3Je9GPPQD7wWeE9wXOD4yGeEa4393DdSAwAAAAAAAMAgSAgAAAAAAAC4OYKEAAAAAAAAgJtz+5yEN5KQkCDx8fGZss5c7+fSpUtum4/P1fvAy8tLPD09nb0bAAAAAAAAmY4g4XVERkbK8ePHxWKx3HRH631okEzv092ShOeWPtB91lLhgYGBzt4VAABynfLly8uRI0eStQ8cOFAmT55sTjKOGDFC5s6dK7GxsdKuXTuZMmWKBAcHO2V/AQAAchuChNeZQagBQn9/fylatOhNB7U0QKaz6HQ2misGyDKDK/eB7vvZs2fNa6JKlSrMKAQAIJNt3LjRjL+sdu7cKW3atJEHHnjAXB82bJj8/PPPMn/+fAkKCpLBgwdLt27dZN26dfwtAAAAMgFBwlRoMEsDQxog9PPzu+mO1vvSpaquGCDLLK7eB/paOHz4sHltsOwYAIDM/5y1995770mlSpWkefPmEhYWJjNnzpQ5c+ZIy5Ytze2zZs2SGjVqyF9//SW33XZbivepMw71YhUeHm5+6soGvWQV631n5WO4AvqBfuC1wHuC4wLHRz4jcsZnZVrv322DhLpsRS/2Z6xT4orBLGQNXgsAAGSPuLg4mT17tgwfPtx8/v7zzz/mJF3r1q1t21SvXl3Kli0r69evTzVI+O6778ro0aOTtZ85c0ZiYmIkq+kqBNAPvB54T3Bs4PjI5wSflSnReNSuXbvkwoULUqhQIalZs2aWTUiKiIhI03ZuGyQcNGiQuegZZV2yAgAAgJxh0aJFEhoaKn369DHXT58+Ld7e3lKgQAGH7TQfod6WmlGjRplAo5WO+8qUKSPFihWT/PnzZ+nZeg0Q6uxIVyzWllnoB/qB1wLvCY4LHB/5jEjZqlWrZOzYsebEpZWOTzT/cosWLSSzpXWFrNsGCQEAAHKSpTtPyYQV++Xg2UipWDRQnm1dRdrXLiHuSJcWd+jQQUqWLHlT9+Pj42MuSWngLjuCd9n1ODkd/UA/8FrgPcFxgeMjnxHX/Pbbb/LCCy9IUhow1PYxY8bY0qtklrSORxi1uDg9S60vHj0bbk3snZrXX39dunTpYruuy3e2bt2apfunj6ePCwAAUrdkxynpP3uz7D0dIXEJFvNTr2vg0N1oheMVK1bIE088YWsrXry4WYKsswvthYSEmNsAAABcZYnxhx9+eN1tdIbhjVLjZRVmEmbxbIBD56KkQpEAGdqqirSqVjjTH2f69OlmzboOmjlTDQBAzi7gdS4yTg6fj5LD56Ku/ow2Y4Xdp64U1LBYtzUn80QmrtzvdrMJtSCJLre59957bW0NGzY0hc9Wrlwp3bt3N2179+6Vo0ePStOmTZ24twAAAGm3ZcsWhyXGKdGToLpdo0aNJLsRJMyiAKGe/c9zdZCvswEGfL1ZJj1UR+6rWzpTH+vQoUNSq1YtpwcINZm4Dt5z+n0CAJDVgcDzUXFy5HyUHDoXbYKBh85HmesaEIyMvZyO+xI5eDZK3C2HnQYJe/fuLXnzXhumav7oxx9/3OQX1MTeuoJiyJAhJkCYWtESAACAnObAgQNp2u7cuXPiDAQJ06HjpLVyNiL2htudi4xNNhtADft2h7z9y740PVbRfD6yeEiz626jy4s1sbcuG/7000+lYMGC5ky7tllpgm+9fvfdd0tGqt9o0szFixeb6507dzbTXgMCAuTw4cNSoUIF+eyzz+Ttt98222q0+7vvvjNr6HUZdI8ePeTyZccvQxoNf/HFF2Xbtm1mkK/bPvnkk+Y2XZa8adMmk1B83rx50rdvX/N4AADktEDgxeh4MwPQNiPw/JWAoF4i0hEIVJ4eecQjj0h8gnXEcIXOJKxYNEDciS4z1tmB/fr1S3bb+PHjzUlRnUkYGxsr7dq1kylTpjhlPwEAANIjJiZGvvzyS/n888/TtH2RIkXEGQgSpoMGCE+HX8pwZ19OtNzU7yc1f/58U/VPA4ETJkwwQbbMzDE4dOhQEwzcuXOn+UJ0//33y7Bhw2TGjBm2bX788UcT2NOKg/v27ZNevXrJggULTLJxDVwOHjzYNkVWqw/ec889ZkCv97V7925p27atVKxYUVq1amW2Wbp0qfm9SZMmmdxDAAA4y8WoODML8Eog8GoQ0MwQjJKIS+kPBJYu6CflCwdI+cL+Ur5IgLlUKBwgpQr6ycrdIVdWIeS5MoPQ+nNoq6riTnRcoGOOlPj6+srkyZPNBQAAwFVWSSxbtszEOG60zNgqODhY6tevL85AkDAddHZfWuhMQg0IJutsjzxSJNAnUx8rK1/IX3/9taxZs0YKF76SS/Gdd94xRVKmTZtm2+61114zQUqls/802NexY0dzvX///jJx4kTbtl999ZU0a9bMzDDU2Y+1a9c2swXnzJljCxJqmwY+lf0yIwAAskJodJwJ+h05fyU34LV8gdESFhOfrvvS2YClrgYCNR+xCQgW8Tc/Sxf0F++8qacG0byD0x5pIBNX7JcDZyOlUtFAGdq6qrSvTVEOAAAAV7Rjxw6zOlInXllpTYnbb79d/vjjj1R/T1d06nbOQBQmHW60/DdZTsIkswEm9Kgj99YtZQJkOU1gYKDt30uWLJGqVauamXzly5e3teuMP13eY782vmzZsrZ/nzx5UsqVK+dwv/bXdVaizhTUZdFWWrHnzjvvTPH+AADIDGHR8VeXBF+ZBWibGXg+SkKj0xcI1I/wUgX8bEHAcoX9r/y7SICUuUEg8EY0UNi2ZrBJ36FnkJ2dbxgAAADpp2O5jz/+2MRW7Gns49lnnzVxkt9++81UObafXajjPw0Q6uQsZyFImAVsswFW7jcJxzWfkFY3blk186sbJw30RUdH265HRUVJePiVaok3EhkZmWwmoS4h1sCevlCV/tvHx8esjdd8Qcr+C0zJkiVl/fr1Dvej21kTimuuQc1rqDMOUwuU8oUIAJAR4ZfirxQJMUHA6CuFQ67OCtT8gemhH1Elg64EAm1BQDMrMEDKFPITn7zOObMLAACAnOvSpUsm7+AXX3xhJljZT7jS4mv2xdY0ENi8eXPZvHmzKWZSqVIladCggdNmEFoRJMzCQKFerDS/jlbrzUr6gho9erTs2bPHzAD83//+l+FZixqs0/yCL730knz77bdm//X+Hn300VQDebqM+K233pKff/7ZJBPX6oSap9BKf1eTjmtxEw0Wqn///df0S+PGjTP4rAEA7iLCBAKjr1QLPnctCKhLhbWicEYCgRoEtOYGNHkCC/tLmUL+4utFIBAAAAA3pvESa95BnUVoFRQUJE8//bR069YtxXRqGhDU4rOlS5fOMatICBLmIhqJ1hegrm/39/c3hUzy5cuX4fvTfIIa7a5Zs6a53qlTp+tWG65WrZrJO/jMM8+YJclafbl9+/a220uVKiU//fSTvPzyyyZfoc5WrFGjhrzxxhsZ3kcAQO4SFZcgO0+EydGLMVdnBl5ZFqwzA89Fpr+gVYkgX9ssQGvBEJ0ZWJZAIAAAAG7Szp07TZxE8w/aB/90EtWTTz4p+fPnd6k+Jkjo4pKWz9Y17XqxeuKJJ2z/1qChvdSqB1rpi1krDadEZyqm9Pv6RtBLaurVq2ci7CnNcEy6fwCA3Ckq9vLVAiFXAoBXCodc+ZmRQGDx/L62AiFXgoHXAoF+3swIBAAAQOY6c+aMmTmYNO+gFmvVvIP29R1cCUFCAACQ6aLjLtuCgLaKwVeXCp+NuJajJa2C8/tIOQ3+XQ0EVijib67rcmF/b4YzAAAAyJ68g1999ZXJO6j/tqpQoYJZidm0aVOX/jO47ah68uTJ5qLVdQEAQPrFxCXYlgKbZcFX8wTq9ZDw9AcCi+bzkVL5vKRKiQJSoejVYiEmKEggEAAAADkz7+BTTz0l3bt3TzHvoKtx/WeQQYMGDTIXrf6rf1QAAJDcpfgEUxjEVA22zgi8ulT4dPi1s6dpVSTQxzYL8FrV4CvX/b08zKArpyRuBgAAAHbu3Cnjxo2T7du3O+Qd1DoMmncwN8WU3DZICAAArgUCj164Ggg0QcBoWzDwVFhGAoHeJuh3JTfglWIh5a8uDc7n65Xq72lBKwAAACCn5B38+OOP5ZdffnFov+OOO0zeQV1inNsQJAQAwE0CgceuBgLNzEBbnsAoORV+SW5QyyqZQgHetmrB1oIhmi+wXBF/yX+dQCAAAACQk126dElmz55tCsXa5x3UYiSad/D222+X3IogIQAAuUTs5SuBQPuqwdalwSfDYtIdCCzo73UtCHh1WbAuEdZZgkF+BAIBAACQu/IO/vrrr/LRRx/J6dOnbe358+c3eQfvv//+XJF38Hpy97MDACCXibucKMcuXi0ScjUIaM0ZeDI0RhLTGQgs4O91tWrwlVmB1iCgzgoM8icQCAAAgNzv33//lbFjxybLO6iBQc07WKBAAXEHBAkBAMhh4hMSr8wIvFo1+Er14CsBwRMX0x8IzO+b90qRELtCIVfyBQZIAX/vrHoaAAAAQI529uxZmTx5svz0008O7bqkeNiwYbky7+D1ECTMBS/oBx98UDZt2iTt2rWT+fPnZ/ljtm/fXp555hm55557MuX+tm7dKvXr1zdTe7ObJhsNDQ01uQYiIiKkXr16smHDBilSpEi27wsA9wsEHr8Yc3VGYKTsOnZOzkQfNUVDToTGSEI6I4H5rgYC7WcFWoOCumw4T548WfZcAAAAAFeiuQa//vprEwuIiYlxyDs4bNgwU5zEHREkzCIrjqyQqdumyuGww1I+qLz0r9NfmpdsnumPM336dDMFVgNdHh4eGb6fw4cPmwj5xYsXrzuNdtWqVSYwmVkBwpwkX7588thjj8nbb78t48ePd/buAMgFLlsDgdYiIVeXBevMwGMX0x8IDPTJ6zAL0H5WoBYSIRAIAAAApE4nJ61YscLkHTx16pRDPODpp592i7yD1+O+zzyLA4TDfh8meSSPWMQi+y/ul+Grh8sHzT6QdhXbZepjHTp0SGrVqnVTAcL00Gm4ffv2TfX2+Ph48fJy3RxWvXv3NrMJNVDo7+/v7N0B4CKBwJOhl2zVgq1BQA0I6pLhy+kMBAZ4e9pmAZpqwYWvFAvR64UJBAIAAAAZsnv3bpN3UFczWumkq+7du5vCJAXcJO/g9RAkTIcHf3pQzsWcu+F252POm58aILT/OWrdKPlg8wdpeqwifkVk3n3zrrvNAw88IIsWLTIzRz799FOZOHGitGnTRh5//HHzor98+bJZR6+BPZ0yq7RSz4gRI0xwUYNg3bp1k6lTp0qTJk3M7aVLl7bNUHz44YeTBQCXLl0qb7zxhq1Np+ZOmDBBOnfubH5Hp+R+9913MnfuXHnnnXfk6NGjUqVKFbNvTZs2Nb+jsx41Qr98+XIpXry4DBw40OFxdF/1Prt06WKu63PUZcE621GFh4fL//73P1m8eLGZ+VitWjX5/vvvpUyZMhIZGSkvvvii/Pjjj2b6sC6NnjRpkgQFBZnfXbNmjQwaNMg8/7Zt20rBggWTPXbhwoVl9erV0qFDhzT9rQDkfjrjT4uC2FcLts4O1CIi8QnpCwT6ayDQbhagBgKD8sRK/SqlpVh+X2YEAgAAAJnk3LlztryD9mnObrvtNrO0uFKlSvT1VQQJ00EDhGeiz0hGXbZcvqnfT0rzD/bp08dEuzWopjSQNnz4cGnRooXExcWZgKFW4tHgoHWm3Pvvvy+PPvqoREVFybZt20z733//bZYbHz9+PNXo+f79+yU6OtoE5ezt3LnTRN41IKiByV9++UVGjhxpAnU6K0+DfB07dpS9e/ea0uFDhw41gULdV72/Tp06pet563PW31u/fr0JMupz8PPzM7f169fPTA3WikQ6o/GJJ56QwYMHy1dffWUCivpY+vy1X5YsWWKmEvfs2dPh/mvWrGmCrAQJAfcMBJrg3/kr1YPNzEBdGnwh/YFAPy9Ph1mA9rMCi+bzcQgEJiYmSkhISLJ2AAAAABkTGxtr8g7OmjXLIe9g2bJlTdxEJzkx9nZEkDAddHZfWuhMwgRLQrL2vHnySiG/Qpn6WEnpTDjrrEFfX1956aWXTHRcv4DqkmQNnP33338mr2DRokXNTMO00iCbzj7U6bj2dJaePo7ev7e3t4nQP/fcc9KgQQNzu85W1Cm9GjzU2Y/z5s2TP/74wwQj9aLbPvTQQ2naB/0SvXDhQjly5IiULFnStGnRE6XPSWcx6lkCa6BTZz3qcmyd8ahnDfR3dBaj0sBly5Ytkz2GBjL1uQJwLUt3npIJK/ab2X4aiHu2dRVpX7uEwzaJiRY5FX7JtizYmidQA4NHz0dLXEJiuh7T18vjyoxAMytQf14pGKKPX4yAHwAAAJDtdLbgypUrzYrGpHkHdVmxxiXcOe/g9dAr6XCj5b+p5SS0/nz3jndNTsKsjFRroExn6mkQLiwszBY918q9GszTAJvm29PZgOXKlZNRo0ZJjx490nTfujRXZ/AlJCQ4BApLlSrlkBNRZwjqcuDXXnvNYanyiRMnTABPZzjqY1vZ//tGNDjo4+NjIv9J6eNqMDRpiXLdt9OnT8vJkyeTPZZe12XJ9nQ5c+3atdO8TwByRoCw/+zNokdXne+393SEud6rSRkJ9PWyBQSPXIiWuMvpCwT65PVwWBpsrRhsDQR6eDDzDwAAAMgJ9uzZYyYpbdmyxSEmoKsfdcIQeQevjyBhFmhdrrWMv3u8TNs2TQ6FHZIKQRWkf93+cleJuySradBPA3mbN282MwV12azOtLOuu9fZfTrbToNpugxYA4TNmzdPU+ETzS2oMwl12bAuybVK+ruaG3DIkCHSv39/h3bdBw3I6WxGDfYFBwebdl2mbC8wMNA8Byv7yL8G9TToeezYMfM4SR9X90WDgSkVHdFZhPq49vSxixUr5tC2a9euZHkSAeRMOjPwwNlIGb14l7luXRBs/Tnn72Npuh/vvB5SrtC1WYD2QcHi+X0JBAIAAAAumHfw1ltvNXkHK1eu7NT9cxUECbMwUKgXK32R6my6rKaz4DRAptHx8+fPy+jRo2236Qw+Xep73333mVmB1gi6TrPVJbYaYDtw4IA0bNgwxfvW4F67du1k1apVDkHCpLQwiBYaady4sQlK6tr/P//808xe1MCgBiZfffVVU9xEg4EffOBYzEV/55tvvjHLlDXgp290K/19LZKiAUgt1qLXNSehzizU/IRa7ERzEI4ZM0aKFCliZhBq7sKuXbvKvffea2775JNPTIXmZcuWyW+//eaw1FmDiHpwueuurA/oAki/kPBLsuVoqGw7HirbjoXK9uNhEhl7OU2/6+3pIWV1ObCZBegv5a7OBtTAYAkCgQAAAIDL0UlEc+bMMXkH7ScbaYxAg4PNmjUj72A6ECTMZTQoqMVJNAiolYo1GafOGLTSN48G8DRgqG8ava7VfJUuD9ZiHXrblClTpFevXikGAF944QXzMzWa609nDGrBlIMHD5rlwVo9+eOPPza3f/TRRyYPgM4KLFGihJm1t2nTJtvvv/XWW6ayss6E1HyCjz32mNkfqy+++MLsQ6NGjcwy6ho1apjZkUpzD+rz0AClBkk1iPjggw+aIGGhQoXkhx9+MIFCPVhoJWh9HF0+bfXll1+awigBAQE3/bcAcHMiLsXLjuNhsvVqQHDbsTA5He6YHuBGdBbghw/UNbMCSwT5iSdLgwEAAIBck3dQ4ws6uch+ZaLGInRykk50QvrksdjPw3RDOvNOc/Vp/j6dTWelQa5Dhw6Z/HZaAORmWWcS6ovU1avn6GxCDTSmt/pvTu8DDTjq0mydeagByqQy6zVhrWKqAcy0LPPOjegD+iHpa6FQkaKyLyTKLiAYKv+djZQbfUJpELBumSBTSXjR1pO2nIR6iNHfnfZIQ2lfu7jkdLwn6Ifsfi2kNv7J7bLrefOeph94PfCe4NjA8ZHPiaz7rNS8g+PGjTNp1qz0PnQ1oq46dKW8g4k5bOzHTEKkmy7TzY200pFWfgaQtfSEwdEL0bL1WKhsOXpRNh06J/vObr1hQZFAn7xSp3SQ1C1TQOqWLiD1yhSQ4kHXAvYaDJy4cr8cPBslFYsGyNBWVV0iQAgAAADgxjQ1mK4yXLx4sUPeQV25qKsoyTt48wgSAgCy1PnIWJNDcOuxsCuzBI+HSmj09XO05vXIIzVK5DezBK0BwUpFA69bQKR97RLmAgAAACB35R3UugWadzAqKsrWrinUdJXjnXfemSNXK7oitw0SajEMvdjnowMA3JyYuAT592SYmSW49WpA8NiFmBv+XrnC/iYQqAFBnSlYq2R+8fXy5M8BAAAAuCmdLaiFUydOnCgnTpxwyDv4xBNPmPoD5B3MXG4bJNTCG3qxrstOjZunbIQdXguAo4REi/x3JlK2HrtomyW4NyTCtF9PoQBvW0CwTun8UsI7TqqWL+W2OToBAAAAONq7d6+MHTs2Wd5BLUqqeQe1WCsyn9sGCW/E0/PKDBat9Ovn5+fs3UEOoK8F+9cG4G5B8lNhl0wg0DpLcMeJMImOu/5sbF8vD7mlVJBthqAGB0sX9LMtB7Am6gUAAACA8+fPy9SpU+WHH35wmKjTuHFjk3ewSpUqdFIWIkiYWsfkzSv+/v5y9uxZM331Zme4WCv76vJmd10r78p9oIEMfS3oa0JfG0BuFxYTLzuOh12bJXg8VM5GxF73dzRdYNXgfLaAoOYTrBacT/J6MkMQAAAAwPUn5Wjewc8++8wh72Dp0qVN3sHmzZu7XBzBFRHtSIW++EqUKCGHDh2SI0eOZEqATANNGmx01xe2q/eB7rcmRnXFfQeuJ/Zyguw5FXGluMjRUNl6PNRUCL6RkkG+Uq/stTyCtUsFmQrEAAAAAJCevIMTJkxwyDsYEBAgTz75pPTo0UO8vb3pzGzCt7nr0BeiTmW1LjO9GRoc03LdRYoUcdu8W67eB/p6cMX9BuwlJlrk8PkoExDcdixMthwLld0nwyUuIfG6HZXPN69DYZG6pYOkWH5fOhcAAABAhuikrNdff13++ecfW5t+5+7SpYvJO1ioUCF6NpsRJLwBfYH6+vpmSoBMly3rfblroIk+ALKfLhHWPIJmlqD+PBYq4ZcuX/d3vD09pEbJ/FKvdNDVZcMFpELhAPHQ9cQAAAAAcBMuXLggU6ZMSZZ3sFGjRjJixAjyDjoRQUIAyCWiYi/LzhNhdgHBMDkRGnPD36tYJODKLMGrlxol8olPXgr0AAAAAMg8ukpz7ty5MnPmTPIO5lAECQHABV1OSJR9IZFXlw1fCQruC4mQxGsn4lJUJNDHBATrlbkyS7BOqQIS5O+VXbsNAAAAwM3obMHVq1ebvIPHjx+3tfv5+ckTTzwhPXv2JO9gDkGQEABc4EP1+MUYh4DgjhNhcin++nkE/bw85ZbSQVLfbpagFhuh+A4AAACA7LB//34ZN26cbNy40dam30c6d+4s3bp1k+rVq7ttSraciCAhAOQwodFxsu142JVcglfzCZ6LvH4BJU+PPFI1OJ/DLMHKRQMlrycfuABch1Y1fOGFF2TJkiUSHR0tlStXllmzZpkcRdaTJq+99pp88sknEhoaKnfccYdMnTqV3EUAAOTAvIP6Ga15B7U+gVXDhg1l+PDh5rM7JCTEqfuI5AgSAoATxV5OlC1HL8r2E+FXA4Jhcuhc1A1/r3RBPxMItM4SrFUyv/h7c0gH4LouXrxogn4tWrQwQcKiRYua2QcFCxa0bTNmzBj56KOP5IsvvpAKFSrIK6+8Iu3atZNdu3ZlSqE5AABwc+Lj42XevHnmhF5U1LXvNaVKlZJnn31W7r77bjOT0D5wiJyDb5QAkE0SEy1y8FykbD0WZls2vPtUuFy+QSLBID8vEwi0zhKsU7qAyS0IALnJ+++/L2XKlDEzB600EGilswg1l9HLL79sliipL7/8UoKDg2XRokXy0EMPJbvP2NhYc7EKDw83P/WLSVZ+ObHet7t/AaIf6AdeC7wnOC64z/FRP6fXrFljTuYdO3bM1h4QECB9+/Y1n9Pe3t5mO73kxj7IiOzqh7TeP0FCAMgiIeGXrlYZvrJkePuxMImIvXzd3/HO62FmBV4JCBaQuqULSLnC/uQRBJDr/fjjj2ZW4AMPPGCSm+uMg4EDB8qTTz5pbj906JCcPn1aWrdubfudoKAgufXWW2X9+vUpBgnfffddGT16dLL2M2fOSEzMjau/36yzZ89m+WO4AvqBfuC1wHuC40LuPj4ePnxYPv30U9m2bZutTWcLtmnTRh555BGzKkBXDOTmPrhZWd0PERERadqOICEAZILI2Muy3RQWuZpL8HionAq7dN3fyZNHpFxBX2lYvrDUK1tQ6pUuINWK5zOBQgBwNwcPHjS5izRP0f/+9z+T4PyZZ54xsw569+5tAoRKZw7a0+vW25IaNWqUuT/7mYQ6W7FYsWKSP3/+LD1br4N9XTLtzsnY6Qf6gdcC7wmOC7n7+KiBv+nTp5sZ/fYz1Ro0aGA+f6tWrZrr++BmZVc/aCXptCBICADpFJ+QKHtPRzjMEtx/JlIs1181LMXy+VyZIVi2gAkI1iyZT2LCLpgvuO78wQgA1kGyFih55513zPX69evLzp07Zdq0aSZImBE+Pj7mkpQec7PjuJtdj5PT0Q/0A68F3hMcF3LX8dGad1BnD0ZGRtradRXA0KFDTX5hnUmYm/sgs2V1P6T1vgkSAsB1aL6MoxeirwYEw0xAcOeJMFNw5HoCvD1N7sBruQQLSPEg32RfiGPC6H4AUCVKlJCaNWs6dEaNGjXku+++M/8uXry4+amVEHVbK71er149OhEAgGzKO6g5gu3zDvr7+0u/fv2kZ8+eKZ6cg+sgSAgAdi5ExdmKimhAUP99MTr++gdSjzxSvUQ+kz/QGhCsWDRQPD3SdvYMACCmsvHevXsdumLfvn1Srlw5WxETDRSuXLnSFhTU5cMbNmyQAQMG0IUAAGSh//77T8aNGyd///23rU1nC3bq1Ml8DhcpUoT+zwUIEgJwW5fiE8yswCsBwSu5BHXW4I1oIRFrQFBnCmqhEV8vz2zZZwDIrYYNGya33367WW7co0cP8yVkxowZ5mL9IvLss8/KW2+9JVWqVDFBw1deeUVKliwpXbp0cfbuAwCQK4WGhprUH99//71D3kFNCzJixAipXr26U/cPmYsgIQC3kJBokQNnI2Xr0VDZenWG4J7TEab9egoFeEvd0kG2ZcMaHCwY4J1t+w0A7qJx48aycOFCU2zkjTfeMEFAXc708MMP27Z5/vnnJSoqSp566inzpaVZs2aydOlS8fV1TOcAAABuPu/gt99+K5988olD3kE9OaeFxVq1apXmvINwHQQJAeTKXBmnwy85BAR3HA+TqLiE6/6eT14PuaXUtYCgXkoX9OPDDwCyyX333WcuqdEvIxpA1AsAAMia71Jr166V8ePHy9GjRx3yDvbt21d69epF3sFcjCAhAJcXfiletl8tKmKtOHwmIva6v6MnvaoWyyd1y1wLClYNzidenlTWAgAAAOB+Dhw4YPIOar5f+xN0HTt2lIEDB5J30A0QJATgUuIuJ8qe0+EmGGgNCB44G3XD3ysZ5GuCgdaAYO1SQRLowyEQAAAAgHvTFB7Tp083eQcTEq6tvtJCYZp3sEaNGk7dP2QfviEDyNFT3Q+fj7ZVG9bLrpPhEpdwLWFuSvL55jW5A3WWYL0yBU1OwWL5yVcFAAAAAFaXL1+25R2MiIiwtZcoUUKGDh1K3kE3RJAQQLZauvOUTFixXw6ejZSKRQPl2dZVpH3tEua2sxGxsv3qkmG9bD8eJmEx8de9Py/PPFKzRP4rswS14nDZAlKhcIB4eJBEFwAAAABSmoyxbt06k3fwyJEjtnY/Pz+Td1CLhvn4+NBxboggIYBsDRD2n71ZNHynNYX3no4w1xuULSAh4bFyIjTmhvdRsUjA1YBgkNQrW1BqlMgnPnk9s2X/AQAAAMCVHTx40AQH169f79CueQcHDRpE3kE3R5AQQLb5YNle81MDhPY/Nx8NTXH7IoHeJn+gdYZgnVIFJMjfK5v2FgAAAAByT97BGTNmyHfffeeQd7Bu3boycuRI8g7CIEgIIEsdOBspS3eelmX/nr5ugRE/L0+5RWcH2gUFtdiIVtMCAAAAAGQs7+CCBQtMgDA8PNzWXrx4cZN3sHXr1nzngg1BQgCZnt/i35PhJiiowcH9ZyKvu72GACsUCZDlw+6SvJ4e/DUAAAAAIBOsXbs2xbyDffr0MXkHfX0p7ghHbhsknDx5srnYT7MFkDGJiRbZfPSiCQou/fe0HL+Ycm7BMgX95NjFGFtOQp0kaLGIPN++OgFCAAAAAMikvIMTJkyQP//806H9vvvuM3kHixYtSj8jRW4bJNQ3hl50um1QUJCzdwdwOfEJifLXwfMmMLh8V4ipTJyUBgEbli0o7WsXl3a1ikuZQv6meMnEFfvNMuRKRQNlaOuq5nYAAAAAQMaFhYWZZcW6vNh+QlSdOnVkxIgRUqtWLboX1+W2QUIA6XcpPkHW7DtrZguu3H1GwmLik22T1yOPNK1U2AQF29YMlmL5Haewt69dwrSHhIRIcHCweHiwxBgAAAAAsiLv4DPPPCNt2rQh7yDShCAhgOuKuBQvq/aeNTMAf997VqLjki/R98nrIXdVLSrtaxWXVjWKSQF/b3oVAAAAALKYLinWvIOHDh2ytWmuQc07+Mgjj5B3EOlCkBBAMhei4mTFrhAzY3Dt/nMSl5CYbJtAn7zSsnoxs1S4edWiEuDD4QQAAAAAssPhw4dl3LhxyfIO3nvvvSa1WrFixfhDIN34Vg/AOBUWI8v/DTE5BjccOi+JWlkkiYL+XtKmZrB0qF1Cbq9cWHzyetJ7AAAAAJCNeQc/+eQTmT9/frK8g8OHD5fatWvzt0CGESQE3Njhc1FmtqAGBrceC01xm+L5faVdrWBpV7u4NClfiCrEAACb2NhY2bBhgxw5ckSio6NNtcT69etLhQoV6CUAADI57+B3331n8g5qoNBK87wPGTJE2rVrR95B3DSChIAbsVgssud0hAkKLvv3tPl3SsoV9jfLiDXHYN3SBcTDI0+27ysAIOdat26dTJw4URYvXizx8fESFBQkfn5+cuHCBRM4rFixojz11FPSv39/yZcvn7N3FwAAl7Z+/XqTd/DgwYMOeQd79+4tjz76KHkHkWkIEgK5XGKiRbYeD5VlO0+bWYNHzkenuF314vmuBAZrF5dqwfk4CwUASFGnTp1k8+bN0qtXL1m+fLk0atTIBAit9AvMH3/8Id98843JlfTll1+aqooAACD9eQc1OKgn5+x16NBBBg8ebGYRApmJICGQC11OSJS/D10wQUGdMRgSHpvidvXLFjCzBdvVKi7liwRk+34CAFyPJkTX5U5eXl4p3q6zCPWisxt27dolp06dyvZ9BADAlYWHh5u8g99++61D3sFbbrlFRowYQd5BZBmChEAuEXs5Qdb9d84sJf51V4hcjI5Pto2uGr6tYmEzW7BtzeJSPMjXKfsKAHBdTz/9dJq3rVmzprkAAIC05R38/vvvZfr06eQdhFMQJARcWFTsZfl971kzY3DVnjMSGXs52Tbenh7SrEoRM2Owdc1gKRTg7ZR9BQAAAACk7K+//jJpOuzzDvr4+JiZ+Y899hh5B5EtCBICLiY0Ok5W7D5jZgyu2X9W4i4nJtvG39tTWlQrZioSt6hWVPL5prwkDACA9CpYsGCa89ZqIRMAAHD9vINaDEzz+doj7yCcgSAh4ALOhF+SZbtCTPGR9QfPS0KiJdk2QX5e0rpGsFlKfGeVIuLr5emUfQUA5G4TJkyw/fv8+fPy1ltvSbt27aRp06a2CozLli2TV155xYl7CQBAzhYZGSlz5sxJlnewdu3aJu+g5h8EshtBQiCHOnYh2swW1KXEm49eFEvyuKAUzecj7WoFS/taJeTWioXEy9PDGbsKAHAjuuzJqnv37vLGG2+YCotWzzzzjHz88ceyYsUKGTZsmJP2EgCAnJ13cOrUqRIREWFrL1asmAwZMsScePPw4HsdnIMgIZBDWCwW+e9MpCzRwODO07LrVHiK25Uu6Ccdahc3MwbrlykoHlqNBAAAJ9AZg++//36y9vbt28uLL77I3wQAADsbNmwweQcPHDjgkHdQcw7qxc/Pj/6CUxEkBJwcGNxxIsw2Y/Dg2agUt6tSLNAEBjXHYM0S+dOcCwoAgKxUuHBh+eGHH8yyKHvaprcBAACRo0ePyvjx45PlHdRZgzp7sHjx4nQTcgSChEA203yCW05EyN8bz8vyf0PkZNilFLerUzpI2tUqbi6ViwXydwIA5DijR4+WJ554Qn7//Xe59dZbbbMkli5dKp988omzdw8AAKfS5cSffvqpzJs3zywztqpVq5b06dNHmjdvztJi5CgECYFsoBWI/zxwTpb9e9oEBs9HxSXbRicHNi5fSNprYLB2cSlVgKnmAICcTb/g1KhRQz766COTX0np9bVr19qChgAAuBsNCC5atEimTZsmoaGhtvaiRYuaPL46g/Ds2bNO3UcgJQQJgSwSHXdZ1uw7a5YSr9xzRiIuXTtzZOXlmUdur1TE5BfUysRaiAQAAFeiwcCvv/7a2bsBAECO8Pfff8vYsWOvm3cwMTHRqfsIpIYgIZCJwmLi5bc9ISYwuHrfWbkUn/zg7+vlIbeVzS+dGpaVVjWKS5CfF38DAIDL0i9Bs2bNkoMHD8qECRNMdcYlS5ZI2bJlzXIqAADcJe+gfg6uWbPGoZ28g3AlBAmBm3QuMlZ+3XUlMKhLiuMTLMm2yeeb18wU1PyCd1YuLOEXz0lwcDD5JwAALm316tXSoUMHueOOO8yXorfeessECbdt2yYzZ86UBQsWOHsXAQDIUpGRkSbv4Ny5cx3yDtasWdMU9qpbty5/AbgMgoRABpwIjZFlVysSbzx8QSzJ44JSOMBb2tYKlva1S0jTioXFO6+Hadep5eH0OgAgF3jxxRdNYHD48OGSL18+W3vLli3l448/duq+AQCQlRISEmx5By9evGhrL1KkiMk7eM899zApBC6HICGQRgfORprZglp8ZPvxsBS3KRnka4qOaPGRRuULiadHHvoXAJBr7dixQ+bMmZOsXWcTnjt3zin7BABAVtu4caOMGzdO9u/f75B38JFHHpHevXuLv78/fwS4JIKEQCosFov8ezLcBAU1OLj/TGSK21UsEmAKj+jlllJBkkfLFAMA4AYKFCggp06dkgoVKji0b9myRUqVKuW0/QIAICscO3ZMJk6cKL///rtDe5s2beSZZ56REiVK0PFwaQQJATuJiRbZfPSiCQrqUuLjF2NS7J9aJfOb2YIaGKxcLJDAIADALT300EPywgsvyPz5881noabUWLdunYwcOdJUcAQAILfkHdRcu998841D3sEaNWqYvIP16tVz6v4BmYUgIdxefEKibDh4QZb+e0qW/RsiZyNik/WJTg5sWLagCQpq8ZEyhZg+DgDAO++8I4MGDZIyZcqY3EyapF1/9urVS15++WU6CADg0vQz7YcffpCpU6eSdxBugSAh3NKl+AT5Y/85M2Nwxe4QCYuJT7aN5hO8vVJhExRsWzNYiuX3dcq+AgCQU3l7e8snn3wir776qslPqDMt6tevL1WqVHH2rgEAcFM2bdokY8eOdcg7qJ97mnewT58+5B1ErkSQELna0p2nZMKK/XLoXJSUK+wvd1YpKqfDLsmqvWckOi4h2fZagfiuKkXNjMHWNYpJAX9vp+w3AACuRGcS6gUAAFd3/Phxk3dw1apVyfIODhkyREqWLOm0fQOyGkFC5OoAYf/Zm23X94VEmktSAd6e0rJGsMkxeHe1ohLgw9sCAIC06N69uzRp0sTkJbQ3ZswYU/lRcxWm1euvvy6jR492aKtWrZrs2bPH/PvSpUsm79PcuXMlNjZW2rVrJ1OmTJHg4GD+WACAm6az4T/77DOTdzA+/tpKM/IOwp0QDUGu9f7SvaneVtDfS9rUDDYzBm+vVER8vTyzdd8AAMgN1qxZY4J7SXXo0MEs0UqvWrVqyYoVK2zX8+a9NlQdNmyY/PzzzybwGBQUJIMHD5Zu3bqZQikAANxM3sHFixebE08XLlywtRcuXNh81tx7773i4eFBB8MtECRErpOQaJGZaw+aJcYp8fLMIxtfai15PTnQAwBws7MuND9Tss9aLy8JDw9P9/1pULB48eLJ2sPCwkxVyTlz5kjLli1N26xZs8zsjr/++ktuu+22FO9PZxzqxcq6T1qFWS9ZxXrfWfkYroB+oB94LfCeyOnHhc2bN8u4ceNk3759tjb9XNMCXL1795aAgIAs2d+c1g/OQB9kbz+k9f4JEiJXOXwuSkbO3yabjlxM8XatUly5WCABQgAAMsEtt9wi8+bNM4VL7OmSYK10nF6aHF5zPfn6+krTpk3l3XfflbJly8o///xjln61bt3atm316tXNbevXr081SKi/n3QJszpz5ozExMRIVjt79myWP4YroB/oB14LvCdy2nHh9OnTZmmxfobYu+OOO0xREj1hpSfC9JKb+yEnoA+ypx8iIiLStB1BQuQKiYkWmb3hiLz7yx6JiU9IFhi0WK79HNqqqtP2EwCA3OSVV14xS34PHDhgm+G3cuVKk88pPfkI1a233iqff/65yUN46tQpE9y78847ZefOnebLnM7sKFCggMPvaD5CvS01o0aNkuHDhzvMJNQCK8WKFZP8+fNLVp6t18F+0aJF3XqJGv1AP/Ba4D2R044LGvTTz5qkeQf1s0fTWjRo0MAt+iEnoA+ytx/8/PzStB1BQri8E6Ex8vyCbbLuv/O2trKF/OXDB+rKhahYmbhyvxw8GyUViwaYAKHmIQQAADevY8eOsmjRInnnnXdkwYIFZgBap04dk1ewefPm6bovzWNopfehQcNy5crJt99+m+aBbVI+Pj7mkpQOwrPjS1l2PU5ORz/QD7wWeE84+7hgzTs4depUOX/+vEPewYEDB8p9990nnp7Zn6ee4yN9kF2vhbTeN0FCuCyLxSLzNx2XN37aJZGxl23tj95WTl7sUN1Wpbh97RJO3EsAAHI3Teiul8ymswarVq0q//33n7Rp00bi4uIkNDTUYTZhSEhIijkMAQCwzzuoxbT27t3rkDv34Ycflr59+9ryDgIgSAgXFRJ+SV78brus2ntt3X6JIF8Zc38dubNKUafuGwAAyJwlYbqM+dFHH5WGDRuaL3S6lLl79+7mdv2yd/ToUZO7EACApE6cOCEfffSR+eywp+kxhg4dKqVKlaLTgCSYSQiXmz3447aT8uoP/0pYzLUcEg80LC2vdKwp+X29nLp/AAC4E126kkeT/l5neVdajRw50ixf1iXGJ0+elNdee80s/erZs6cEBQXJ448/bvILFipUyOQTHDJkiAkQpla0BACQu+lnzJYtW+TcuXNSpEgRqV+/vvnciIqKklmzZsnXX3/tkHdQZ6ePGDHCnHgCkDKChHAZ5yNj5eVFO2XJzmsJyovm85H3ut0irWoEO3XfAABwRwsXLnS4rl/G9AvbF198kWJV4es5fvy4CQhqrihN3t2sWTP566+/zL/V+PHjTVBSZxLGxsZKu3btZMqUKZn6fAAAruG3336TDz/80FSrt9KiVHfddZesWrXKIe+gnlzSvIN6IsoZeQcBV0KQEC5h6c7T8tLCHXI+Ks7W1rFuSXmjUy0pGODt1H0DAMBdde7cOVnb/fffL7Vq1ZJ58+aZ2X9pNXfu3Ove7uvrK5MnTzYXAIB7Bwiff/75ZO0aMNQiWlaapqJXr14m72BgYGA27yXgmggSIkcLi46X1xf/Kwu3nLC1FfT3kre63CL31qEgCQAAOZEuAX7qqaecvRsAgFy4xFhnEN7I3XffLc8++6yULl06W/YLyC0IEiLHWrX3jClOEhIea2trUzNY3ul6i1lmDAAAcp6YmBiTKJ6E8ACAzKYpLeyXGKfmoYceIkAIZABBQuQ4EZfi5e2fd8vcjcdsbfl888roTrWka/1S102QDgAAsk/BggUdPpe1wFhERIT4+/vL7Nmz+VMAADJVWgKESouZAHDTIGHXrl3l999/l1atWjnkIIDr+fPAOXlu/nY5ERpja7uralF5v/stUiLIz6n7BgAAHGkxEfsgoRYW0UIjt956qwkgAgCQWf755x/55JNP0rStVjsG4KZBwqFDh0q/fv1MJT24ppi4BHl/6R75/M/DtrYAb0956d6a0rNJGWYPAgCQA/Xp08fZuwAAyOWOHz9u0lhowZK0CA4Olvr162f5fgG5Ua4IEmpSUp1JCNf0z5ELMnL+djl0LsrWdmuFQvLhA3WlTCF/p+4bAABwtH379jR3SZ06deg+AECGREZGyqxZs2TOnDkSHx9vay9RooScOnUq1d8bMWKEeHp60uuAKwYJ16xZIx988IGZOqxv9IULF0qXLl0ctpk8ebLZ5vTp01K3bl2ZNGmSNGnSxGn7jMxxKT5Bxq/YJ5+sOSiJlittPnk95IX21aXP7eXFw4PcgwAA5DT16tUzM/w1/6C6Xq5grUIJAEB66GfH4sWLZcqUKXLhwgVbe+HChWXAgAHSsWNHWb16talybJ+jUGcQaoCwZcuWdDjgqkHCqKgoE/jT5cLdunVLdvu8efNk+PDhMm3aNJPfZsKECdKuXTvZu3evFCtWLN2PFxsbay5W4eHh5mdiYqK5ZBXrfWflY+R09n2w40SYmT24/0yk7fb6ZQrIB/ffIhWLBmrqc0m0Rg5zGV4L9AGvBd4PHBc4Pjr7M+Jm7v/QoUMOVSZHjhwpzz33nDRt2tS0rV+/XsaOHStjxozJlH0FALiPTZs2ybhx42Tfvn22Nm9vb3n44YdNiouAgADTpoHA5s2bm88hLVKiOQh1iTEzCAEXDxJ26NDBXFKjB4gnn3xS+vbta65rsPDnn3+Wzz77TF588cV0P967774ro0ePTtauZyBiYq4Vy8gqZ8+eFXd2OcEi7y7eJp//fVoSrsYAvTzzyBO3lpBeDYMlb2KUhIRcW3acm7n7a0HRB/QDrwXeExwbnHN81ArEGVWuXDnbvx944AGTJ+qee+5xWGJcpkwZeeWVV5KtDgEAICXHjh2TiRMnJksj1rp1axkyZIiUKlUq2e9oQLBRo0Z0KJCbgoTXExcXZ5Yhjxo1yqFqnh4o9Cx1Ruh96cxE+5mEOpDVWYn58+eXrDxjrwN+rfinz8Ed7T4ZJsPmbZF9Z68FY2uVzC8f3l9HqhXPJ+6C1wJ9wGuB9wPHBY6Pzv6M8PPzy5T72bFjh1SoUCFZu7bt2rUrUx4DAJC78w7OnDlTvvnmG7l8+bKtvXr16mbpMAVIgOyVo4OEOm1Y8xFobgF7en3Pnj226xo03LZtm1m6XLp0aZk/f75tyUtSPj4+5pKUDsSzI3iXXY+TkyQkWmTGmoMy/te9End1+qCnRx4Z3KKyDG5ZWbw83as/3Pm1kBR9QD/wWuA9wbHBOcfHzLrvGjVqmFUan376qVkOZj3Jq216GwAAKdGAoOYd1JWCFy9etLXrsuFBgwbJvffe6/bflQBnyNFBwrRasWKFs3cBqTh4NlJGzN8mW46G2tqqFguUsT3qyS2lg+g3AABcmH650wTyepLWWslYqx9rMRP98gcAQFJbt26VL774Qv777z9bm07keeSRR6R3797i7+9PpwFOkqODhHoWQfMMhISEOLTr9eLFizttv3BjWnTk8z8Py5hle+RS/JXk6Fr88OEGwfJSp7ri5+NFNwIA4OKaNGkiBw8elK+//tq2yuPBBx+UXr162ZLLAwCgjhw5YgqR/vHHHw4d0rZtW5N3sESJEnQU4GQ5Okioy1YaNmwoK1eutCW+1lw9en3w4MHO3j2k4tiFaHluwTb56+C1cvXlC/vLB/fXkdK+ceLj5UnfAQCQS2gw8KmnnnL2bgAAciitA6BpKebNm2fSiVnVrFnT1AuoV6+eU/cPQA4KEmqiUvtpxocOHTLTjwsVKiRly5Y1Bw2dcqxVi/RstZ550NyD1mrHyDksFot88/cxefvnXRIVd+3g3+f28vJ8+2rim9cj2axQAADg2r766iuZPn26mVGoheW0+vH48eOlYsWK0rlzZ2fvHgDAiXkHFy5caFJThIWF2dr1u/4zzzwj99xzD3kHgRzG6UHCTZs2SYsWLWzXrZWHNTD4+eefmyUrWuXv1VdfldOnT5uzDEuXLk1WzCS9Jk+ebC72ZzKQcafCYuSF73bImn1nbW2lCviZ2YO3Vy5imwUKAAByj6lTp5ox2rPPPitvvfWWbVxVsGBBc2KXICEAuKe//vpLxo0bZ04g2ecdfPTRR6Vdu3bmhJK7F3EEciKnBwnvvvtuMwPtenRpcWYvL9aKSXrRqc9BQRTQyCj9232/+YS8vvhfibh0rWR9zyZl5H/31JB8vuQeBAAgt5o0aZJ88sknJi3Me++9Z2vXFSAjR4506r4BALLf4cOHzUmitWvXOrR36NDBfP8uVqwYq8uAHMzpQUK4rrMRsfK/hTvk113XlhAH5/eR97rXkRbVijl13wAAQNbTNDH169dP1q6zRTQ9DADAPehyYj1pNH/+fIfVerVr15YRI0bILbfcYq6zugzI2QgSIkN+3n5KXl60Qy5Gx9vautYvJa93rCVB/sweBADAHVSoUMHkktZlY/Y0NUyNGjWctl8AgOzLO7hgwQKZMWOGWaVnpenBtGKxVi5mWTHgJkHCY8eOmZ9lypTJrP1BDncxKk5e/fFfWbztpK2tcIC3vN31Fmlfu7hT9w0AAGQvzSWty8cuXbpkUpD8/fff8s0338i7775rKlkCAHKvdevWmUJVusTYytfX19QX0NyD+m8AuTxIqGcKRo8eLR999JGpTKwCAwPNWYLXXntNvLyYRZZbrdwdIi9+v8MsM7bqULu4vNWlthQO9HHqvgEAgOz3xBNPiJ+fn7z88ssSHR0tvXr1kpIlS8rEiRPloYce4k8CALmQFiPRvIN//vmnQ7tWK9YTRzdbZBSACwUJNRj4/fffy5gxY6Rp06ambf369fL666/L+fPnTZU75C7hl+LljcW7ZME/x21tQX5e8kbnWtKpbknJkyePU/cPAABkPz1xPGfOHFOl8uGHHzZBQj2BrEnpAQC5T2hoqEyfPt3EA+zzDtapU8fMLNf8gwDcLEiog8G5c+ea6kT2BwVdctyzZ0+ChLnM2v3n5PkF2+Rk2CVbW4tqRU1xkuD8TB8HAMBd5c2bV/r37y+7d+821/39/c0FAJC7xMfHm4IkWpgkIiLC1l68eHF55plnpE2bNkwcAdw1SKjV6sqXL59i4mpvb29xFZMnTzYX+zMguCYq9rK8u2S3zP7rqK0t0CevvHpfTXmgUWk+BAAAgDRp0kS2bNmSrHAJAMD1aa7ZtWvXmryDR49e+16oaSb69OljZpGTdxBw8yDh4MGD5c0335RZs2aZgKGKjY2Vt99+29zmKjRXgl60AlNQUJCzdydH+fvQBRk5f5scvRBta7u9UmEZc38dKV2QGQIAAOCKgQMHyogRI+T48ePSsGFDCQgIcOgaXW0CAHA9//33nwkObtiwwaG9Y8eO5thftGhRp+0bgBwUJNSzxStXrpTSpUtL3bp1Tdu2bdskLi5OWrVqJd26dbNtq7kK4DouxSfIh8v2ysx1h8RiudLm5+Upo+6pLo/cWk48PMg9CAAArrEWJ9HlZlaaq1hnn+hPVmwAgGu5ePGiLe9gYmKirb1evXom72DNmjWdun8AcliQsECBAtK9e3eHNs1HCNe29ViojPh2qxw4G2Vra1SuoHz4QF0pX8RxVgAAAIA6dOgQHQEAuSTv4Lx58+TTTz81RaistGK9ngjSCUEUrARyv3QHCXWZMXKPuMuJ8tHK/TJ19QFJSLwyfdA7r4eMbFtVHm9WUTyZPQgAAFJBLkIAcG0683v16tUyceJEOXbsmK1dC1H169fPFCe1phkDkPulO0ioLl++LL///rscOHBAevXqJfny5ZOTJ09K/vz5JTAwMPP3Elli18lwGf7tVtlz+lqFqjqlg2TsA3WlSnA+eh0AAAAAcqn9+/fLuHHjZOPGjbY2nS3YqVMnGTBggBQpUsSp+wfABYKER44ckfbt25vqRlqwRMuda5Dw/fffN9enTZuWNXuKTHM5IVGmrT4gE1ful/iEK7MH83rkkaGtqkj/uyuJl6cHvQ0AAAAAudCFCxdk6tSp8sMPPzjkHWzQoIHJO1i9enWn7h8AFwoSDh06VBo1amSKlRQuXNjW3rVrV3nyyScze/+Qyf47EyEjvt0m246H2dqqF88nY3vUlVolqfIMAAAAALmRFhudO3euzJw5U6KiruWiL1WqlPme36JFC/IOAm4u3UHCP/74Q/7880/x9vZ2aC9fvrycOHEiM/cNmUjzDX629pB8sHyvyUOoNN3ggLsryTOtqohPXk/6GwAApMnBgwelYsWK9BYAuEjeQU0XNmHCBIfv7AEBAba8g0m/3wNwT+kOEup05ISEhGTtx48fN8uOkfMcOR8lI+dvk42HL9raKhYNMLkH65ct6NR9AwAArqdOnTrmBLHmrercubPceuutzt4lAEAK9uzZY/IObt682SHvYJcuXaR///4OqwMBIN1BwrZt25ozEDNmzLAdYLRE+muvvSb33HOPy/To5MmTzSWlgGdukZhoka83HJF3ftkjMfFXnmeePCL97qggz7WrJr5ezB4EAADpd+7cOfn1119NPisNEup48L777jNBQ81X7evrS7cCgJOP05p38McffzQzCa00dZjmHaxatSp/HwA3HyQcO3astGvXTmrWrCmXLl0y1Y21KpJWPvrmm2/EVQwaNMhcwsPDJSgo9+XiOxEaIy8s2C5r/ztnaytTyE8+vL+u3FqRs0UAACDjNAjYsWNHc9Evn+vXrzdfRF944QWzbK1169YmYKi3Fy1alK4GgGyixUT1e/lnn30m0dHRtvbSpUvLs88+K82bNyfvIIDMCxLqwUWLlsybN8/81FmEjz/+uDz88MPi5+eX3rtDJtOB+vx/jsubi3dJROxlW/vDt5aV/91TQwJ80v0nBwAASJXOIrz99tvN5b333jMnjzVg+Pnnn8uAAQPMMjc9MQsAyNrvgStXrpSPPvpITp486ZB3UAuM9ujRg7yDAG4o3RGjNWvWmEGgBgX1YnX58mVz21133ZXeu0QmORN+SUZ9v0NW7jljaysR5Cvvd68jd1XlLD4AAMh6VapUkREjRpjL+fPn5cKFC3Q7AGSh3bt3mxMyW7ZssbV5eHhI165dTd7BggXJQw8gbTwknbQsekqDvbCwMHMbnHPW6IetJ6TN+DUOAcL7G5aWpc/eRYAQAAA4hSbE16BhRuisRJ2lqMvjrDTVjc5K1PsNDAyU7t27S0hISCbuMQC4Vt7B0aNHy2OPPeYQIGzSpInMmTNHRo0aRYAQQNbOJNSAlA7YktIzxTqVGdnrfGSsvPLDTvllx2lbW5FAH3m32y3SpmYwfw4AAOByNm7cKNOnTzdVlO0NGzZMfv75Z5k/f77JKT148GDp1q2brFu3zmn7CgDZTU+YfP311yatQ0xMjK29bNmy5jjZrFkz8g4CyNogoQ7AlAYI+/TpIz4+PrbbtELw9u3bzTJkZJ9l/56WlxbukHORcba2++qUkDc715aCAd78KQAAgMvRfNea0uaTTz6Rt956y2HVysyZM83smJYtW5q2WbNmSY0aNeSvv/6S2267LcUE/nqx0oJ1KjEx0VyyivW+s/IxXAH9QD/wWsjc94RO2FmxYoVMmjRJTp++NkkkX7588sQTT8j9998vXl5eZjv7isY5CccF+oHXgnPeE2m9/zQHCa0VgPVgowch+yIl3t7eZmCmCVGR9cKi42X04n/l+y0nbG0F/b3kzS615b46JfkTAAAAl6XLie+9915TIdk+SPjPP/9IfHy8abeqXr26mTmj1ZVTChK+++67ZileUmfOnHGYfZNVzp49m+WP4QroB/qB18LNvyf27dsnn376qck/aJ93sEOHDqaqvH5fd6UcsBwX6AdeC9n7noiIiMjcIKGeqVXly5eXkSNHsrTYSX7fe0Ze+G67hIRfOyveukawvNOtthTL5+us3QIAALhpc+fOlc2bN5vlxknprBk9MV2gQAGH9uDgYIcZNfY0H9fw4cMdZhKWKVNGihUrJvnz58/Ss/U62C9atKj5Eu+u6Af6gdfCzb8n9KTGlClT5JdffnFov/XWW03O1kqVKokr4bhAP/BacM57wn6iX6bmJHz++ecdpi4fOXJEFi5cKDVr1pS2bdum9+6QRpGxl+Xtn3fLN38ftbXl880rr3WsJd0blCLnBAAAyBZaJTOl/NQpSc+slmPHjsnQoUPl119/FV/fzDnxqelx7FPkWOkgPDuCd9n1ODkd/UA/8FpI/3tC8w7Onj3b5B3Uf1uVK1fO5B284447XPo7IMcF+oHXQva+J9J63+kOEnbu3NnkJ9RS6qGhoaZykp7V1cpKWnZ9wIAB4gomT55sLppPMadbf+C8PLdgmxy/eG1ZzJ1Visj73etIyQJpiwYDAABkhgkTJmRJR+pyYp0x06BBA1ubjtPWrFkjH3/8sSxbtkzi4uLM+NN+NqFWNy5evHiW7BMAZDedkKPHO807aF+9XWc/P/XUUybvYN686f4aDwBpku6jiy4BGT9+vPn3ggULzKBMy61/99138uqrr7pMkFDz3ehFl51Y8y3mNDFxCfL+0j3y+Z+HbW3+3p7y0r01pFeTsi595ggAALim3r17Z8n9tmrVSnbs2OHQ1rdvX5N38IUXXjDLhDUh/8qVK6V79+7m9r1798rRo0eladOmWbJPAJCd9BioE2/sj4Wenp7ywAMPmPz/OfV7KwA3DhJGR0ebwiVq+fLlZlahTlvUZNG69BiZ458jF2Xk/G1y6FyUra1JhULy4f11pWxhf7oZAAA4hbVCcFqkJ++fji9r167t0BYQECCFCxe2tT/++OMmx2ChQoXMfQ8ZMsQECFMqWgIArkLzquoqtyVLlji065JizTtYoUIFp+0bAPeS7iBh5cqVZdGiRdK1a1czDVrzIShdHpKVCaDdRezlBBn/636ZseaAJF5N/eiT10Oeb19d+t5eXjw8mD0IAACcR5f63mg1gy6X020yO62LrmbRk9M6kzA2NlbatWtnEvoDgCvSKutffvmluegxzUqDgvo9+/bbb3fq/gFwP+kOEuqS4l69epmDli4LsS7v0FmF9evXz4p9dBs7T4TJ8G+3yr6QSFtbvTIFZGyPulKpaKBT9w0AAECtWrUq2zri999/d7iuBU2seaUBwJWrmeqsQT2W6WQbK11O/PTTT5vVeuQdBOASQUJNlNqsWTM5deqU1K1b19auAUOdXYj0i09IlMmr/pOPf/tPLl+dPujlmUeGtakqT91ZUfJ6UhUPAADkDM2bN3f2LgCAy9q9e7fJs7pr1y6HvIMPPvigPPHEE6zOA+BUGSqLpMVKklaR0yrHSL+9pyNkxPytsvPEtfw+NUvkN7MHa5Rg+TYAAMi5tPLw9dx1113Zti8AkJPpJButWKwr8OzdeeedJu9guXLlnLZvAGBF7XQnSUi0yIw1B2X8r/skLiHRtHl65JFBLSrL4BaVxTsvswcBAEDOdvfddydrs89XmNk5CQHA1Wjhzy+++EJmz57tkHewUqVKJoUXhZcA5CQECZ3g4NlIU7l489FQW1vlYoEyrkddqVO6gDN2CQAAIN0uXrzocD0+Pl62bNkir7zyirz99tv0KAC3zjv4yy+/yMcffyznzp2ztWuxzwEDBphUXeQdBJDTECTMBkt3npIJK/ab4GABf2+5GB0n8QlXcg/qyXbNO6j5B329PLNjdwAAADKFJtlPqk2bNuLt7S3Dhw+Xf/75h54G4Ha2bt0qY8eONfkHrTQgqHkH77vvPjOLUCu1A4DLBgm1qnHnzp2lYcOGWbtHuTBA2H/2ZtGFNxoWPBNxbYp5+cL+8uEDdaVR+UJO3UcAAIDMFBwcLHv37qVTAbiVEydOmLyDK1asSJaa4ZlnnpHSpUtLSEiI0/YPADItSHj8+HHp0KGDOTPcsWNH6dSpk6lorNeROp1BqK7MG7ymoL+X/DL0TvH3ZjInAABwTdu3b3e4brFYTHL+9957T+rVq+e0/QKA7BQVFSWff/65fP311xIXF2drr1KliplV3bhxY9sSZADIydIcofrss8/MQW3dunWyePFiU4FJB4G6pERnGOq06UKFmBGX1MFzUSn2Z3RcAgFCAADg0jQQqIVKNDhoTxPx69gRAHIzLc70008/yZQpU+T8+fO29oIFC8rAgQPNxBpPT1JKAXAd6ZrGpnkTtES7XsaMGWNyLGjAcPr06fLUU09JkyZNzIGwZ8+eUqpUKcnJJk+ebC5ZXXWvYpEA2Xs6wmEmoeYhrFg0IEsfFwAAIKsdOnQo2VixaNGi4uvrS+cDyNU05+q4ceMcUit4eXmZ78L9+vWTwMBAp+4fAGTETWVLrVGjhjz//PNmduGxY8ekd+/e8scff8g333wjOd2gQYNk165dsnHjxix9nGdbVzEBQs1JaA0Q6sn2oa2qZunjAgAAZAVdOWKt1Dl69GhzvVy5cuZSpkwZAoQAcjVNw/Xcc8/J008/7RAgbNGihcyfP9/kHiRACMBVZVpCPD1r/Pjjj5sLrmlfu4RMe6SBTFyxXw6cjZRKRQNlaOuq0r52cboJAAC4HM23FR4eLkWKFJEvvvhC3n//fcmXL5+zdwsAslRkZKRJo6ATYuLj423tVatWNXkHGzVqxF8AgMujakY2BQrb1gw2lay02h/l7gEAgKtq2rSpdOnSRRo2bGhyEeqsGT8/vxS3JS8hAFen6al+/PFHmTp1qly4cMHWXrhwYZN3UHPzk3cQQG5BkBAAAABpNnv2bBk/frwcOHDAFC0JCwuTS5cu0YMAcp1NmzbJ2LFjZf/+/bY2b29vefjhh6VPnz4SEECeeQC5C0FCAAAApJmuinjvvffMvytUqCBfffWVmVEDALmF5tufOHGi/P777w7trVu3liFDhuT4Ip0AkFEECQEAAJAp1Y0BwJVFRETIzJkzZe7cuXL58mVbe/Xq1WXEiBFSv359p+4fAOS46saaoPrnn3+2XdfqxgUKFJDbb79djhw5ktn7BwAAgBxEvzynZzbOunXrsnR/AOBmaUBwwYIF0rVrV5NSwRog1AJNr732mnz55ZcECAG4hXQHCd955x1bcur169fL5MmTZcyYMeYAOmzYsKzYRwAAAOQQmry/Ro0aZvy3e/fuZLdrjsJffvlFevXqJQ0aNJDz5887ZT8BIC02bNhgcgxqGoXQ0FDT5uPjI48//rh8//330rFjRwpPAnAb6V5urGeEK1eubP69aNEi6d69uzz11FNyxx13yN13350V+wgAAIAcYvXq1abS56RJk2TUqFEmcb/mKfT19ZWLFy/K6dOnzcljTeq/c+dOcxsA5DS6Cm7ChAnyxx9/OLS3bdvW5B0sUaKE0/YNAFwmSBgYGGjOCJctW1aWL18uw4cPN+06MIyJicmKfQQAAEAO0qlTJ3M5d+6crF271nzZ1nGgBgc1Z5dePDzSvWAFALJceHi4fPrppzJv3jxJSEiwtdesWdN8t61Xrx5/BQBuK91BwjZt2sgTTzxhBn/79u2Te+65x7T/+++/Ur58+azYRwAAAORAGhTs0qWLs3cDAG5I8wzq8uHp06ebtAhWRYsWlcGDB0uHDh04uQHA7aU7SKg5CF9++WWz7Pi7776TwoULm/Z//vlHevbs6fYdCgAAAADIOTSX/vjx4+XgwYO2Ns07+Nhjj5mLNec+ALi7dAcJtZLxxx9/nKx99OjRmbVPAAAAAADclMOHD5vgYNIq6zprcNCgQVK8eHF6GABuJkioLl26JNu3b5czZ85IYmKirT1Pnjym+hMAAAAAAM6gy4k/+eQTmT9/vkPewdq1a8uIESPklltu4Q8DAJkRJFy6dKk8+uijpnhJUhoktD8IAwAAAACQXXkHFyxYIDNmzDAFSqy0yrpWLNbKxRRVAoBMDBLqwbVHjx7y6quvmoOtq9LcinohqAkAAAAArk2XFOvSYl1ibOXr6yu9e/c2k1z03wCATA4ShoSEmNLwrhwgVJqDQi96hikoKMjZuwMAAAAASCctRqLBQS1OYu/ee+813/eKFStGnwJAVgUJ77//fvn999+lUqVK6f1VAAAAAABuWmhoqEyfPl2+//57h9VhderUMZNaNP8gACCLg4Ra2fiBBx6QP/74wyR89fLycrj9mWeeSe9dAgAAAABwQ/Hx8aYgiRYmiYiIsLVrpWL9LtqmTRuTKx8AkA1Bwm+++UaWL19ucjrojEL7A7D+myAhAAAAACAzWSwWWbt2rVlafPToUVu7n5+f9OnTRx5++GHyDgJAdgcJX3rpJRk9erS8+OKLVIYCAABwcxcvXpSZM2fK7t27zfUaNWpIv379pFChQs7eNQC5xH///WeCgxs2bHBo79ixowwcOFCKFi3qtH0DALcOEsbFxcmDDz5IgBAAAMDNrVmzRjp16iT58+eXRo0ambZJkybJm2++KYsXL5a77rrL2bsIwMVPQljzDiYmJtra69WrJyNGjDAnJQAATgwSagn5efPmyf/+979M3A0AAAC4Gq0c2qNHD5k6dap4enqaNi0goDN79LYdO3Y4excBuGjeQf3O+emnn0pkZKStvWTJkia9VatWrcg7CAA5IUioA78xY8bIsmXLTOWopIVLxo0bl5n7BwAAgBy8BHDBggW2AKHSf2tl0S+//NKp+wbANfMOrl69WiZOnCjHjh2ztfv7+5s0Bj179hQfHx+n7iMA5GbpDhLqGeH69eubf+/cudPhNqpIAQAAuI8GDRqYXITVqlVzaNe2unXrOm2/ALie/fv3mwknGzdudPh+qSkNBgwYIEWKFHHq/gGAO0h3kHDVqlVZsycAAABwKbrsb+jQoWZG4W233Wba/vrrL5k8ebK89957sn37dtu2ugIFAJK6cOGCSVnwww8/OOQd1JMQOiu5evXqdBoA5NQgIQAAAKB06Z96/vnnU7xNZwHp8kH9qSlrALgffe9v3rxZDhw4IJUqVTLBP01LoAUxv/nmG/nss88kKirKtn2pUqXMyYcWLVqwUg0AshlBQgAAAGTIoUOH6DkAqfrtt9/kww8/lDNnztjaihUrJu3atTO3nThxwtYeEBBgyzvo7e1NrwKAExAkBAAAQIaUK1cu03pOlxvq5fDhw+Z6rVq15NVXX5UOHTqY65cuXZIRI0bI3LlzJTY21gQZpkyZIsHBwfz1gBxIg4ApzTLWgOFXX31lu64zjbt06SL9+/eXwoULZ/NeAgDsESQEAABAhtyogvFjjz2W5vsqXbq0yWNYpUoVs0T5iy++kM6dO8uWLVtMwHDYsGHy888/y/z58yUoKEgGDx4s3bp1k3Xr1vHXA3LgEmOdQXgjDRs2NMH/qlWrZst+AQCujyAhAAAAMkTzhtmLj4+X6Ohos1TQ398/XUHCjh07Olx/++23zcxCLYSiAcSZM2fKnDlzpGXLlub2WbNmSY0aNczt1qIpSemMQ71YhYeHm59aHMG+QEJms953Vj6GK6Af3LcfNAeh/RLj1Dz++ONSuXJlt+kbd3wtJEUf0A+8Fpzznkjr/RMkBAAAQIZcvHgxWdv+/ftlwIAB8txzz93ULCSdMajFDJo2bSr//POPCUC2bt3ato1WPC1btqysX78+1SDhu+++K6NHj07WrsGLmJgYyWpnz57N8sdwBfSD+/XD7t2707TdwYMHpUyZMuJu3Om1kBr6gH7gtZC974mIiIg0bUeQEAAAAJlGlwvrsuFHHnlE9uzZk67f3bFjhwkKav7BwMBAWbhwodSsWVO2bt1qZicWKFDAYXvNR3j69OlU72/UqFEyfPhwh5mEGpDQwgn58+eXrDxbr4P9okWLioeHh7gr+sH9+kHfuwsWLDAzf9NCqx27U15Rd3otpIY+oB94LTjnPeHn55em7QgSAgAAIFPlzZtXTp48me7fq1atmgkIhoWFmUBD7969ZfXq1RneDx8fH3NJSgfh2fEFPbseJ6ejH3J/P8TFxcmiRYvks88+k3PnzqXpdzQ42KBBg1zbJ+76Wkgr+oB+4LWQve+JtN43QUIAAABkyI8//uhwXQuOnDp1Sj7++GO544470n1/OltQ85NZCxps3LhRJk6cKA8++KAJQoSGhjrMJgwJCZHixYvz1wOc5PLly/LTTz/Jp59+6jCrVysW16tXzxQeSo0WLPH09MymPQUApIXbBgknT55sLprzBgAAAOnXpUsXh+saGNDlMlpcZOzYsZmyBEcLj2jA0MvLS1auXCndu3c3t+3du1eOHj1qlicDyF76HerXX3+V6dOny7Fjxxxua9GihTz99NMm4P/bb7+ZKsf2RUx0BqEGCK1FiAAAOYfbBgkHDRpkLpqbJigoyNm7AwAA4HIysxKf5g/s0KGDKUaiybW1kvHvv/8uy5YtM2M1rYKq+QULFSpk8gkOGTLEBAhTK1oCIPPpbOFVq1bJtGnTTNERezp7uH///qbquJUGAps3b26qHR84cMDkINQlxswgBICcyW2DhAAAAMg5dKbRY489ZpYra1CwTp06JkDYpk0bc/v48eNNPh2dSaizC9u1aydTpkxx9m4DbhMc/PPPP2Xq1KnJChLpTF+taK7Li1OiAUHdpnTp0mYWobvn4gOAnIwgIQAAANLMvlrwjYwbNy7N296oGqqvr68tXQyA7LNp0yYTkN++fbtD+y233CIDBw6Uxo0b8+cAgFyCICEAAADSLGkhAl1GqMULtDKx2rdvn23mEADXpUFBXVb8999/O7RXrVrVBAd1ebHmIQUA5B4ECQEAAJBmmo/MfqZgvnz55IsvvpCCBQuatosXL0rfvn3lzjvvpFcBF6TLiTU4uHbtWof2ihUrmoIkWpiEJcMAkDsRJAQAAECGaAXj5cuX2wKESv/91ltvSdu2bU0FUwCuQQuRaLVirSJuT3MJPvXUUyYPKAVHACB3I0gIAACADAkPD5ezZ88ma9c2rVAMIOc7duyYzJgxQ5YuXWoKlFhpkZEnn3xS7rvvPsmbl6+NAOAOONoDAAAgQ7p27WqWFuuMwiZNmpi2DRs2yHPPPSfdunWjV4Ec7PTp06Zg0I8//igJCQm29sKFC5v3tb6Hvb29nbqPAIDsRZAQAAAAGaJ5y0aOHCm9evWS+Pj4K4PLvHnl8ccflw8++IBeBXKgc+fOyaxZs+T777+3vW9VUFCQ9O7dW3r06GGqiQMA3A9BQgAAAGSIv7+/TJkyxQQEDxw4YNoqVaokAQEB9CiQw4SGhsqXX34p8+bNk9jYWFu7vl8feeQR6dmzpwQGBjp1HwEAzkWQEAAAADfl1KlT5nLXXXeJn5+fyWuWJ08eehXIASIjI+Xrr7+WOXPmSFRUlK1dZwtqYFADhDqLEAAAgoQAAADIkPPnz5uliatWrTJBwf3790vFihXNcmOtcqy5CgE4R0xMjJk1qLMHtciQleYZ7N69u/Tp08fkHwQAwIogIQAAADJk2LBh4uXlJUePHpUaNWrY2h988EEZPnw4QULACXQp8XfffSeff/65XLhwwdbu6ekpXbp0kX79+pnKxQAAJEWQEAAAABmyfPlyWbZsmZQuXdqhvUqVKnLkyBF6FchGWoREKxVrxeIzZ87Y2j08POSee+6RJ554Itl7FQAAewQJAQAAkCGa30yLlySls5d8fHzoVSAbJCQkyJIlS+STTz6REydOONzWpk0befrpp6V8+fL8LQAAN0SQEAAAABly5513mnxnb775prmueQkTExNlzJgx0qJFC3oVyEL6Xlu5cqVMnz5dDh8+7HCbFhHq37+/VK1alb8BACDNCBICAAAgQzQY2KpVK9m0aZPExcXJ888/L//++6+ZSbhu3Tp6FcgCWj18zZo1Mm3aNFMsyN5tt91mgoO1a9em7wEA6UaQEAAAABmigYh9+/bJpEmTJF++fBIZGSndunWTQYMGSYkSJehVIJODgxs2bJCpU6eaYLy9+vXry4ABA6RBgwb0OQAgwwgSAgAAIMOCgoLk5ZdfpgeBLLRlyxYTHNy8ebNDe82aNWXgwIFy6623muX+AADcDIKEAAAAyLA//vjD5EQ7ePCgzJ8/X0qVKiVfffWVVKhQQZo1a0bPAjdBZwzqsuL169c7tFeuXNksK27evDnBQQBApiFICAAAgAz57rvv5NFHH5WHH37YzHCKjY017WFhYfLOO+/IL7/8Qs8CGaC5BjU4uHr1aof2smXLmuBg69atxcPDg74FAGQqgoQAAADIkLfeessEMh577DGZO3eurf2OO+4wtwFIH61SPGPGDPn1119NDkKrkiVLypNPPikdOnSQvHn5CgcAyBp8wgAAACBD9u7dK3fddVeKeQpDQ0PpVSCNTpw4IZ9++qn8/PPPkpiYaGsvWrSoPP7449K5c2fx8vKiPwEAWYogIQAAADKkePHi8t9//0n58uUd2teuXSsVK1akV4EbOHPmjHz22WeyaNEiuXz5sq29YMGC0qdPH+nevbv4+vrSjwCAbEGQEAAAABmiyx+HDh1qghxaWfXkyZOmwMLIkSPllVdeoVeBVFy4cEE+//xzWbBggcTFxdna8+XLZ5bvP/jgg+Lv70//AQCyFUFCAAAAZMiLL75olka2atVKoqOjzdJjHx8fEyQcMmQIvQokER4ebqp/aw7PmJgYW7sGBHv16mWKAGmgEAAAZyBICAAAgAzR2YMvvfSSPPfcc2bZcWRkpNSsWVMCAwPpUcCOBtFnzpwpX3/9tXmfWGlQvUePHtK7d28pUKAAfQYAcCq3DRJOnjzZXBISEpy9KwAAAC7N29vbzH7SCwFC4JpLly7Jt99+K7NmzZKIiAhbu1Yo7tatm/Tr10+KFClClwEAcgQPcVODBg2SXbt2ycaNG529KwAAAC5JCy1o7kGtZqzFS/Si/3755ZclPj7e2bsHOI3mGZw3b56pSvzRRx/ZAoSenp7SpUsXWbhwoTz//PMECAEAOYrbziQEAADAzdG8g99//72MGTNGmjZtatq0cMnrr78u58+fl6lTp9LFcLvA+U8//SSffvqpnD592mFpfrt27eTpp5+WMmXKOHUfAQBIDUFCAAAAZMicOXNMAYYOHTrY2urUqWOCID179iRICLehKYyWL18uM2bMkGPHjjnc1rJlS7O0uEmTJuLh4bYLuQAALoAgIQAAADJEiy7oEuOkKlSoYPIUArmdxWKRVatWybRp0+TgwYMOt91xxx0yYMAAqVq1qoSEhDhtHwEASCuChAAAAMiQwYMHy5tvvmmKMmjAUMXGxsrbb79tbgNyc3Bw3bp1Jji4Z88eh9saN24s/fv3l7p165rriYmJTtpLAADShyAhAAAAMmTLli2ycuVKKV26tC0gsm3bNlO0oVWrVmaJpZXmLgRyg02bNsmUKVNk+/btDu261F5nDmqQEAAAV0SQEAAAABlSoEAB6d69u0MbRRmQW2lQUIvxbNy40aG9WrVqJjioy4u1QAkAAK6KICEAAAAyRJcZA7mdLifWZcVr1651aK9YsaJZVtyiRQuCgwCAXIEgIQAAADIkJibG5Gbz9/c3148cOSILFy6UmjVrStu2belVuDQtRKLBwd9++82hXZfXP/300+Y17unp6bT9AwAgsxEkBAAAQIZ07tzZ5B3U2VShoaHSpEkTU9X43LlzMm7cOLMEE3A1x44dkxkzZsjSpUtNENyqePHi8sQTT8h9990nefPyNQoAkPvw6QYAAIAM2bx5s4wfP978e8GCBSaIosVMvvvuO3n11VcJEsKlnD59Wj799FNZvHixJCQk2NoLFy4s/fr1k65du5ogOAAAuRVBQgAAAGRIdHS05MuXz/x7+fLlZlahh4eH3HbbbWbpMeAKdOar5tfUCtzx8fG29qCgIOnTp4888MAD4uvr69R9BAAgO3hky6MAAAAg16lcubIsWrTILM9ctmyZLQ/hmTNnJH/+/Om6r3fffVcaN25sgo7FihWTLl26yN69ex22uXTpkgwaNMjM7AoMDDSVlUNCQjL1OcF96BL5iRMnmmXz8+bNswUIAwICzBL6H3/8UR599FEChAAAt0GQEAAAABmiS4pHjhwp5cuXl1tvvVWaNm1qm1VYv379dN3X6tWrTQDwr7/+kl9//dUEbDToGBUVZdtm2LBhZino/PnzzfYnT540sxeB9IiMjJTp06eb4OBXX30lsbGxpt3Pz0/69u1rgoOae1CDhQAAuBOWGwMAACBD7r//fmnWrJmcOnVK6tata2tv1aqVyd+WHlokwt7nn39uZhT+888/ctddd0lYWJjMnDlT5syZIy1btjTb6BLRGjVqmMCiLnFOSoM/1gCQCg8PNz8TExPNJatY7zsrH8MV5LR+0GrcOmNw9uzZtteC0jyDOiu1d+/eUqhQoUzf55zWD85AH9APvBZ4T3BccO7xMa33T5AQAAAAGabFSvRiT6sc3ywNCipr0EaDhTq7sHXr1rZtqlevLmXLlpX169enGCTUJcyjR49O1q7LoTVglNXOnj2b5Y/hCrKrH7TYyK5du+TChQvmdVOzZk3x9PSUuLg4WbJkiSmuo0uMrfQ2na364IMPmiXs+vrKyuXrvB7oA14LvB84LnBsdNZnRERERJq2I0gIAACAHEXPdj/77LNyxx13SO3atW2VZ3XGV4ECBRy2DQ4ONrelZNSoUTJ8+HDbdZ09VqZMGTNDMb05E9O7/zrYL1q0qCnk4q6ysx9WrVolY8eONQFgK31cfQ39+eefDu26L/fcc488/vjjUqpUKclqvB7oA14LvB84LnBsdPZnhKbUSAuChAAAAMhRNDfhzp07Ze3atTd1Pz4+PuaSlA7CsyN4l12Pk9NldT/89ttv8sILLyRr1y9dWljHns4cfOqpp0wezezG64E+4LXA+4HjAsdGZ31GpPW+CRICAAAgxxg8eLD89NNPsmbNGildurStXZc067JRXS5qP5tQl4cmXe4M96FLjD/88MMbbqd5LQcMGCBVqlTJlv0CAMAVcWoTAAAATmexWEyAcOHChWZmWIUKFRxub9iwoXh5ecnKlSttbXv37pWjR4/aqirD/WzZssVhKXFqevXqRYAQAIAbYCYhAAAAcsQSY61c/MMPP0i+fPlseQaDgoJMHh39qTnkNMegFqXQnIJDhgwxAcKUipbAPRw7dixN2507dy7L9wUAAFdHkBAAAABON3XqVPPz7rvvdmifNWuW9OnTx/x7/PjxJqdO9+7dJTY2Vtq1aydTpkxxyv7C+cuMddbppEmT0rR9kSJFsnyfAABwdQQJAQAAkCOWG9+Ir6+vTJ482VzgvjZt2mQqGe/fvz9N22sF7Pr162f5fgEA4OoIEgIAAADI8U6dOiUTJkxwyEupNACouQlTM2LECPH09MyGPQQAwLURJAQAAACQY126dEm++OIL+fLLL80yc6saNWrIyJEjpW7duqbYjVY5ti9iojMINUDYsmVLJ+05AACuhSAhAAAAgBy5BP3XX3+ViRMnSkhIiK1dC9dooZuOHTuaHJVKA4HNmzc3Mwq1SInmINQZhswgBAAg7QgSAgAAAMhR9uzZY/IO2i8j1oBfz5495YknnpDAwMBkv6O3N2rUKJv3FACA3IMgIQAAAIAc4eLFi6Zi9aJFixyK2dx+++0yfPhwKV++vFP3DwCA3IwgIQAAAACnunz5ssyfP19mzJghERERtvayZcua4GCzZs2cun8AALgDgoQAAAAAnOavv/6ScePGycGDB21tAQEBZlnxQw89JF5eXvx1AADIBgQJAQAAAGS748ePy/jx42X16tUO7VqQRAuTaPERAACQfQgSAgAAAMg20dHR8tlnn8nXX38t8fHxtvZbbrlFRo4cKbVq1eKvAQCAExAkBAAAAJDlEhMTZcmSJTJp0iQ5d+6crV1nDA4ZMkQ6dOggHh4e/CUAAHASgoQAAAAAstS///4rH374oezYscPWprkGH3nkEenbt6/4+/vzFwAAwMkIEgIAAADIEhcvXjQVixcvXuzQ3rx5cxk2bJiULl2angcAIIcgSAgAAAAgU2muwTlz5sinn34qMTExtvYKFSrIiBEj5LbbbqPHAQDIYQgSAgAAAMg0a9eulXHjxsnRo0dtbYGBgfL000/LAw88IHnz8hUEAICciE9oAAAAADft8OHDJjj4559/2try5MkjXbp0kYEDB0rBggXpZQAAcjCChAAAAAAyLDIy0iwr/uabbyQhIcHWXq9ePVOUpGnTplQtBgDABRAkBAAAAJBuiYmJ8uOPP8qUKVPkwoULtvbg4GAZOnSotGrVSs6cOUPPAgDgIggSAgAAAEiXbdu2yYcffii7d++2tfn4+Mhjjz0mvXv3Fl9fXxNEBAAAroMgIQAAAIA00ZmBkyZNkiVLlji0t27d2sweLFGiBD0JAICL8pBc4KeffpJq1apJlSpVTD4UAAAAAJknNjZWZs6cKd26dXMIEOr4e9q0afLee+8RIAQAwMW5/EzCy5cvy/Dhw2XVqlUSFBQkDRs2lK5du0rhwoWdvWsAAACAS7NYLPL777/LhAkT5MSJE7Z2HXcPGDDAVC7Om9flv1IAAIDcECT8+++/pVatWlKqVClzvUOHDrJ8+XLp2bOns3cNAAAAcFkHDhyQsWPHmvG2laenp3Tv3l2efvppEygEAAC5h9OXG69Zs0Y6duwoJUuWlDx58siiRYuSbTN58mQpX768SYB86623OgxUTp48aQsQKv23/VlOAAAAAGkXHh4uH3zwgfTq1cth3N24cWP5+uuv5fnnnydACABALuT0mYRRUVFSt25d6devn8lxktS8efPMcmLNdaIBQl3q0K5dO9m7d68UK1YsQ/lU9GI/CFJafS0rK7BZ79udq7zRB/QDrwXeExwXOD7yOZEzPivdeTyC1CUkJMjChQtl6tSpEhYWZmvXk/nDhg2Tu+++25zUBwAAuZPTg4S6PFgvqRk3bpw8+eST0rdvX3Ndg4U///yzfPbZZ/Liiy+aQYv9zEH9d5MmTVK9v3fffVdGjx6dYqW2mJgYyWpnz54Vd0cf0A+8FnhPcFzg+MjnhHM/KyMiIrL0/uF6/vnnH/nwww9l//79tjZdxaNj8EceeUR8fHycun8AAMANgoTXExcXZwYso0aNsrV5eHhI69atZf369ea6BgR37txpgoOaF0Wrrb3yyiup3qfel85MtJ9JWKZMGTMrMX/+/Fl6xl4H/EWLFjXPwR3RB/QDrwXeExwXOD7yOZEzPiv9/Pyy7L7hWk6dOiUTJ06UFStWOLS3b99ehgwZIsHBwU7bNwAAkL1ydJDw3LlzZtlD0sGJXt+zZ4/5t1ZT04TKLVq0MANrzZFyvcrGehY0pTOhOhDPjuBddj1OTkYf0A+8FnhPcFzg+MjnhHM/K919LAKRS5cuyRdffCFffvmlQyqe6tWry8iRI6VevXp0EwAAbiZHBwnTqlOnTuYCAAAAIHUWi8XMGtQ83yEhIbb2ggULyqBBg0xBQa1gDAAA3E+ODhIWKVLEDFLsBzBKrxcvXtxp+wUAAAC4Gi38pytwNm/ebGvTsfZDDz1kcoAHBgY6df8AAIBz5eggobe3tzRs2FBWrlwpXbp0MW26pFivDx482Nm7BwAAAOR4oaGhMmXKFFm0aJFDZevbb7/d5OouX768U/cPAADkDE4PEkZGRsp///1nu37o0CHZunWrFCpUSMqWLWsGLr1795ZGjRqZIiW6NCIqKspW7RgAAABAcpcvX5YFCxbI9OnTHSpaa9E+HWM3a9ZM8uTJQ9cBAICcESTctGmTKTpiZa08rIHBzz//XB588EFT5e/VV1+V06dPmyTKS5cuvelKa5MnTzYXLYwCAAAA5CYbNmwwS4sPHjxoa/P395fHH39cevbsaVbsAAAA5Kgg4d13320SKF+PLi3O7OXFmphZL+Hh4RIUFJSp9w0AAAA4w/Hjx2X8+PGyevVqh3YtSKJjX835DQAAkCODhAAAAABuTnR0tMyaNUtmz54t8fHxtvbatWvLyJEjzU8AAIDrIUgIAAAAuChdkbNkyRKZNGmSSdFjVbhwYRkyZIjcc8894uHh4dR9BAAAroEgIQAAAOCCdu3aJR9++KFs377d1ubl5SW9evWSfv36SUBAgFP3DwAAuBZOKwIAAMDp1qxZY/LmlSxZ0lTcXbRoUbIZc1rIrkSJEuLn5yetW7eW/fv3izs6f/68vPHGG6bQn32A8K677pJvv/3WzCAkQAgAANKLICEAAACcLioqSurWrSuTJ09O8fYxY8bIRx99JNOmTTOVezUI1q5dO7l06ZK4C801+NVXX0m3bt3kxx9/tBX/q1ChglluPG7cOClTpoyzdxMAALgolhsDAADA6Tp06GAuKdFg2IQJE+Tll1+Wzp07m7Yvv/xSgoODzYzDhx56SHK7tWvXmiDg0aNHbW2BgYHy1FNPSY8ePSRvXob1AADg5jCaAAAAQI526NAhOX36tFlibBUUFCS33nqrrF+/PtUgYWxsrLlYhYeHm5+JiYnmklWs950Zj3HkyBETIF23bp2tTZdja7B0wIABUrBgwUx7rJzcD66MfqAPeC3wfuC4wLHR2Z8Rab1/tw0S6lIWvSQkJDh7VwAAAHAdGiBUOnPQnl633paSd999V0aPHp2s/cyZMxITE5PlfW5fbTi9oqOjZe7cubJ48WK5fPmyrb1mzZry5JNPSuXKlSUuLk5CQkIkp7uZfshN6Af6gNcC7weOCxwbnfUZERERkabt3DZIOGjQIHPRM8p6JhoAAAC5y6hRo2T48OG26zru05x9xYoVk/z582fp2Xod7BctWlQ8PDzS/bs//fSTTJkyRS5cuGBr131+5plnpE2bNmYmoSu4mX7ITegH+oDXAu8HjgscG539GaFF39LCbYOEAAAAcA3Fixc3P3XWnFY3ttLr9erVS/X3fHx8zCUpHYRnR9AqvY+jlYo//PBD2bVrl61N9//RRx81lYzTOsDPabKrv3M6+oE+4LXA+4HjAsdGZ31GpPW+CRICAAAgR9PqvRooXLlypS0oqLMCtcqx5uVzdTqDQCs3L1myxKG9ZcuW8uyzz0rJkiWdtm8AAMB9ECQEAACA00VGRsp///3nUKxk69atUqhQISlbtqwJlr311ltSpUoVEzR85ZVXTPCsS5cu4qq0qMqcOXPks88+c8iRWKlSJXnuueekUaNGTt0/AADgXggSAgAAwOk2bdokLVq0sF235hLUZbaff/65PP/88xIVFSVPPfWUhIaGSrNmzWTp0qXi6+srrsZiscjq1atl/PjxcuLECVu75snu37+/dO3aVfLmZZgOAACyF6MPAAAAON3dd99tgmep0WIdb7zxhrm4soMHD8rYsWPNUmn7PEHdu3c3AUIK6gEAAGchSAgAAABkkoSEBNm8ebMcOHDALBtu0KCBeHp6mhyKM2bMkPnz55ttrHRJ8ciRI6Vy5cr8DQAAgFMRJASQrVYcWSFTtk6Rw+GHpXz+8jKw3kBpXa41fwUAgMv77bffTHXiM2fO2NqKFStmlkZr0ZWwsDBbu1Zp1jyLWpxEZ0kCAAA4m9sGCSdPnmwu9mdyswpBEefRvp+6baocDjss5YPKy4C6A9w2IJVoSTQXXcqVYEm48m+x2NqTXvQ223ZXf0d/6nXz76u/m/T+EhKv3Wb/O3rZcmaLzNw507ZP+0P3y7Dfh0nfWn2lbrG6kufqfx55PGxfmExLnjwp/7zebWnYRpnHsm/LI+IhVx7f+t+V/+Vx2Na6byntq/X3TZvd/ei29m3aN/GJ8ebiKZ7J9hHpw/sdro7xgusHCDVvYlIaMPz+++9t1zWHYp8+feSRRx5xyXyKAAAg98pjuV7yFzegSz8094ue2c2fP3+WDPg1CJJUqcBSkt87f4qBjKQBCZVS4CPZdlev32i7FIMvad0uleBKWrZT0dHREhAQkDwwY90uhf1JSwAopcffe2GvLPxvYbK+71Spk1QpUMUh0JU0mGUu+l/ilZ8OgTDdTlIJrFmDade7X0uiXIq9JHm98iYP0qV2v9cJvl0vmGd/n3A9GQmImiCn9f2VUkDT/r2VNKBpbUvp2GN9nOvtz3W2ud57PD4uXnx8fJIFVJMeB1ILyOo2IVEhsvnM5mR92Ci4kZQMLCk5nkUk5lKM+Pn6meeTHtb+zW5ZEcjWY9elmEvi6+eb6v074/lmR9D+ROQJ2XDqWo46q/F3j8+Sk1tZPf7JqbLqeetJ544dOzrMIExJ27ZtZejQoRIcHCzuQMdRISEh5vlq3kV3RT/QB7wWeD9wXODY6OzPiLSOgdx2JmF20VlsqX0Z0P+Q/X488CPdDpegAWCH8zhufUonYzaFbBIJcfZeABmjAdFp26a57Qx4V7Jly5YbBghVt27d3CZACAAAXA9Bwiymy1xT7fw8V2aSWQMB+h/cl86Sss7q8szjaZvRZbvItX9bt7HOtnLYzu6it9m2S+l+r96nw23iIZ4enmm63+vtY0r3O3fPXDl36Vyy517Er4g8XONh82/re8E6S/LK/669T6yzIu3fN8l+Xu+2JNso62PZvw+TtV3dD9t+3WBfk7Wl8vixcbFmVqn9/dm2u96+prL/1r6xzixN674mfXzTlqTv09rXADKHvp8OhR2iO13AuXPnMnU7AAAAZyBImMU0D97+i/sdvjhrgKVqwaqyoNOCFH/H/gu3LhVN7Yu62dZuO1tbCl/ere0pBTiuFwhJMciRSsDiRtvpNNrzF85LwUIFTdAoLQGdpP1h/1jW55za40/4Z4Kcjj7t0Lfa9yUCSsjIxiNTDLpZg17m3zcIktm2S2sw7+q/LYkWOXf2nBQPLi55PfMmW7adm1UqUMksv9fnq38j68+Xb31ZWpVrJe4mty4/Su0963BsuXps0GX02gdFihWxvQfsj2kpBTkdHuNq26AVg0wxnKTHWi2O81HLj7L+Od9kcNQcF86fkyKFi0gej7QfC5wWlM2ih9W/5flz56VwkcK2HJ7Ofr7ZlZXl2d+flaPhR5O9hisEVciWx8fNKVKkSKZuBwAA4AwECbOYFspIKSii7amxz02oxQxyVUBEsi8g4uPpk2LfP9/4eacGpLQfdN+8Pb1zVWAoLXTJnObX0mX4OjtGv/wOrDvQLQOEuZn9MSwt74cArwCTo/Vm3g/PNHgmxff70AZDzcmanE77we+SnwQH5a6AcXqZ18OlALfsh2cbPJvu8QJyjvr165sqxtdbcqzjH90OAAAgp3KvEbgTgyJVClYRLw8v83PC3RMIimRj3+usTW8Pb/OTvs8Zf5f5982Xn1v/bH4SIERmva54v8OVMV5wbZ6enjJy5MjrbjNixAizHQAAQE7FTMJsGvi3LNMyVy4rdIW+J+E74B54v8PVMV5wbS1btpQxY8bIhx9+6DCjUMd+GiDU2wEAAHIygoQAAABAJtBAYPPmzWXz5s1y4MABqVSpkjRo0IAZhAAAwCUQJAQAAAAyiS4pbtiwoZQuXZoVJAAAwKWw7hUAAAAAAABwc24bJJw8ebLUrFlTGjdu7OxdAQAAAAAAAJzKbYOEgwYNkl27dsnGjRudvSsAAAAAAACAU7ltkBAAAAAAAADAFQQJAQAAAAAAADdHkBAAAAAAAABwcwQJAQAAAAAAADeXV9ycxWIxP8PDw7P0cRITEyUiIkL8/PzEw8M9Y7P0Af3Aa4H3BMcFjo98TuSMz0rruMc6DnIXjPuyF2M/+oHXAu8JjgscH/mMcK2xn9sHCfWPocqUKZNlfwwAAICcOg4KCgoSd8G4DwAAuLOIG4z98ljc7RRyClHbkydPSr58+SRPnjxZGrXVQOSxY8ckf/784o7oA/qB1wLvCY4LHB/5nMgZn5U6/NNBYsmSJd1qhQPjvuzF2I9+4LXAe4LjAsdHPiNca+zn9jMJtXNKly4t2UX/6O4aJLSiD+gHXgu8JzgucHzkc8L5n5XuNIPQinGfczD2ox94LfCe4LjA8ZHPCNcY+7nPqWMAAAAAAAAAKSJICAAAAAAAALg5goTZxMfHR1577TXz013RB/QDrwXeExwXOD7yOcFnpTtgzEM/8HrgPcGxgeMjnxN8VrrimMHtC5cAAAAAAAAA7o6ZhAAAAAAAAICbI0gIAAAAAAAAuDmChAAAAAAAAICbI0gIAAAAAAAAuDmChDfh3XfflcaNG0u+fPmkWLFi0qVLF9m7d6/DNpcuXZJBgwZJ4cKFJTAwULp37y4hISEO2xw9elTuvfde8ff3N/fz3HPPyeXLl8UVvffee5InTx559tln3a4PTpw4IY888oh5nn5+fnLLLbfIpk2bbLdbLBZ59dVXpUSJEub21q1by/79+x3u48KFC/Lwww9L/vz5pUCBAvL4449LZGSkuIKEhAR55ZVXpEKFCub5VapUSd58803zvHNzH6xZs0Y6duwoJUuWNK/9RYsWOdyeWc95+/btcuedd4qvr6+UKVNGxowZI67QB/Hx8fLCCy+Y90NAQIDZ5rHHHpOTJ0/mqj5Iy2vBXv/+/c02EyZMcMt+2L17t3Tq1EmCgoLM60I/S/VzILd8btyoD/RvOnjwYCldurQ5LtSsWVOmTZvmsI2r90FuxLgvZe469nP3cZ9i7MfYj7EfYz/Gfblw7GdBhrVr184ya9Ysy86dOy1bt2613HPPPZayZctaIiMjbdv079/fUqZMGcvKlSstmzZtstx2222W22+/3Xb75cuXLbVr17a0bt3asmXLFssvv/xiKVKkiGXUqFEu95f5+++/LeXLl7fUqVPHMnToULfqgwsXLljKlStn6dOnj2XDhg2WgwcPWpYtW2b577//bNu89957lqCgIMuiRYss27Zts3Tq1MlSoUIFS0xMjG2b9u3bW+rWrWv566+/LH/88YelcuXKlp49e1pcwdtvv20pXLiw5aeffrIcOnTIMn/+fEtgYKBl4sSJuboP9PX60ksvWb7//nuNhloWLlzocHtmPOewsDBLcHCw5eGHHzbHm2+++cbi5+dnmT59uiWn90FoaKh5b8+bN8+yZ88ey/r16y1NmjSxNGzY0OE+XL0P0vJasNLb9bmWLFnSMn78eLfrBz0uFipUyPLcc89ZNm/ebK7/8MMPlpCQkFzzuXGjPnjyySctlSpVsqxatcocL/Xv5+npafoht/RBbsS4Lzl3Hfsx7ruCsR9jP8Z+jP0Y9+W+sR9Bwkx05swZ84JYvXq17cuxl5eXCZZY7d6922yjX5SV/mE9PDwsp0+ftm0zdepUS/78+S2xsbEWVxEREWGpUqWK5ddff7U0b97cNlB0lz544YUXLM2aNUv19sTEREvx4sUtH3zwga1N+8bHx8d8yVe7du0y/bJx40bbNkuWLLHkyZPHcuLECUtOd++991r69evn0NatWzcTzHCXPkj6gZBZz3nKlCmWggULOrwf9DVXrVo1S05zveCY/ZdK3e7IkSO5sg+u1w/Hjx+3lCpVygT49MSCfZDQXfrhwQcftDzyyCOp/k5u+9xIqQ9q1apleeONNxzaGjRoYAaXubEPcit3Hve5+9iPcd8VjP0Y+ynGflcw9mPcl1vGfiw3zkRhYWHmZ6FChczPf/75xyy10+UFVtWrV5eyZcvK+vXrzXX9qcsTgoODbdu0a9dOwsPD5d9//xVXodNiddqr/XN1pz748ccfpVGjRvLAAw+Yab/169eXTz75xHb7oUOH5PTp0w79oMvsbr31Vod+0KUmej9Wur2Hh4ds2LBBcrrbb79dVq5cKfv27TPXt23bJmvXrpUOHTq4TR8klVnPWbe56667xNvb2+E9oukNLl68KK54rNRp+Pq83akPEhMT5dFHHzXLAmrVqpXsdnfoB+2Dn3/+WapWrWr2W4+X+n6wX5LhDp8berzUzw1drqhjyVWrVpljZ9u2bd2mD3IDdx73ufvYj3HfFYz9kmPslzLGfu459mPc55pjP4KEmfgG0Fwsd9xxh9SuXdu0aXBA38zWL8JW+kfV26zb2P+Rrbdbb3MFc+fOlc2bN5tcPUm5Sx8cPHhQpk6dKlWqVJFly5bJgAED5JlnnpEvvvjC4Xmk9Dzt+0G/MNvLmzev+fLhCv3w4osvykMPPWQOZl5eXiZQqu8JzbXjLn2QVGY959zwHrHPtaE5Cnv27GlyMLlTH7z//vvmeemxISXu0A9nzpwxOVk0h1n79u1l+fLl0rVrV+nWrZusXr3abT43Jk2aZHLRaF4afa7aF5MnTzZfAtylD1ydO4/7lLuP/Rj3XcHYLznGfskx9nPfsR/jPtcc++XNtHtyc3o2defOnWbmlDs5duyYDB06VH799VeTRN+dvyzoGaB33nnHXNcAmb4eNBlp7969xR18++238vXXX8ucOXPMmbKtW7eaL1CavNVd+gDXp2fHevToYc6eaVDdnejZwYkTJ5ov1TqL0p2Plapz584ybNgw8+969erJn3/+aY6XzZs3F3egA8W//vrLnFEuV66cSXat4wg9XiadlYWcyV3HfYqxH+M+K8Z+uBHGfu499mPc55pjP2YSZgKtUvPTTz+ZKaMaGbYqXry4xMXFSWhoqMP2WqFGb7Nuk7RijfW6dZuc/sVXzxA0aNDAnPXQi84G+eijj8y/NbKd2/tAaeU6PTNgr0aNGrZqndbnkdLztO8H7Ut7WqlIK9+5Qj/oEkrrGWWdBq3LKjUIYJ1l4A59kFRmPefc8B6xDhKPHDliTipYZxG6Sx/88ccf5jnqkgHrsVL7YsSIEVK+fHm36YciRYqY536j42Vu/tyIiYmR//3vfzJu3DhTBa9OnTpmHPHggw/Khx9+6BZ94OrcedynGPsx7rNi7JccY79rGPsx9mPc55pjP4KEN0Fnw+gfd+HChfLbb79JhQoVHG5v2LChWXapedqsNHeAfhFq2rSpua4/d+zY4fDF0PoFOumXqJyoVatWZv911pj1ojPqdImp9d+5vQ+ULjfS52VPcwzoWQKlrw1949r3g+YO0FwT9v2gBwUdfFvp60rPwGjOrpwuOjra5M+w5+npaTuD5A59kFRmPWfdRs826WDL/j1SrVo1KViwoLjKIHH//v2yYsUKKVy4sMPt7tAHGjTfvn27w7FSzxzqFyxNUeAu/aDLKBo3bnzd42Vu/+zUv51erne8zO194KoY913B2I9xnxVjv+QY+13B2I+xn2Lc56Jjv0wrgeKGBgwYYAkKCrL8/vvvllOnTtku0dHRDmWsy5Yta/ntt99MGeumTZuaS9Iy1m3btrVs3brVsnTpUkvRokUzvYx1drKvcOcufaDVWvPmzWt5++23Lfv377d8/fXXFn9/f8vs2bNt27z33nuWAgUKmDLn27dvt3Tu3NlSoUIFS0xMjG2b9u3bW+rXr2/ZsGGDZe3ataZqYM+ePS2uoHfv3qZq608//WTKumv5dy3J/vzzz+fqPtDqjlqCXi96SB03bpz5t7Vyb2Y8Z612FRwcbHn00UdNVdy5c+ea19f06dMtOb0P4uLiLJ06dbKULl3avL/tj5X2VbhcvQ/S8lpIKml1Y3fpBz02aPW2GTNmmOPlpEmTLJ6enpY//vgj13xu3KgP9HNSq9ytWrXKcvDgQcusWbMsvr6+pnp1bumD3IhxX+rcbezHuO8Kxn6M/Rj7MfZj3Jf7xn4ECW+m80RSvOgf3EoDAQMHDrQULFjQfJHr2rWr+XJs7/Dhw5YOHTpY/Pz8TFBlxIgRlvj4eEtuGSi6Sx8sXrzYvGl9fHws1atXN1+A7SUmJlpeeeUV8wVft2nVqpVl7969DtucP3/eBAQCAwNNKfO+ffuaA44rCA8PN393PbDpAa9ixYqmpLt9ICg39oEe6FM6DujAOTOf87Zt2yzNmjUz96HBWA0+ukIfaMA4tWOl/l5u6YO0vBbSEiR0l36YOXOmpXLlyuZYUbduXcuiRYsc7sPVPzdu1Af6XPr06WMpWbKk6YNq1apZxo4da44XuaUPciPGfalzx7Gfu4/7FGM/xn6M/Rj7Me7LfWO/PPp/mTcvEQAAAAAAAICrISchAAAAAAAA4OYIEgIAAAAAAABujiAhAAAAAAAA4OYIEgIAAAAAAABujiAhAAAAAAAA4OYIEgIAAAAAAABujiAhAAAAAAAA4OYIEgIAAAAAAABujiAhADjZ3XffLc8++6yzdwMAAABZjHEfgJwsj8VisTh7JwDAHfz+++/SokULuXjxohQoUMDWfuHCBfHy8pJ8+fI5df8AAACQORj3AXBFeZ29AwDgDuLj41O9rVChQtm6LwAAAMg6jPsAuCqWGwNAEjNmzJCSJUtKYmKiQ3vnzp2lX79+5t8//PCDNGjQQHx9faVixYoyevRouXz5sm3bPHnyyNSpU6VTp04SEBAgTz75pJlFqAoWLGhu79OnT4rLTmJjY+WFF16QMmXKiI+Pj1SuXFlmzpxpu33nzp3SoUMHCQwMlODgYHn00Ufl3LlzttsXLFggt9xyi/j5+UnhwoWldevWEhUVxd8ZAACAcR8ApIogIQAk8cADD8j58+dl1apVDkuCly5dKg8//LD88ccf8thjj8nQoUNl165dMn36dPn888/l7bffdrif119/Xbp27So7duwwQcTvvvvOtO/du1dOnTolEydOTLHv9b6/+eYb+eijj2T37t3m/jUgqEJDQ6Vly5ZSv3592bRpk9mnkJAQ6dGjh7ld77dnz54mmKm/q0tdunXrJmSWAAAASI5xHwBcQ05CAEhBly5dzCw86ww+nV2ogb5jx45J27ZtpVWrVjJq1Cjb9rNnz5bnn39eTp48eeXgmiePmR04fvz4G+am0ZmE9erVkwkTJsi+ffukWrVq8uuvv5oZgEm99dZbJki5bNkyW9vx48fNrEMNPkZGRkrDhg3l8OHDUq5cOf62AAAAN8C4DwCuYCYhAKRAZwzqzD9d+qu+/vpreeihh8TDw0O2bdsmb7zxhpndZ73ocmKdxRcdHW27j0aNGqW7b7du3Sqenp7SvHnzFG/Xx9YZjvaPXb16dXPbgQMHpG7duiaAqcuN9cz4J598YoKSAAAASBnjPgC4gsIlAJCCjh07miW6P//8szRu3NjM3rPOCtTZejqrUJfxJqU5Cq00F2F6aR7B69HH1n17//33k91WokQJE2DUWYh//vmnLF++XCZNmiQvvfSSbNiwQSpUqMDfGgAAIAnGfQBwBUFCAEiBBvs0CKgzCP/77z+zBFgLlSj9qUt7taBIenh7e5ufCQkJqW6jMwC1YMrq1atTXG6sj60zHMuXLy9586Z8CNelznfccYe5vPrqq2bZ8cKFC2X48OH8rQEAAJJg3AcAV7DcGACus/REZxJ+9tln5t9WGnj78ssvzWzCf//91xQImTt3rrz88svX7UsN1mkA76effpKzZ8+aWYFJafCvd+/epvDIokWL5NChQyaX4bfffmtuHzRokCmiosVJNm7caJYYa37Cvn37muCjzhh85513TFGTo0ePyvfff28eq0aNGvydAQAAGPcBQKoIEgJAKrSKcKFChcyswV69etna27VrZwJ9upxXlyLfdtttZinyjQqFlCpVygQWX3zxRQkODpbBgwenuN3UqVPl/vvvl4EDB5p8g5rvMCoqytxWsmRJWbdunQkIagEVnXmoBVK0EIrmS8yfP7+sWbNG7rnnHqlataoJXI4dO1Y6dOjA3xkAAIBxHwCkiurGAAAAAAAAgJtjJiEAAAAAAADg5ggSAgAAAAAAAG6OICEAAAAAAADg5ggSAgAAAAAAAG6OICEAAAAAAADg5ggSAgAAAAAAAG6OICEAAAAAAADg5ggSAgAAAAAAAG6OICEAAAAAAADg5ggSAgAAAAAAAG6OICEAAAAAAAAg7u3/GT2Bbd/yaxQAAAAASUVORK5CYII=", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "def bench(solver, warmup=2, n=12):\n", " \"\"\"Median ms/step for one solver under the metronome rig (after warmup).\"\"\"\n", " J, p0 = solver.J, solver.p0\n", " s = np.zeros((solver.state_dim, 1)); s_dot = np.zeros_like(s)\n", " p = p0.copy(); p_dot = np.zeros_like(p); ts = []\n", " for i in range(warmup + n):\n", " p_next = handle(p0, i)\n", " t0 = time.perf_counter()\n", " s_new = solver.step(s, s_dot, p, p_dot, p_next)\n", " dt = time.perf_counter() - t0\n", " s_dot = (s_new - s) / DT; p_dot = (p_next - p) / DT; s, p = s_new, p_next.copy()\n", " if i >= warmup: ts.append(dt)\n", " return 1000 * np.median(ts)\n", "\n", "res_list = [(40, 5), (70, 7), (100, 8), (140, 9), (180, 10)]\n", "verts, t_fulls, t_reds = [], [], []\n", "for (w, h) in res_list:\n", " Xr, Tr = beam(w, h)\n", " tf = bench(FullCD(Xr, Tr))\n", " tr = bench(ReducedCD(Xr, Tr))\n", " verts.append(Xr.shape[0]); t_fulls.append(tf); t_reds.append(tr)\n", " print(f\"{Xr.shape[0]:5d} verts / {Tr.shape[0]:5d} tris: full {tf:6.1f} ms reduced {tr:5.2f} ms ({tf/tr:5.1f}x)\")\n", "\n", "fig, (axL, axR) = plt.subplots(1, 2, figsize=(13, 4.2))\n", "utils.line_plot(verts, {\"full-order\": t_fulls, \"fast (reduced)\": t_reds},\n", " xlabel=\"vertices\", ylabel=\"ms / step\", markers=True, logy=True,\n", " colors={\"full-order\": FULL_C, \"fast (reduced)\": RED_C}, ax=axL)\n", "axL.set_title(\"per-step cost\")\n", "axR.plot(verts, np.array(t_fulls) / np.array(t_reds), \"-o\", color=\"0.2\", lw=2)\n", "axR.set_xlabel(\"vertices\"); axR.set_ylabel(\"speedup (full / reduced)\")\n", "axR.set_title(\"speedup grows with resolution\"); axR.grid(True, color=\"0.9\")\n", "fig.tight_layout(); plt.show()" ] }, { "cell_type": "markdown", "id": "b2804ffd", "metadata": {}, "source": [ "## Takeaways\n", "\n", "* **Complementary dynamics = ordinary elastodynamics + one constraint.** Forcing the secondary\n", " motion to be rig-orthogonal, $J^{\\top} M D\\,u_c = 0$, is the whole idea. The momentum-leaking field\n", " $D$ chooses where tissue rides with the rig versus resists it.\n", "* **Full-order makes the constraint explicit** through a KKT global step in an ARAP local/global\n", " solve — the same block coordinate descent as tutorials 21 & 22, with one extra constraint block.\n", " Correct, but its cost scales with the mesh.\n", "* **Fast complementary dynamics moves both costs off the mesh.** A *rig-orthogonal skinning-eigenmode*\n", " subspace makes the constraint structural (so the global solve is tiny and unconstrained), and\n", " *rotation clustering* fixes the number of local rotations. Per-step cost becomes essentially\n", " independent of resolution — the result reported by Benchekroun et al. (2023).\n", "* **Everything here is one call away in `simkit`:** `lbs_jacobian`, `diffuse_field`,\n", " `lbs_weight_space_constraint`, `skinning_eigenmodes`, `spectral_cubature`,\n", " `cluster_grouping_matrices`, and `solvers.block_coord`." ] } ], "metadata": { "kernelspec": { "display_name": "simkit", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.13.9" } }, "nbformat": 4, "nbformat_minor": 5 }