diff --git a/icerm20_pymor_talk.ipynb b/icerm20_pymor_talk.ipynb new file mode 100644 index 0000000000000000000000000000000000000000..d0cda9628276f458bed220f3c453f7092761bbd2 --- /dev/null +++ b/icerm20_pymor_talk.ipynb @@ -0,0 +1,2169 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "<h1>pyMOR: Generic Interfaces for Model Order Reduction</h1>\n", + "<h3>René Fritze, <u>Petar Mlinarić</u>, Stephan Rave, Felix Schindler</h3>\n", + "<h3>ICERM Semester Program “Model and dimension reduction in uncertain and dynamic systemsâ€<br>\n", + " April 15, 2020</h3>" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Overview\n", + "\n", + "- 27k lines of Python code, 5k single commits.\n", + "- Reduced basis and system-theoretic MOR methods.\n", + "- Integration with external PDE solver packages.\n", + "- Support for MPI distributed computing.\n", + "- Permissive open source license (BSD-2-clause).\n", + "- Copyright holders are \"pyMOR developers and contributors\"." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# History\n", + "- Oct 2012,\n", + " development start at WWU Münster\n", + " (reduced basis methods)\n", + "- Apr 2013,\n", + " version 0.1\n", + "- May 2015,\n", + " first contributions from MPI Magdeburg\n", + "- Oct 2016,\n", + " pyMOR paper published:\n", + " R. Milk, S. Rave, F. Schindler,\n", + " \"pyMOR - Generic Algorithms and Interfaces for Model Order Reduction\"\n", + "- Jan 2019,\n", + " DFG project \"pyMOR - Sustainable Software for Model Order Reduction\"\n", + "- Jan 2019,\n", + " version 0.5\n", + " (the first version to include system-theoretic methods)\n", + "- Dec 2019,\n", + " version 2019.2" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Vectors (`pymor.vectorarrays`)\n", + "\n", + "## `VectorArrays` (`pymor.vectorarrays.interface.VectorArrayInterface`)\n", + "- a short sequence of high-dimensional vectors\n", + "- all vectors belong to the same `VectorSpace`\n", + "- mutable\n", + "- sub-classes:\n", + " - `NumpyVectorArray` (`NumPy`-based)\n", + " - `ListVectorArray` (typically for external solvers implementing `Vector` interface)\n", + " - `BlockVectorArray` (for use with `BlockOperators`)\n", + " - `MPIVectorArray` (for MPI distributed vectors)\n", + "\n", + "## `VectorSpaces` (`pymor.vectorarrays.interface.VectorSpaceInterface`)\n", + "- information to construct new `VectorArrays` (e.g., dimension, discrete function space, ...)\n", + "- sub-classes similar as for `VectorArrays` (`NumpyVectorSpace`, `ListVectorSpace`, ...)" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "from pymor.basic import NumpyVectorSpace" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [], + "source": [ + "NumpyVectorSpace?" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [], + "source": [ + "V = NumpyVectorSpace(5).ones(2)" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "NumpyVectorArray(\n", + " [[1. 1. 1. 1. 1.]\n", + " [1. 1. 1. 1. 1.]],\n", + " NumpyVectorSpace(5))" + ] + }, + "execution_count": 4, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "V" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "NumpyVectorArray(\n", + " [[2. 2. 2. 2. 2.]\n", + " [2. 2. 2. 2. 2.]],\n", + " NumpyVectorSpace(5))" + ] + }, + "execution_count": 5, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "2 * V" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [], + "source": [ + "w = NumpyVectorSpace(5).from_numpy(np.linspace(0, 1, 5))" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "NumpyVectorArray([[0. 0.25 0.5 0.75 1. ]], NumpyVectorSpace(5))" + ] + }, + "execution_count": 7, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "w" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "NumpyVectorArray(\n", + " [[1. 1.25 1.5 1.75 2. ]\n", + " [1. 1.25 1.5 1.75 2. ]],\n", + " NumpyVectorSpace(5))" + ] + }, + "execution_count": 8, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "V + w" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "NumpyVectorArray(\n", + " [[2. 2. 2. 2. 2.]\n", + " [2. 2. 2. 2. 2.]],\n", + " NumpyVectorSpace(5))" + ] + }, + "execution_count": 9, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "V + V" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([2.23606798, 2.23606798])" + ] + }, + "execution_count": 10, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "V.l2_norm()" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": {}, + "outputs": [], + "source": [ + "V.append(w)" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "NumpyVectorArray(\n", + " [[1. 1. 1. 1. 1. ]\n", + " [1. 1. 1. 1. 1. ]\n", + " [0. 0.25 0.5 0.75 1. ]],\n", + " NumpyVectorSpace(5))" + ] + }, + "execution_count": 12, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "V" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([[1. , 1. ],\n", + " [1. , 1. ],\n", + " [0. , 0.25]])" + ] + }, + "execution_count": 13, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "V.dofs([0, 1])" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Operators (`pymor.operators`)\n", + "\n", + "## `Operators` (`pymor.operators.interface.OperatorInterface`)\n", + "- mapping between two `VectorSpaces`\n", + "- can be nonlinear\n", + "- can be parametric\n", + "- immutable\n", + "- $x \\mapsto A_{\\mu}(x)$, $f(t, x) \\rightsquigarrow f_t(x)$\n", + "- sub-classes:\n", + " - `NumpyMatrixOperator` (`NumPy`/`SciPy`-based)\n", + " - `ZeroOperator`, `IdentityOperator`,\n", + " `BlockOperator`, `LincombOperator` (constructions)\n", + " - `FenicsMatrixOperator`, `NGSolveMatrixOperator` (for\n", + " external solvers)\n", + " - ..." + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "metadata": {}, + "outputs": [], + "source": [ + "from pymor.basic import NumpyMatrixOperator" + ] + }, + { + "cell_type": "code", + "execution_count": 15, + "metadata": {}, + "outputs": [], + "source": [ + "A = NumpyMatrixOperator(np.arange(1, 26).reshape(5, 5))" + ] + }, + { + "cell_type": "code", + "execution_count": 16, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "NumpyMatrixOperator(<5x5 dense>)" + ] + }, + "execution_count": 16, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "A" + ] + }, + { + "cell_type": "code", + "execution_count": 17, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "NumpyMatrixOperator: R^5 --> R^5 (parameter type: None, class: NumpyMatrixOperator)\n" + ] + } + ], + "source": [ + "print(A)" + ] + }, + { + "cell_type": "code", + "execution_count": 18, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([[ 1, 2, 3, 4, 5],\n", + " [ 6, 7, 8, 9, 10],\n", + " [11, 12, 13, 14, 15],\n", + " [16, 17, 18, 19, 20],\n", + " [21, 22, 23, 24, 25]])" + ] + }, + "execution_count": 18, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "A.matrix" + ] + }, + { + "cell_type": "code", + "execution_count": 19, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "LincombOperator((NumpyMatrixOperator(<5x5 dense>)), (2))" + ] + }, + "execution_count": 19, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "2 * A" + ] + }, + { + "cell_type": "code", + "execution_count": 20, + "metadata": {}, + "outputs": [], + "source": [ + "B = NumpyMatrixOperator(np.eye(5))" + ] + }, + { + "cell_type": "code", + "execution_count": 21, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([[1., 0., 0., 0., 0.],\n", + " [0., 1., 0., 0., 0.],\n", + " [0., 0., 1., 0., 0.],\n", + " [0., 0., 0., 1., 0.],\n", + " [0., 0., 0., 0., 1.]])" + ] + }, + "execution_count": 21, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "B.matrix" + ] + }, + { + "cell_type": "code", + "execution_count": 22, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "LincombOperator((NumpyMatrixOperator(<5x5 dense>), NumpyMatrixOperator(<5x5 dense>)), (1.0, 1.0))" + ] + }, + "execution_count": 22, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "A + B" + ] + }, + { + "cell_type": "code", + "execution_count": 23, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "NumpyVectorArray(\n", + " [[1. 1. 1. 1. 1. ]\n", + " [1. 1. 1. 1. 1. ]\n", + " [0. 0.25 0.5 0.75 1. ]],\n", + " NumpyVectorSpace(5))" + ] + }, + "execution_count": 23, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "V" + ] + }, + { + "cell_type": "code", + "execution_count": 24, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "NumpyVectorArray(\n", + " [[ 15. 40. 65. 90. 115. ]\n", + " [ 15. 40. 65. 90. 115. ]\n", + " [ 10. 22.5 35. 47.5 60. ]],\n", + " NumpyVectorSpace(5))" + ] + }, + "execution_count": 24, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "A.apply(V)" + ] + }, + { + "cell_type": "code", + "execution_count": 25, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "NumpyVectorArray(\n", + " [[55. 60. 65. 70. 75. ]\n", + " [55. 60. 65. 70. 75. ]\n", + " [40. 42.5 45. 47.5 50. ]],\n", + " NumpyVectorSpace(5))" + ] + }, + "execution_count": 25, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "A.apply_adjoint(V)" + ] + }, + { + "cell_type": "code", + "execution_count": 26, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "NumpyVectorArray(\n", + " [[-0.27717391 -0.14130435 -0.00543478 0.13043478 0.26630435]\n", + " [-0.27717391 -0.14130435 -0.00543478 0.13043478 0.26630435]\n", + " [ 0.05434783 0.0326087 0.01086957 -0.01086957 -0.0326087 ]],\n", + " NumpyVectorSpace(5))" + ] + }, + "execution_count": 26, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "(A + B).apply_inverse(V)" + ] + }, + { + "cell_type": "code", + "execution_count": 27, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "NumpyVectorArray(\n", + " [[-0.05978261 -0.0326087 -0.00543478 0.02173913 0.04891304]\n", + " [-0.05978261 -0.0326087 -0.00543478 0.02173913 0.04891304]\n", + " [ 0.2173913 0.14130435 0.06521739 -0.01086957 -0.08695652]],\n", + " NumpyVectorSpace(5))" + ] + }, + "execution_count": 27, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "(A + B).apply_inverse_adjoint(V)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Models (`pymor.models`)\n", + "\n", + "## `Models` (`pymor.models.interface.ModelInterface`)\n", + "- a structured collection of `Operators`\n", + "- can be nonlinear\n", + "- can be parametric\n", + "- immutable\n", + "- sub-classes:\n", + " - `StationaryModel`, `InstationaryModel`\n", + " (\"models for reduced basis methods\")\n", + " - `LTIModel`, `SecondOrderModel`, `LinearDelayModel`,\n", + " `LinearStochasticModel`, `BilinearModel`, `TransferFunction`\n", + " (\"models for system-theoretic methods\")" + ] + }, + { + "cell_type": "code", + "execution_count": 28, + "metadata": {}, + "outputs": [], + "source": [ + "from pymor.basic import StationaryModel, InstationaryModel, LTIModel" + ] + }, + { + "cell_type": "code", + "execution_count": 29, + "metadata": {}, + "outputs": [], + "source": [ + "StationaryModel?" + ] + }, + { + "cell_type": "code", + "execution_count": 30, + "metadata": {}, + "outputs": [], + "source": [ + "InstationaryModel?" + ] + }, + { + "cell_type": "code", + "execution_count": 31, + "metadata": {}, + "outputs": [], + "source": [ + "LTIModel?" + ] + }, + { + "cell_type": "code", + "execution_count": 32, + "metadata": {}, + "outputs": [], + "source": [ + "L = NumpyMatrixOperator(np.diag(np.arange(1, 6)))\n", + "F = L.source.ones(1)\n", + "m = StationaryModel(L, F)" + ] + }, + { + "cell_type": "code", + "execution_count": 33, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "StationaryModel(\n", + " NumpyMatrixOperator(<5x5 dense>),\n", + " VectorOperator(NumpyVectorArray([[1. 1. 1. 1. 1.]], NumpyVectorSpace(5)), name='rhs'),\n", + " products={})" + ] + }, + "execution_count": 33, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "m" + ] + }, + { + "cell_type": "code", + "execution_count": 34, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "00:28 StationaryModel: Solving StationaryModel for {} ...\n" + ] + }, + { + "data": { + "text/plain": [ + "NumpyVectorArray([[1. 0.5 0.33333333 0.25 0.2 ]], NumpyVectorSpace(5))" + ] + }, + "execution_count": 34, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "m.solve()" + ] + }, + { + "cell_type": "code", + "execution_count": 35, + "metadata": {}, + "outputs": [], + "source": [ + "from pymor.algorithms.timestepping import ExplicitEulerTimeStepper" + ] + }, + { + "cell_type": "code", + "execution_count": 36, + "metadata": {}, + "outputs": [], + "source": [ + "F = L.source.zeros(1)\n", + "u0 = L.source.ones(1)\n", + "m2 = InstationaryModel(1, u0, L, F, time_stepper=ExplicitEulerTimeStepper(10))" + ] + }, + { + "cell_type": "code", + "execution_count": 37, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "00:32 InstationaryModel: Solving InstationaryModel for {} ...\n" + ] + }, + { + "data": { + "text/plain": [ + "NumpyVectorArray(\n", + " [[1.00000000e+00 1.00000000e+00 1.00000000e+00 1.00000000e+00\n", + " 1.00000000e+00]\n", + " [9.00000000e-01 8.00000000e-01 7.00000000e-01 6.00000000e-01\n", + " 5.00000000e-01]\n", + " [8.10000000e-01 6.40000000e-01 4.90000000e-01 3.60000000e-01\n", + " 2.50000000e-01]\n", + " [7.29000000e-01 5.12000000e-01 3.43000000e-01 2.16000000e-01\n", + " 1.25000000e-01]\n", + " [6.56100000e-01 4.09600000e-01 2.40100000e-01 1.29600000e-01\n", + " 6.25000000e-02]\n", + " [5.90490000e-01 3.27680000e-01 1.68070000e-01 7.77600000e-02\n", + " 3.12500000e-02]\n", + " [5.31441000e-01 2.62144000e-01 1.17649000e-01 4.66560000e-02\n", + " 1.56250000e-02]\n", + " [4.78296900e-01 2.09715200e-01 8.23543000e-02 2.79936000e-02\n", + " 7.81250000e-03]\n", + " [4.30467210e-01 1.67772160e-01 5.76480100e-02 1.67961600e-02\n", + " 3.90625000e-03]\n", + " [3.87420489e-01 1.34217728e-01 4.03536070e-02 1.00776960e-02\n", + " 1.95312500e-03]\n", + " [3.48678440e-01 1.07374182e-01 2.82475249e-02 6.04661760e-03\n", + " 9.76562500e-04]],\n", + " NumpyVectorSpace(5))" + ] + }, + "execution_count": 37, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "m2.solve()" + ] + }, + { + "cell_type": "code", + "execution_count": 38, + "metadata": {}, + "outputs": [], + "source": [ + "A = np.diag(-np.arange(1, 6))\n", + "B = np.ones((5, 1))\n", + "C = np.ones((1, 5))\n", + "lti = LTIModel.from_matrices(A, B, C)" + ] + }, + { + "cell_type": "code", + "execution_count": 39, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "LTIModel(\n", + " NumpyMatrixOperator(<5x5 dense>, source_id='STATE', range_id='STATE'),\n", + " NumpyMatrixOperator(<5x1 dense>, range_id='STATE'),\n", + " NumpyMatrixOperator(<1x5 dense>, source_id='STATE'),\n", + " D=ZeroOperator(NumpyVectorSpace(1), NumpyVectorSpace(1)),\n", + " E=IdentityOperator(NumpyVectorSpace(5, id='STATE')))" + ] + }, + "execution_count": 39, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "lti" + ] + }, + { + "cell_type": "code", + "execution_count": 40, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "LTIModel\n", + " class: LTIModel\n", + " number of equations: 5\n", + " number of inputs: 1\n", + " number of outputs: 1\n", + " continuous-time\n", + " linear time-invariant\n", + " parameter_space: None\n", + " solution_space: NumpyVectorSpace(5, id='STATE')\n" + ] + } + ], + "source": [ + "print(lti)" + ] + }, + { + "cell_type": "code", + "execution_count": 41, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([[1.62760181-0.89728507j]])" + ] + }, + "execution_count": 41, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "lti.eval_tf(1j)" + ] + }, + { + "cell_type": "code", + "execution_count": 42, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEaCAYAAAAPGBBTAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3dd5hU5fnG8e8zu+xSFhbpvXekKNhjxy5WRE00mpigSYxRYmKLRmPUxJZfTOzGGBsCVixRY9QQRZQivRiKFEFAEFjawu4+vz/mEIfN7rBlZs/s2ftzXedi5p1T7ndndx5OmfOauyMiIlKeWNgBREQks6lQiIhIUioUIiKSlAqFiIgkpUIhIiJJqVCIiEhSKhQi1WRmc83sqBSu73MzG5aq9SWs9wkz+22q1yvRp0IhtVrwobrTzFqUav/UzNzMuqQ7g7v3d/f3g+3ebGZPp3ub6Rb87HqEnUMygwqFRMFS4PzdT8xsANAwvDgi0aJCIVHwFPDdhOcXAU8mzmBmpwR7GZvNbIWZ3Vzq9e+a2TIzW29mNyYe/gn2EsaZ2ZNmVhAcahqasOznZjbMzE4ErgfONbMtZjYz8fWE+ffY6zCzCxO2fUOpXDEzu9bMFgevjzOzZmX9EMzsKDNbaWbXm9lXwXa/U94Pzcx+aGaLzGyDmU0ws3ZB+8RglplBP84tbx1SN6hQSBRMBpqYWV8zywLOA0of/tlKvJg0BU4BfmRmZwCYWT/gAeA7QFsgH2hfavnTgOeC5ScAfy4dwt3fBG4Hxrp7nrsP2lvwYNsPAhcC7YDmQIeEWX4KnAEcGbz+NXB/klW2AVoE+S8CHjGz3mVs9xjgDmAk8T4vC/qHux8RzDYo6MfYvfVDok2FQqJi917FccB84IvEF939fXef7e4l7j4LGEP8wxdgBPCqu3/g7juBm4DSN0H7wN3fcPfiYFt7LQIVNAJ4zd0nunshcCNQkvD6ZcAN7r4yeP1mYISZZSdZ543uXuju/wJeJ14MSvsO8Li7Tw/Wex1wSE2c05HaJ9kvm0ht8hQwEehKqcNOAGZ2EPA7YF8gB8gFxgcvtwNW7J7X3beZ2fpSq/gy4fE2oL6ZZbt7UTVzl9721lLb7gy8ZGaJxaMYaE2pYhj42t23JjxfFmyjrO1OT9julmC77YHPK9sJiTbtUUgkuPsy4ie1TwZeLGOWZ4kfMuro7vnAQ4AFr60m4XCPmTUgfgioSlHKaNvKnifX2yQ8Xg10TNh2w1LbXgGc5O5NE6b67l5WkQDYx8waJTzvBKwqY75VxIvQ7u02CrZb3nqlDlOhkCi5BDim1P+od2sMbHD3HWZ2IPDthNeeB4ab2aFmlkP88I6VsY6KWAN0MbPEv60ZwHlmVi84CT6i1LZPNbNvBdv+DXv+XT4E3GZmnQHMrKWZnb6XDLeYWY6ZHQ6cyjd7TonGAN8zs8Fmlkv83MrH7v55Qj+6VaTDEn0qFBIZ7r7Y3aeW8/KPgd+YWQHxcxDjEpabS/yk8XPE/4e/BVgLFFYhxu4P5fVmtvvQzo1Ad+Inom8hvneTuO2fBG2rg3lWJqzvj8T3hN4Osk8GDkqy/S+DdawCngEuc/cFpWdy93eCXC8E2+1O/CKA3W4G/mZmG82srHMcUoeYBi4S2ZOZ5QEbgZ7uvjTsPBUVfDv8aXfvsLd5RSpDexQigJkNN7OGwbH6u4HZ6KSuCKBCIbLb6cQP16wCegLnuXa3RQAdehIRkb3QHoWIiCSlQiEiIklF8pvZLVq08C5duoQdQ0Sk1pg2bdpX7t6yrNciVSjMbDgwvEePHkydWt7l9CIiUpqZLSvvtUgdenL3V919VH5+fthRREQiI1KFQkREUk+FQkREklKhEBGRpFQoREQkKRUKERFJKlKXx1bXpEVfsXVncdgxMk55AzOYlfN49xJ7/oMFM1kwv2HBv/HGmMWXjMUsmMeIGWTFjJjFp/hjyM6KkR0zsrPibfViMbKyjJysGLnZsf9uS0SqT4UiwY2vzGHxurLGvJHaJic7Rm5WjNx6MXKzs6hfL0Ze/Xrk5WbRKCebvPrZ5OVm06R+PZrn5dA8L5cWeTm0yMuleaMcmjbMISumYiMCKhR7eOiCIRQWlex9xjqkvHtGesKIn4nz+H/bvNTzhCU93r67rcR3twX/erxtd3uJO8UlHrRBUYlTUuLsKi6huMTZVeIUF5dQVOIUFpVQWFTCzqISCouK4893lbB9VxFbCovZWljEVwXb2FJYxNadRWzevouSMvqYHTM67NOAzs0b0aV5w/i/LYJ/mzdSEZE6RYUiQc/WjcOOIDWspMTZuH0X67cU8tWWnazfWshXBYWsKShk+fptLNuwlenLvqagsOi/y+TlZjOgfT6DOzVlcMem7NexKa2a1A+xFyLppUIhdVosZjRrlEOzRjn0bF32PO7Ohq07WbZhG4vXbmHWyk3MWLGRRycuoSjYHWmXX5/DerTg+P5tOLxnC+rXy6rBXoikVyTHoxg6dKjrXk+Sbjt2FTN31SY+Xb6RT5dvZOJ/1lGwo4j69WIc0bMlx/dvwzF9WtGsUU7YUUX2ysymufvQsl6L1B5F4k0BRdKtfr0shnRuxpDOzQDYWVTCx0vX8495a3h77hrenreGmMFhPVpw4cGdObZva53bkFpJexQiaeDuzPliM2/N/ZLnp63ky807aN+0Ad8+qBPnHdCR5nm5YUcU2UOyPQoVCpE0Kyou4Z35a3jyo2VMWryenKwYpwxsy0WHdmFwx6ZhxxMBVChEMsaitQU89dEyXpj+BVsKixjWtxU/P743fds2CTua1HEqFCIZZkthEX+b9DkP/2sxBYVFDB/YjquO60XXFo3CjiZ1lAqFSIbatG0XD09czF8//JydxSWMHNqBK47tSdv8BmFHkzpGhUIkw60t2MED7y3mmY+XETPjimN78sPDu5GTrft2Ss1IVij0WyiSAVo1rs/Np/Xn3Z8fxTF9WnHXWwsZ/qcPmLbs67CjiahQiGSSjs0a8uAFQ3jsu0Mp2LGLEQ9N4lcvz2bzjl1hR5M6TIVCJAMN69eat0cfycWHduHZj5cz7J5/8ffZq4nioWLJfCoUIhkqLzebXw/vz0s/Pozmebn86JnpXDl2BlsSblAoUhNUKEQy3KCOTXn18sO4algvXp25ilPv+zdzvtgUdiypQ1QoRGqB7KwYPxvWkzE/PJgdu0o464FJPPHhUh2KkhqhQiFSixzUrTlv/OxwvtWzBTe/Oo9Ln5rGpm060S3ppUIhUss0a5TDXy4ayq9O6ct7C9dy8n3/5tPluoxW0keFQqQWMjN+cHg3nr/sUGIxOO+Rybw+a3XYsSSiVChEarFBHZvy8o8PY9/2+fzk2ek88P4inbeQlFOhEKnlmufl8swPDuK0Qe24882FXPvCbHYVl4QdSyIkUiPcidRV9etl8cfzBtOleUPue3cRK77exoMXDCG/Qb2wo0kEZPwehZk1MrO/mdmjZvadsPOIZCozY/Txvbn7nEFM+XwDZz84iRUbtoUdSyIglEJhZo+b2Vozm1Oq/UQzW2hmi8zs2qD5LOB5d/8hcFqNhxWpZUYM6cCT3z+IdQWFnPXgJBatLQg7ktRyYe1RPAGcmNhgZlnA/cBJQD/gfDPrB3QAVgSzFddgRpFa65DuzXn+skNwj18RtfBLFQupulAKhbtPBDaUaj4QWOTuS9x9J/AccDqwknixgCR5zWyUmU01s6nr1q1LR2yRWqVn68Y8N+pgYmac/+hk5q3aHHYkqaUy6RxFe77Zc4B4gWgPvAicbWYPAq+Wt7C7P+LuQ919aMuWLdObVKSW6NEqj7GXHkJudoxvPzZZ94iSKsmkQlEmd9/q7t9z9x+5+zNh5xGpbbq2aMTYUYfQKCebbz86mZkrNoYdSWqZTCoUXwAdE553CNoqzMyGm9kjmzbpf00iiTo1b8hzow4mv2E9LnjsY6brlh9SCZlUKKYAPc2sq5nlAOcBEyqzAnd/1d1H5efnpyWgSG3WsVlDxo46hOZ5OXz3L5/oMJRUWFiXx44BPgJ6m9lKM7vE3YuAy4G3gPnAOHefG0Y+kahq17QBY0YdTH6Delz81yksW7817EhSC1gU7wszdOhQnzp1atgxRDLWorVbGPHQJPIb1OP5yw6lZePcsCNJyMxsmrsPLeu1TDr0VG06RyFSMT1a5fHXiw9g7eZCLv7rJxTs0JgWUr5IFQqdoxCpuP067cMDF+zPgi8LuPSpaRQW6fusUrZIFQoRqZyje7fizrMHMmnxekaPnUlxSfQORUv16e6xInXc2UM6sH5rIbe/sYDmeTncclp/zCzsWJJBVChEhFFHdGddQSGP/nspnZo15AeHdws7kmSQSB160slskaq77qS+nNC/Nbe/MZ/3Fq4NO45kkEgVCp3MFqm6WMy4d+RgerdpwhXPfsqitVvCjiQZIlKFQkSqp1FuNo9+dwg52TF+8LcpbNy2M+xIkgFUKERkDx32acjDFw7hi43bufzZTynS+Nt1ngqFiPyPoV2acdsZA/hg0Vf89vX5YceRkEWqUOhktkjqjDygI5d8qytPTPqcMZ8sDzuOhChShUIns0VS67qT+nBEr5bc+PIcPllaelBKqSsiVShEJLWys2L86fz96NisIT95djprC3aEHUlCoEIhIknlN6jHgxfsT8GOXfxUJ7frJBUKEdmrPm2acPuZA/h46QbufvuzsONIDYtUodDJbJH0OWv/Dnz7oE489K/FvD33y7DjSA2KVKHQyWyR9Lrp1H4MaJ/Pz8fP1Oh4dUikCoWIpFf9elk88J39iZlx2dPT2bFLY1jUBSoUIlIpHZs15A/nDmL+6s3c9MqcsONIDVChEJFKO6ZPay4/ugfjpq5k3JQVYceRNFOhEJEqueq4XhzavTk3TZjDf9YUhB1H0kiFQkSqJCtm/N+5g2mUk83lz36q8xURpkIhIlXWqkl97hk5iIVrCrj1tXlhx5E0iVSh0PcoRGreUb1bcekR3Xjm4+W8MXt12HEkDSJVKPQ9CpFwXH1CbwZ3bMo1L8xixYZtYceRFItUoRCRcNQLbh6Iw0/HfMou3Q8qUlQoRCQlOjZryO/OHsiMFRu5R/eDihQVChFJmVMGtuX8A+P3g5r42bqw40iKqFCISErddGo/erXOY/S4GXy1pTDsOJICKhQiklINcrL40/n7s3lHEb8YPxN3DzuSVJMKhYikXO82jbn+pD68t3AdT01eFnYcqSYVChFJi4sO7cJRvVty2+vz+Uy3+KjVIlUo9IU7kcxhZtw1YhCN62dzxRjd4qM2i1Sh0BfuRDJLy8a53DViEAu+LODONxeGHUeqKFKFQkQyz9F9WnHRIZ15/MOlvL9wbdhxpApUKEQk7a47uS+9Wzfm6vGzWK9LZmsdFQoRSbv69bL44/mD2bxjF798fpYuma1lVChEpEb0adOEa0/swz8XrOXZT5aHHUcqQYVCRGrMxYd24fCeLfjta/NZsm5L2HGkglQoRKTGxGLG3ecMIrdejKvGztBdZmsJFQoRqVGtm9Tn9jMHMHPlJv707qKw40gFVKhQWNwFZnZT8LyTmR2Y3mgiElUnD2jLWfu35/73FjFt2ddhx5G9qOgexQPAIcD5wfMC4P60JBKROuGW0/rTNr8+o8fNYGthUdhxJImKFoqD3P0nwA4Ad/8ayElbKhGJvMb163HvyMEs37CNW1+bF3YcSaKihWKXmWUBDmBmLQGdhRKRajmwazMuO7I7z01Zwdtzvww7jpSjooXiPuAloJWZ3QZ8ANyetlRVpJsCitQ+Vw3rRf92Tbj2xdmsLdgRdhwpQ4UKhbs/A/wSuANYDZzh7uPTGawqdFNAkdonJzvG/507mK2FRVz3wmx9azsDJS0UZtZs9wSsBcYAzwJrgjYRkWrr2box1wTf2h47ZUXYcaSU7L28Po34eQkDOgFfB4+bAsuBrmlNJyJ1xsWHduGfC9bwm9fmcUj35nRu3ijsSBJIukfh7l3dvRvwDjDc3Vu4e3PgVODtmggoInVDLBYf6CgrZvx83EyKS3QIKlNU9GT2we7+xu4n7v534ND0RBKRuqpd0wbcevq+TF32NQ9PXBx2HAlUtFCsMrNfmVmXYLoBWJXOYCJSN50+uB2nDGjLH/7xGXNX6QrGTFDRQnE+0JL4JbIvAa345lvaIiIpY2b89ox92adhDleNnaGxtjNARS+P3eDuP3P3/YLpZ+6+Id3hRKRu2qdRDneOGMhna7Zwz9saaztse7vqCQAze4/gW9mJ3P2YlCcSEQGO6t2KCw7uxGMfLOWYPq05pHvzsCPVWRUqFMDVCY/rA2cDuouXiKTV9Sf35cNF67l6/EzevPJwGtevF3akOqmih56mJUwfuvto4Kj0RhORuq5hTjb3jBzE6k3b+c2runFgWCo6HkWzhKmFmZ0A6D4ZIpJ2+3fahx8f1YPx01bqxoEhqeihp8RvaBcBS4FL0hVKRCTRFcf25N0Fa7nuxdns33kfWuTlhh2pTqno5bF93b1b8E3tnu5+PDAlncFERHbLyY7xh3MHU7CjiOte1I0Da1pFC8WkMto+SmUQEZFkerdpzNUn9OIf89bwwvQvwo5TpyQ99GRmbYD2QAMz24/4oSeAJkDDNGcTEdnDJd/qxjvz13LzhLkc3K0ZHfbRx1BN2NsexQnA3UAH4F7gnmAaDVyf3mgiInvKihn3nDMId+fq8TMp0Y0Da8Te7h77N3c/GrjY3Y9OmE5z9xdrKKOIyH91bNaQm4b3Y/KSDTz+4dKw49QJezv0dIG7Pw10MbPRpV9393vTluybDN2AG4B8dx+R7u2JSOYbObQjb89dw51vLeTIXi3p2bpx2JEibW+HnnaPHJIHNC5jSsrMHjeztWY2p1T7iWa20MwWmdm1ydbh7kvcXZfiish/mRl3nD2ARjlZjB43k13FJWFHirSkexTu/nDw7y1VXP8TwJ+BJ3c3mFkWcD9wHLASmGJmE4As4mNyJ/q+u6+t4rZFJMJaNa7PbWcO4MfPTOf+9xZx5bBeYUeKrIreFLAl8EOgS+Iy7v79ZMu5+0Qz61Kq+UBgkbsvCdb9HHC6u99BfOS8KjGzUcAogE6dOlV1NSJSi5w8oC1nDG7Hn95dxDF9WjGwQ9OwI0VSRb9H8QrxW3a8A7yeMFVFeyBx9PSVQVuZzKy5mT0E7Gdm15U3n7s/4u5D3X1oy5YtqxhNRGqbW07bl5Z5uYweN1NjV6RJRW/h0dDdr0lrknK4+3rgsjC2LSKZL79hPe4cMZDvPv4Jd721kBtP7Rd2pMip6B7Fa2Z2coq2+QXQMeF5h6BNRKRKjujVkgsP7szjHy7lo8Xrw44TORUtFD8jXiy2m9lmMysws81V3OYUoKeZdTWzHOA8YEIV17UHMxtuZo9s2qRxdkXqmutO7kPnZg25evxMCnbsCjtOpFR0PIrG7h5z9wbu3iR43mRvy5nZGOL3hOptZivN7BJ3LwIuB94C5gPj3H1udTqRkPNVdx+Vn687oIvUNYljV9z6msauSKWKXvW0fxnNm4BlwQd/mdz9/HLa3wDeqFBCEZEKGtK5GZce2Z0H31/M8f3aMKxf67AjRUJFDz09AEwGHg2mycB4YKGZHZ+mbCIilXblsJ70adOYa1+czYatO8OOEwkVLRSrgP3cfYi7DwEGA0uIf2nuznSFqyydoxCR3Ows/nDuYDZt38mNL8/R2BUpUNFC0SvxPIK7zwP67P7SXKbQOQoRAejbtglXDuvF67NXM2HmqrDj1HoVLRRzzexBMzsymB4A5plZLqDLC0Qk41x6RDf269SUG1+ew5ebdoQdp1araKG4GFgEXBlMS4K2XcDR6QgmIlId2Vkx7h05mJ3FJVzzwiwdgqqGil4eu93d73H3M4Ppbnff5u4l7r4l3SErSucoRCRR1xaNuP7kvvzrs3U8+8nysOPUWhUqFGbW08yeN7N5ZrZk95TucJWlcxQiUtoFB3XmWz1acNvr81m2fmvYcWqlih56+ivwIFBE/FDTk8DT6QolIpIqsZhx54iBZMWMq8fPpFjDp1ZaRQtFA3f/J2DuvszdbwZOSV8sEZHUade0ATcP78+Uz7/msX9n3MGQjFfRQlFoZjHgP2Z2uZmdSXzUOxGRWuGs/dtzfL/W3PP2Z3y2piDsOLVKZW4K2BC4AhgCXAhclK5QVaWT2SJSHjPj9rMG0Lh+NqPHzdDwqZVQ0aueprj7Fndf6e7fc/ez3H1yusNVlk5mi0gyLfJyue3MAcz5YjN/endR2HFqjaQ3BQzGsi6Xu5+W2jgiIul14r5tOGu/9tz/3iKO7dOKQR01fOre7O3usYcQH7Z0DPAxYGlPJCKSZr8+rT8fLVnP6HEzeP2Kw6lfLyvsSBltb4ee2gDXA/sCfyR+E8Cv3P1f7v6vdIcTEUmH/Abx4VMXr9vKXW8tDDtOxktaKNy92N3fdPeLgIOJ38bjfTO7vEbSiYikyeE9vxk+dfISDZ+azF5PZptZrpmdRfwLdj8B7gNeSnewqtBVTyJSGRo+tWKSFgoze5L4UKb7A7e4+wHufqu7f1Ej6SpJVz2JSGXsHj511cbt3Pb6/LDjZKy97VFcAPQk/j2KSWa2OZgKzGxz+uOJiKTXkM7NGHVEd56bsoJ3F6wJO05G2ts5ipi7Nw6mJglTY3dvUlMhRUTS6arj4sOnXvPCbL7W8Kn/o6LfzBYRiazc7CzuGTmIjdt2cuMrc8KOk3FUKEREgP7t8vnZsT15bdZqXtXwqXtQoRARCVx2ZHcGdWzKja/MYe1mDZ+6W6QKhS6PFZHqyM6Kcc85g9i+s1jDpyaIVKHQ5bEiUl09WuVxzYl9eG/hOsZOWRF2nIwQqUIhIpIKFx/ahUO6NefW1+axYsO2sOOEToVCRKSUWMy465yBmMWHTy2p48OnqlCIiJShwz4NuenUfny8dAOPf7g07DihUqEQESnHOUM7cGyfVtz51kIWra27w6eqUIiIlMPMuOPsATTMyeLn42ZSVEeHT1WhEBFJolXj+vz2jH2ZuXITD7y/OOw4oVChEBHZi1MHtmP4oHbc98//MOeLuvc9rUgVCn3hTkTS5dbT+9OsUQ6jx81gx67isOPUqEgVCn3hTkTSpWnDHH5/9kA+W7OFP7zzWdhxalSkCoWISDod3acV5x/YkUcmLmHq5xvCjlNjVChERCrhhlP60WGfBvx8/Ey2FhaFHadGqFCIiFRCXm42d48YxPIN27jj73Vj+FQVChGRSjqoW3MuOawrT09ezsTP1oUdJ+1UKEREquDqE3rTo1Uev3x+Fpu27wo7TlqpUIiIVEH9elncO3IQ67YUcsuEuWHHSSsVChGRKhrYoSk/OboHL376BW/O+TLsOGmjQiEiUg0/PaYH/ds14YaXZvPVlsKw46SFCoWISDXUy4px78jBFOwo4voXZ0dy+FQVChGRaurdpjE/P74Xb89bw0uffhF2nJRToRARSYEfHN6NA7rsw68nzGXVxu1hx0mpSBUK3RRQRMKSFTPuPmcQRcXONS/MitQhqEgVCt0UUETC1Ll5I64/pS///s9XPP3x8rDjpEykCoWISNguOKgTh/dswe2vz+fzr7aGHSclVChERFLIzLhzxECys4yrx8+kuKT2H4JSoRARSbG2+Q34zen9mbrsax7795Kw41SbCoWISBqcMbg9J/Zvwz1vf8bCLwvCjlMtKhQiImlgZtx25r40rp/N6HEz2FlUEnakKlOhEBFJk+Z5udx+1gDmrtrMn99bFHacKlOhEBFJoxP6t+Gs/dtz/3uLmLliY9hxqkSFQkQkzX49vD+tGucyetwMduwqDjtOpalQiIikWX6Detw5YiCL123lrrcWhh2n0lQoRERqwOE9W3LhwZ15/MOlTF6yPuw4laJCISJSQ647uQ+dmjXk6vEz2VJYFHacClOhEBGpIQ1zsrnnnEF8sXE7t70+L+w4FaZCISJSg4Z2acaoI7ox5pMVvLdgbdhxKkSFQkSkho0+rhe9WudxzQuz2LhtZ9hx9kqFQkSkhuVmZ3HvyMFs2LqTm16ZG3acvVKhEBEJwb7t87ni2J5MmLmK12etDjtOUioUIiIh+fFR3RnUIZ9fvTybtQU7wo5TrowvFGZ2hpk9amZjzez4sPOIiKRKdlaMe0YOZtvOYq57YXbGDp+a1kJhZo+b2Vozm1Oq/UQzW2hmi8zs2mTrcPeX3f2HwGXAuenMKyJS03q0yuMXJ/TmnwvWMn7ayrDjlCndexRPACcmNphZFnA/cBLQDzjfzPqZ2QAze63U1Cph0V8Fy4mIRMr3D+vKQV2b8ZtX57Hy621hx/kfaS0U7j4R2FCq+UBgkbsvcfedwHPA6e4+291PLTWttbjfA3939+npzCsiEoZYzLj7nEG4O798fhYlGTZ8ahjnKNoDKxKerwzayvNTYBgwwswuK28mMxtlZlPNbOq6detSk1REpIZ0bNaQG0/tx6TF63nyo8/DjrOHjD+Z7e73ufsQd7/M3R9KMt8j7j7U3Ye2bNmyJiOKiKTEuQd05KjeLfndmwtYsm5L2HH+K4xC8QXQMeF5h6BNRKROMzN+f/ZAcrOzGD1uJkXFmTF8ahiFYgrQ08y6mlkOcB4wIRUrNrPhZvbIpk2bUrE6EZEa17pJfW49Y19mrNjIwxOXhB0HSP/lsWOAj4DeZrbSzC5x9yLgcuAtYD4wzt1T8h12d3/V3Ufl5+enYnUiIqE4bVA7ThnYlv975zPmrdocdhwsU7/gUR1Dhw71qVOnhh1DRKTKvt66k+P+MJEWeTm8cvlh5GZnpXV7ZjbN3YeW9VrGn8wWEamL9mmUw+/PHsCCLwu475//CTVLpAqFzlGISJQc27c1I4d24MH3FzN9+deh5YhUodA5ChGJmhtP7Ufb/AZcPW4m23cWh5IhUoVCRCRqGtevx13nDGTJV1v5/ZsLQsmgQiEikuEO7d6Ciw/twhOTPmfSoq9qfPuRKhQ6RyEiUXXNiX3o1qIRv3h+Fpt37KrRbUeqUOgchYhEVYOcLO4eOYjVm7Zz66vzanTbkSoUIiJRtn+nffjRUd0ZP20l78xbU2PbVaEQEalFfnZsL/q2bd7vCwkAAAjLSURBVMK1L85mw9adNbJNFQoRkVokJzvGvSMHsWn7Tm58eU6NDJ8aqUKhk9kiUhf0bduEK4f14vXZq5kwc1XatxepQqGT2SJSV1x6RDf269SUm16Zy5rNO9K6rUgVChGRuiI7K8a9IwdTWFTMNS/MSushKBUKEZFaqmuLRlx3Ul/eX7iO56as2PsCVaRCISJSi114cGcO7d6c3742jxUbtqVlGyoUIiK1WCxm3HXOIGJm/Hz8TEpKUn8IKlKFQlc9iUhd1L5pA24a3o/1WwpZt6Uw5evXCHciIhHg7uwsLqnySHga4U5EJOLMLG3DpapQiIhIUioUIiKSlAqFiIgkpUIhIiJJRapQ6PJYEZHUi1Sh0E0BRURSL1KFQkREUi+SX7gzs3XAsoSmfGBTBR+3AL6q4qYT11fZecpqL92W7Pnux4lttbEvqX5PkuWsyDyV7Uum/n6V91pl+1Lbf78SH9fGvqTz96uzu7cs8xV3j/wEPFLRx8DUVGynsvOU1V66LdnzhPyJbbWuL6l+T2q6L5n6+5WqvtT236/a3pd0/n4lm+rKoadXK/k4Fdup7DxltZduS/b81XLmqaqw+pLq96Si60lVXzL196u81yrbl9r4npR+Xpv7ks7fr3JF8tBTdZjZVC/nfie1TVT6EpV+gPqSqaLSl3T1o67sUVTGI2EHSKGo9CUq/QD1JVNFpS9p6Yf2KEREJCntUYiISFIqFCIikpQKhYiIJKVCUQlm1tfMHjKz583sR2HnqSozO8PMHjWzsWZ2fNh5qsPMupnZX8zs+bCzVIWZNTKzvwXvx3fCzlMdtf292C1ifx+p+cxKx5czMnECHgfWAnNKtZ8ILAQWAddWcF0x4OkI9GMf4C8ReU+eD/t3rCr9Ai4EhgePx4adPRXvUSa9F9XsR6h/HynuS7U+s0LvdA3+cI8A9k/84QJZwGKgG5ADzAT6AQOA10pNrYJlTgP+Dny7NvcjWO4eYP/a/p4Ey2XMh1Ml+3UdMDiY59mws1enL5n4XlSzH6H+faSqL6n4zMqmjnD3iWbWpVTzgcAid18CYGbPAae7+x3AqeWsZwIwwcxeB55NX+KypaIfZmbA74C/u/v09CYuX6rek0xTmX4BK4EOwAwy8FBwJfsyr2bTVVxl+mFm88mAv4/yVPY9ScVnVsb9Ytaw9sCKhOcrg7YymdlRZnafmT0MvJHucJVQqX4APwWGASPM7LJ0BquCyr4nzc3sIWA/M7su3eGqobx+vQicbWYPkubbMKRQmX2pRe/FbuW9J5n891Ge8t6TlHxm1Zk9ilRw9/eB90OOUW3ufh9wX9g5UsHd1wO15Y/5f7j7VuB7YedIhdr+XuwWsb+P90nBZ1Zd36P4AuiY8LxD0FbbRKUfEK2+JIpSv6LSl6j0A9Lcl7peKKYAPc2sq5nlAOcBE0LOVBVR6QdEqy+JotSvqPQlKv2AdPcl7DP4NXilwBhgNbCL+PG7S4L2k4HPiF8xcEPYOetKP6LWl6j2Kyp9iUo/wuqLbgooIiJJ1fVDTyIishcqFCIikpQKhYiIJKVCISIiSalQiIhIUioUIiKSlAqF1HpmVmxmMxKmLmFnShUz28/M/lLNdTxhZiMSnp9nZjeUM29LM3uzOtuT6NG9niQKtrv74LJeCO6Ua+5eUsOZUuV64LelG80s292LqrjOkyjnXkbuvs7MVpvZYe7+YRXXLxGjPQqJHDPrYmYLzexJYA7Q0cx+YWZTzGyWmd2SMO8NZvaZmX1gZmPM7Oqg/X0zGxo8bmFmnwePs8zsroR1XRq0HxUs87yZLTCzZ4IihZkdYGaTzGymmX1iZo3NbKKZDU7I8YGZDSrVj8bAQHefGTy/2cyeMrMPgaeCfv7bzKYH06HBfGZmfw5+Bu8ArRLWacBgYLqZHZmwF/ZpsD2Al4FaPdqepJb2KCQKGpjZjODxUuAqoCdwkbtPtvhwlj2J37PfiN+b/whgK/F74gwm/rcwHZi2l21dAmxy9wPMLBf40MzeDl7bD+gPrAI+BA4zs0+AscC57j7FzJoA24G/ABcDV5pZL6D+7oKQYCjxQpeoH/Atd99uZg2B49x9h5n1JH5rh6HAmUDvYN7WxMeJeDwh40x396Ao/sTdPzSzPGBHMM9UytiLkbpLhUKiYI9DT8E5imXuPjloOj6YPg2e5xEvHI2Bl9x9W7BcRW6idjwwMOGYf36wrp3AJ+6+MljXDKALsAlY7e5TANx9c/D6eOBGM/sF8H3giTK21RZYV6ptgrtvDx7XA/4c7JkUA72C9iOAMe5eDKwys3cTlj+R+GhnEC9m95rZM8CLu7MTH2azXQV+FlJHqFBIVG1NeGzAHe7+cOIMZnZlkuWL+ObQbP1S6/qpu79Val1HAYUJTcUk+fty921m9g/io5CNBIaUMdv2UtuGPft1FbAGGBRk3cHeHQ+cHWT4ncVHPTuZ+J7RCe6+INjm9iTrkDpG5yikLngL+H5weAUza29mrYCJwBlm1iA4Pj88YZnP+ebDe0Spdf3IzOoF6+plZo2SbHsh0NbMDgjmb2xmuwvIY8RPKk9x96/LWHY+0CPJuvOJ762UABcSHzeZoF/nBudT2gJHB9vOB7I9PsAQZtbd3We7+++J36a6T7B8L/73kJfUYdqjkMhz97fNrC/wUXB+eQtwgbtPN7OxxAeiX0v8w3K3u4FxZjYKeD2h/THih5SmByeG1wFnJNn2TjM7F/iTmTUg/j/1YcAWd59mZpuBv5az7AIzyzezxu5eUMYsDwAvmNl3gTf5Zm/jJeAY4ucmlgMfBe3HAe8kLH+lmR0NlABz+eaQ1NGl+ix1nG4zLhIws5uJf4DfXUPba0d8mMo+5V2+a2ZXAQXu/lgKtvcY8FjCuZvy5psInF7OXo7UQTr0JBKCYC/gY+IDzCT7jseD7Hnuo8rc/QcVKBItgXtVJCSR9ihERCQp7VGIiEhSKhQiIpKUCoWIiCSlQiEiIkmpUIiISFIqFCIiktT/A/4Q4lcyzkVYAAAAAElFTkSuQmCC\n", + "text/plain": [ + "<Figure size 432x288 with 1 Axes>" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "_ = lti.mag_plot(np.logspace(-3, 3, 50))" + ] + }, + { + "cell_type": "code", + "execution_count": 43, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([1.05594884e+00, 8.22854101e-02, 3.35759106e-03, 7.41367413e-05,\n", + " 6.87894883e-07])" + ] + }, + "execution_count": 43, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "lti.hsv()" + ] + }, + { + "cell_type": "code", + "execution_count": 44, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "['A',\n", + " 'B',\n", + " 'C',\n", + " 'D',\n", + " 'E',\n", + " '_BasicInterface__auto_init',\n", + " '_BasicInterface__implementors',\n", + " '_CacheableInterface__auto_init',\n", + " '_CacheableInterface__implementors',\n", + " '_ImmutableInterface__auto_init',\n", + " '_ImmutableInterface__implementors',\n", + " '_InputOutputModel__auto_init',\n", + " '_InputOutputModel__implementors',\n", + " '_InputStateOutputModel__auto_init',\n", + " '_InputStateOutputModel__implementors',\n", + " '_LTIModel__auto_init',\n", + " '_ModelBase__auto_init',\n", + " '_ModelBase__implementors',\n", + " '_ModelInterface__auto_init',\n", + " '_ModelInterface__implementors',\n", + " '_Parametric__implementors',\n", + " '__abstractmethods__',\n", + " '__add__',\n", + " '__class__',\n", + " '__copy__',\n", + " '__deepcopy__',\n", + " '__delattr__',\n", + " '__dict__',\n", + " '__dir__',\n", + " '__doc__',\n", + " '__eq__',\n", + " '__format__',\n", + " '__ge__',\n", + " '__getattribute__',\n", + " '__gt__',\n", + " '__hash__',\n", + " '__init__',\n", + " '__init_subclass__',\n", + " '__le__',\n", + " '__lt__',\n", + " '__module__',\n", + " '__mul__',\n", + " '__ne__',\n", + " '__neg__',\n", + " '__new__',\n", + " '__reduce__',\n", + " '__reduce_ex__',\n", + " '__repr__',\n", + " '__setattr__',\n", + " '__signature__',\n", + " '__sizeof__',\n", + " '__str__',\n", + " '__sub__',\n", + " '__subclasshook__',\n", + " '__weakref__',\n", + " '_abc_cache',\n", + " '_abc_negative_cache',\n", + " '_abc_negative_cache_version',\n", + " '_abc_registry',\n", + " '_cached_method_call',\n", + " '_format_repr',\n", + " '_generate_sid',\n", + " '_hsv_U_V',\n", + " '_implements_reduce',\n", + " '_init_arguments',\n", + " '_locked',\n", + " '_logger',\n", + " '_name',\n", + " '_parameter_space',\n", + " '_solve',\n", + " '_uid',\n", + " 'build_parameter_type',\n", + " 'cache_region',\n", + " 'cached_method_call',\n", + " 'cont_time',\n", + " 'disable_caching',\n", + " 'disable_logging',\n", + " 'enable_caching',\n", + " 'enable_logging',\n", + " 'estimate',\n", + " 'estimator',\n", + " 'eval_dtf',\n", + " 'eval_tf',\n", + " 'freq_resp',\n", + " 'from_abcde_files',\n", + " 'from_files',\n", + " 'from_mat_file',\n", + " 'from_matrices',\n", + " 'generate_sid',\n", + " 'gramian',\n", + " 'h2_norm',\n", + " 'hankel_norm',\n", + " 'has_interface_name',\n", + " 'hinf_norm',\n", + " 'hsv',\n", + " 'implementor_names',\n", + " 'implementors',\n", + " 'input_dim',\n", + " 'input_space',\n", + " 'linear',\n", + " 'logger',\n", + " 'logging_disabled',\n", + " 'mag_plot',\n", + " 'name',\n", + " 'order',\n", + " 'output',\n", + " 'output_dim',\n", + " 'output_space',\n", + " 'parameter_space',\n", + " 'parameter_type',\n", + " 'parametric',\n", + " 'parse_parameter',\n", + " 'poles',\n", + " 'products',\n", + " 'sid_ignore',\n", + " 'solution_space',\n", + " 'solve',\n", + " 'solver_options',\n", + " 'strip_parameter',\n", + " 'uid',\n", + " 'visualize',\n", + " 'visualizer',\n", + " 'with_']" + ] + }, + "execution_count": 44, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "dir(lti)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Parameters (`pymor.parameters`)\n", + "\n", + "## `Parameters`\n", + "- a dictionary mapping strings to `NumPy` arrays (parameter name $\\mapsto$ parameter value)\n", + "\n", + "## `ParameterTypes`\n", + "- a dictionary mapping strings to tuples (parameter name $\\mapsto$ parameter shape)\n", + "\n", + "## `ParameterSpaces`\n", + "- the set of possible parameter values\n", + "- `CubicParameterSpace` is the only sub-class\n", + "\n", + "## `ParameterFunctionals`\n", + "- a function mapping a `Parameter` to a number\n", + "- sub-classes: `ProjectionParameterFunctional`, `GenericParameterFunctional`, `ExpressionParameterFunctional`, ..." + ] + }, + { + "cell_type": "code", + "execution_count": 45, + "metadata": {}, + "outputs": [], + "source": [ + "from pymor.basic import Parameter" + ] + }, + { + "cell_type": "code", + "execution_count": 46, + "metadata": {}, + "outputs": [], + "source": [ + "mu = Parameter({'diffusion': [[1, 2], [3, 4]], 'p': 1})" + ] + }, + { + "cell_type": "code", + "execution_count": 47, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "{'diffusion': array([[1, 2],\n", + " [3, 4]]),\n", + " 'p': array(1)}" + ] + }, + "execution_count": 47, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "mu" + ] + }, + { + "cell_type": "code", + "execution_count": 48, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "ParameterType({'diffusion': (2, 2), 'p': ()})" + ] + }, + "execution_count": 48, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "mu.parameter_type" + ] + }, + { + "cell_type": "code", + "execution_count": 49, + "metadata": {}, + "outputs": [], + "source": [ + "from pymor.basic import ProjectionParameterFunctional" + ] + }, + { + "cell_type": "code", + "execution_count": 50, + "metadata": {}, + "outputs": [], + "source": [ + "theta = ProjectionParameterFunctional('diffusion', (2, 2), (0, 1))" + ] + }, + { + "cell_type": "code", + "execution_count": 51, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "ProjectionParameterFunctional('diffusion', (2, 2), index=(0, 1))" + ] + }, + "execution_count": 51, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "theta" + ] + }, + { + "cell_type": "code", + "execution_count": 52, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "2" + ] + }, + "execution_count": 52, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "theta.evaluate(mu)" + ] + }, + { + "cell_type": "code", + "execution_count": 53, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "ParameterType({'diffusion': (2, 2)})" + ] + }, + "execution_count": 53, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "theta.parameter_type" + ] + }, + { + "cell_type": "code", + "execution_count": 54, + "metadata": {}, + "outputs": [], + "source": [ + "L1 = NumpyMatrixOperator(np.diag(np.arange(1, 6)))\n", + "L2 = NumpyMatrixOperator(np.eye(5))\n", + "L = L1 + ProjectionParameterFunctional('p', ()) * L2\n", + "F = L.source.ones(1)\n", + "m = StationaryModel(L, F)" + ] + }, + { + "cell_type": "code", + "execution_count": 55, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "LincombOperator(\n", + " (NumpyMatrixOperator(<5x5 dense>), NumpyMatrixOperator(<5x5 dense>)),\n", + " (1.0, ProjectionParameterFunctional('p', ())))" + ] + }, + "execution_count": 55, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "L" + ] + }, + { + "cell_type": "code", + "execution_count": 56, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "StationaryModel(\n", + " LincombOperator(\n", + " (NumpyMatrixOperator(<5x5 dense>), NumpyMatrixOperator(<5x5 dense>)),\n", + " (1.0, ProjectionParameterFunctional('p', ()))),\n", + " VectorOperator(NumpyVectorArray([[1. 1. 1. 1. 1.]], NumpyVectorSpace(5)), name='rhs'),\n", + " products={})" + ] + }, + "execution_count": 56, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "m" + ] + }, + { + "cell_type": "code", + "execution_count": 57, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "NumpyVectorSpace(5)" + ] + }, + "execution_count": 57, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "m.solution_space" + ] + }, + { + "cell_type": "code", + "execution_count": 58, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "00:49 StationaryModel: Solving StationaryModel for {p: 1} ...\n" + ] + }, + { + "data": { + "text/plain": [ + "NumpyVectorArray([[0.5 0.33333333 0.25 0.2 0.16666667]], NumpyVectorSpace(5))" + ] + }, + "execution_count": 58, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "m.solve(mu=Parameter({'p': 1}))" + ] + }, + { + "cell_type": "code", + "execution_count": 59, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "00:51 StationaryModel: Solving StationaryModel for {p: 1} ...\n" + ] + }, + { + "data": { + "text/plain": [ + "NumpyVectorArray([[0.5 0.33333333 0.25 0.2 0.16666667]], NumpyVectorSpace(5))" + ] + }, + "execution_count": 59, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "m.solve(1)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Algorithms (`pymor.algorithms`)\n", + "- `gram_schmidt`\n", + "- `genericsolvers`, `newton`\n", + "- `pod`, `hapod`\n", + "- `ei`, `greedy`, `adaptivegreedy`\n", + "- `lyapunov`, `lradi`, `riccati`, `lrradi`, `sylvester`\n", + "- `timestepping`\n", + "- `to_matrix`\n", + "- ..." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Modified Gram-Schmidt using NumPy:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "for i in range(A.shape[1]):\n", + " for j in range(i):\n", + " p = A[:, i] @ P @ A[:, j]\n", + " A[:, i] -= p * A[:, j]\n", + " A[:, i] /= np.sqrt(A[:, i] @ P @ A[:, i])" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Modified Gram-Schmidt using pyMOR:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "for i in range(len(A)):\n", + " for j in range(i):\n", + " p = A[j].inner(A[i], product=P)[0, 0]\n", + " A[i].axpy(-p, A[j])\n", + " A[i].scal(1 / A[i].norm(product=P))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Reductors (`pymor.reductors`)\n", + "- object with a `reduce` method which returns a reduced-order model\n", + "- other reduction data, e.g. basis matrices, are attributes of the reductor\n", + " object\n", + "- modules:\n", + " - `basic`\n", + " - `coercive`, `parabolic`, `residual` (reduced basis)\n", + " - `bt`, `h2`, `interpolation`, `sobt`, `sor_irka` (system theory)" + ] + }, + { + "cell_type": "code", + "execution_count": 60, + "metadata": {}, + "outputs": [], + "source": [ + "from pymor.basic import StationaryRBReductor, gram_schmidt" + ] + }, + { + "cell_type": "code", + "execution_count": 61, + "metadata": {}, + "outputs": [], + "source": [ + "StationaryRBReductor?" + ] + }, + { + "cell_type": "code", + "execution_count": 62, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "StationaryModel(\n", + " LincombOperator(\n", + " (NumpyMatrixOperator(<5x5 dense>), NumpyMatrixOperator(<5x5 dense>)),\n", + " (1.0, ProjectionParameterFunctional('p', ()))),\n", + " VectorOperator(NumpyVectorArray([[1. 1. 1. 1. 1.]], NumpyVectorSpace(5)), name='rhs'),\n", + " products={})" + ] + }, + "execution_count": 62, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "m" + ] + }, + { + "cell_type": "code", + "execution_count": 63, + "metadata": {}, + "outputs": [], + "source": [ + "RB = gram_schmidt(m.solution_space.random(2))" + ] + }, + { + "cell_type": "code", + "execution_count": 64, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "NumpyVectorArray(\n", + " [[ 0.26734547 0.67861666 0.5224948 0.4273204 0.11136558]\n", + " [-0.0812647 -0.56548911 0.43790783 0.2440862 0.64982827]],\n", + " NumpyVectorSpace(5))" + ] + }, + "execution_count": 64, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "RB" + ] + }, + { + "cell_type": "code", + "execution_count": 65, + "metadata": {}, + "outputs": [], + "source": [ + "reductor = StationaryRBReductor(m, RB)" + ] + }, + { + "cell_type": "code", + "execution_count": 66, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "01:00 StationaryRBReductor: Operator projection ...\n", + "01:00 StationaryRBReductor: Building ROM ...\n" + ] + }, + { + "data": { + "text/plain": [ + "StationaryModel(\n", + " LincombOperator(\n", + " (NumpyMatrixOperator(<2x2 dense>), NumpyMatrixOperator(<2x2 dense>)),\n", + " (1.0, ProjectionParameterFunctional('p', ()))),\n", + " NumpyMatrixOperator(<2x1 dense>, name='rhs'),\n", + " products={},\n", + " name='StationaryModel_reduced')" + ] + }, + "execution_count": 66, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "reductor.reduce()" + ] + }, + { + "cell_type": "code", + "execution_count": 67, + "metadata": {}, + "outputs": [], + "source": [ + "from pymor.basic import rb_greedy" + ] + }, + { + "cell_type": "code", + "execution_count": 68, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "01:00 weak_greedy: Started greedy search on training set of size 10.\n", + "01:00 weak_greedy: Estimating errors ...\n", + "01:00 | RBSurrogate: Reducing ...\n", + "01:00 | | StationaryRBReductor: Operator projection ...\n", + "01:00 | | StationaryRBReductor: Building ROM ...\n", + "01:00 | StationaryModel: Solving StationaryModel for {p: 0.0} ...\n", + "01:00 | StationaryModel: Solving StationaryModel for {p: 0.1111111111111111} ...\n", + "01:00 | StationaryModel: Solving StationaryModel for {p: 0.2222222222222222} ...\n", + "01:00 | StationaryModel: Solving StationaryModel for {p: 0.3333333333333333} ...\n", + "01:00 | StationaryModel: Solving StationaryModel for {p: 0.4444444444444444} ...\n", + "01:00 | StationaryModel: Solving StationaryModel for {p: 0.5555555555555556} ...\n", + "01:00 | StationaryModel: Solving StationaryModel for {p: 0.6666666666666666} ...\n", + "01:00 | StationaryModel: Solving StationaryModel for {p: 0.7777777777777777} ...\n", + "01:00 | StationaryModel: Solving StationaryModel for {p: 0.8888888888888888} ...\n", + "01:00 | StationaryModel: Solving StationaryModel for {p: 1.0} ...\n", + "01:00 weak_greedy: Maximum error after 0 extensions: 1.209797962930634 (mu = 0.0)\n", + "01:00 weak_greedy: Extending surrogate for mu = 0.0 ...\n", + "01:00 | RBSurrogate: Computing solution snapshot for mu = 0.0 ...\n", + "01:00 | | StationaryModel: Solving StationaryModel for {p: 0.0} ...\n", + "01:00 | RBSurrogate: Extending basis with solution snapshot ...\n", + "01:00 | RBSurrogate: Reducing ...\n", + "01:00 | | StationaryRBReductor: Operator projection ...\n", + "01:00 | | StationaryRBReductor: Building ROM ...\n", + " \n", + "01:00 weak_greedy: Estimating errors ...\n", + "01:00 | StationaryModel: Solving StationaryModel for {p: 0.0} ...\n", + "01:00 | StationaryModel: Solving StationaryModel for {p: 0.1111111111111111} ...\n", + "01:00 | StationaryModel: Solving StationaryModel for {p: 0.2222222222222222} ...\n", + "01:00 | StationaryModel: Solving StationaryModel for {p: 0.3333333333333333} ...\n", + "01:00 | StationaryModel: Solving StationaryModel for {p: 0.4444444444444444} ...\n", + "01:00 | StationaryModel: Solving StationaryModel for {p: 0.5555555555555556} ...\n", + "01:00 | StationaryModel: Solving StationaryModel for {p: 0.6666666666666666} ...\n", + "01:00 | StationaryModel: Solving StationaryModel for {p: 0.7777777777777777} ...\n", + "01:00 | StationaryModel: Solving StationaryModel for {p: 0.8888888888888888} ...\n", + "01:00 | StationaryModel: Solving StationaryModel for {p: 1.0} ...\n", + "01:00 weak_greedy: Maximum error after 1 extensions: 0.13877814174930383 (mu = 1.0)\n", + "01:00 weak_greedy: Extending surrogate for mu = 1.0 ...\n", + "01:00 | RBSurrogate: Computing solution snapshot for mu = 1.0 ...\n", + "01:00 | | StationaryModel: Solving StationaryModel for {p: 1.0} ...\n", + "01:00 | RBSurrogate: Extending basis with solution snapshot ...\n", + "01:00 | RBSurrogate: Reducing ...\n", + "01:00 | | StationaryRBReductor: Operator projection ...\n", + "01:00 | | StationaryRBReductor: Building ROM ...\n", + " \n", + "01:00 weak_greedy: Maximum number of 2 extensions reached.\n", + "01:00 weak_greedy: Greedy search took 0.03968214988708496 seconds\n" + ] + }, + { + "data": { + "text/plain": [ + "{'max_errs': [1.209797962930634, 0.13877814174930383],\n", + " 'max_err_mus': [0.0, 1.0],\n", + " 'extensions': 2,\n", + " 'time': 0.03968214988708496,\n", + " 'rom': StationaryModel(\n", + " LincombOperator(\n", + " (NumpyMatrixOperator(<2x2 dense>), NumpyMatrixOperator(<2x2 dense>)),\n", + " (1.0, ProjectionParameterFunctional('p', ()))),\n", + " NumpyMatrixOperator(<2x1 dense>, name='rhs'),\n", + " products={},\n", + " name='StationaryModel_reduced')}" + ] + }, + "execution_count": 68, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "rb_greedy(\n", + " m,\n", + " StationaryRBReductor(m),\n", + " np.linspace(0, 1, 10),\n", + " use_estimator=False,\n", + " max_extensions=2,\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": 69, + "metadata": {}, + "outputs": [], + "source": [ + "from pymor.basic import BTReductor" + ] + }, + { + "cell_type": "code", + "execution_count": 70, + "metadata": {}, + "outputs": [], + "source": [ + "bt = BTReductor(lti)" + ] + }, + { + "cell_type": "code", + "execution_count": 71, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "01:03 LTIPGReductor: Operator projection ...\n", + "01:03 LTIPGReductor: Building ROM ...\n" + ] + } + ], + "source": [ + "lti_rom = bt.reduce(2)" + ] + }, + { + "cell_type": "code", + "execution_count": 72, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "LTIModel(\n", + " NumpyMatrixOperator(<2x2 dense>),\n", + " NumpyMatrixOperator(<2x1 dense>),\n", + " NumpyMatrixOperator(<1x2 dense>),\n", + " D=ZeroOperator(NumpyVectorSpace(1), NumpyVectorSpace(1)),\n", + " E=NumpyMatrixOperator(<2x2 dense>, name='IdentityOperator'),\n", + " name='LTIModel_reduced')" + ] + }, + "execution_count": 72, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "lti_rom" + ] + }, + { + "cell_type": "code", + "execution_count": 73, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "NumpyVectorArray(\n", + " [[ 0.66177775 0.48450457 0.38528025 0.32091271 0.27546268]\n", + " [-0.67242346 0.08147675 0.34253034 0.44237423 0.47768843]],\n", + " NumpyVectorSpace(5, id='STATE'))" + ] + }, + "execution_count": 73, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "bt.V" + ] + }, + { + "cell_type": "code", + "execution_count": 74, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "0.003612034902573163" + ] + }, + "execution_count": 74, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "(lti - lti_rom).h2_norm() / lti.h2_norm()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Bindings (`pymor.bindings`)\n", + "- PDE solvers:\n", + " `FEniCS`, `NGSolve`\n", + "- linear equations solvers:\n", + " `SciPy`, `PyAMG`\n", + "- matrix equations solvers:\n", + " `SciPy`, `Slycot`, `Py-M.E.S.S.`" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Minimal C++ demo (`pymor/src/pymordemos/minimal_cpp_demo`)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "lines_to_next_cell": 2 + }, + "outputs": [], + "source": [ + "# %load minimal_cpp_demo/model.hh\n", + "#ifndef DISCRETIZATION_HH\n", + "#define DISCRETIZATION_HH\n", + "\n", + "#include <vector>\n", + "\n", + "\n", + "class Vector {\n", + " friend class DiffusionOperator;\n", + "public:\n", + " Vector(int dim, double value);\n", + " Vector(const Vector& other);\n", + " const int dim;\n", + " void scal(double val);\n", + " void axpy(double a, const Vector& x);\n", + " double dot(const Vector& other) const;\n", + " double* data();\n", + "private:\n", + " std::vector<double> _data;\n", + "};\n", + "\n", + "\n", + "class DiffusionOperator {\n", + "public:\n", + " DiffusionOperator(int n, double left, double right);\n", + " const int dim_source;\n", + " const int dim_range;\n", + " void apply(const Vector& U, Vector& R) const;\n", + "private:\n", + " const double h;\n", + " const double left;\n", + " const double right;\n", + "};\n", + "\n", + "\n", + "#endif // DISCRETIZATION_HH" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "lines_to_next_cell": 2 + }, + "outputs": [], + "source": [ + "# %load minimal_cpp_demo/model.cc\n", + "#include \"model.hh\"\n", + "\n", + "#include <assert.h>\n", + "#include <iostream>\n", + "#include <cmath>\n", + "\n", + "Vector::Vector(int dim, double value) : _data(dim, value), dim(dim) {}\n", + "\n", + "Vector::Vector(const Vector& other) : _data(other._data), dim(other.dim) {}\n", + "\n", + "void Vector::scal(double val) {\n", + " for (int i = 0; i < dim; i++) {\n", + " _data[i] *= val;\n", + " }\n", + "}\n", + "\n", + "void Vector::axpy(double a, const Vector& x) {\n", + " assert(x.dim == dim);\n", + " for (int i = 0; i < dim; i++) {\n", + " _data[i] += a * x._data[i];\n", + " }\n", + "}\n", + "\n", + "double Vector::dot(const Vector& other) const {\n", + " assert(other.dim == dim);\n", + " double result = 0;\n", + " for (int i = 0; i < dim; i++) {\n", + " result += _data[i] * other._data[i];\n", + " }\n", + " return result;\n", + "}\n", + "\n", + "double* Vector::data() {\n", + " return _data.data();\n", + "}\n", + "\n", + "// -------------------------------------------------------\n", + "\n", + "DiffusionOperator::DiffusionOperator(int n, double left, double right)\n", + " : dim_source(n+1), dim_range(n+1), h(1./n), left(left), right(right) {}\n", + "\n", + "void DiffusionOperator::apply(const Vector& U, Vector& R) const {\n", + " const double one_over_h2 = 1. / (h * h);\n", + " for (int i = 1; i < dim_range - 1; i++) {\n", + " if ((i * h > left) && (i * h <= right)) {\n", + " R._data[i] = -(U._data[i-1] - 2*U._data[i] + U._data[i+1]) * one_over_h2;\n", + " }\n", + " }\n", + "}\n", + "\n", + "#include <pybind11/functional.h>\n", + "#include <pybind11/numpy.h>\n", + "#include <pybind11/operators.h>\n", + "#include <pybind11/pybind11.h>\n", + "\n", + "PYBIND11_MODULE(model, m)\n", + "{\n", + " namespace py = pybind11;\n", + "\n", + " py::class_<DiffusionOperator> op(m, \"DiffusionOperator\",\n", + " \"DiffusionOperator Docstring\");\n", + " op.def(py::init([](int n, double left, double right) {\n", + " return std::make_unique<DiffusionOperator>(n, left, right);\n", + " }));\n", + " op.def_readonly(\"dim_source\", &DiffusionOperator::dim_source);\n", + " op.def_readonly(\"dim_range\", &DiffusionOperator::dim_range);\n", + " op.def(\"apply\", &DiffusionOperator::apply);\n", + "\n", + " py::class_<Vector> vec(m, \"Vector\", \"Vector Docstring\",\n", + " py::buffer_protocol());\n", + " vec.def(py::init([](int size, double value) {\n", + " return std::make_unique<Vector>(size, value); }));\n", + " vec.def(py::init([](const Vector& other) {\n", + " return std::make_unique<Vector>(other); }));\n", + "\n", + " vec.def_readonly(\"dim\", &Vector::dim);\n", + " vec.def(\"scal\", &Vector::scal);\n", + " vec.def(\"axpy\", &Vector::axpy);\n", + " vec.def(\"dot\", &Vector::dot);\n", + " vec.def(\"data\", &Vector::data);\n", + "\n", + " vec.def_buffer([](Vector& vec) -> py::buffer_info {\n", + " return py::buffer_info(\n", + " vec.data(), sizeof(double), py::format_descriptor<double>::format(),\n", + " 1, {vec.dim}, {sizeof(double)});\n", + " });\n", + "}" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# %load minimal_cpp_demo/wrapper.py\n", + "# This file is part of the pyMOR project (http://www.pymor.org).\n", + "# Copyright 2013-2019 pyMOR developers and contributors. All rights reserved.\n", + "# License: BSD 2-Clause License (http://opensource.org/licenses/BSD-2-Clause)\n", + "\n", + "from pymor.operators.basic import OperatorBase\n", + "from pymor.vectorarrays.list import CopyOnWriteVector, ListVectorSpace\n", + "\n", + "import numpy as np\n", + "import math\n", + "from model import Vector, DiffusionOperator\n", + "\n", + "\n", + "class WrappedVector(CopyOnWriteVector):\n", + "\n", + " def __init__(self, vector):\n", + " assert isinstance(vector, Vector)\n", + " self._impl = vector\n", + "\n", + " @classmethod\n", + " def from_instance(cls, instance):\n", + " return cls(instance._impl)\n", + "\n", + " def to_numpy(self, ensure_copy=False):\n", + " result = np.frombuffer(self._impl, dtype=np.float)\n", + " if ensure_copy:\n", + " result = result.copy()\n", + " return result\n", + "\n", + " def _copy_data(self):\n", + " self._impl = Vector(self._impl)\n", + "\n", + " def _scal(self, alpha):\n", + " self._impl.scal(alpha)\n", + "\n", + " def _axpy(self, alpha, x):\n", + " self._impl.axpy(alpha, x._impl)\n", + "\n", + " def dot(self, other):\n", + " return self._impl.dot(other._impl)\n", + "\n", + " def l1_norm(self):\n", + " raise NotImplementedError\n", + "\n", + " def l2_norm(self):\n", + " return math.sqrt(self.dot(self))\n", + "\n", + " def l2_norm2(self):\n", + " return self.dot(self)\n", + "\n", + " def sup_norm(self):\n", + " raise NotImplementedError\n", + "\n", + " def dofs(self, dof_indices):\n", + " raise NotImplementedError\n", + "\n", + " def amax(self):\n", + " raise NotImplementedError\n", + "\n", + "\n", + "class WrappedVectorSpace(ListVectorSpace):\n", + "\n", + " def __init__(self, dim):\n", + " self.dim = dim\n", + "\n", + " def zero_vector(self):\n", + " return WrappedVector(Vector(self.dim, 0))\n", + "\n", + " def make_vector(self, obj):\n", + " return WrappedVector(obj)\n", + "\n", + " def __eq__(self, other):\n", + " return type(other) is WrappedVectorSpace and self.dim == other.dim\n", + "\n", + "\n", + "class WrappedDiffusionOperator(OperatorBase):\n", + " def __init__(self, op):\n", + " assert isinstance(op, DiffusionOperator)\n", + " self.op = op\n", + " self.source = WrappedVectorSpace(op.dim_source)\n", + " self.range = WrappedVectorSpace(op.dim_range)\n", + " self.linear = True\n", + "\n", + " @classmethod\n", + " def create(cls, n, left, right):\n", + " return cls(DiffusionOperator(n, left, right))\n", + "\n", + " def apply(self, U, mu=None):\n", + " assert U in self.source\n", + "\n", + " def apply_one_vector(u):\n", + " v = Vector(self.range.dim, 0)\n", + " self.op.apply(u._impl, v)\n", + " return v\n", + "\n", + " return self.range.make_array([apply_one_vector(u) for u in U._list])\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "lines_to_next_cell": 2 + }, + "outputs": [], + "source": [ + "# %load minimal_cpp_demo/demo.py\n", + "# This file is part of the pyMOR project (http://www.pymor.org).\n", + "# Copyright 2013-2019 pyMOR developers and contributors. All rights reserved.\n", + "# License: BSD 2-Clause License (http://opensource.org/licenses/BSD-2-Clause)\n", + "\n", + "import numpy as np\n", + "\n", + "from pymor.algorithms.pod import pod\n", + "from pymor.algorithms.timestepping import ExplicitEulerTimeStepper\n", + "from pymor.models.basic import InstationaryModel\n", + "from pymor.grids.oned import OnedGrid\n", + "from pymor.gui.visualizers import OnedVisualizer\n", + "from pymor.operators.constructions import VectorOperator, LincombOperator\n", + "from pymor.parameters.functionals import ProjectionParameterFunctional\n", + "from pymor.parameters.spaces import CubicParameterSpace\n", + "from pymor.reductors.basic import InstationaryRBReductor\n", + "\n", + "# import wrapped classes\n", + "from pymordemos.minimal_cpp_demo.wrapper import WrappedDiffusionOperator\n", + "\n", + "\n", + "def discretize(n, nt, blocks):\n", + " h = 1. / blocks\n", + " ops = [WrappedDiffusionOperator.create(n, h * i, h * (i + 1))\n", + " for i in range(blocks)]\n", + " pfs = [ProjectionParameterFunctional('diffusion_coefficients',\n", + " (blocks,),\n", + " (i,))\n", + " for i in range(blocks)]\n", + " operator = LincombOperator(ops, pfs)\n", + "\n", + " initial_data = operator.source.zeros()\n", + "\n", + " # use data property of WrappedVector to setup rhs\n", + " # note that we cannot use the data property of ListVectorArray,\n", + " # since ListVectorArray will always return a copy\n", + " rhs_vec = operator.range.zeros()\n", + " rhs_data = rhs_vec._list[0].to_numpy()\n", + " rhs_data[:] = np.ones(len(rhs_data))\n", + " rhs_data[0] = 0\n", + " rhs_data[len(rhs_data) - 1] = 0\n", + " rhs = VectorOperator(rhs_vec)\n", + "\n", + " # hack together a visualizer ...\n", + " grid = OnedGrid(domain=(0, 1), num_intervals=n)\n", + " visualizer = OnedVisualizer(grid)\n", + "\n", + " time_stepper = ExplicitEulerTimeStepper(nt)\n", + " parameter_space = CubicParameterSpace(operator.parameter_type, 0.1, 1)\n", + "\n", + " fom = InstationaryModel(T=1e-0, operator=operator, rhs=rhs,\n", + " initial_data=initial_data,\n", + " time_stepper=time_stepper, num_values=20,\n", + " parameter_space=parameter_space,\n", + " visualizer=visualizer, name='C++-Model',\n", + " cache_region=None)\n", + " return fom\n", + "\n", + "\n", + "# discretize\n", + "fom = discretize(50, 10000, 4)\n", + "\n", + "# generate solution snapshots\n", + "snapshots = fom.solution_space.empty()\n", + "for mu in fom.parameter_space.sample_uniformly(2):\n", + " snapshots.append(fom.solve(mu))\n", + "\n", + "# apply POD\n", + "reduced_basis = pod(snapshots, modes=4)[0]\n", + "\n", + "# reduce the model\n", + "reductor = InstationaryRBReductor(fom, reduced_basis,\n", + " check_orthonormality=True)\n", + "rom = reductor.reduce()\n", + "\n", + "# stochastic error estimation\n", + "mu_max = None\n", + "err_max = -1.\n", + "for mu in fom.parameter_space.sample_randomly(10):\n", + " U_RB = (reductor.reconstruct(rom.solve(mu)))\n", + " U = fom.solve(mu)\n", + " err = np.max((U_RB-U).l2_norm())\n", + " if err > err_max:\n", + " err_max = err\n", + " mu_max = mu\n", + "\n", + "# visualize maximum error solution\n", + "U_RB = (reductor.reconstruct(rom.solve(mu_max)))\n", + "U = fom.solve(mu_max)\n", + "\n", + "# avoid blocking in CI run\n", + "if __name__ == '__main__':\n", + " fom.visualize((U_RB, U), title=f'mu = {mu}', legend=('reduced',\n", + " 'detailed'))\n", + " fom.visualize((U-U_RB), title=f'mu = {mu}', legend=('error'))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Getting Started with pyMOR\n", + "- home page:\n", + " https://pymor.org\n", + "- GitHub page:\n", + " https://github.com/pymor/pymor\n", + "- documentation:\n", + " https://docs.pymor.org\n", + "- pyMOR School material:\n", + " https://zivgitlab.uni-muenster.de/pymor/school-material\n", + "- pyMOR School 2020 (Sep 28 - Oct 2):\n", + " https://school.pymor.org" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Contributing to pyMOR\n", + "- for questions, discussions, bug reports, feature requests, etc.:\n", + " - GitHub issues:\n", + " https://github.com/pymor/pymor/issues\n", + " - pyMOR-dev mailing list:\n", + " http://listserv.uni-muenster.de/mailman/listinfo/pymor-dev\n", + "- for code:\n", + " - development guide\n", + " https://github.com/pymor/pymor/blob/master/README.md#setting-up-an-environment-for-pymor-development\n", + " - contribution guide\n", + " https://github.com/pymor/pymor/blob/master/CONTRIBUTING.md\n", + " - fork the project, work on a new branch, and submit a pull request\n", + " https://github.com/pymor/pymor/pulls\n", + " - by contributing, you agree to your work being published under pyMOR's\n", + " license (and confirming that you have copyright over your code or\n", + " permission from the copyright holder)\n", + " https://github.com/pymor/pymor/blob/master/CONTRIBUTING.md#license\n", + " - attribution\n", + " https://github.com/pymor/pymor/blob/master/CONTRIBUTING.md#attribution\n", + " - becoming a main developer (i.e., owner of the GitHub pyMOR organization)\n", + " https://github.com/pymor/pymor/blob/master/CONTRIBUTING.md#becoming-a-main-developer" + ] + } + ], + "metadata": { + "jupytext": { + "encoding": "# -*- coding: utf-8 -*-", + "formats": "ipynb,py:percent" + }, + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.6.8" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +}