diff --git a/rb-intro/Sheet3.ipynb b/rb-intro/Sheet3.ipynb
new file mode 100644
index 0000000000000000000000000000000000000000..db101dd2bad013f2acc45f56deace1ec9ac85cee
--- /dev/null
+++ b/rb-intro/Sheet3.ipynb
@@ -0,0 +1,347 @@
+{
+ "cells": [
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "# pyMOR RB-Intro\n",
+ "\n",
+ "# Exercise Sheet 3 - Solutions"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 1,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "from pymor.basic import *\n",
+ "set_log_levels({'pymor': 'WARN',\n",
+ " 'pymor.functions.bitmap': 'CRITICAL'})"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "We will continue to work with the models defined on exercise sheet 1. To simplify working with these models, you can find their definition in `problems.py` in the same directory as this notebook:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 2,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "from pymor.basic import *\r\n",
+ "\r\n",
+ "\r\n",
+ "def problem_1_1_a():\r\n",
+ " return StationaryProblem(\r\n",
+ " domain=RectDomain([[0., 0.], [1., 1.]]),\r\n",
+ "\r\n",
+ " diffusion=ConstantFunction(1, 2),\r\n",
+ "\r\n",
+ " rhs=ExpressionFunction('(sqrt( (x[...,0]-0.5)**2 + (x[...,1]-0.5)**2) <= 0.3) * 1.', 2, ()),\r\n",
+ "\r\n",
+ " name='Sheet1_Problem_1a'\r\n",
+ " )\r\n",
+ "\r\n",
+ "\r\n",
+ "def problem_1_1_b():\r\n",
+ " return StationaryProblem(\r\n",
+ " domain=RectDomain(bottom='neumann'),\r\n",
+ "\r\n",
+ " neumann_data=ConstantFunction(-1., 2),\r\n",
+ "\r\n",
+ " diffusion=ExpressionFunction('1. - (sqrt( (x[...,0]-0.5)**2 + (x[...,1]-0.5)**2) <= 0.3) * 0.999', 2, ()),\r\n",
+ "\r\n",
+ " name='Sheet1_Problem_1b'\r\n",
+ " )\r\n",
+ "\r\n",
+ "\r\n",
+ "def problem_1_1_c():\r\n",
+ " return StationaryProblem(\r\n",
+ " domain=RectDomain(bottom='neumann'),\r\n",
+ "\r\n",
+ " neumann_data=ConstantFunction(-1., 2),\r\n",
+ "\r\n",
+ " diffusion=ExpressionFunction(\r\n",
+ " '1. - (sqrt( (np.mod(x[...,0],1./K)-0.5/K)**2 + (np.mod(x[...,1],1./K)-0.5/K)**2) <= 0.3/K) * 0.999',\r\n",
+ " 2, (),\r\n",
+ " values={'K': 10}\r\n",
+ " ),\r\n",
+ "\r\n",
+ " name='Sheet1_Problem_1c'\r\n",
+ " )\r\n",
+ "\r\n",
+ "\r\n",
+ "def problem_1_1_d():\r\n",
+ " return StationaryProblem(\r\n",
+ " domain=RectDomain(bottom='neumann'),\r\n",
+ "\r\n",
+ " neumann_data=ConstantFunction(-1., 2),\r\n",
+ "\r\n",
+ " diffusion=BitmapFunction('RB.png', range=[0.001, 1]),\r\n",
+ "\r\n",
+ " name='Sheet1_Problem_1d'\r\n",
+ " )\r\n",
+ "\r\n",
+ "\r\n",
+ "def problem_1_2_a():\r\n",
+ " return StationaryProblem(\r\n",
+ " domain=RectDomain(bottom='neumann'),\r\n",
+ "\r\n",
+ " neumann_data=ExpressionFunction('-cos(pi*x[...,0])**2*neum', 2, (), parameter_type={'neum': ()}),\r\n",
+ "\r\n",
+ " diffusion=ExpressionFunction(\r\n",
+ " '1. - (sqrt( (np.mod(x[...,0],1./K)-0.5/K)**2 + (np.mod(x[...,1],1./K)-0.5/K)**2) <= 0.3/K) * 0.999',\r\n",
+ " 2, (),\r\n",
+ " values={'K': 10}\r\n",
+ " ),\r\n",
+ "\r\n",
+ " parameter_space=CubicParameterSpace({'neum': ()}, -10, 10),\r\n",
+ "\r\n",
+ " name='Sheet1_Problem_2a'\r\n",
+ " )\r\n",
+ "\r\n",
+ "\r\n",
+ "def problem_1_2_b():\r\n",
+ " return StationaryProblem(\r\n",
+ " domain=RectDomain(bottom='neumann'),\r\n",
+ "\r\n",
+ " neumann_data=ExpressionFunction('-cos(pi*x[...,0])**2*neum', 2, (), parameter_type={'neum': ()}),\r\n",
+ "\r\n",
+ " diffusion=ExpressionFunction(\r\n",
+ " '1. - (sqrt( (np.mod(x[...,0],1./K)-0.5/K)**2 + (np.mod(x[...,1],1./K)-0.5/K)**2) <= 0.3/K) * (1 - diffu)',\r\n",
+ " 2, (),\r\n",
+ " values={'K': 10},\r\n",
+ " parameter_type={'diffu': ()}\r\n",
+ " ),\r\n",
+ "\r\n",
+ " parameter_space=CubicParameterSpace(\r\n",
+ " {'neum': (), 'diffu': ()},\r\n",
+ " ranges={'diffu': [0.0001, 1.], 'neum': [-10, 10]}\r\n",
+ " ),\r\n",
+ "\r\n",
+ " name='Sheet1_Problem_2b'\r\n",
+ " )\r\n",
+ "\r\n",
+ "\r\n",
+ "def problem_1_2_c():\r\n",
+ " return StationaryProblem(\r\n",
+ " domain=RectDomain(bottom='neumann'),\r\n",
+ "\r\n",
+ " neumann_data=ConstantFunction(-1., 2),\r\n",
+ "\r\n",
+ " diffusion=LincombFunction(\r\n",
+ " [ConstantFunction(1., 2),\r\n",
+ " BitmapFunction('R.png', range=[1, 0]),\r\n",
+ " BitmapFunction('B.png', range=[1, 0])],\r\n",
+ " [1.,\r\n",
+ " ExpressionParameterFunctional('-(1 - R)', {'R': ()}),\r\n",
+ " ExpressionParameterFunctional('-(1 - B)', {'B': ()})]\r\n",
+ " ),\r\n",
+ "\r\n",
+ " parameter_space=CubicParameterSpace({'R': (), 'B': ()}, 0.0001, 1.),\r\n",
+ "\r\n",
+ " name='Sheet1_Problem_2c'\r\n",
+ " )\r\n"
+ ]
+ }
+ ],
+ "source": [
+ "%cat problems.py"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "The file can be imported as usual:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 3,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import problems"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "We have four non-parametric problems:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "scrolled": false
+ },
+ "outputs": [],
+ "source": [
+ "for f in [problems.problem_1_1_a, problems.problem_1_1_b,\n",
+ " problems.problem_1_1_c, problems.problem_1_1_d]:\n",
+ " p = f()\n",
+ " m, _ = discretize_stationary_cg(p, diameter=1/100)\n",
+ " m.visualize(m.solve(), legend=m.name)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "The parametric problems have already been assigned a `ParameterSpace` which is inherited by the generated models:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "for f in [problems.problem_1_2_a, problems.problem_1_2_b,\n",
+ " problems.problem_1_2_c]:\n",
+ " p = f()\n",
+ " m, _ = discretize_stationary_cg(p, diameter=1/100)\n",
+ " mu = m.parameter_space.sample_randomly(1)[0]\n",
+ " m.visualize(m.solve(mu), legend=str(mu))"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Problem 1"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Problem 1 a)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Write a method\n",
+ "\n",
+ "```python\n",
+ "def orthogonal_projection(U, basis, product=None,\n",
+ " basis_is_orthonormal=False,\n",
+ " return_projection_matrix=False):\n",
+ " ...\n",
+ "```\n",
+ "that encapsulates the orthogonal projection algorithm you implemented in Sheet 2, Problem 2."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Problem 1 b)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Repeat the error calculations from Sheet 2, Problem 2 with an orthonormal basis obtained from `U`\n",
+ "using `pymor.algorithms.gram_schmidt`. Compare your results."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Problem 2"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Problem 2 a)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Extend the `greedy` function from Sheet 2, Problem 3 to a function\n",
+ "\n",
+ "```python\n",
+ " def greedy(U, basis_size, product=None,\n",
+ " rtol=0,\n",
+ " orthonormalize=False):\n",
+ " ...\n",
+ "```\n",
+ "\n",
+ "The algorithm should stop, when the relative approximation error goes below the prescribed value `rtol`. When `orthonormalize` is set to `True`, the algorithm shall return an orthonormal basis w.r.t. to `product` by orthonormalizing each new snapshot vector w.r.t. `product`.\n",
+ "\n",
+ "**Hint:** When using `gram_schmidt` for orthonormalization, the `offset` parameter can be used to avoid re-orthonormalizing the old basis. "
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Problem 2 b)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Repeat the experiments from Sheet 2, Problem 3, using `strong_greedy` with `orthonormalize=True`."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Bonus Problem"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Write a function `word_problem` which, given an arbitrary string, generates a `StationaryProblem` similar to Sheet 1, Problem 2 (c). Each letter should have a corresponding parameter component.\n",
+ "\n",
+ "**Hint:** Use the modules `PIL.ImageDraw` und `PIL.ImageFont`."
+ ]
+ }
+ ],
+ "metadata": {
+ "kernelspec": {
+ "display_name": "Python 3",
+ "language": "python",
+ "name": "python3"
+ },
+ "language_info": {
+ "codemirror_mode": {
+ "name": "ipython",
+ "version": 3
+ },
+ "file_extension": ".py",
+ "mimetype": "text/x-python",
+ "name": "python",
+ "nbconvert_exporter": "python",
+ "pygments_lexer": "ipython3",
+ "version": "3.7.2"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}
diff --git a/rb-intro/problems.py b/rb-intro/problems.py
new file mode 100644
index 0000000000000000000000000000000000000000..79b65e5e3f73b9fcb9aec99851d02074b8cece42
--- /dev/null
+++ b/rb-intro/problems.py
@@ -0,0 +1,114 @@
+from pymor.basic import *
+
+
+def problem_1_1_a():
+ return StationaryProblem(
+ domain=RectDomain([[0., 0.], [1., 1.]]),
+
+ diffusion=ConstantFunction(1, 2),
+
+ rhs=ExpressionFunction('(sqrt( (x[...,0]-0.5)**2 + (x[...,1]-0.5)**2) <= 0.3) * 1.', 2, ()),
+
+ name='Sheet1_Problem_1a'
+ )
+
+
+def problem_1_1_b():
+ return StationaryProblem(
+ domain=RectDomain(bottom='neumann'),
+
+ neumann_data=ConstantFunction(-1., 2),
+
+ diffusion=ExpressionFunction('1. - (sqrt( (x[...,0]-0.5)**2 + (x[...,1]-0.5)**2) <= 0.3) * 0.999', 2, ()),
+
+ name='Sheet1_Problem_1b'
+ )
+
+
+def problem_1_1_c():
+ return StationaryProblem(
+ domain=RectDomain(bottom='neumann'),
+
+ neumann_data=ConstantFunction(-1., 2),
+
+ diffusion=ExpressionFunction(
+ '1. - (sqrt( (np.mod(x[...,0],1./K)-0.5/K)**2 + (np.mod(x[...,1],1./K)-0.5/K)**2) <= 0.3/K) * 0.999',
+ 2, (),
+ values={'K': 10}
+ ),
+
+ name='Sheet1_Problem_1c'
+ )
+
+
+def problem_1_1_d():
+ return StationaryProblem(
+ domain=RectDomain(bottom='neumann'),
+
+ neumann_data=ConstantFunction(-1., 2),
+
+ diffusion=BitmapFunction('RB.png', range=[0.001, 1]),
+
+ name='Sheet1_Problem_1d'
+ )
+
+
+def problem_1_2_a():
+ return StationaryProblem(
+ domain=RectDomain(bottom='neumann'),
+
+ neumann_data=ExpressionFunction('-cos(pi*x[...,0])**2*neum', 2, (), parameter_type={'neum': ()}),
+
+ diffusion=ExpressionFunction(
+ '1. - (sqrt( (np.mod(x[...,0],1./K)-0.5/K)**2 + (np.mod(x[...,1],1./K)-0.5/K)**2) <= 0.3/K) * 0.999',
+ 2, (),
+ values={'K': 10}
+ ),
+
+ parameter_space=CubicParameterSpace({'neum': ()}, -10, 10),
+
+ name='Sheet1_Problem_2a'
+ )
+
+
+def problem_1_2_b():
+ return StationaryProblem(
+ domain=RectDomain(bottom='neumann'),
+
+ neumann_data=ExpressionFunction('-cos(pi*x[...,0])**2*neum', 2, (), parameter_type={'neum': ()}),
+
+ diffusion=ExpressionFunction(
+ '1. - (sqrt( (np.mod(x[...,0],1./K)-0.5/K)**2 + (np.mod(x[...,1],1./K)-0.5/K)**2) <= 0.3/K) * (1 - diffu)',
+ 2, (),
+ values={'K': 10},
+ parameter_type={'diffu': ()}
+ ),
+
+ parameter_space=CubicParameterSpace(
+ {'neum': (), 'diffu': ()},
+ ranges={'diffu': [0.0001, 1.], 'neum': [-10, 10]}
+ ),
+
+ name='Sheet1_Problem_2b'
+ )
+
+
+def problem_1_2_c():
+ return StationaryProblem(
+ domain=RectDomain(bottom='neumann'),
+
+ neumann_data=ConstantFunction(-1., 2),
+
+ diffusion=LincombFunction(
+ [ConstantFunction(1., 2),
+ BitmapFunction('R.png', range=[1, 0]),
+ BitmapFunction('B.png', range=[1, 0])],
+ [1.,
+ ExpressionParameterFunctional('-(1 - R)', {'R': ()}),
+ ExpressionParameterFunctional('-(1 - B)', {'B': ()})]
+ ),
+
+ parameter_space=CubicParameterSpace({'R': (), 'B': ()}, 0.0001, 1.),
+
+ name='Sheet1_Problem_2c'
+ )
diff --git a/rb-intro/solutions/Sheet3_Solutions.ipynb b/rb-intro/solutions/Sheet3_Solutions.ipynb
new file mode 100644
index 0000000000000000000000000000000000000000..b6eb495be2aef96220d4c34527ba2007e0f8cbf8
--- /dev/null
+++ b/rb-intro/solutions/Sheet3_Solutions.ipynb
@@ -0,0 +1,858 @@
+{
+ "cells": [
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "# pyMOR RB-Intro\n",
+ "\n",
+ "# Exercise Sheet 3 - Solutions"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 1,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "from pymor.basic import *\n",
+ "set_log_levels({'pymor': 'WARN',\n",
+ " 'pymor.functions.bitmap': 'CRITICAL'})"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "We will continue to work with the models defined on exercise sheet 1. To simplify working with these models, you can find their definition in `problems.py` in the same directory as this notebook:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 2,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "from pymor.basic import *\r\n",
+ "\r\n",
+ "\r\n",
+ "def problem_1_1_a():\r\n",
+ " return StationaryProblem(\r\n",
+ " domain=RectDomain([[0., 0.], [1., 1.]]),\r\n",
+ "\r\n",
+ " diffusion=ConstantFunction(1, 2),\r\n",
+ "\r\n",
+ " rhs=ExpressionFunction('(sqrt( (x[...,0]-0.5)**2 + (x[...,1]-0.5)**2) <= 0.3) * 1.', 2, ()),\r\n",
+ "\r\n",
+ " name='Sheet1_Problem_1a'\r\n",
+ " )\r\n",
+ "\r\n",
+ "\r\n",
+ "def problem_1_1_b():\r\n",
+ " return StationaryProblem(\r\n",
+ " domain=RectDomain(bottom='neumann'),\r\n",
+ "\r\n",
+ " neumann_data=ConstantFunction(-1., 2),\r\n",
+ "\r\n",
+ " diffusion=ExpressionFunction('1. - (sqrt( (x[...,0]-0.5)**2 + (x[...,1]-0.5)**2) <= 0.3) * 0.999', 2, ()),\r\n",
+ "\r\n",
+ " name='Sheet1_Problem_1b'\r\n",
+ " )\r\n",
+ "\r\n",
+ "\r\n",
+ "def problem_1_1_c():\r\n",
+ " return StationaryProblem(\r\n",
+ " domain=RectDomain(bottom='neumann'),\r\n",
+ "\r\n",
+ " neumann_data=ConstantFunction(-1., 2),\r\n",
+ "\r\n",
+ " diffusion=ExpressionFunction(\r\n",
+ " '1. - (sqrt( (np.mod(x[...,0],1./K)-0.5/K)**2 + (np.mod(x[...,1],1./K)-0.5/K)**2) <= 0.3/K) * 0.999',\r\n",
+ " 2, (),\r\n",
+ " values={'K': 10}\r\n",
+ " ),\r\n",
+ "\r\n",
+ " name='Sheet1_Problem_1c'\r\n",
+ " )\r\n",
+ "\r\n",
+ "\r\n",
+ "def problem_1_1_d():\r\n",
+ " return StationaryProblem(\r\n",
+ " domain=RectDomain(bottom='neumann'),\r\n",
+ "\r\n",
+ " neumann_data=ConstantFunction(-1., 2),\r\n",
+ "\r\n",
+ " diffusion=BitmapFunction('RB.png', range=[0.001, 1]),\r\n",
+ "\r\n",
+ " name='Sheet1_Problem_1d'\r\n",
+ " )\r\n",
+ "\r\n",
+ "\r\n",
+ "def problem_1_2_a():\r\n",
+ " return StationaryProblem(\r\n",
+ " domain=RectDomain(bottom='neumann'),\r\n",
+ "\r\n",
+ " neumann_data=ExpressionFunction('-cos(pi*x[...,0])**2*neum', 2, (), parameter_type={'neum': ()}),\r\n",
+ "\r\n",
+ " diffusion=ExpressionFunction(\r\n",
+ " '1. - (sqrt( (np.mod(x[...,0],1./K)-0.5/K)**2 + (np.mod(x[...,1],1./K)-0.5/K)**2) <= 0.3/K) * 0.999',\r\n",
+ " 2, (),\r\n",
+ " values={'K': 10}\r\n",
+ " ),\r\n",
+ "\r\n",
+ " parameter_space=CubicParameterSpace({'neum': ()}, -10, 10),\r\n",
+ "\r\n",
+ " name='Sheet1_Problem_2a'\r\n",
+ " )\r\n",
+ "\r\n",
+ "\r\n",
+ "def problem_1_2_b():\r\n",
+ " return StationaryProblem(\r\n",
+ " domain=RectDomain(bottom='neumann'),\r\n",
+ "\r\n",
+ " neumann_data=ExpressionFunction('-cos(pi*x[...,0])**2*neum', 2, (), parameter_type={'neum': ()}),\r\n",
+ "\r\n",
+ " diffusion=ExpressionFunction(\r\n",
+ " '1. - (sqrt( (np.mod(x[...,0],1./K)-0.5/K)**2 + (np.mod(x[...,1],1./K)-0.5/K)**2) <= 0.3/K) * (1 - diffu)',\r\n",
+ " 2, (),\r\n",
+ " values={'K': 10},\r\n",
+ " parameter_type={'diffu': ()}\r\n",
+ " ),\r\n",
+ "\r\n",
+ " parameter_space=CubicParameterSpace(\r\n",
+ " {'neum': (), 'diffu': ()},\r\n",
+ " ranges={'diffu': [0.0001, 1.], 'neum': [-10, 10]}\r\n",
+ " ),\r\n",
+ "\r\n",
+ " name='Sheet1_Problem_2b'\r\n",
+ " )\r\n",
+ "\r\n",
+ "\r\n",
+ "def problem_1_2_c():\r\n",
+ " return StationaryProblem(\r\n",
+ " domain=RectDomain(bottom='neumann'),\r\n",
+ "\r\n",
+ " neumann_data=ConstantFunction(-1., 2),\r\n",
+ "\r\n",
+ " diffusion=LincombFunction(\r\n",
+ " [ConstantFunction(1., 2),\r\n",
+ " BitmapFunction('R.png', range=[1, 0]),\r\n",
+ " BitmapFunction('B.png', range=[1, 0])],\r\n",
+ " [1.,\r\n",
+ " ExpressionParameterFunctional('-(1 - R)', {'R': ()}),\r\n",
+ " ExpressionParameterFunctional('-(1 - B)', {'B': ()})]\r\n",
+ " ),\r\n",
+ "\r\n",
+ " parameter_space=CubicParameterSpace({'R': (), 'B': ()}, 0.0001, 1.),\r\n",
+ "\r\n",
+ " name='Sheet1_Problem_2c'\r\n",
+ " )\r\n"
+ ]
+ }
+ ],
+ "source": [
+ "%cat problems.py"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "The file can be imported as usual:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 3,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import problems"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "We have four non-parametric problems:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 4,
+ "metadata": {
+ "scrolled": false
+ },
+ "outputs": [
+ {
+ "data": {
+ "application/vnd.jupyter.widget-view+json": {
+ "model_id": "a1341fdeaaab4341b98c5d06e11771cf",
+ "version_major": 2,
+ "version_minor": 0
+ },
+ "text/plain": [
+ "ThreeJSPlot(children=(HBox(children=(Renderer(children=(Renderer(camera=PerspectiveCamera(position=(0.0, 0.0, …"
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ },
+ {
+ "data": {
+ "application/vnd.jupyter.widget-view+json": {
+ "model_id": "883de4ed34f84071ac9d0cb5f1cf72f2",
+ "version_major": 2,
+ "version_minor": 0
+ },
+ "text/plain": [
+ "ThreeJSPlot(children=(HBox(children=(Renderer(children=(Renderer(camera=PerspectiveCamera(position=(0.0, 0.0, …"
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ },
+ {
+ "data": {
+ "application/vnd.jupyter.widget-view+json": {
+ "model_id": "b176afd6bf0040609eb8d7f6238ec7e4",
+ "version_major": 2,
+ "version_minor": 0
+ },
+ "text/plain": [
+ "ThreeJSPlot(children=(HBox(children=(Renderer(children=(Renderer(camera=PerspectiveCamera(position=(0.0, 0.0, …"
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ },
+ {
+ "data": {
+ "application/vnd.jupyter.widget-view+json": {
+ "model_id": "ccdeb44fb8b84b4c86d73e90629ac04a",
+ "version_major": 2,
+ "version_minor": 0
+ },
+ "text/plain": [
+ "ThreeJSPlot(children=(HBox(children=(Renderer(children=(Renderer(camera=PerspectiveCamera(position=(0.0, 0.0, …"
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "for f in [problems.problem_1_1_a, problems.problem_1_1_b,\n",
+ " problems.problem_1_1_c, problems.problem_1_1_d]:\n",
+ " p = f()\n",
+ " m, _ = discretize_stationary_cg(p, diameter=1/100)\n",
+ " m.visualize(m.solve(), legend=m.name)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "The parametric problems have already been assigned a `ParameterSpace` which is inherited by the generated models:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 5,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "application/vnd.jupyter.widget-view+json": {
+ "model_id": "7d87b8468d4f41d2af49d7ceeac79d96",
+ "version_major": 2,
+ "version_minor": 0
+ },
+ "text/plain": [
+ "ThreeJSPlot(children=(HBox(children=(Renderer(children=(Renderer(camera=PerspectiveCamera(position=(0.0, 0.0, …"
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ },
+ {
+ "data": {
+ "application/vnd.jupyter.widget-view+json": {
+ "model_id": "1649d8e9e93246878262f2f2b1d3deee",
+ "version_major": 2,
+ "version_minor": 0
+ },
+ "text/plain": [
+ "ThreeJSPlot(children=(HBox(children=(Renderer(children=(Renderer(camera=PerspectiveCamera(position=(0.0, 0.0, …"
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ },
+ {
+ "data": {
+ "application/vnd.jupyter.widget-view+json": {
+ "model_id": "de720dc55ebe4d31813af2a551bb20b4",
+ "version_major": 2,
+ "version_minor": 0
+ },
+ "text/plain": [
+ "ThreeJSPlot(children=(HBox(children=(Renderer(children=(Renderer(camera=PerspectiveCamera(position=(0.0, 0.0, …"
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "for f in [problems.problem_1_2_a, problems.problem_1_2_b,\n",
+ " problems.problem_1_2_c]:\n",
+ " p = f()\n",
+ " m, _ = discretize_stationary_cg(p, diameter=1/100)\n",
+ " mu = m.parameter_space.sample_randomly(1)[0]\n",
+ " m.visualize(m.solve(mu), legend=str(mu))"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Problem 1"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Problem 1 a)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Write a method\n",
+ "\n",
+ "```python\n",
+ "def orthogonal_projection(U, basis, product=None,\n",
+ " basis_is_orthonormal=False,\n",
+ " return_projection_matrix=False):\n",
+ " ...\n",
+ "```\n",
+ "that encapsulates the orthogonal projection algorithm you implemented in Sheet 2, Problem 2."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Solution"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 6,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "def orthogonal_projection(U, basis, product=None,\n",
+ " basis_is_orthonormal=False,\n",
+ " return_projection_matrix=False):\n",
+ "\n",
+ " rhs = basis.inner(U, product)\n",
+ " if basis_is_orthonormal:\n",
+ " M = np.eye(len(basis))\n",
+ " v = rhs\n",
+ " else:\n",
+ " M = basis.gramian(product)\n",
+ " try:\n",
+ " v = np.linalg.solve(M, rhs)\n",
+ " except np.linalg.LinAlgError:\n",
+ " v = np.zeros(rhs.shape)\n",
+ "\n",
+ " U_proj = basis.lincomb(v.T)\n",
+ "\n",
+ " if return_projection_matrix:\n",
+ " return U_proj, M\n",
+ " else:\n",
+ " return U_proj"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Problem 1 b)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Repeat the error calculations from Sheet 2, Problem 2 with an orthonormal basis obtained from `U`\n",
+ "using `pymor.algorithms.gram_schmidt`. Compare your results."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Solution"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "In the following we use `problem_1_2_b`, for which the effect of different basis generation methods can be observed best:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 7,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "m, _ = discretize_stationary_cg(problems.problem_1_2_b(), diameter=1/100)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "First, we compute snapshots for the approximation space (`U`) and for the calculation of the approximation error (`V`):"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 8,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "U: ..............................\n",
+ "V: ...................................................................................................."
+ ]
+ }
+ ],
+ "source": [
+ "U = m.solution_space.empty()\n",
+ "print('U: ', end='')\n",
+ "for mu in m.parameter_space.sample_randomly(30):\n",
+ " U.append(m.solve(mu))\n",
+ " print('.', end='', flush=True)\n",
+ " \n",
+ "V = m.solution_space.empty()\n",
+ "print('\\nV: ', end='')\n",
+ "for mu in m.parameter_space.sample_randomly(100):\n",
+ " V.append(m.solve(mu))\n",
+ " print('.', end='', flush=True)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Now we compute the approximation error using our new `orthogonal_projection` method:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 9,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "errors = []\n",
+ "conds = []\n",
+ "for N in range(len(U)):\n",
+ " V_proj, M = orthogonal_projection(V, U[:N], product=m.h1_0_semi_product, return_projection_matrix=True)\n",
+ " errors.append(np.max(m.h1_0_semi_norm(V - V_proj)))\n",
+ " if N > 0:\n",
+ " conds.append(np.linalg.cond(M))"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "To compute an orthonormal basis for the space spanned by the vectors in `U`, we use `pymor.algorithms.gram_schmidt.gram_schmidt`-Method. Note that, unless the parameter `copy` is set to `False`, the input `VectorArray` will remain untouched:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 10,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "basis = gram_schmidt(U, product=m.h1_0_semi_product)\n",
+ "errors_orth = []\n",
+ "conds_orth = []\n",
+ "for N in range(len(U)):\n",
+ " V_proj, M = orthogonal_projection(V, basis[:N], product=m.h1_0_semi_product, basis_is_orthonormal=True,\n",
+ " return_projection_matrix=True)\n",
+ " errors_orth.append(np.max(m.h1_0_semi_norm(V - V_proj)))\n",
+ " if N > 0:\n",
+ " conds_orth.append(np.linalg.cond(M))"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Finally we plot the results. We observe the improved numerical stability due to orthonormalization:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 11,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "<Figure size 640x480 with 1 Axes>"
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "from matplotlib import pyplot as plt\n",
+ "plt.semilogy(errors_orth, marker='o', label='w. orth')\n",
+ "plt.semilogy(errors, label='wo. orth')\n",
+ "plt.legend()\n",
+ "plt.show()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 12,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stderr",
+ "output_type": "stream",
+ "text": [
+ "/home/stephan/.virtualenvs/pymor37/lib/python3.7/site-packages/matplotlib/axes/_base.py:3507: UserWarning: Attempting to set identical bottom==top results\n",
+ "in singular transformations; automatically expanding.\n",
+ "bottom=1.0, top=1.0\n",
+ " self.set_ylim(upper, lower, auto=None)\n"
+ ]
+ },
+ {
+ "data": {
+ "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAD8CAYAAABw1c+bAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xl81NW9//HXhyxEWQQhKhKEgIBERCkRBMUitgpVxKtoBaTAxSK9V221i9h7+7Pt7/bWW7WtWq9clAh1AbmoiMqibVVWlYAbi2BElgASFtmXEHLuH2fQELJMMpN8Z3k/Hw8emTmZme/nm3nwnjPne77na845REQkeTQIugAREalfCn4RkSSj4BcRSTIKfhGRJKPgFxFJMgp+EZEko+AXEUkyCn4RkSSj4BcRSTKpQRdQkZYtW7p27doFXYaISFxZtmzZDudcZnWPi8ngb9euHfn5+UGXISISV8xsQziP01CPiEiSUfCLiCQZBb+ISJKJyTH+ihw9epTCwkIOHz4cdClxISMjg6ysLNLS0oIuRURiTNwEf2FhIU2aNKFdu3aYWdDlxDTnHDt37qSwsJDs7OygyxGRGBM3Qz2HDx+mRYsWCv0wmBktWrTQtyMRqVDcBD+g0K8B/a1EpDJxFfwiIjHhwA74ZAbE6aVrFfwx5O2332bx4sVf3x81ahQzZswIsCIROYlzMPNH8OIYWPFi0NXUioI/RpSUlJwU/CISg9bMgc/egLRG8MavoPhA0BXVmII/TA8++CCPPvooAHfffTf9+/cH4B//+AfDhw+v8rnr16+nf//+dOvWjSuvvJKNGzcCvkc/btw4evXqxc0338yECRP405/+xEUXXcSCBQsAmD9/Pn369KF9+/bq/YsE7eghmHsvZJ4Hw6fDvi2w4I9BV1VjcTOds6zfvLqSVVv2RvU1c85uyv2Dzq/093379uXhhx/mrrvuIj8/nyNHjnD06FEWLFjA5ZdfXuVr33nnnYwcOZKRI0eSl5fHXXfdxcyZMwE/TXXx4sWkpKTw61//msaNG/Ozn/0MgEmTJrF161YWLlzIp59+ynXXXceQIUOit9MiseSr9XBqC2jYJOhKKrfwz7B7I4x8FdpdBt2+D4sfhe7D4fT2QVcXNvX4w9SjRw+WLVvG3r17adiwIb179yY/P58FCxbQt2/fKp+7ZMkShg0bBsCIESNYuHDh17+76aabSElJqfS5119/PQ0aNCAnJ4dt27ZFZ2dEYknxAZj3b/Dot+DJ/rBnc9AVVWzXF7DwT3D+DZAd6ux95zfQIA3m/XuwtdVQXPb4q+qZ15W0tDSys7OZPHkyffr0oVu3brz11lsUFBTQpUuXWr9uo0aNqvx9w4YNv77t4nQGgUilPnsTXrsH9myErkNg7TzIuxp+8Aq06BB0dSea90tokApX/cc3bU1bwbd/Dn/7NRT8Dc79TmDl1YR6/DXQt29fHnroIS6//HL69u3LhAkT6N69e7Vz5vv06cO0adMAeO655yr9htCkSRP27dsX9bpFYs7+Ipjxz/DcEEg7BUbPgSGTYNRrcPSgD/+tHwdd5TfWvgFrZsO3fwGntT7xd5f8ix/mmTMeSoqDqa+Goh78ZtbezCaZ2YwybTlmNt3MnjCzuB2k7tu3L1u3bqV3796ceeaZZGRknBDit912W4XXEXjsscd4+umn6datG8888wyPPPJIha8/aNAgXn755RMO7ooklNJSWDYF/pILq1+Ffr+EcQugbR//+7MvgtFzIaUhTL4WNiwJtl6Ao4dhzi+gZScf8uWlNoQBD8DOz+D9ifVfXy1YOMMHZpYHXAsUOee6lmkfADwCpABPOeceKPO7Gc65IaHbPwXed84tMLNZzrnrqtpebm6uKx+gq1evjmhIJRnpbyYxZfsaePUnsHExtL0MBv0ZWnas+LG7N8Ez1/vx/u8/Cx0DHEJ550F46z9gxEzocEXlj3vuJv9BdecyaHJm/dVXhpktc87lVve4cHv8k4EB5TaQAjwODARygKFmllPJ858BbjGzB4EWYW5TRBJByRF46/fwxKVQtAqu+4sf0qks9AGatfE9/5bnwtRbgjtRavdGWPAw5AyuOvQBrv49lByGv/+2fmqLQFjB75ybD+wq19wTKHDOrXPOFQPTgMGVPL/IOfevwHhgR0WPMbOxZpZvZvnbt28PewdEJIZtfM8H/jsPwPnXwx358K0REM5aUo0zYdTrkJULM8ZA/tN1X295c+/ztV71u+of2/Jc6P0v8OGzULis7muLQCRj/K2BTWXuFwKtzayFmU0AupvZfQBm1s7MJgJ/BR6s6MWccxOdc7nOudzMzGqvFSwS+w7v8WPayWrFSzD5Gjh2BG59EW58yod5TWScBre+5GfLvPYTP52yvhT8DT59Dfr+1H8DCcflP4fGZ8Kcn8f0ex/1g7vOuZ3OuXHOuQ7Oud+H2tY758Y654Y75xZW9xoica/g7/BQZ3juRjiShDO13pvoZ+1k5cLt8yOb5ph+KtzyPHS90U+bfPP+ul8creQIzP4FnN4B+twZ/vMaNoHv/hY2L4OPp9VdfRGKZB7/ZqDsx2BWqE0kua2ZA9N/AE1bw7p3/OyU4TNq3tuNR87BW7+D+Q9C5+/BkDw/XTNSqelww5P+G8CiP8Ph3dBzLBzeC0f2+m9Xh/eUuV2mveSI/+DpPgIahXmIcclfYNfnMPxFP2unJi64GZY+5T+gzrsWMprWfH/rWCTBvxToaGbZ+MC/BRgWlapE4tXKmX7VxrO6wYiXYNP7MH0kTPquvx9Hp/XX2LESeP0eWD4Fut8K1z4CKVE8R7RBClzzR8hoBgv/CMsmV/y4lHT/AdGwqQ/d0mPwt/vhrf/0xxlyx0CbnpUfZ9i9CeY/5EO7NrOJGjSAgf8FT14J8/9w4glfMSKsd8XMpgL9gJZmVgjc75ybZGZ3APPw0znznHMr66zSJPb222+Tnp5Onz5+rvOoUaO49tprY3fdng2L4aXbof+/w4XfD7qa+vPxdHj5dmjTC4ZN96HT6Wq/rsvzN8Gkq3zP/+yLgq40+o4e9h94x8fE+/8qvAO4NWUG37kf2n8bDu32f+OGp/mgz2jqwz4t4+TnFa2G/Dz4aBp8/AKc2RVy/xm63Xzy2kBv/Bu4Urj6P2tfZ+se/sPv3Seg+w8gs1PtX6sOhDurZ6hzrpVzLs05l+WcmxRqn+2c6xQazw/jsLfUVNwt17xnsx/m2FvoQ3DppKArqh/L/wovjfULd9364olf79tcDP/8BqRm+IOdn78VXJ114dBuePYGH/oD/guu/H91E/plte/ne+8d+kNWDz+jpvEZFYc+wBld4HsPwj2rYdAjvr7X74GHu8DrP4Vtq/zjPn8LVr3iP7yat42sxivvh7RTYe74mLtgi5ZsCFM4yzJPnTqVCy64gK5du3LvvfdW+5oJt1zz0cMwfYRfunbs2763+/o9fkXDRPb+kzDrTj+OPGw6pFew/lJmJxjzJjRr60/0+SRG3rNI7d0KT3/PD2ndOAkuGRd0RVVr2Bh6jILbF8CYv0GXa2H5M/BEb8gbAK/dDc2zoc9dkW+rcSb0uw8+/zusnRv560VRXC7Sxpzx8OUn0X3Nsy6AgQ9U+uvqlmXesmUL9957L8uWLaN58+ZcddVVzJw5k+uvv77S10yo5Zqdg9k/9bMZvv8stLrQ/3z5dj++emSfH/pJtGsBL3oU3vyVHw8eklf1gcCmrWD0bJg2zA+L7N8Gvf+1/mqNth0F8Ow/wYGdfm36Dv2Drih8Zv6bWJuL/ZDOB8/6oaCvvoBh/1v5N4ea6vlDfyzipduh01V+Vc/sy6F5u+i8fi2pxx+m6pZlXrp0Kf369SMzM5PU1FSGDx/O/Pnzq3zNhFquOT/P/+e5/OfQZZBvS0nzMzG+9QNY8BDMuTem5zbXiHPwzh986He9EW6aHN7sj1Oa+XnpXa7zqz2+8avY+Js457+phTsksXkZ5F3ll1Qe9Vp8hX55p54Ol94Fdy6Hu1f5gI6WlDS4eYo/SLzuHf/N8JEL4c/d4JU74OP/hX1fRm97YYrPHn8VPfO6Ut2yzJ999llUtxdXyzVvfNeHeser/FfbshqkwKBH/UG3JX+B4v3+fjRne9Q35/xp+Qv/CBcOg8F/8fsZrrQM/0Ex5xf+Ih77i/xrpKTVWckVcg62fODHtFfPgl3rAIP0xn64qmHoZ3qT0M9QW2oGfPCcv2jKiJf9+HoiaNDg5JU3o+GMLv7boHN+vaIv5sMX7/hF6j54xj+mZedvvg20u8x/GNWhOP7fV/+OL8ucl5fHBRdcwD333EOPHj0wM3r27Mldd93Fjh07aN68OVOnTuXOO6s+8eP4cs0jRoyodrnmvXuje8WxqNm7BV4Y4c9svOHJigPQzE9pa9gU3v5PP+xz41M1nx9dV3Z9Afu2+tkdDZv4Ohs2qTiInfM99Xf/G3qM9tMLG9Tii3ODFPjeQ9DkLPjHf8CB7f6bUXpjf8JSeiN/O+3Ub0K3Jh8ulSkthcKlPuhXzfLr4FuKD5wLh8GxYv/hXLwfjuz3Pfri/bD/yxPvn9HFD+U1bRV5TcnCDM44z//rNdZPM/3yk9AHwXz48HlY+qQ/xvDjD+u0FAV/DfTt25ff/e539O7dm0aNGp2wLHOrVq144IEHuOKKK3DOcc011zB4sF+66LbbbmPcuHHk5p64aN5jjz3G6NGjefDBB8nMzOTppytei2TQoEEMGTKEV155hccee6xud7ImSo74GTzFB/yFM05pVvljzaDfvb7HOO+XMHWoD470U+uv3rJ2fAarZvrebmXHi1IzynwYhD4QSo5A4fvQ60cw4PeRHbMw++YU/1d/4g8CViU1I9TrbgJNs6DZOSf/a9r65G9Tpcdg45JQz/5V/yGXkg7tr4B+46HzwDrvYUoFGqT4qb1nX+SHmo4dhc3L4dBXdb7psJZlrm9aljk66vxv9uqP/YGrm6b4qXXhWjbFP/ec3jBsmp+DXdecg+2f+vBb9YpfJRL8nPucwXBGTqiXuy/0b2+Z2/tCZ4Lug+J9/vF9fxbdA9UHdvghn+M96qMHv7ldfACKD35z+8he2FPoV47cuwUo83/YUnz4H/8gaNDAX9XqwHb/wXHud3z9na6un7+71Ktwl2VWj19qJ/9pH/qX3V2z0AfoMdL3/F8aC1Ou8wc7wz2Vviacg20rvgn7HWsB8xf9GPgHfxC66dnR325tNGrp/9VUSbE/Z2L3xpP/ffGOH54590rIuQ7O/a7/u0vSU/BLzW16H2b/3M/k6P+r2r1G1xv9GPYLI+Cp/tD2Un8CTqMz/M/GZ/ghkMZn+FP0K+pdHzvqe7L7i8r8LIL92/1UyS3L/QFLa+APmPW6Hc4bFNhFMupEarpfBiKRl4KQqIur4HfOVXt9W/HqbAhv35c+rJue7U/YieSAY6er/Vmuf/8NrHvbB3fp0ZMfl5Ie+kDIhLRGcDA0LHKo/CUiQtIa+Q+MFh3g0h/7Ofa16U2LJKi4Cf6MjAx27txJixYtFP7VcM6xc+dOMjKidBLKcSXFfsGxI3vh1jejc0Awuy/c9jd/2zl/YOvrnnuR77nvL/qmrfiAv3LT8W8IZb8lNMr0Pys6c1ZEvhY3wZ+VlUVhYSG6Old4MjIyyMrKis6LlRTDzgJY8jhsetf39M/qWv3zasrMf5icejpwXvRfX0SAOAr+4ydQSR0qKfZrkBet9jNgtn8KRZ/6ttIS/5hLfwwXxOiqoCISlrgJfomy49Mb18yGrR/5Mwp3FnwT8Bicng2ZXeC8a/wJO2fkwJnnB1q2iEROwZ9MSkv9GiufvgqrX/M9efBnCp7RxZ/Ik9nFn1nYslN0rpwkIjFHwZ/ojh2F9Qv9GZtrZvuzNhuk+umNl/zI9+ZjZS67iNQLBX8iKj7oT/9f/ZpfB/zwbkg9xZ/I02WQn0Z5SvOgqxSRgCj4E81X6/3ZsLs3+BOfOg/089g79A9uXRwRiSkK/kSyowCmDPLrvAx9wffw63upXxGJeVEPfjNrD/wbcJpzbkiorQHw/4GmQL5zbkq0t5v0ilb7nr4rhVGv1808exFJCGEtJG5meWZWZGYryrUPMLM1ZlZgZuMBnHPrnHNjyr3EYCALOAoURqNwKWPrx/4i3tbAX9pPoS8iVQj3ChKTgQFlG8wsBXgcGAjkAEPNLKeS53cGFjvn7gF+VLtSpUKFy2DKtf7g7ejZkNk56IpEJMaFFfzOuflA+RWxegIFoR5+MTAN37OvSCFw/OoCx2pTqFRgwxL462A/Q2f0bL8omYhINSK52HprYFOZ+4VAazNrYWYTgO5mdvwCrC8BV5vZY0CFVyA3s7Fmlm9m+VqPJwzr3oFnb/BLDI+eA83bBl2RiMSJqB/cdc7tBMaVazsIlB/3L/+8icBE8FfginZdCeWzN+GFW/0a7D94xa9IKSISpkh6/JuBNmXuZ4XapC6tfs1fr7ZlJxj5mkJfRGoskuBfCnQ0s2wzSwduAWZFpyyp0IoX/cXNW10II1+tm8sVikjCC3c651RgCdDZzArNbIxzrgS4A5gHrAamO+dW1l2pSW7Fi/Dibf7i4D+YCac0C7oiEYlTYY3xO+eGVtI+G5gd1YrkZCXFMPc+OLs73DpDV5gSkYhoyYZ4sPJlfwnCwf+t0BeRiEUyxi/1wTl47wl/MLdD/6CrEZEEoOCPdZvegy0fQK/boYHeLhGJnJIk1r3735BxGlxY4WEWEZEaU/DHst0b/ZWzeozS2L6IRI2CP5a9/yRgcPEPg65ERBKIgj9WFR+A5VP8pRKbtan+8SIiYVLwx6qPpsLhPXDJvwRdiYgkGAV/LCothXcnwNnfgjY9g65GRBKMgj8Wff532PkZXPIjMAu6GhFJMAr+WPTuE9D4LMi5PuhKRCQBKfhjTdGnvsff8zZITQ+6GhFJQAr+WPPeBEhpCD1GB12JiCQoBX8sObgLPpoG3W6GRi2DrkZEEpSCP5YsnwIlh/xBXRGROqLgjxXHjvozdbMvhzPPD7oaEUlgCv5YsfpV2LtZJ2yJSJ1T8MeKd5+A5tnQ8eqgKxGRBKfgjwWFy6Dwfeg1Tmvui0idU8rEgveegIZNofvwoCsRkSQQ9eA3s/ZmNsnMZpRp62JmE8xshplpykpZe7f4a+p2HwENmwRdjYgkgbCC38zyzKzIzFaUax9gZmvMrMDMxgM459Y558aUfZxzbrVzbhxwM3BptIpPCEsnQekx6Kk190WkfoTb458MDCjbYGYpwOPAQCAHGGpmOZW9gJldB7wOzK5VpYno6CHIz4PzroHTs4OuRkSSRFjB75ybD+wq19wTKAj18IuBacDgKl5jlnNuIKCB7OM+ng6HdumELRGpV5GM8bcGNpW5Xwi0NrMWZjYB6G5m9wGYWT8ze9TM/odKevxmNtbM8s0sf/v27RGUFSeOHYWFf4RWF0JbjX6JSP1JjfYLOud2AuPKtb0NvF3N8yYCEwFyc3NdtOuKOR9Ng6/Ww9AXtOa+iNSrSHr8m4GyF4PNCrVJdUqKYf4f/BW2OumELRGpX5EE/1Kgo5llm1k6cAswKzplJbiPnofdG+GKX6q3LyL1LtzpnFOBJUBnMys0szHOuRLgDmAesBqY7pxbWXelJoiSYpj/ELTOhXO/E3Q1IpKEwhrjd84NraR9NpqeWTMfPAN7NsGgR9TbF5FAaMmG+lRyBBY8DG0ugQ79g65GRJKUgr8+Lf+rX3r5ivvU2xeRwCj468vRQ7633/ZSyP520NWISBKL+jx+qcSyybBvK9zwpHr7IhIo9fjrQ/FBWPgnaNcXsvsGXY2IJDn1+OtDfh7s3wY3TQ66EhER9fjrXPEBWPRnaN8P2vYJuhoREQV/nVv6FBzYDv1+GXQlIiKAgr9uHdkPix6BDlfCOb2CrkZEBFDw1633J8LBnX5NHhGRGKHgryuH98LiR6HjVZCVG3Q1IiJfU/DXlff/Bw59Bf3GB12JiMgJFPx14fAeWPwX6DQQWvcIuhoRkRMo+OvCuxPg8G719kUkJin4o+3QbljyOJx3LZx9UdDViIicRMEfbR8+D0f2wLfvDboSEZEKKfij7Yt34PQO0Kpb0JWIiFRIwR9NpcdgwxJod2nQlYiIVErBH03bVvhhnraXBV2JiEilFPzRtH6R/6kev4jEsKgHv5m1N7NJZjajTFsjM5tiZk+a2fBobzNmbFgEzdrCaVlBVyIiUqmwgt/M8sysyMxWlGsfYGZrzKzAzMYDOOfWOefGlHuJG4AZzrkfAtdFpfJYU1rqg7+dhnlEJLaF2+OfDAwo22BmKcDjwEAgBxhqZjmVPD8L2BS6fazmZcaB7av9Eg1tNcwjIrEtrOB3zs0HdpVr7gkUhHr4xcA0YHAlL1GID/9Kt2lmY80s38zyt2/fHk5ZsUXj+yISJyIZ42/NN7148OHe2sxamNkEoLuZ3Rf63UvAjWb2BPBqRS/mnJvonMt1zuVmZmZGUFZANiyEpll+jF9EJIZF/Zq7zrmdwLhybQeA0dHeVsxwDjYshg79wSzoakREqhRJj38z0KbM/axQW/LZsdZfXlHj+yISByIJ/qVARzPLNrN04BZgVnTKijPrF/qfmtEjInEg3OmcU4ElQGczKzSzMc65EuAOYB6wGpjunFtZd6XGsA2LoPFZcHr7oCsREalWWGP8zrmhlbTPBmZHtaJ445yf0dPuUo3vi0hc0JINkdq1DvZ/qfF9EYkbCv5IaXxfROKMgj9SGxZBo0xo2SnoSkREwqLgj8Tx8f22fTS+LyJxQ8Efid0bYW+h1t8Xkbii4I/EBq3PIyLxR8EfifWL4JTmkNkl6EpERMKm4I/EhoV+GmcD/RlFJH4osWprz2b4ar3m74tI3FHw15bG90UkTin4a2v9Qmh4GpzZNehKRERqRMFfWxsWQdve0CAl6EpERGpEwV8b+76EnQUa3xeRuKTgrw2N74tIHFPw18b6RZDeBM66MOhKRERqTMFfGxsWwTm9ICXqlywWEalzCv6aOrADtn+q8X0RiVsK/pr6enxfC7OJSHxS8NfU+kWQdiqc3T3oSkREaqVeBqnNrC8wPLS9HOdcn/rYbp3YsAja9ISUtKArERGplVr3+M0sz8yKzGxFufYBZrbGzArMbDyAc26Bc24c8BowJbKSA3RwF2xbqfX3RSSuRTLUMxkYULbBzFKAx4GBQA4w1MxyyjxkGPB8BNsM1sYlgNP8fRGJa7UOfufcfGBXueaeQIFzbp1zrhiYBgwGMLNzgD3OuX213Wbg1i+C1Axo3SPoSkREai3aB3dbA5vK3C8MtQGMAZ6u7IlmNtbM8s0sf/v27VEuK0o2LISsiyG1YdCViIjUWr3N6nHO3e+cW1zF7yc653Kdc7mZmZn1VVb4Du+BLz/R/H0RiXvRDv7NQJsy97NCbfFv47vgSjW+LyJxL9rBvxToaGbZZpYO3ALMivI2grF+IaSk+6EeEZE4Fsl0zqnAEqCzmRWa2RjnXAlwBzAPWA1Md86tjE6pAduwyB/UTTsl6EpERCJS6xO4nHNDK2mfDcyudUWx6Mg+2PIhXHZ30JWIiERMSzaEY+N74I5pfF9EEoKCPxxr50LqKXBO76ArERGJmIK/Os754O9whcb3RSQhKPirs20l7NkEnQZU/1gRkTig4K/O2jn+Z6erg61DRCRKFPzVWTMXzv4WNDkr6EpERKJCwV+V/UWweRl0/l7QlYiIRI2Cvypr5wEOOmt8X0QSh4K/KmvnQtMsOLNr0JWIiESNgr8yRw/D5//wvX2zoKsREYkaBX9lvpgPRw9Cp4FBVyIiElUK/sqsnQNpjaCdrq8rIolFwV8R5/yB3Q5XQFpG0NWIiESVgr8iX34MezdDZw3ziEjiUfBXZM1cwKCjztYVkcSj4K/I2jmQlQuNY/DavyIiEVLwl7d3K2z5QMM8IpKwFPzlfTbP/9Q0ThFJUAr+8tbMgWbnwBldgq5ERKROKPjLKj4I6972vX2drSsiCaregt/MGplZvpldW1/brLEv3oGSw1qUTUQSWq2D38zyzKzIzFaUax9gZmvMrMDMxpf51b3A9Npur16smQPpTaCtztYVkcQVSY9/MnBC19jMUoDHgYFADjDUzHLM7LvAKqAogu3VrdJSf7buuf0hNT3oakRE6kxqbZ/onJtvZu3KNfcECpxz6wDMbBowGGgMNMJ/GBwys9nOudKyTzSzscBYgHPOOae2ZdXe1g9h/5eazSMiCa/WwV+J1sCmMvcLgV7OuTsAzGwUsKN86AM45yYCEwFyc3NdlOuq3tq5YA2g41X1vmkRkfoU7eCvknNucn1ur0bWzIGsntCoRdCViIjUqWjP6tkMtClzPyvUFtv2bPYLs2k2j4gkgWgH/1Kgo5llm1k6cAswK8rbiL61c/xPXVRdRJJAJNM5pwJLgM5mVmhmY5xzJcAdwDxgNTDdObcyOqXWoTVzoXk2tOwUdCUiInUuklk9Qytpnw3MrnVF9a34gL/M4sVjdLauiCQFLdnw+Vtw7Ah00vi+iCQHBf/aOdDwNGjbJ+hKRETqRXIHf2kprH0Dzr0SUtKCrkZEpF4kd/BvWQ4HinTRFRFJKskd/GvmgKXAud8JuhIRkXqT3MG/di6ccwmcenrQlYiI1JvkDf5dX8C2FZrNIyJJJ3mDf/kUvyjb+f8UdCUiIvUqOYP/6GFY/le/BHOzNtU/XkQkgSRn8K96BQ7uhJ63BV2JiEi9S87gX/oknN4BsvsFXYmISL1LvuDf8iEULoWLb4MGybf7IiLJl3xLn4LUU+CiYUFXIiISiOQK/kNfwSczoNtNcEqzoKsREQlEcgX/B89BySG4+IdBVyIiEpjkCf7SUsifBG16QatuQVcjIhKY5An+df+AXev8QV0RkSSWPMH//lNwakvIGRx0JSIigUqO4P9qg1+QrcdISG0YdDUiIoGql+A3s/ZmNsnMZtTH9k6y7Gl/Pd0eowPZvIhILKl18JtZnpkVmdmKcu0DzGyNmRWY2XgA59w659yYSIutlZIjWpdHRKSMSHr8k4ET1jQ2sxTgcWAgkAMMNbOcCLYRuZUztS6PiEgZtQ5+59x8YFe55p5AQaiHXwxMA4J8mDsvAAAFPklEQVQ9mqp1eUREThDtMf7WwKYy9wuB1mbWwswmAN3N7L6KnmhmY80s38zyt2/fHp1qtC6PiMhJUutjI865ncC4ah4zEZgIkJub66KyYa3LIyJykmh3gzcDZY+gZoXa6p/W5RERqVC0g38p0NHMss0sHbgFmBXlbYTnw+e1Lo+ISAUimc45FVgCdDazQjMb45wrAe4A5gGrgenOuZXRKbUGSkv9MI/W5REROUmtx/idc0MraZ8NzK51RdFwfF2efhUeRxYRSWqJOdVl6SStyyMiUonEC/7dG7Uuj4hIFRIv+PPz/E+tyyMiUqF6mcdfb0qOcODdp1mV3ouHXijEnz8mIhI/cs5uyv2Dzq/TbSRW8B/ey+dNezKXK4KuREQkZiVW8DfOpNtd/4smcIqIVC7xxvhFRKRKCn4RkSSj4BcRSTIKfhGRJKPgFxFJMgp+EZEko+AXEUkyCn4RkSRjzkXnKofRZGbbgQ3lmlsCOwIop64l6n5B4u6b9iv+JOq+ld+vts65zOqeFJPBXxEzy3fO5QZdR7Ql6n5B4u6b9iv+JOq+1Xa/NNQjIpJkFPwiIkkmnoJ/YtAF1JFE3S9I3H3TfsWfRN23Wu1X3Izxi4hIdMRTj19ERKIg5oPfzAaY2RozKzCz8UHXE01mtt7MPjGzD80sP+h6asvM8sysyMxWlGk73czeNLPPQj+bB1ljbVWyb782s82h9+1DM/tekDXWhpm1MbO3zGyVma00sx+H2uP6fativ+L6PTOzDDN738w+Cu3Xb0Lt2Wb2XigfXzCz9LBeL5aHeswsBVgLfBd/HcWlwFDn3KpAC4sSM1sP5Drn4np+sZldDuwH/uqc6xpq+wOwyzn3QOgDu7lz7t4g66yNSvbt18B+59xDQdYWCTNrBbRyzi03sybAMuB6YBRx/L5VsV83E8fvmZkZ0Mg5t9/M0oCFwI+Be4CXnHPTzGwC8JFz7onqXi/We/w9gQLn3DrnXDEwDRgccE1SjnNuPrCrXPNgYEro9hT8f764U8m+xT3n3Fbn3PLQ7X3AaqA1cf6+VbFfcc15+0N300L/HNAfmBFqD/v9ivXgbw1sKnO/kAR4E8twwBtmtszMxgZdTJSd6ZzbGrr9JXBmkMXUgTvM7OPQUFBcDYeUZ2btgO7AeyTQ+1ZuvyDO3zMzSzGzD4Ei4E3gc2C3c64k9JCw8zHWgz/RXeac+xYwEPjX0LBCwnF+PDF2xxRr7gmgA3ARsBV4ONhyas/MGgMvAj9xzu0t+7t4ft8q2K+4f8+cc8eccxcBWfjRkPNq+1qxHvybgTZl7meF2hKCc25z6GcR8DL+zUwU20LjrcfHXYsCridqnHPbQv8JS4EnidP3LTRW/CLwnHPupVBz3L9vFe1XorxnAM653cBbQG+gmZmlhn4Vdj7GevAvBTqGjlynA7cAswKuKSrMrFHo4BNm1gi4ClhR9bPiyixgZOj2SOCVAGuJquPBGPJPxOH7FjpYOAlY7Zz7Y5lfxfX7Vtl+xft7ZmaZZtYsdPsU/ISX1fgPgCGhh4X9fsX0rB6A0LSrPwMpQJ5z7ncBlxQVZtYe38sHSAWej9d9M7OpQD/8SoHbgPuBmcB04Bz8Sqs3O+fi7iBpJfvWDz9k4ID1wO1lxsXjgpldBiwAPgFKQ82/xI+Hx+37VsV+DSWO3zMz64Y/eJuC77BPd879NpQj04DTgQ+AW51zR6p9vVgPfhERia5YH+oREZEoU/CLiCQZBb+ISJJR8IuIJBkFv4hIklHwi4gkGQW/iEiSUfCLiCSZ/wOXj+AgN1lnBQAAAABJRU5ErkJggg==\n",
+ "text/plain": [
+ "<Figure size 432x288 with 1 Axes>"
+ ]
+ },
+ "metadata": {
+ "needs_background": "light"
+ },
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "plt.semilogy(range(1, len(U)), conds_orth, label='w. orth')\n",
+ "plt.semilogy(range(1, len(U)), conds, label='wo. orth')\n",
+ "plt.legend()\n",
+ "plt.show()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Problem 2"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Problem 2 a)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Extend the `greedy` function from Sheet 2, Problem 3 to a function\n",
+ "\n",
+ "```python\n",
+ " def greedy(U, basis_size, product=None,\n",
+ " rtol=0,\n",
+ " orthonormalize=False):\n",
+ " ...\n",
+ "```\n",
+ "\n",
+ "The algorithm should stop, when the relative approximation error goes below the prescribed value `rtol`. When `orthonormalize` is set to `True`, the algorithm shall return an orthonormal basis w.r.t. to `product` by orthonormalizing each new snapshot vector w.r.t. `product`.\n",
+ "\n",
+ "**Hint:** When using `gram_schmidt` for orthonormalization, the `offset` parameter can be used to avoid re-orthonormalizing the old basis. "
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Solution"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 13,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "def greedy(U, basis_size, rtol=0, product=None, orthonormalize=False):\n",
+ " # initialization\n",
+ " basis = U.space.empty()\n",
+ " max_errors = []\n",
+ "\n",
+ " # compute norms of vectors in U for relative approximation error computation\n",
+ " norms = U.norm(product)\n",
+ "\n",
+ " # main loop\n",
+ " for n in range(basis_size):\n",
+ "\n",
+ " # compute projections\n",
+ " # since the greedy basis is hierarchical, calling orthogonal_projection is\n",
+ " # suboptimal since previous results could be reused\n",
+ " U_proj = orthogonal_projection(U, basis, product=product, basis_is_orthonormal=orthonormalize)\n",
+ "\n",
+ " # compute projection errors\n",
+ " U_err = U - U_proj\n",
+ " errors = U_err.norm(product)\n",
+ "\n",
+ " if np.all(errors / norms < rtol):\n",
+ " break # relative approximation error threshold reached\n",
+ "\n",
+ " # select vector for basis extension\n",
+ " max_errors.append(np.max(errors))\n",
+ " basis.append(U[np.argmax(errors)])\n",
+ "\n",
+ " if orthonormalize:\n",
+ " gram_schmidt(basis, offset=n, product=product, copy=False)\n",
+ " if len(basis) == n:\n",
+ " break # Gram-Schmidt has removed newly added vector since it is numerically linear dependent\n",
+ "\n",
+ " return basis, max_errors"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Problem 2 b)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Repeat the experiments from Sheet 2, Problem 3, using `strong_greedy` with `orthonormalize=True`."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Solution"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 14,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "greedy_basis, _ = greedy(U, 30, product=m.h1_0_semi_product, orthonormalize=False)\n",
+ "greedy_orth_basis, _ = greedy(U, 30, product=m.h1_0_semi_product, orthonormalize=True)\n",
+ "\n",
+ "errors_greedy = []\n",
+ "errors_greedy_orth = []\n",
+ "for N in range(len(U)):\n",
+ " V_proj = orthogonal_projection(V, greedy_basis[:N], product=m.h1_0_semi_product)\n",
+ " errors_greedy.append(np.max(m.h1_0_semi_norm(V - V_proj)))\n",
+ " \n",
+ " V_proj = orthogonal_projection(V, greedy_orth_basis[:N], product=m.h1_0_semi_product,\n",
+ " basis_is_orthonormal=True)\n",
+ " errors_greedy_orth.append(np.max(m.h1_0_semi_norm(V - V_proj)))"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 15,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAD8CAYAAAB3u9PLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xd8VMX6+PHP7GbTIQkJNQECCBGQ0PFKE7GAooCNKmDlWtCrXr1gQ/h6KQp6xfJTvNJUpKiIqChXAaWISpBeRWoSIJCwaaRsduf3x9mEhBRSNtls8rxfr7zYPXt2zpwsOc/OnJlnlNYaIYQQwuTuCgghhKgeJCAIIYQAJCAIIYRwkoAghBACkIAghBDCSQKCEEIIQAKCEEIIJwkIQgghAAkIQgghnLzcXYGyCAsL05GRke6uhhBCeJRt27ad01rXv9x+HhUQIiMjiYmJcXc1hBDCoyiljpdmP+kyEkIIAUhAEEII4SQBQQghBOBh9xCKYrPZiI2NJTMz091VERXg6+tLREQEFovF3VURotby+IAQGxtLnTp1iIyMRCnl7uqIctBak5iYSGxsLC1atHB3dYSotTy+yygzM5PQ0FAJBh5MKUVoaKi08oRwM48PCIAEgxpAPkMh3K9GBAQhhPA0jqwsrF+soDotYywBwQP89NNP/PLLL3nP7733Xj7//HM31kgIUVFp63/i1AsvkLlzp7urksfjbyqX1crtccxac5B4awZNgv14dkAUQzuHu7taxcrJyeGnn34iMDCQnj17urs6QggXsZ9PAiA7Ng6/Tp3cXBtDrWohrNwex3MrdhNnzUADcdYMnluxm5Xb48pd5qxZs3jrrbcAeOqpp+jfvz8A69atY/To0SW+99ixY/Tv35/o6Giuv/56Tpw4ARgtgIcffpirr76aYcOG8f777/Of//yHTp06sXHjRgA2bNhAz549admypbQWhPBA9uRkAGxx5b/+uFqNaiFM/Xov++JTin19+wkr2XZHgW0ZNjv/+nwXS34/UeR72jWpy8u3tS+2zD59+vD666/zxBNPEBMTQ1ZWFjabjY0bN9K3b98S6/v4448zbtw4xo0bx/z583niiSdYuXIlYAyn/eWXXzCbzUyZMoXAwECeeeYZAObNm8epU6fYtGkTBw4cYPDgwdx1110lHksIUb3YrdUvINSqFsKlweBy20uja9eubNu2jZSUFHx8fLjmmmuIiYlh48aN9OnTp8T3btmyhVGjRgEwZswYNm3alPfa3XffjdlsLva9Q4cOxWQy0a5dO86cOVPu+gsh3MNutQLVKyDUqBZCSd/kAXrNXEecNaPQ9vBgP5b9/ZpyHdNisdCiRQsWLlxIz549iY6OZv369Rw+fJi2bduWq0yAgICAEl/38fHJe1ydRikIIUqnOnYZ1aoWwrMDovCzFPzW7Wcx8+yAqAqV26dPH2bPnk3fvn3p06cP77//Pp07d77s2PqePXuydOlSABYvXlxsi6JOnTqkpqZWqI5CiOolr4UQH492lL+XwpXcGhCUUs2UUiuVUvOVUpMq+3hDO4cz444OhAf7oTBaBjPu6FDhUUZ9+vTh1KlTXHPNNTRs2BBfX98CF/cHH3ywyHUc3n77bRYsWEB0dDQff/wxc+bMKbL82267jS+//LLATWUhhGfLbSHo7Gxyzp1zc20MytXdDUqp+cCtQILW+qp82wcCcwAz8KHWeqZSahAQorX+RCm1TGs9vKSyu3Xrpi+9sO7fv79CXTOi+pDPUtQmh3r1Rnl7k3PqFM2XfIp/586Vdiyl1DatdbfL7VcZLYSFwMBLKmMG3gVuBtoBI5VS7YBfgQeUUuuA7yuhLkIIUe1orbEnJ+Pbvh0Atrh4N9fI4PKAoLXeACRdsrkHcFhrfURrnQ0sBYYA9wEva637A4NcXRchhKiOHOnpkJODr7NFXF1uLFfVPYRw4GS+57HObd8DTyil3geOFfVGpdR4pVSMUirm7NmzlV5RIYSobLk3lC2Nm2CuV6/aBAS3DjvVWu8BSpxRpbX+APgAjHsIVVEvIYSoTLmT0szBQVjCw7HFxrq5RoaqaiHEAU3zPY9wbhNCiFrHnmy0EMxBzoBQTVoIVRUQtgKtlVItlFLewAhgVRUdWwghqpXcLiNzcDCW8CbVZi6CywOCUmoJsAWIUkrFKqUe0FrnABOANcB+YLnWeq+rj13bSFpsITxT7hyE3BaCttnIOev+uQguv4egtR5ZzPbVwGpXH6+2krTYQnguR25AqFsX73BjYqwtLg5LwwburFbtSl1RGUqT/nrJkiV06NCBq666iokTJ162TEmLLUTNZrdaMQUEoLy9sUREANVj6GmNSm7Hd5Pg9G7XltmoA9w8s9iXL5f+Oj4+nokTJ7Jt2zZCQkK46aabWLlyJUOHDi22TEmLLUTNZrcmYw4KAsDSpAlQPQKCtBAq6HLpr7du3Uq/fv2oX78+Xl5ejB49mg0bNpRYpqTFFqJms1utmIODATD5+WEODa0WAaFmtRBK+CZfWS6X/vrPP/906fEkLbYQns+enIw5OCjveXUZeiotBBcoKf11jx49+Pnnnzl37hx2u50lS5Zw7bXXlliepMUWomazJydjCsofEJpIQKgpSkp/3bhxY2bOnMl1111Hx44d6dq1K0OGDAEkLbYQtVX+LiMA7/DwajEXweXpryuTpL+u2eSzFLWBdjg4cFUHQsc/RIMnnwTg/JIlnJ76f1zx88+VMvTUnemvhRBCFMORng4OB+agiy0ES765CO4kAUEIIapQ/rQVuSQgCCFELZSX6TT/TeW8uQjuzXoqAUEIIapQUS2E6jIXQQKCEEJUobzEdvnmIUD1mIsgAUEIIapQ/rUQ8rOENyFbAoIorX79+hU5b8FdduzYwerVFxPYTpkyhdmzZ7uxRkJUf3ldRpcEBO+ICGzxp9w6F6Fmpa64jH7L+pGYmVhoe6hvKD8N/6lSj52Tk4OXV835defk5LBjxw5iYmK45ZZb3F0dITyGPTkZU2Ag6pLrgSU8HGw2cs6exdKwoVvqVqtaCEUFg5K2l9Yrr7xCVFQUvXv3ZuTIkXnfkvv168eTTz5Jt27dmDNnDmfPnuXOO++ke/fudO/enc2bNwOQnp7O/fffT48ePejcuTNfffUVABkZGYwYMYK2bdty++23k5GRAcD8+fN50jmhBeC///0vTz31VIE6ffbZZzz99NMAzJkzh5YtWwJw5MgRevXqVeL5JCUlMXToUKKjo/nb3/7Grl27AKMFMGbMGHr16sWYMWOYPHkyy5Yto1OnTixbtgyAffv20a9fP1q2bJmXFlwIcZEjOblQ6wCqx9DTmvOVFXj191c5kHSgXO+97/v7itx+Zb0rmdij+DUMtm7dyhdffMHOnTux2Wx06dKFrl275r2enZ2d180zatQonnrqKXr37s2JEycYMGAA+/fvZ9q0afTv35/58+djtVrp0aMHN9xwA3PnzsXf35/9+/eza9cuunTpAsCwYcOYNm0as2bNwmKxsGDBAubOnVugXn369OG1114DYOPGjYSGhhIXF5eXlrskL7/8Mp07d2blypWsW7eOsWPHsmPHDsC44G/atAk/Pz8WLlxITEwM77zzDmAEjAMHDrB+/XpSU1OJiorikUcewWKxlHg8IWqTnEvSVuQqEBCcf+tVrUYFBHfYvHkzQ4YMwdfXF19fX2677bYCrw8fPjzv8Y8//si+ffvynqekpJCWlsb//vc/Vq1aldeyyMzM5MSJE2zYsIEnnngCgOjoaKKjowEIDAykf//+fPPNN7Rt2xabzUaHDh0KHLdRo0akpaWRmprKyZMnGTVqFBs2bGDjxo3ccccdJZ7Tpk2b+OKLLwDo378/iYmJpKSkADB48GD8/PyKfe+gQYPw8fHBx8eHBg0acObMGSKcC4AIIcBhLaaFUA3WRahRAaGkb/IAHRZ1KPa1BQMXuLo6QMF01Q6Hg19//RVfX98C+2it+eKLL4iKiip1uQ8++CDTp0/nyiuv5L77im7d9OzZkwULFhAVFUWfPn2YP38+W7Zs4fXXXy/fyVC29Ntms5mcnJxyH0uImshutea1BvIz+fpiDgtza0CoVfcQKkOvXr34+uuvyczMJC0tjW+++abYfW+66SbefvvtvOe53TADBgzg7bffzlu/YPv27QD07duXTz/9FIA9e/bk9eUDXH311Zw8eZJPP/2UkSOLXMa6QFruzp07s379enx8fAgq4tvJpe9bvHgxAD/99BNhYWHUrVu30H6SfluIsrt0LYT83J0G260BQSnVTym1USn1vlKqX2UfL9Q3tEzbS6N79+4MHjyY6Ohobr75Zjp06FDsBfett94iJiaG6Oho2rVrx/vvvw/ASy+9hM1mIzo6mvbt2/PSSy8B8Mgjj5CWlkbbtm2ZPHlygXsTYNxL6NWrFyEhIUUer0+fPpw8eZK+fftiNptp2rQpvXv3znt98uTJrFq1qtD7pkyZwrZt24iOjmbSpEksWrSoyPKvu+469u3bV+CmshCieNrhwJ6SUmAthPy8w8PdOxdBa+3SH2A+kADsuWT7QOAgcBiY5Nx2LfAdsBC44nJld+3aVV9q3759hbZVtdTUVK211unp6bpr165627ZtVXLcQYMG6R9//LFKjlUVqsNnKURlyrFa9b6oK3XiwoVFvn5m9my976oO2mG3u/S4QIwuxfW7MloIC50X/zxKKTPwLnAz0A4YqZRqB2zUWt8MTASmVkJdqsT48ePp1KkTXbp04c4778wbDVRZrFYrbdq0wc/Pj+uvv75SjyWEcJ3ctBXFtRDy5iIkJFRltfK4/Kay1nqDUiryks09gMNa6yMASqmlwBCtde6Qm/OAD0VQSo0HxgM0a9bM1dV1idx+/qoSHBzMoUOHqvSYQoiKy8tjVFJAwBhpZGnUqMrqlauq7iGEAyfzPY8FwpVSdyil5gIfA+8U9Uat9Qda625a627169evgqoKIUTlKCrTaX7unpzm1mGnWusVwAp31kEIIarKxbUQigkIbp6LUFUthDigab7nEc5tQghRa+S1EEKKDggmX1/M9cPcNtKoqgLCVqC1UqqFUsobGAEUHu8ohBA1WN49hDp1it3Hu4n71kVweUBQSi0BtgBRSqlYpdQDWuscYAKwBtgPLNda73X1sWu66pb+ujiSFluIotmTkzHVqVMo02l+xkI58VVYq4sqY5RRkdNmtdargdVFvVYb1LT018WRtNhCFM9eTGK7/Czh4aT88APabkeZzVVUM4OkrnABT01/vXbtWjp37kyHDh24//77ycrKKvE8JS22EBVjT7YWO+Q0V/51EapajfrKenr6dLL2ly/9dXF82l5Jo+efL/Z1T01/nZmZyb333svatWtp06YNY8eO5b333isQaC4labGFqBh7MWsh5OfOuQjSQqig/Omv69Spc9n01xMmTKBTp04MHjy4QPrrmTNn0qlTJ/r161cg/fU999wDFJ/++sCBA2VOf92nTx8OHjxIixYtaNOmDQDjxo1jw4YNJZ7rpk2bGDNmDFD+tNhhYWF5abGFqG1K22UE7hl6WqNaCCV9k3eX6pr++tixY+U6n+JIWmwhLq+4tRDyszRpDLgnIEgLoYI8Nf11VFQUx44d4/DhwwB8/PHHXHvttSWeq6TFFqL8tN2OPSXlsi0Ed85FkIBQQZ6a/trX15cFCxZw991306FDB0wmEw8//DAgabGFqAyO1FTQuti1EPJz21yE0qRErS4/kv66IEl/LYTnyDp2TO+LulJbV6687L6xTz2t/7zhRpcdGzemv651JP21EOJyLpfYLj9LeDi2U6fQdntlV6uAGnVT2V0k/bUQ4nIul/o6P0t4OOTkkJOQgKVx48quWp4a0ULQzpuxwnPJZyhqusstjpOfJSICqPqRRh4fEHx9fUlMTJQLigfTWpOYmFhoOK4QNYn9fFm6jNyTBtvju4wiIiKIjY3lrBumeQvX8fX1JcL5rUiImsienAxKYS5iqPalctdFqOqhpx4fECwWCy1atHB3NYQQokR2qxVT3bqlSlhn8vHBq3596TISQoiaqDR5jPJzRxpsCQhCCFEFyhcQpIUghBA1TmkS2+XnjrkIEhCEEKIKlKeFkDsXoapIQBBCiCpgt15+cZz83JEGWwKCEEJUMm234yhFptP83DEXwa3DTpVSQ4FBQF1gntb6f+6sjxBCVAa7cyGpMrUQ3DAXweUtBKXUfKVUglJqzyXbByqlDiqlDiulJgForVdqrR8CHgaGF1WeEEJ4uouJ7UofEPLmIsR6cEAAFgID829QSpmBd4GbgXbASKVUu3y7vOh8XQghahxHbmK7MnQZQdUPPXV5QNBabwCSLtncAzistT6itc4GlgJDlOFV4Dut9R+urosQQlQHZcl0mp8lIsKzA0IxwoGT+Z7HOrc9DtwA3KWUerioNyqlxiulYpRSMZKvSAjhicqyFkJ+lvBwbKdPo6toDXK33lTWWr8FvHWZfT4APgDo1q2bpDQVQniccrcQwptcXBfBeZO5MlVVCyEOaJrveYRzmxBC1Hh2qxWUwlSnTpneV9VzEaoqIGwFWiulWiilvIERQOFV3IUQogayW5MxlzLTaX7ezoBQVUNPK2PY6RJgCxCllIpVSj2gtc4BJgBrgP3Acq31XlcfWwghqiN7cjKmMgw5zeXVpGonp7n8HoLWemQx21cDq119PCGEqO6MtBVlu6EMYPL2xqtBgypLgy2pK4QQopLZk5PLNCktv6qciyABQQghKpmR6bTsLQSQgCCEEDVKWddCyK8q5yJIQBBCiEqkc3JwpKaWeQ5CrvxzESqbBAQhhKhE5cl0ml9VzkVw60zlqtBvWT8SMxMLbQ/1DeWn4T9VfYWEELWK3Vq+xHa58uYixMbh3727y+pVlBrfQsgLBloXvV0IISqRPbnsqa/z82rSBJSqkhZCjQ8IAIEXNP/+yE6How53V0UIUcvkJbYrZ5fRxbkIEhBcIiQdAjPgpaUOHv7WTkCG5MgTQlQNeznXQsjPv0cPvMJCXVWlYtX4ewgAJ+srnn3AzF2bHQz+VdP5LzvzBphw2O2YyphbRAghysJRzkyn+YXPes1V1SlRrWghANgsiiX9zDx3r5nzgfDMCgdfD45m+89furtqQogaLMdqBZOpzJlO3aHGB4RQ34LNrGONFC+MM7O0nxeRxxzYn3iepU/eyIX0VDfVUAhRkzmSnZlOTdX/clvju4yKHVp6P+zcsJKE/3uRjt/H8s2BHswdaCYhRBXYTYanCiEqwkhsV/7uoqpU/UNWJerYdygD1uxg7x2daBUPs+fZGfS7A+W4eNNZhqcKISrCbk2u0A3lqlSrAwKA2ezFXdOX8PRDZvY2V4xb6+ClJQ78smQkkhCi4sq7FoI71PqAkCupruLVu0y8O8jElbGaF5ba8c+UoCCEqBjpMvJUSvFztIk3bjfR4gxM/tRO4AUJCkKI8jPWQpAuI48V08bEa3eaiEiElz+1c/LwTndXSQjhgbTNhiMtTVoInubS4ak7W5mYebeJhlb4a+xIYv/c4aaaCSE8VV6mUw9pIbh12KlSqiXwAhCktb7LnXUpbmjplz730fzTXzkyZhT6o09o2qZL1VZMCOGx8tJWlHO1tKpW7haCUmq+UipBKbXnku0DlVIHlVKHlVKTSipDa31Ea/1AeetQFW6ftIATI7oQkKE5fs89HN/3m7urJITwEBVNbFfVKtJltBAYmH+DUsoMvAvcDLQDRiql2imlOiilvrnkp0EFjl2lhkz6mLg72+GTrYkbdx9Hdm9yd5WEEB6gomshVLVydxlprTcopSIv2dwDOKy1PgKglFoKDNFazwBuLe+x3E2ZTNz24mesdtxK2KqjnLn3IR4baSa2vsxqFkIU72Km05rfQihKOHAy3/NY57YiKaVClVLvA52VUs8Vs894pVSMUirm7Nmzrq1tGSiTiYEvrsJ6axMcCqYsttP8jCy6I4QoXm3qMqowrXWi1vphrXUrZyuiqH0+0Fp301p3q1+/flVXsQCzlxfXvfAtL482k20xhqRGnpZ5CkKIotmTrWA2e0SmU3B9QIgDmuZ7HuHcVmN4+/hypp7i5dFmLvjAi0vthJ+ToCCEKMyem+lUqcvvXA24OiBsBVorpVoopbyBEcAqFx+jWjgbrPj3CDMOkxEU6lslKAghCrJbrR5zQxkqNux0CbAFiFJKxSqlHtBa5wATgDXAfmC51nqva6pa/Zyup3hlhBkfG7y0xE5wmgQFIcRFjuRkj7l/ABUICFrrkVrrxlpri9Y6Qms9z7l9tda6jfO+wDTXVbX6CPQKyXt8soFi+jAzwelGUEg6fdyNNRNCVCc5HpTYDiR1RblsGb2BV6K/p+6pOaTtn4mV1zhxQwANz8PWkbeRak1wdxWFENWAw4PWQoBasGJaZRnaOZyhnS+OqM3MuJbv6U+rNSlsHHYj/Vb8jH+g5/xHEEK4npHpVFoItY6vXwA3z1jPof5BND+RzY/D+5GdccHd1RJCuInOzsaRno5JuoxqJx9ffwa/8TN7+9Wl9V9ZfDuyNzm2bHdXSwjhBp6W6RSky8jlLN4+3PH2Rr54pDfRm1J5/55OzL3ZBPnGIUuKCyFqvouZTqWFUKt5Wby5e+4WvuipuH6nZuxaB+iLQ1IlxYUQNV9e2gppIQiT2cyyvib8sxzculWT6q/5sqdnzFYUQlScp62FANJCqFxKsfBGE5vbKu7e6KBJokxcE6K2sJ/PbSFIl5Fw0kqx4EYT2RYY+6PD3dURQlSRi6mvpYUg8kkJUHzW20SXI5rOh42g8PucUWSkp7q5ZkKIymJPTjYynQYEuLsqpSYBoRKF+obmPf6+qyK+Hoxb6yA4x4tuSas58/o1HNsf48YaCiEqi92ZtsJTMp2C3FSuVJcOLf016F2CXniHFxN6svf6yTRe9yQBS2/m9+jn6X77P1Amic9C1BTGLGXP6S4CaSFUqb/d+RjH2tUjbOk6GrXqAA9v5LBve3rsnsIf/7mTzzfvpdfMdbSY9C29Zq5j5fYatZSEELWK3cMS24EEhCrX+uXpeNvgt6n/IKxRM9r/ay2/Rj5Gx5Sf6LbmdkKS96KBOGsGz63YLUFBCA8lLQRxWVd0vJbjN7WnxYYj7N/yLSazmb/dO52/e03FW9n4wnsK15m2A5BhszNrzUE311gIUR72ZGkhiFLo/eIc0v0VJ1+ZisNhjDpal96KW7JmcFA35R3LW7RXRwGIt2a4s6pCiHKyWz1rcRyQgOAWwWHhpNx7G02PpLLpk1cBaBLsh5U6PJD9DOepw3zvWTQmkSbBfm6urRCirBzZ2egLFzCHSJeRKIV+j7zCqca+mN9dTHpqEs8OiMLPYuYsIdyX/S/8yGKB92s81aeBu6sqhCgjhwcmtgMJCG7jZfEmaOLT1Eu28/NrTzO0czgz7uhAeLAfh3UEE83P0krFE7XhcWzZWe6urhCiDPIS23lYQJB5CG7UdeAYvu0yn/AvfyP+3l0M7RydbxW2Qfz+pYUeO1/i9/fup/vjH8s8BSE8hCemrYBq0EJQSgUopWKUUre6uy7uED1lNgr4Y+rThV7rcfsTbIm4nx7nv+HXj1+s+soJIcolt4XgSaulQQUCglJqvlIqQSm155LtA5VSB5VSh5VSk0pR1ERgeXnr4ematulK3OButPo9jh0/LCn0+t/uf52YujdwzdF3ifnmv26ooRCirHJbCF4e1kKoSJfRQuAd4KPcDUopM/AucCMQC2xVSq0CzMCMS95/P9AR2Af4VqAeHq/vxDfZ8UNvkqf9H2PjpqEvWV1tzaOr2ffGTURvncT++s1oe/UAN9ZWCHE5dqsREEwetBYCVKCFoLXeACRdsrkHcFhrfURrnQ0sBYZorXdrrW+95CcB6Af8DRgFPKSUKlQfpdR4Z5dSzNmzZ8tb3WotMCiUT64z0fI09NtVcM2ExMxEfHz9Cf/7F5w2N6Txd/dz8s+dbqqpEKI07FYreHlhCvB3d1XKxNX3EMKBk/mexzq3FUlr/YLW+kngU+C/WutCCwZorT/QWnfTWnerX7++i6tbfWxupzgQDiN/dhCQUXghnaDQhpjv+RwHJlh8NwOnr5CcR0JUU7lpKzwp0ylUg5vKAFrrhVrrb9xdD7dSioU3mgnMgKmL7YSmFA4K4S3bsSb6TerrJP6dOQNvsiXnkRDVkD3Z82Ypg+sDQhzQNN/zCOc2UQpHGitmDDMRlgz//shO8zOFg8I7h0J40vYYXdSfzLG8i4UcyXkkRDXjiZlOwfUBYSvQWinVQinlDYwAVrn4GDXa7hYmJo8xAzD1EzsdjhbsRYu3ZvC9owdTc8Yy0LyV9yz/wYdsyXkkRDXiiZlOoWLDTpcAW4AopVSsUuoBrXUOMAFYA+wHlmut97qmqjVb/tXVTjRQvDDWTEIwPLfcwbr3X8p7LTe30SL7AF6w3c8N5u18aJlNyyDP6qsUoibz1BZCuYedaq1HFrN9NbC63DWqpS5dXQ0gZfRpfrnvdpq/+TmrY08y8P/m8+yAKJ5bsZsMm53F9hvI0N7Mssxlsd9rkNkbfOtWfeWFEAXUuhaCqHx16zWi//K1/NWzOS0+/43VD97CLe1C8nIeKWBz4I087XiCMOsuHIsGw4VLRwILIaqSIysLnZFRu1oIomp4+/hzy4er+X7yfbT6/HdWDu3O7NsVGY0VgY3hArAO6JHTkt9j96AWDkKN/QoCJUuqqGHitoElABpc6e6alCh3Upo52PMCgrQQPIDJZOKWfy/i1D/upO1xB1M/sROSWnAEks0ri3FZz5Bz7ggsuAWSZXCXqCG0ht/moufegJ53MyTHurtGJbInOzOdSpeRqEz9H/k3M4aZaGCFaYvsNE0oGBQCrryB0VkTyUk5BQtuhvPH3FNRIVzFboNv/4lj1USO/tSc49964fj0XsjJdnfNiuWpayGABASPkzssVQFTFtsJP3cxKLw+rCOJ9bpwn/1FHBnJRkvh3GH3VVaIisiwwuK7IWYep0/+jayz2WQkWji18k/0D5PdXbti5XjoWgggAcEjnWigmDLajN0MLyy9OKu5jq+FuWO6sd3ekif9/o3OyTIDSBBJAAAgAElEQVRaCmdk5K/wMIl/wbwb4dgmrMHjSf79OGGPPEL9f/yDlOP+JC74GPZ95e5aFsnhoWshgAQEj3UmRDFtuBn/bCMo1LlgBIUrGgTyxrCOrDpdj1lN/oM2mY2g8OkI+GoC/DgVtrwLu5bDX+vg9G5IPW00zYWoDo5tgg+vh/SzZF33HqcXrcW/Rw/CHnuU0L+Pp+7NAzm7qy6pbz9hBI5qxlNXSwMZZeRxQn1DScxMBOB4Q8Wrd5l5Yamd5z5zcGF0Ev5163FT+0Y80f8K3lp3mJwO/6H3X28QdmAfDUwp1FOpmHRO0YX7BoN3ICgFOCe6KfI9zrfdZAafuuAXAn7Bxr++wQWf524LaQ7eAa75BWRYYfdnkHQEzN75fizg5WP8a/YGs/OxbxC07GfUV1R/2z+Br5+Eei1w3LGIuIefw+TrS5NZs1Bm4zNsPH062Uf/In7jISLfG4XPxJ/A4ufeeudjT05GWSwo/8tkOrXnQOopCIpw/m25nwQED1PUBLafW82i5Svz2TRuMP2X/IiXry9P3tCGH/ad4YPdqXzAk3n7+llMzL61OYNaWSD9LKSfK/iv7YIxqgMAXfxjhx2yUiDjvHHzOuM8ZFqhcMJa8PKDqJshehi0uh68vMt20lrDiV/hj0Ww90vIyTSGHzpywF6K9aavuBHummcEh2JkHjqEUgqf1q3LVjdRpJQffsDk509g716le4PDDj9OgV/egpbXwd0LOTPzTbIOHaLpfz/A0vDiMGqTnx8R73/A0duHcPLLJCJb/gOvUR9UzomUg92ajCk4qHCmU3sOnN4JRzcaraATWyA7DYKawZWD4MpboFlPMLvvsiwBoQa4dtSzfJeUQOQ737B2/BBunP8tJi8vrBmFu4EybA6mrz/NoKv7Q5iLL34OB2SnGsEhw+r89zwc3wx7VsDeFUarod1QIzg0/RuUtE50eiLsXAJ/fATnDoJ3Heg0CrqMgyadjH20dgaGbOePDXKyLj4++jOseR4+vAFGLoXQVoUOk/z1N5x6/nmUxUKzjz/Cr3171/5e3MFug+O/QPNeVXqBsaelcXrq/5Hy9dcAhD0+gbBHHil5PfCsNFgxHg5+C90fhIEzSVnzA9Zlywh96EEC+/Qp9BZLo0Y0fW8ux+8ZTdz7P9Ks9SJU93GVdVplkpe2wp4Dp3fBMWcAOL7F+PsACIuCjiMgtDUc+Qli5sNv7xl/H20GGgGiVX/XtaxLSWldOKNmddWtWzcdExPj7mpUW19NvY82S34l9qYO3DBnGS2fW01Rn64Cjs4cVLWVs9uMexa7P4MD3xotkaCmcNWdRnBo6LwIOxzGH9Afi2D/18aFPaI7dL0X2t9evj+Qoxth+VjQdrh7ofGHBmitOffee5x76238unUlJ/4UjuxsIpd8infTpiWXWZ3ZMuGze+HQd8a53jXfuNBUsoxdu4j75zPY4uIIe/RRbCdPkvzVV9S58UaazJyBKaCIzy7piPHZnNkLA1+Fq8eTfeIER2+/A5/WrWn+8Ucoi6XYY1q/XMGp514gpE0mjT74GhpdVYlneAmtjS88ea1s4+f4vxejM1OJvC7BaEUDhLWByN4Q2cf499KJo1lpxt/HgW/h0PdGa9vL1/j8rhxkBImAsHJXVSm1TWvd7bL7SUCoObTWfPnkUNquOcTpEdfyou8w4orIghoe7MfmSf3dUEOnrDQ4+B3sXg6H1xoX6gbtoUVf44/h/FGje6fjSOgy9mKwqIiko7B0FJw9CAOmozvfx6nJL5P81VcEDRlMo1dewRYby/GRozAFBxH56ad4hYZevtzqJisVlow0gmqn0cbggeBmRuuofptKOaR2OEj8cB5n33oLrwb1CZ89G/8uXdBac/6jjzjz6mv4tGpFxP9792Kg1Rp2LIbvJoIyG0Gr9Q04srM5PnIU2bGxtFzxBZbwYtfXynNm2lSSPl5Ko2sthMzZVDn5vI5tgh2fGn3+ac6L/4VzRuvUyZGjOLcnkMSDgQRFedHkgRsuBoE6DUt/rNzW3YFvjZ+UWFAm4/Mc8k65qi8BoZay2W18/cAA2v56iu0jbuEV+41k2OwF9unTOowF93bHy1wNBpmlnzPuC+z+DE7+ZnRxdBkH7Qa7/kZhViqsGI9913fE7mzLhSNWo0vj0Ufz+nsvbN/OifvuN76dLlxQ9Lfa6upCkjFuP347DH0POg43uimW3WO0tO6aD61vdOkhbWcSiJ80kQtbfqXOwIE0njql0OiatM2biXv6nyggfM6bBES3ga//AftXGRfLoe9BsBEoTk+bzvmPPybi3Xeoc/31paqDtts5ee8I0mN20+zeKwn415euu0l74ldYPw2ObjBaWfVaQUB949t6YAPn4/qk7o7j9NzPyTlzjqA7bqfBs8/iFeKCVpnWcGqnERjqNILuD5SrGAkItdiFzFR+HH0TrfZamTPExJZ2BS/8jpxAOus3eXdUF4L8i2+OV7mcLGOkUCXKPnaMk2OHYTuXQuObwwia+iUEFlyaNXXdOmInPE5A7140fffdErssipNz/jxn33gDe1oa9caMwb9LlwrXPefcOaxfrMBctw7Bw4bljboxKn0GPr4dEv+EuxZA21svvmY9CUtHwuk9cOP/Qc/H8y6Y9uRkEt58E/u5c/hfcw2BvXvj3axZqeqTun49p55/AUdmJo1eeJ6gO+8sdsnI7OPHOfnYY2QfOUrDq3MIiUxE3fASXDMhbwRY6o8/EjvhcULGjqHR88+X6XdjT03l2G03YU9KJHLaA3jfNrFM7y/k5Fb4abrRjRPQAHo/Bd3uK/QlxRYfz+lp00lbuxaf1lfQaMoU/Lt2rdixK4EEhFrOmnyGdXf2o3U8zLzbxK6WBYNC5qFXiQjx58Nx3WhVP9BNtaxaF/74g9jHJoDDQcQ/78L/0CzjG97IJdCoQ4F9zy9fzunJLxM0dCiNZ0wv09q4qevWc2ryZOzJyZj8/XEkJ+PXqRP1HrifOv37F7yQl0Lm/v0kLfqIlG+/RduMgQK+0dE0mfZvY1SU9QR8NMQICiMWQ6vrCheSnQ4rHzEmc0WPgNvmkLr5V06/NJmc8+exNGiALT4eAEvTpgT06klg7974X3015jp1ChTlyMoiYdZszn/yCT5XXkn4G6/j07JlySdhy8T+7YvEv/05aXF+BN18HY1efROTtzHizBYXx5Hb78C7aVOaL/k0b3tZZB89ytHbb8Pim0XzRfMwR/UtcxnEbYP1M+DwD+AfBr2fhG4PgHfBIaTaZiPpo484+867ANSf8Bj1xo4t15eHqiABQXD13KuYsthOo/OwoqeJHzsp0vyNC9v8a3/m7x9vw2Z38O6oLvRtU/8ypXm25G+/5dRzz+PVuBHN5s7FOzLS+ONfOhoyk+H2uUY3VT5n33mXc++8Q+jf/06Dp54suuB87KmpnJkxk+QVK/CJiqLJqzPxbtYM64ovSVq4EFtsLN7Nm1PvvvsIGjoEk69vsWVpu5209etJWvQRF7ZuRfn7E3z77YTcM5qs/fs5/cq/caSlEXbvcEL1J6icdBj9OTTtUXwFtYYNs7Cvmc6Zgy1J3puBT5s2NJk5A5+2bbEdP07a5s2kb/6FC7/+iuPCBTCb8evYMS9AKD8/4p/9F1kHDxIydgwN/vlPTD6XadWd3gNfPAhn96O7jefcn405N/e/+HXqRPhbc/AKCeH4PWPIOnyYFl+uKHULpShp69dw8tEnCWymafTALVgi20Jwc2MuTHCz4rsh43fATzOMe1h+9aDXE9D9IfAp/GXpwh9/cPrlKWT9+SeB/fvT6IXnS3Wvw50kIAg6LOpAcJrmsa8ddDymyfKCDVcpvutm4rtn9nAy6QIPfRTDoTOpDO3UhN+OJhFvzaRJsB/PDohiaOeK/Se3JyeTsXsPpgB//Dp1KtO3bFfRWpM4dy5n35yDX7euRLz9dsG+3dTTxs3muG0QPRzCuxk3sRu2R/sGcfrlKViXL6fhiy9S757RxR4nfcsW4l94gZzTZwgd/xD1H30Ule9brs7JIfWHH0icN5/MPXsw16tHyD2jCRk5skB97GlpJH/xBUkff4ItNhZLkyaE3HMPwXfdibnuxZulOUlJnHnhGVLWb8GnnqbJq7Pw7XPbZX8faRs3cWrSM+QkWQntCPVnLkJFXl3495adTcbOnUaA2LSZzL178+ahmENCaDxjOnX69Sv5YA4H/Pr/YO1Uo/99yLt59zBSvl9D/HPPYa5TB/8ePUj55hvC//MGdW+++bLncDlJ78zkzDuLAPAJshHQKIuAxln418/CFNTIGRycQSKoKfz5PzjwjTGJsufjcPXfwadOoXJzuwGtn32OV+PGNHrxhVLf53A3CQiCDosudoM0Pau5OcZB3z0a7xwI6NmTeuPGontcw4j//sbe+JQC7/WzmJlxR4dSBwVts5F56BCZu3aRsWMnGbt2kX30aN7rPq2vIHj4CIKGDC7UBVEZ7MnJZB0+jHX5cpK/WkXd226j8bR/F90VYcuE7yfBvpXGMMJcQU3R9dsR+1USaXviCZ/6LHXvHFNgXL/jwgUSXn+D84sX4x0ZSZNXJuPXtpUxrNaW6bzpeHG0ktaaC1u3kjRvPmk//4zy9SX4jjuoO+gWUtasIfmLFTjS0/Hr2pV6Y8dS5/r+KK8i5hGc+BUWDyMlPpDT24Kwp6QRNv4hQh9+uMhztKelkfDqa1g/+wzvK1rR5JkH8dv+EqSdMUauRA8r8feZc/48F7ZsIevIUYKH3Y2lwWXW27CehK8eM+aBRA2CwW8VGjaZeeAAsY8+hi0+nuDhw2k8dUrJZZZB5sFDpG/4mbSf15OxczfaloOymPFvEURAhCIwzIq3KQ6FA3yCoOcEIxD4BmFPTcUWG0t2bCy22DhsJ0+SHRdL5o6dxj2he8dR/9FHPWrAgQQEUSAg5KpzQXPDDs2Y/WHkJCTgHRnJ/Ibd+aJ+JzIvuaFb3PBUnZ2NLeEsmXt2k7FzFxm7dpG5dy86MxMAc2goftHR+HXsiF/HaGzx8ZxfspTMPXtQfn7UHXQLIcNH4Neh4mPG7VYrWYcPk3X4L+Pfvw6Tdfgw9rPn8vYJe+wxwiY8dvkWitZGi+HMHufPXjizF8fpQ5xYG0zmeQtN+6cQ0K45aM2F2HTi19qwpSjqRV2gfodkTF5F/D017mjM0G7VH5penTdTO+vPP0lcsJDkr78Gmw28vKh7y83UGzO25N/N4bXGyKE6jWHsV9ipY3RVffUVPq2voPG0afhFR+ftnv7LL8S/+KLRenngfsImTDC6edLPGXMAjm82WkfRw42hv+Zy9oPnZMOfa2D7YuNbt5cPDJxpDB0u5nefk5RE6v/+R9DQoSV2oVWE48IFLmzdStqmzaRv2pT3RcWrUUMCunfEK7Qh2acSsMXGYouNxe5MTpfLVKcOlogIfFq0IPTvf8c3qnKG71YmjwgISqk+wGiMGdPttNY9S9pfAkLZ9FvWLy/vUX71fOvx0x0/krLmfyR99BGZu3aR5uXL95FXsyvsCupmpxOUlUZQdjoPRYdgTzqPPSmJnPPGv460tLyylLc3vu3a4dcxGt/oaPw6dsIS3qTIi2/G7j1Yly8j+Ztv0RkZ+LZvT8jIEdS95RZMJeR90Tk52OLiyD5+nOxjx8g+doysI0eNC/+5ixd+k78/3ldcgU+rVvhccQU+V7TCJyoKS6NGFftF5mSR89c2jj8ykZzzyTS7pwWp+1NI3HwaS7AvjUd1JaBtU6N/2uJnTCiy+IPFFxKPwF9r4eTvxnwLSwC06HMxQIS2wpZwlvRffiGgZ88CKRoAo9sl/awxFj05zpix/fNrxkSnMV8WmOCU9vPPnHp5CjkJCdS7715CH3iAs2+/jXXJUrwjI2k8Yzr+nTtfcm7ZRpfOtoVGGgW/esYIpfa3Q2Tf0s1yPr3bCAK7l8OFRAhsZMzC7Xa/0S1Tzdji4vK6wtK3bEFnZmIJD8cSEYElIhzvpk2xhEdgiYjAu2mERyapu1SlBwSl1HzgViBBa31Vvu0DgTmAGfhQaz2zFGUNBRpqreeWtJ8EhIrZm7iXe1bfw7UR1/Kffv/Ju2iPeXoevXb+SO/4XZjz5SLKMZnxDQvFXK8eXvVCMIfUu/g4NBTftm3xjYoq0FdeGvbUVJJXrcK6dClZfx7GFBhI0JAhBA0ZjONCRt5FPy8AxMYa36CdTHXq4N2ihfOi77zwX3EFXo0bV+p9CltcHMdGjiInIQGA4LvvpsHEiZgDS9F1kJlijGX/a50RIHIXLwpuZgSGZj2NWa0pccaFPyXOWBks9ZQxhyC/Zj1h5KdFzj62p6aS8NosrJ99BmYzOBzUGzeO+k/+o+Rv4LYMo+Wxb6UxaTA7DfxD4crc4NCnYHBITzTmjuxYbKRnMHtD1C3Q+R4jF5Eb8/GUhXYY/99LTK1RA1RFQOgLpAEf5QYEpZQZOATcCMQCW4GRGMFhxiVF3K+1TnC+bznwgNY6taRjSkCouEV7FzE7ZjYv/e0lhkUZ/cYrt8fx3Ird+KUk0SDjPMnegST7BNCqeQNWPNYbSyVNYNNak7F9O+eXLCX1++/zhlQCKB8fvJs3xzsy8uK/LSLxjozEHBLilhvUAJkHD5Iwazb1xtxD4LXXlr+gpCNGcDi8zggUuTluTBao2xjqRkBQONR1/uQ+DoowLtSXOf/0X37h/JIl1Bs3Dv9ul70OFJQbHPZ+aYy6yQ0ObW8z8k8dXG0EDYcNGncyZtB2uAv865XzlyEqW5V0GSmlIoFv8gWEa4ApWusBzufPAWitLw0G+ctoBryktX7ocseTgFBxDu3g0bWPEnM6hiWDltA6xEhwt3J7HLPWHCTemkGTYD+ubhHCiu3x3BrdmDeHd6r0Wc0558+TvmkTXqGheEdG4tWoUY3/1pbHboPEw0Z3TUD9khP+VTVbBhz+EfY6Ww62dGN8fvRwI9FgVeYOEuXmroBwFzBQa/2g8/kY4Gqt9YQSypgKrNFa/1LM6+OB8QDNmjXrevz48XLXVxjOZZzjrlV3EeIbwpJBS/D1Kror4YMNfzF99QFu7xzO7Ls7YjZVj5ztwk1sGZCwDxp2KHsKc+FWpQ0Ibv8qorV+ubhg4Hz9A611N611t/r1a/bkqaoS5hfG9N7TOWw9zOyY2cXuN75vK565qQ1fbo/juRW7cDg8Z0SaqAQWPwjvKsGgBnN1QIgD8ucNjnBuE9VMz/Ce3Nf+PpYdXMba42uL3W9C/9Y8cX1rlsfE8tJXe/CkYcpCiLJxdUDYCrRWSrVQSnkDI4BVLj6GcJHHOz9O+9D2TP5lMqfTTxe731M3tObha1ux+LcTTP16nwQFIWqoco8NU0otAfoBYUqpWOBlrfU8pdQEYA3GyKL5Wuu9LqmpcDmL2cJrfV/j7q/vZuKGicwfMB9zEWsPK6WYODAKm93BvE1HOZ6YzqEzqS5NcyGEcD+ZqSy45tNrSLOlFdoe6htaYA1nrTVj5v3GpsMFJ7uVNc2FEKJqecxNZeF+RQUDoNAsZ6UUR86lF9ovw2Zn1pqDlVI3IUTVkYAgyuSUNbPI7fFFLNUphPAsEhBEmTQJLjqffHHbhRCeQwKCKFFyVsHMj88OiMLPUvDGs9mkeHZAVFVWSwhRCSQgiBIN/2Y4B5IO5D0f2jmcGXd0IDzYDwUEeJuxOzTB1WltZiFEuUhAEIT6hha5Pcg7CJvDxj2r72HVXxenkwztHM7mSf05OnMQ2166kaiGdXjms52cTc2qqioLISqBDDsVJUrMSOTZDc+y9fRWhkcNZ2L3iVguWUDl4OlUbntnEz1bhTJ/XHdMkvNIiGrFIxbIKSsJCO6R48jhrT/eYsHeBXgpL3J0TqF9/M3BnNkzicm3tuP+3i3cUEshRHFkHoJwGS+TF093e5rXr329yGAAcMFu5Ya2DZj53QH2XbI+sxDCM0hAEKV2U+RNJb7+2l0dCfa38PiSP8jItldRrYQQriIBQbhMvQBv/jO8E0fOpfPKt/vcXR0hRBlJQBAu1euKMMb3bcmnv53g+z2n3F0dIUQZeMZK2MKj/PPGKLb8lchTy3YwZdVezqRkSVZUITyAtBBEmRQ3ZyHQEpj32NvLxG3RTciwOTidkoUG4qwZPLdiNyu3y3pJQlRX0kIQZZI/HTaAzWHjwTUPsj9pP39Z/6JVcCsAFv5yrNB7M2x2Zny3n8EdmxSaq7Byexyz1hwk3ppRYmuitPsJIcpO5iGICku4kMCwr4dR16cuSwYtIcASQItJ31Lc/ywfLxORoQE0D/WnRVgA5y9ks3JHPNk5jrx9ilpjYeX2OJ5bsZsMm73E/YQQBcnENFGlfj/1Ow/98BADmg/g1b6v0vvV9cQVkRI7yM/C3V0jOJaYzrHEC5xIvEC23VFEiaCA+nV88Pc242sx89fZNGz2wv9fw4P92Dypv6tPSYgao7QBQbqMhEv0aNyDxzs/zpw/5tC5YWeeHdC3yG/zUwe3L/Bt3u7QXPH86iJbExq4LqoBGTY7F7LtHDidWuSxZS0GIVxDbioLl7n/qvu5NuJaXtv6Gi0jEgtkRQ0P9iuya8dsUsWupRAe7Merd0Xz1sjOfDiuG+GyFoMQlarKAoJSqqVSap5S6vN82wKUUouUUv9VSo2uqrqIymFSJqb1nkZD/4b88+d/0q9tQF5W1M2T+hfbz1/UGgt+FnOhNRaK2g+ge2SI605CiFqsVPcQlFLzgVuBBK31Vfm2DwTmAGbgQ631zFKU9bnW+i7n4zGAVWv9tVJqmdZ6eEnvlXsInmHvub2M+HZEka+F+oYWGqkE5Rtl1DjYl7BAb3bFpjB1cHvG9Yx08ZkIUTO4+h7CQuAd4KN8BzAD7wI3ArHAVqXUKozgMOOS99+vtU4ootwIYLfzsSS/qSHah7Uv9rXEzMQitw/tHF6qkUKX7mezO3h08R+8vGovXmbF6Kubl73CQgiglAFBa71BKRV5yeYewGGt9REApdRSYIjWegZGa6I0YjGCwg7kfoYoB4vZxDujOvPIJ3/wwpd78DIphndv5u5qCeGRKnIRDgdO5nse69xWJKVUqFLqfaCzUuo55+YVwJ1KqfeAr4t533ilVIxSKubs2bMVqK6oqXy8zPy/0V3o26Y+k1bs5vNtse6ukhAeqcqGnWqtE4GHL9mWDtx3mfd9AHwAxj2ESqug8Gi+FjMfjOnKg4tiePbznVjMiiGdZLKaEGVRkYAQBzTN9zzCuU2IEq0+sppbWt7i8nJ9LWb+O7Yb9y38naeW7eCP4+f5cX+CpLkQopQq0mW0FWitlGqhlPIGRgCrLvMeUUsUlwTPS3kxceNE/rXhX6Rku35lNT9vM/PGdad5qD+LthwnzpohyfWEKKVStRCUUkuAfkCYUioWeFlrPU8pNQFYgzGyaL7Wem+l1VR4lKKGloKxPvOHuz/k/Z3vsz1hO9N6TaNH4x4uPXaAjxeZtsLpMDJsdmatOSitBCGKIbmMhFvsPrub5zY9x4mUE/iYfci0Zxba59I5C/2W9Sty2GpRcxvaz7sGk1daoX0dOYHsfWBLucos7b6u3k/KlDJLs29JSjsPQYZ6CrfoUL8Dy29dzt1t7i4yGEDhOQvFzWEoantRwSB3+/5TF7uqylJmafd19X5SppRZ2u0VJcnthNv4W/x56ZqXWH5oebH73LnqTnzMPlhMlhLL+teGf5X6uHd89iiRYQFc2aiOy8os7b6u3k/KrL1lVgYJCKJaCw8MJ9ueTbYju8T99iXuK3WZwSFnOJGeQ+wRoIQ4U5YyS7uvq/eTMmtvmZVB7iEIt+uwqEOxr+0et7vM+5Vm34OnU5n69V52ez/ksjIraz8pU8oszb4lkfUQhChBVKM6LH7waqI/Kn6fG9/4ueCGokfSFt7X1ftJmVJmFZGAINwu1De02JEU5dmvtPsqpXDkBBY7Gql1w8AC2+IzS7evq/eTMqXMovatDNJlJGq1XjPXFbnUZ1HLcpZ2X1fvJ2VKmaXZtyQy7FSIUijt4jxl2dfV+0mZUmZp9nUF85QpUyql4MrwwQcfTBk/fry7qyFqkCsb1yUixI/dccmkZeYQHuzH5NvaFTmbubT7uno/KVPKLM2+JZk6deqpKVOmfHC5/aTLSAghajjpMhJCCFEmEhCEEEIAEhCEEEI4SUAQQggBSEAQQgjh5FGjjJRSZ4HjFSgiDDjnoupUBzXtfKDmnZOcT/VX086pqPNprrWuf7k3elRAqCilVExphl55ipp2PlDzzknOp/qraedUkfORLiMhhBCABAQhhBBOtS0gXHbqtoepaecDNe+c5Hyqv5p2TuU+n1p1D0EIIUTxalsLQQghRDFqRUBQSg1USh1USh1WSk1yd31cQSl1TCm1Wym1QynlcRn/lFLzlVIJSqk9+bbVU0r9oJT60/lviDvrWFbFnNMUpVSc83PaoZS6xZ11LAulVFOl1Hql1D6l1F6l1D+c2z3ycyrhfDz5M/JVSv2ulNrpPKepzu0tlFK/Oa95y5RS3qUqr6Z3GSmlzMAh4EYgFtgKjNRau3c16wpSSh0DummtPXL8tFKqL5AGfKS1vsq57TUgSWs90xm4Q7TWE91Zz7Io5pymAGla69nurFt5KKUaA4211n8opeoA24ChwL144OdUwvkMw3M/IwUEaK3TlFIWYBPwD+BpYIXWeqlS6n1gp9b6vcuVVxtaCD2Aw1rrI1rrbGApMMTNdar1tNYbgKRLNg8BFjkfL8L4Y/UYxZyTx9Jan9Ja/+F8nArsB8Lx0M+phPPxWNqQu8amxfmjgf7A587tpf6MakNACAdO5nsei4f/J3DSwP+UUtuUUjVl1aCGWutTzsengYburIwLTVBK7XJ2KXlE98qllFKRQGfgN2rA53TJ+YAHf0ZKKbNSageQAPwA/AVYtdY5zl1Kfc2rDQGhpuqtte4C3Aw85uyuqDG00ZdZE/oz3/c7RO8AAAGsSURBVANaAZ2AU8Dr7q1O2SmlAoEvgCe11in5X/PEz6mI8/Hoz0hrbddadwIiMHpErixvWbUhIMQBTfM9j3Bu82ha6zjnvwnAlxj/ETzdGWc/b25/b4Kb61NhWuszzj9YB/BfPOxzcvZLfwEs1lqvcG722M+pqPPx9M8ol9baCqwHrgGClVJezpdKfc2rDQFhK9DaedfdGxgBrHJznSpEKRXgvCmGUioAuAnYU/K7PMIqYJzz8TjgKzfWxSVyL5xOt+NBn5PzhuU8YL/W+o18L3nk51Tc+Xj4Z1RfKRXsfOyHMXhmP0ZguMu5W6k/oxo/ygjAOYzsTcAMzNdaT3NzlSpEKdUSo1UA4AV86mnnpJRaAvTDyMx4BngZWAksB5phZLUdprX2mJu0xZxTP4yuCA0cA/6er/+9WlNK9QY2ArsBh3Pz8xj97h73OZVwPiPx3M8oGuOmsRnjC/5yrfX/Oa8RS4F6wHbgHq111mXLqw0BQQghxOXVhi4jIYQQpSABQQghBCABQQghhJMEBCGEEIAEBCGEEE4SEIQQQgASEIQQQjhJQBBCCAHA/wcghgSjD0PP1QAAAABJRU5ErkJggg==\n",
+ "text/plain": [
+ "<Figure size 432x288 with 1 Axes>"
+ ]
+ },
+ "metadata": {
+ "needs_background": "light"
+ },
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "plt.semilogy(errors_orth, marker='o', label='w. orth')\n",
+ "plt.semilogy(errors, label='wo. orth')\n",
+ "plt.semilogy(errors_greedy_orth, marker='s', label='greedy w. orth')\n",
+ "plt.semilogy(errors_greedy, label='greedy wo. orth')\n",
+ "plt.legend()\n",
+ "plt.show()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Bonus Problem"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Write a function `word_problem` which, given an arbitrary string, generates a `StationaryProblem` similar to Sheet 1, Problem 2 (c). Each letter should have a corresponding parameter component.\n",
+ "\n",
+ "**Hint:** Use the modules `PIL.ImageDraw` und `PIL.ImageFont`."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Solution"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 16,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "from PIL import Image, ImageDraw, ImageFont"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "We can use `PIL` to programmatically create bitmap graphics. However, as `BitmapFunction` always reads its data from a file, we need to save the generated images to disk. To this end, we use the `tempfile` module of the Python standard library, to create temporary files that are automatically deleted after use."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 17,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "from tempfile import NamedTemporaryFile\n",
+ "\n",
+ "def word_problem(text):\n",
+ " font = ImageFont.truetype('DejaVuSans-Bold.ttf', 64) # load some font from file of given size\n",
+ " size = font.getsize(text) # compute width and height of rendered text\n",
+ " size = (size[0] + 20, size[1] + 20) # add a border of 10 pixels around the text\n",
+ " \n",
+ " def make_bitmap_function(char_num): # we need to genereate a BitmapFunction for each character\n",
+ " img = Image.new('L', size) # create new Image object of given dimensions\n",
+ " d = ImageDraw.Draw(img) # create ImageDraw object for the given Image\n",
+ " \n",
+ " # in order to position the character correctly, we first draw all characters from the first\n",
+ " # up to the wanted character\n",
+ " d.text((10, 10), text[:char_num + 1], font=font, fill=255)\n",
+ " \n",
+ " # next we erase all previous character by drawing a black rectangle\n",
+ " d.rectangle(((0,0), (font.getsize(text[:char_num])[0] + 10, size[1])), fill=0, outline=0)\n",
+ " \n",
+ " # open a new temporary file\n",
+ " with NamedTemporaryFile(suffix='.png') as f: # after leaving this 'with' block, the temporary\n",
+ " # file is automatically deleted\n",
+ " img.save(f, format='png')\n",
+ " return BitmapFunction(f.name, bounding_box=[(0,0), size], range=[0., 1.])\n",
+ " \n",
+ " # create BitmapFunctions for each character\n",
+ " dfs = [make_bitmap_function(n) for n in range(len(text))]\n",
+ " \n",
+ " # create an indicator function for the background\n",
+ " background = ConstantFunction(1., 2) - LincombFunction(dfs, np.ones(len(dfs)))\n",
+ " \n",
+ " # form the linear combination\n",
+ " dfs = [background,] + dfs\n",
+ " coefficients = [1,] + [ProjectionParameterFunctional('chars', (len(text),), (i,)) for i in range(len(text))]\n",
+ " diffusion = LincombFunction(dfs, coefficients)\n",
+ " \n",
+ " return StationaryProblem(\n",
+ " domain=RectDomain(dfs[1].bounding_box, bottom='neumann'),\n",
+ " neumann_data=ConstantFunction(-1., 2),\n",
+ " diffusion=diffusion,\n",
+ " parameter_space=CubicParameterSpace(diffusion.parameter_type, 0.01, 1.)\n",
+ " )"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Lets test the function with the string `'pyMOR'`:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 18,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "application/vnd.jupyter.widget-view+json": {
+ "model_id": "c84bcffc613647ae87ec9fa1a7d33c07",
+ "version_major": 2,
+ "version_minor": 0
+ },
+ "text/plain": [
+ "ThreeJSPlot(children=(HBox(children=(Renderer(children=(Renderer(camera=PerspectiveCamera(position=(0.0, 0.0, …"
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "m, _ = discretize_stationary_cg(word_problem('pyMOR'), diameter=1.)\n",
+ "m.visualize(m.solve([0.01, 0.1, 0.01, 0.1, 0.01]))"
+ ]
+ }
+ ],
+ "metadata": {
+ "kernelspec": {
+ "display_name": "Python 3",
+ "language": "python",
+ "name": "python3"
+ },
+ "language_info": {
+ "codemirror_mode": {
+ "name": "ipython",
+ "version": 3
+ },
+ "file_extension": ".py",
+ "mimetype": "text/x-python",
+ "name": "python",
+ "nbconvert_exporter": "python",
+ "pygments_lexer": "ipython3",
+ "version": "3.7.2"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}
diff --git a/rb-intro/solutions/problems.py b/rb-intro/solutions/problems.py
new file mode 120000
index 0000000000000000000000000000000000000000..6114018a0f95b34b34f7c2548092f2231b20f875
--- /dev/null
+++ b/rb-intro/solutions/problems.py
@@ -0,0 +1 @@
+../problems.py
\ No newline at end of file