From f3ea0a650a107564f10b94a32649311eb82eea44 Mon Sep 17 00:00:00 2001 From: LANZIERI Denise Date: Wed, 14 Apr 2021 19:07:26 +0200 Subject: [PATCH 1/4] test for the Jacobian computed by TensorFlow --- tests/test_3D_jacobian.py | 81 +++++++++++++++++++++++++++++++++++++++ 1 file changed, 81 insertions(+) create mode 100755 tests/test_3D_jacobian.py diff --git a/tests/test_3D_jacobian.py b/tests/test_3D_jacobian.py new file mode 100755 index 0000000..8a614cd --- /dev/null +++ b/tests/test_3D_jacobian.py @@ -0,0 +1,81 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Thu Apr 8 15:19:45 2021 + +@author: Denise Lanzieri +""" + +import tensorflow as tf +import numpy as np +import flowpm +import flowpm.tfpower as tfpower +import flowpm.scipy.interpolate as interpolate +from DifferentiableHOS.pk import pk as pkl +from numpy.testing import assert_allclose + + +#%% +z_source=1. +field=5. +box_size=200. +nc=16 +Omega_c= 0.2589 +sigma8= 0.8159 +nsteps=2 + +#%% +def compute_initial_cond(Omega_c): + cosmology = flowpm.cosmology.Planck15(Omega_c=Omega_c) + k = tf.constant(np.logspace(-4, 1, 256), dtype=tf.float32) + pk = tfpower.linear_matter_power(cosmology, k) + pk_fun = lambda x: tf.cast(tf.reshape(interpolate.interp_tf(tf.reshape(tf.cast(x, tf.float32), [-1]), k, pk), x.shape), tf.complex64) + initial_conditions = flowpm.linear_field([nc, nc, nc], + [box_size, box_size, + box_size], + pk_fun, + batch_size=1) + return initial_conditions + + + +@tf.function +def compute_powerspectrum(Omega_c): + cosmology = flowpm.cosmology.Planck15(Omega_c=Omega_c) + stages = np.linspace(0.1, 1., nsteps, endpoint=True) + state = flowpm.lpt_init(cosmology, initial_conditions, 0.1) + final_state = flowpm.nbody(cosmology, state, stages, [nc, nc, nc]) + final_field = flowpm.cic_paint(tf.zeros_like(initial_conditions), final_state[0]) + final_field=tf.reshape(final_field, [nc, nc, nc]) + k, power_spectrum = pkl(final_field,shape=final_field.shape,boxsize=np.array([box_size, box_size, + box_size]),kmin=0.1,dk=2*np.pi/box_size) + return power_spectrum + + +@tf.function +def Flow_jac(Omega_c): + params = tf.stack([Omega_c]) + with tf.GradientTape() as tape: + tape.watch(params) + cosmology = flowpm.cosmology.Planck15(Omega_c=params[0]) + stages = np.linspace(0.1, 1., nsteps, endpoint=True) + state = flowpm.lpt_init(cosmology, initial_conditions, 0.1) + + final_state = flowpm.nbody(cosmology, state, stages, [nc, nc, nc]) + final_field = flowpm.cic_paint(tf.zeros_like(initial_conditions), final_state[0]) + final_field=tf.reshape(final_field, [nc, nc, nc]) + k, power_spectrum = pkl(final_field,shape=final_field.shape,boxsize=np.array([box_size, box_size, + box_size]),kmin=0.1,dk=2*np.pi/box_size) + return tape.jacobian(power_spectrum, params,experimental_use_pfor=False) + +#%% +initial_conditions=compute_initial_cond(0.2589) + +def test_Nbody_jacobian(): + theoretical, numerical_jac=tf.test.compute_gradient( compute_powerspectrum, [Omega_c], delta=0.01) + FlowPM_jac= Flow_jac(Omega_c) + assert_allclose(numerical_jac[0],FlowPM_jac, rtol=1e-1) + + + + \ No newline at end of file From b58c0324d7a798da9be010d25aec62d8dde00973 Mon Sep 17 00:00:00 2001 From: LANZIERI Denise Date: Wed, 14 Apr 2021 19:19:13 +0200 Subject: [PATCH 2/4] test for the Jacobian computed by TensorFlow --- tests/test_3D_jacobian.py | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/tests/test_3D_jacobian.py b/tests/test_3D_jacobian.py index 8a614cd..c006c81 100755 --- a/tests/test_3D_jacobian.py +++ b/tests/test_3D_jacobian.py @@ -69,9 +69,13 @@ def Flow_jac(Omega_c): return tape.jacobian(power_spectrum, params,experimental_use_pfor=False) #%% -initial_conditions=compute_initial_cond(0.2589) +initial_conditions=compute_initial_cond(Omega_c) +#%% def test_Nbody_jacobian(): + """ This function tests the Automatic differentiation implemented by TensorFlow + computing the Jacobians of the 3D power spectrum both by Automatic differentiation + and finite differences""" theoretical, numerical_jac=tf.test.compute_gradient( compute_powerspectrum, [Omega_c], delta=0.01) FlowPM_jac= Flow_jac(Omega_c) assert_allclose(numerical_jac[0],FlowPM_jac, rtol=1e-1) From 3838b3ee16c76e76c6100d57db80b99f4044e485 Mon Sep 17 00:00:00 2001 From: LANZIERI Denise Date: Wed, 14 Apr 2021 19:25:54 +0200 Subject: [PATCH 3/4] YAPF Formatting --- tests/test_3D_jacobian.py | 103 +++++++++++++++++++++----------------- 1 file changed, 56 insertions(+), 47 deletions(-) diff --git a/tests/test_3D_jacobian.py b/tests/test_3D_jacobian.py index c006c81..6faf6d9 100755 --- a/tests/test_3D_jacobian.py +++ b/tests/test_3D_jacobian.py @@ -14,72 +14,81 @@ from DifferentiableHOS.pk import pk as pkl from numpy.testing import assert_allclose - #%% -z_source=1. -field=5. -box_size=200. -nc=16 -Omega_c= 0.2589 -sigma8= 0.8159 -nsteps=2 +z_source = 1. +field = 5. +box_size = 200. +nc = 16 +Omega_c = 0.2589 +sigma8 = 0.8159 +nsteps = 2 + #%% def compute_initial_cond(Omega_c): - cosmology = flowpm.cosmology.Planck15(Omega_c=Omega_c) - k = tf.constant(np.logspace(-4, 1, 256), dtype=tf.float32) - pk = tfpower.linear_matter_power(cosmology, k) - pk_fun = lambda x: tf.cast(tf.reshape(interpolate.interp_tf(tf.reshape(tf.cast(x, tf.float32), [-1]), k, pk), x.shape), tf.complex64) - initial_conditions = flowpm.linear_field([nc, nc, nc], - [box_size, box_size, - box_size], + cosmology = flowpm.cosmology.Planck15(Omega_c=Omega_c) + k = tf.constant(np.logspace(-4, 1, 256), dtype=tf.float32) + pk = tfpower.linear_matter_power(cosmology, k) + pk_fun = lambda x: tf.cast( + tf.reshape( + interpolate.interp_tf(tf.reshape(tf.cast(x, tf.float32), [-1]), k, pk + ), x.shape), tf.complex64) + initial_conditions = flowpm.linear_field([nc, nc, nc], + [box_size, box_size, box_size], pk_fun, batch_size=1) - return initial_conditions - + return initial_conditions @tf.function def compute_powerspectrum(Omega_c): - cosmology = flowpm.cosmology.Planck15(Omega_c=Omega_c) - stages = np.linspace(0.1, 1., nsteps, endpoint=True) - state = flowpm.lpt_init(cosmology, initial_conditions, 0.1) - final_state = flowpm.nbody(cosmology, state, stages, [nc, nc, nc]) - final_field = flowpm.cic_paint(tf.zeros_like(initial_conditions), final_state[0]) - final_field=tf.reshape(final_field, [nc, nc, nc]) - k, power_spectrum = pkl(final_field,shape=final_field.shape,boxsize=np.array([box_size, box_size, - box_size]),kmin=0.1,dk=2*np.pi/box_size) - return power_spectrum + cosmology = flowpm.cosmology.Planck15(Omega_c=Omega_c) + stages = np.linspace(0.1, 1., nsteps, endpoint=True) + state = flowpm.lpt_init(cosmology, initial_conditions, 0.1) + final_state = flowpm.nbody(cosmology, state, stages, [nc, nc, nc]) + final_field = flowpm.cic_paint(tf.zeros_like(initial_conditions), + final_state[0]) + final_field = tf.reshape(final_field, [nc, nc, nc]) + k, power_spectrum = pkl(final_field, + shape=final_field.shape, + boxsize=np.array([box_size, box_size, box_size]), + kmin=0.1, + dk=2 * np.pi / box_size) + return power_spectrum @tf.function def Flow_jac(Omega_c): - params = tf.stack([Omega_c]) - with tf.GradientTape() as tape: - tape.watch(params) - cosmology = flowpm.cosmology.Planck15(Omega_c=params[0]) - stages = np.linspace(0.1, 1., nsteps, endpoint=True) - state = flowpm.lpt_init(cosmology, initial_conditions, 0.1) + params = tf.stack([Omega_c]) + with tf.GradientTape() as tape: + tape.watch(params) + cosmology = flowpm.cosmology.Planck15(Omega_c=params[0]) + stages = np.linspace(0.1, 1., nsteps, endpoint=True) + state = flowpm.lpt_init(cosmology, initial_conditions, 0.1) + + final_state = flowpm.nbody(cosmology, state, stages, [nc, nc, nc]) + final_field = flowpm.cic_paint(tf.zeros_like(initial_conditions), + final_state[0]) + final_field = tf.reshape(final_field, [nc, nc, nc]) + k, power_spectrum = pkl(final_field, + shape=final_field.shape, + boxsize=np.array([box_size, box_size, box_size]), + kmin=0.1, + dk=2 * np.pi / box_size) + return tape.jacobian(power_spectrum, params, experimental_use_pfor=False) - final_state = flowpm.nbody(cosmology, state, stages, [nc, nc, nc]) - final_field = flowpm.cic_paint(tf.zeros_like(initial_conditions), final_state[0]) - final_field=tf.reshape(final_field, [nc, nc, nc]) - k, power_spectrum = pkl(final_field,shape=final_field.shape,boxsize=np.array([box_size, box_size, - box_size]),kmin=0.1,dk=2*np.pi/box_size) - return tape.jacobian(power_spectrum, params,experimental_use_pfor=False) #%% -initial_conditions=compute_initial_cond(Omega_c) +initial_conditions = compute_initial_cond(Omega_c) + #%% def test_Nbody_jacobian(): - """ This function tests the Automatic differentiation implemented by TensorFlow + """ This function tests the Automatic differentiation implemented by TensorFlow computing the Jacobians of the 3D power spectrum both by Automatic differentiation and finite differences""" - theoretical, numerical_jac=tf.test.compute_gradient( compute_powerspectrum, [Omega_c], delta=0.01) - FlowPM_jac= Flow_jac(Omega_c) - assert_allclose(numerical_jac[0],FlowPM_jac, rtol=1e-1) - - - - \ No newline at end of file + theoretical, numerical_jac = tf.test.compute_gradient(compute_powerspectrum, + [Omega_c], + delta=0.01) + FlowPM_jac = Flow_jac(Omega_c) + assert_allclose(numerical_jac[0], FlowPM_jac, rtol=1e-1) From 4807561e8095eb0a76743516977f44da4c66625e Mon Sep 17 00:00:00 2001 From: LANZIERI Denise Date: Wed, 14 Apr 2021 20:22:29 +0200 Subject: [PATCH 4/4] notebook to compare the autodiff and finitediff --- notebooks/num_theo_jacobian.ipynb | 240 ++++++++++++++++++++++++++++++ 1 file changed, 240 insertions(+) create mode 100644 notebooks/num_theo_jacobian.ipynb diff --git a/notebooks/num_theo_jacobian.ipynb b/notebooks/num_theo_jacobian.ipynb new file mode 100644 index 0000000..4e7fd7e --- /dev/null +++ b/notebooks/num_theo_jacobian.ipynb @@ -0,0 +1,240 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Populating the interactive namespace from numpy and matplotlib\n" + ] + } + ], + "source": [ + "%pylab inline \n", + "%load_ext autoreload\n", + "%autoreload 2\n", + "import numdifftools as nd\n", + "import tensorflow as tf\n", + "import numpy as np\n", + "import flowpm\n", + "import flowpm.tfpower as tfpower\n", + "import flowpm.scipy.interpolate as interpolate\n", + "from DifferentiableHOS.pk import pk as pkl\n", + "import pickle" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [], + "source": [ + "nsteps=11\n", + "nc=32\n", + "box_size=128\n", + "Omega_c=0.2589\n", + "sigma8=0.8159" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [], + "source": [ + "def compute_initial_cond(Omega_c):\n", + " # Instantiates a cosmology with desired parameters\n", + " cosmology = flowpm.cosmology.Planck15(Omega_c=Omega_c)\n", + " stages = np.linspace(0.1, 1., nsteps, endpoint=True)\n", + " # Compute linear matter power spectrum\n", + " k = tf.constant(np.logspace(-4, 1, 256), dtype=tf.float32)\n", + " pk = tfpower.linear_matter_power(cosmology, k)\n", + " pk_fun = lambda x: tf.cast(tf.reshape(interpolate.interp_tf(tf.reshape(tf.cast(x, tf.float32), [-1]), k, pk), x.shape), tf.complex64)\n", + "\n", + " # And initial conditions\n", + " initial_conditions = flowpm.linear_field([nc, nc, nc],\n", + " [box_size, box_size,\n", + " box_size],\n", + " pk_fun,\n", + " batch_size=1)\n", + " return initial_conditions " + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "WARNING:tensorflow:@custom_gradient grad_fn has 'variables' in signature, but no ResourceVariables were used on the forward pass.\n" + ] + } + ], + "source": [ + "initial_conditions=compute_initial_cond(Omega_c).numpy()" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [], + "source": [ + "@tf.function\n", + "def compute_Nbody_deriv(Omega_c):\n", + " \"\"\" Computes a N-body simulation for a given\n", + " set of cosmological parameters\n", + " \"\"\"\n", + " # Instantiates a cosmology with desired parameters\n", + " cosmology = flowpm.cosmology.Planck15(Omega_c=Omega_c)\n", + " stages = np.linspace(0.1, 1., nsteps, endpoint=True)\n", + "\n", + " state = flowpm.lpt_init(cosmology, initial_conditions, 0.1)\n", + "\n", + "\n", + " # Evolve particles from initial state down to a=af\n", + " final_state = flowpm.nbody(cosmology, state, stages, [nc, nc, nc]) \n", + "\n", + " # Retrieve final density field i.e interpolate the particles to the mesh\n", + " final_field = flowpm.cic_paint(tf.zeros_like(initial_conditions), final_state[0])\n", + " final_field=tf.reshape(final_field, [nc, nc, nc])\n", + " params = tf.stack([Omega_c])\n", + " k, power_spectrum = pkl(final_field,shape=final_field.shape,boxsize=np.array([box_size, box_size,\n", + " box_size]),kmin=0.1,dk=2*np.pi/box_size)\n", + " return power_spectrum" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": { + "scrolled": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "WARNING:tensorflow:@custom_gradient grad_fn has 'variables' in signature, but no ResourceVariables were used on the forward pass.\n", + "WARNING:tensorflow:From /Users/dl264294/.local/lib/python3.8/site-packages/tensorflow_probability/python/math/ode/base.py:459: calling while_loop_v2 (from tensorflow.python.ops.control_flow_ops) with back_prop=False is deprecated and will be removed in a future version.\n", + "Instructions for updating:\n", + "back_prop=False is deprecated. Consider using tf.stop_gradient instead.\n", + "Instead of:\n", + "results = tf.while_loop(c, b, vars, back_prop=False)\n", + "Use:\n", + "results = tf.nest.map_structure(tf.stop_gradient, tf.while_loop(c, b, vars))\n", + "WARNING:tensorflow:@custom_gradient grad_fn has 'variables' in signature, but no ResourceVariables were used on the forward pass.\n", + "WARNING:tensorflow:@custom_gradient grad_fn has 'variables' in signature, but no ResourceVariables were used on the forward pass.\n", + "WARNING:tensorflow:@custom_gradient grad_fn has 'variables' in signature, but no ResourceVariables were used on the forward pass.\n" + ] + } + ], + "source": [ + "theoretical, numerical=tf.test.compute_gradient(\n", + " compute_Nbody_deriv, [0.2589], delta=0.01)" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [], + "source": [ + "@tf.function\n", + "def Flow_jac(Omega_c):\n", + " \"\"\" Computes a N-body simulation for a given\n", + " set of cosmological parameters\n", + " \"\"\"\n", + " params = tf.stack([Omega_c])\n", + " with tf.GradientTape() as tape:\n", + " tape.watch(params)\n", + " cosmology = flowpm.cosmology.Planck15(Omega_c=params[0])\n", + " stages = np.linspace(0.1, 1., nsteps, endpoint=True)\n", + " state = flowpm.lpt_init(cosmology, initial_conditions, 0.1)\n", + "\n", + "\n", + " # Evolve particles from initial state down to a=af\n", + " final_state = flowpm.nbody(cosmology, state, stages, [nc, nc, nc]) \n", + "\n", + " # Retrieve final density field i.e interpolate the particles to the mesh\n", + " final_field = flowpm.cic_paint(tf.zeros_like(initial_conditions), final_state[0])\n", + " final_field=tf.reshape(final_field, [nc, nc, nc])\n", + " k, power_spectrum = pkl(final_field,shape=final_field.shape,boxsize=np.array([box_size, box_size,\n", + " box_size]),kmin=0.1,dk=2*np.pi/box_size)\n", + " return tape.jacobian(power_spectrum, params,experimental_use_pfor=False),k\n", + " " + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [], + "source": [ + "jacobian,k =Flow_jac(0.2589)" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZEAAAEGCAYAAACkQqisAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAxbUlEQVR4nO3dd3hUVf7H8fc3HQg1oadSE0JC6NJBURARfiCyqBRBRFHXrrvorrjrqovrqutiw7oqIqA0RVBRqtQEIgFCqAFCJ/SE1Dm/PxIwYMokmcydJN/X88zzZG459zPDJF/OvXfOEWMMSimlVGm4WR1AKaVUxaVFRCmlVKlpEVFKKVVqWkSUUkqVmhYRpZRSpeZhdQBn8/f3NyEhIVbHUEqpCiM2NvaUMaZ+QeuqXBEJCQkhJibG6hhKKVVhiMiBwtbp6SyllFKlpkVEKaVUqWkRUUopVWpV7ppIQbKyskhOTiY9Pd3qKMoBfHx8CAgIwNPT0+ooSlV6WkSA5ORkatasSUhICCJidRxVBsYYUlJSSE5OJjQ01Oo4SlV6ejoLSE9Px8/PTwtIJSAi+Pn5aa9SKSfRIpJHC0jlof+WSjmPns5SSqnKKCud88f2kHJwJ0dTDScbdGdodFOHH0aLiFJKVVC2tDOcTt7FiVMpbPOK5EBKKn22/YXQ1C342VKohaEWcDSnDa/V8tciopRSVYoxZJ87wskjB9jt0ZIDKakExL9NcMpK/DMPU8tcwB+4YGvI05mv4+4mBFevRoZPNJm1gnHza0b1hi1pEBzO0saOLyCgRcRl3XPPPQwfPpzk5GSmTp1Ko0aNOH/+PM899xx33313gfvcd999jB07lh49ehTaXkREBA888AB79+7FGMP48eOZMmVKOb8apVShbDlkpBzgQI4/SSlpeOz4mkaHllDz0iHqZx3Bh0zqGk/GZXyMwY0pXseo7eXDId8+ZNUOxsO/Ob5NWrOyeXua1KmGp/sgp8aXqjY9bqdOncy1Y2clJCQQHh5uUaKCdezYkYULF/Lyyy8TGRnJ/fffz+bNm7nxxhtJSUkpcJ/o6GhiY2Nxd3cvsL358+czbNgwpk6dypAhQ8jIyODOO+9kwIABTJo0qbxfklO54r+pqsJsOWTk2Dh4Op0zu9ZSbddCPM7tp3baQepnH8OTbDqkv8tpajHZfRG3ef5CilcTUqsHYqsbilf95ni3uoGQBrVoUNPb6TePiEisMaZTQeu0J+Iidu3axYQJEzh37hwTJkzg2LFjBAQEEB8fz8iRIwEICAggJyenwP0TEhJo1arVlQJSUHvbt28nJCSEIUOGAODt7c306dPp06dPsUUkKSmJm2++mZ49e7J27VqaNm3KwoULOX78OIMHD2bbtm0AvPrqq1y8eJHnn3+epKQkBg4cSM+ePVm/fj3t2rVj/PjxTJ06lRMnTjBz5ky6dOniqLdQKctlZts4cmA36TuWkHVyNx5n91Mr7SD1s48yPPPvbLeFcLv7Cv7uMZdD0ogDXkEk+PXG1G3GC607ENCoIcF+N1KnuhctrH4xdtIico2/fbOdHUfOO7TNNk1qMfXWiELXZ2dnM3r0aKZPn06XLl144IEHCAsLAyA+Pp6wsDCMMbz55psMHjy4wDaWLFnCwIEDi2wvISGBdu3aXbVf48aNOX/+PJmZmXh5eRX5Onbv3s2sWbN4//33GTlyJF9//TU9e/Yscp89e/Ywd+5cZsyYQefOnfniiy9Ys2YNixYt4qWXXmLBggVF7q+Uq8m6dIGTO9dy/vBOsk7uwePsfmqmHuQN9zHMO9+GHrKVz7z+Sbrx5JA04qBXAAn1enJbaDj3BrQmpF4XMv3/Qasa3la/FIfQIuIC5s2bR3h4+JX/lUdERFCtWjUOHTrExYsXGTBgAJ6ennTp0oW33noLgAkTJvDRRx9daeP777/n448/LrI9d3d3Ll68eNWxjTGkpaXh4VH8RyE0NJTo6Ggg9/RYUlJSsUUkNDSUyMjIKzluuOEGRITIyEiSkpKKf3OUslja8X0krfuaFWcbMPtkEN5n9vCD15M0ATKMJ4doSLJXAIEN/HmoYwua12lFfM2hBAQ2o0UNb1pW8u8taRG5RlE9hvKydetWOnbseOV5bGwsffv2ZevWrdxwww0sXbr0qu3T0tKoXbs2y5cvZ+nSpUydOpWzZ8/SpEmTIttr164dd911F9OmTbtyTvXHH3+kQ4cOuLm5kZ2dzdNPP42IEBwczMMPP3zVcb29f/ufk7u7O5cuXcLDwwObzXZl+bXfFM+/j5ub25Xnl4+nlMsxhlO71nJ0w3zqHFpGYNZ+2gCbGEx488dp1bYbq+RD6gS2pmlQC5r7+tBChOuszm0RLSIuwM/P78o1hdjYWGbNmsWjjz7Kd99997vTTwCbN28mLi6O1q1bM23aNBYvXky/fv2KbS8qKor27dvz3HPP8cILL3D8+HEef/xxZsyYAcA777zD0KFD6dOnj93ZGzZsyIkTJ0hJScHX15dvv/32ymk1pSoKk5nK3oQ4vj1Zn2U7jjH91ETC5STx7uEkBDxM/Y5DuTOqA+PcLw/y0dbSvK5Ei4gLGDNmDIMGDSI6OprWrVtTp04dwsPDmTZtGoMG/f52vU2bNtG1a1dq1KgB5F4PGTFiRLHt/fOf/yQmJobPP/+cfv36MX36dA4cOMADDzzAggUL2Lx5M5MnTy5Rdk9PT5577jm6du1KaGjolWs5Srm69DNHSFo3D7PzO5qd34Sf8WJ65ju0C/JnU+fXsLWLIjowUIfRKYbe4kvFux303nvv5b333uOvf/0rAwcO5JFHHmHDhg1lHvp84cKFLFq0iHr16jFlyhTq1avnoMTOV9H+TZUTGEPKxQx+TjwJ697i9pR3ADhs6pNYpyfuYTcT0eNW/GtVtzio6ynqFl8tIugfnMpI/00VgMnO4Mivy0jZvIhGR3/mofT72WgLo69vMnf576Z2+6FEte+Gj5eelCmKfk9EKVVlZOfY2JKwm+o/TSH4zDqakkY948VWr2iGdw7luS49iWhSS09TOYgWEaVUhWXSz3Fy7xZS9m4h+0g8u7L8+fvp/qReusT33gnE1uyDrdXNhHUfTFd/P7paHbgS0iKilHJ9thwundjDwUMHibW1Zuex84zcPpm2mVtpADQAzpnq7PHsw41tGtI/vCGNWsbT3Fv/xJU3fYeVUi7FZjMkn7lESsxXeO77iRpnd9IoYz/VyKCG8eeZjDfx9fagRc3OHPW7Do/Gkfg170Cz5q0Z5uPJMKtfQBWjRUQpZY2cbFKPJnJ8dyyph+LwPJlA7bQkBmb/m3MZhr97LGKw+zqSPJqxr/ZgbA0i8A1ux+o23WlapxpubgOsfgUKLSJKKSc6cyKZT39NY/vR83Q69BGTsmbSDMgy7iRJE/b4tGZUmzqEBAYQXn863o396OBTtlvXVfnSIqKUcor9a7/G/4eHOJ51B3vqDSGoyQ386NOSmsHtCWwVTQu/2rQUoZfVQVWJaBFRSpUrk5PN9lnP0HbPeyRKM8aMvpuXwqOsjqUcRIuIUqrcpJ87yYH376TtxY2sqn4TbSd9QL06ta2OpRzIrfhNlBXuueceFi9ezHvvvUejRo2Ijo6mWbNmfPLJJ4Xuc9999/HLL78U2V5SUhKDBg2idevWtGrVipdfftnuTO7u7kRHR195XB7K3dfXtyQvTVURh06n8dL7Mwm+sIUfmk2hxxOztYBUQlpEXFRcXBzt2rVj69atPP/888TFxfHVV1/xxBNPFLrPhg0buO66ggekjouLIzIykttuu43777+fxMRE4uPjiYmJuTKKb3GqVatGXFzclUdISEhpXpqqAjZs/IXB/13D/IvhbByynJvG/hl3d/1zUxnpv6qL2LVrFz179iQyMpLXX3/9qulxL48BVdLpca9tr7DpcV999VWHvY7XXnuNtm3b0rZtW9544w0AXnnlFd58800AHnvsMa6//noAfvrpJ0aPHk1SUhJhYWGMGzeOqKgoRowYQVpa2pXlEydOpG3bttx1110sW7aMHj160LJlSzZu3Oiw3MoxbJmX2PbOWDouHkyvGsl881BPeneMtDqWKkd6TaQgH9/y+2UR/wdd7oXMNJh5++/XR98J7e+C1BSYM/bqdeMXF3m4ijI97qVLl67MbBgaGsr8+fOvWh8bG8vHH3/Mhg0bMMbQtWtX+vTpQ+/evfn3v//Nww8/TExMDBkZGWRlZbFmzRp69cq9FycxMZEPP/yQHj16MGHCBN5++21GjBih0+tWIOeP7iXl4z/QNnM3P/qP5l/3jqaaT9GfKVXxaU/EBRQ0nW10dPRV0+N26dKFM2fOXDU9bn7ff//9lSJSWHtlnR43/+msawsIwJo1axg2bBg1atTA19eX4cOHs3r1ajp27EhsbCwXLlzA29ubbt26ERMTw+rVq68UkcDAQHr06AHA6NGjWbNmDfDb9Lpubm46va4LS1q/CPNeb/wzklne4T/0f2i6FpAqQnsiBSmq5+BVvej1NfyK7Xlcq6JMj1ucwqYV8PT0JCQkhI8//pju3bsTFRXF8uXL2bt3L+Hh4Rw4cOB3I6pefq7T67q+r2KTObB4MUM9/Mi47RP6RXawOpJyIpfuiYjIv0Rkp4hsFZH5IlIn37opIrJHRBJFpEKPf1DQdLbR0dHEx8cXOT1uYmIi06ZNY/ny5cVOjxsdHU27du2uTI8LXJke96WXXgJ+mx738qmnkurduzcLFiwgLS2N1NRU5s+ff6Wn0bt3b1599VV69+5Nr169ePfdd4mOjr5SLA4ePMi6desAmDVrFj179izx8ZVzZVw8zdszZ/Pk3F+JDRhHnYdXEaEFpMpx6SIC/Ai0NcZEAbuAKQAi0gYYBUQAA4G3RcTdspRlNGbMGOLi4oiOjuaVV165Mp1tfHw8UVG//1JWQdPj5p/XvLD2Lk+P+49//IOff/6ZyZMnX5keNykpic2bN185pVQaHTp04O6776ZLly507dqViRMn0r59ewB69erF0aNH6datGw0bNsTHx+dKgQEIDw/nf//7H1FRUZw+fbrE0/Qq5zq+O4bTr3VnxK6neahXUz6deB3+detYHUtZwRhTIR7AMGBm3s9TgCn51n0PdLOnnY4dO5pr7dix43fLXNnEiRNNTk6OeeaZZ8yqVatM+/btTWZmZpnbXbBggZkwYYJ58sknTUpKigOS2mf//v0mIiLCoW1WtH/TimTn0hkmbaq/OfZcsFm3/Fur4ygnAGJMIX9TK9I1kQnA7LyfmwLr861LzltWJbz//vsAvPjii0Du6S1HGDp0KEOHDnVIW6rysWVnsf2jyUQemcuv7m2pPfZzrgsOtTqWspjlRURElgGNClj1rDFmYd42zwLZwMzLuxWwfaGTxYvIJGASQFBQUJnyKscLCQm5cg1Huabz6Vk8MTuOWw8lc9z/D3Sb9CY1qvlYHUu5AMuLiDGmf1HrRWQcMBi4Ia9bBbk9j8B8mwUAR4o4xgxgBkCnTp0KLTZKqd87GLuEKctOs+FcXbrd/BbjezbT+cnVFS59YV1EBgJ/AoYYY9LyrVoEjBIRbxEJBVoC+vVlpRzJGLbP/TtNF93B+PTPmDXpOib0aq4FRF3F8p5IMaYD3sCPeR/c9caY+40x20VkDrCD3NNcDxpjCh4PxE7GGP3lqCR+67Cq0spMPcveGWOJOLeSddV60W7iJ9T3r2d1LOWCXLqIGGNaFLHuReBFRxzHx8eHlJQU/Pz8tJBUcMYYUlJS8PHR8/WldSQpEfPZMFpmH+WHwIfpd/fzeHpU2DvoVTlz6SLiLAEBASQnJ3Py5EmroygH8PHxISAgwOoYFUbWhVMkbfyG3ckn+O/Z7uw7eor3vetzoM/L3HS93q2niqZFhNxhOUJD9VZFVXWc2LOZoxu+xvfgz4SkJ9BSDDZbILWa9uaRgZGERi0hsF51q2OqCkCLiFJVQMaFFPZtWMy89A6s2JXC2NNvMsZjGTukBSsajqN6xM1EdO7L7Op6GlCVjBYRpSojm41juzZyNOYbfA8tp1n6DsLFsDn7RRqEdMK97ZPsbfkq4SEhtNHrgKoMtIgoVUmkn08hNukkyw7kkLr9B15Jf55GQIJbC1Y1Godv25v5tFM/alTzLrYtpeylRUSpCsrYbBzeuZHjm7/B99AKmqfvIDZnCF8wil4hnVhV+wWCu9xKWFAI4drbUOVEi4hSFcy5S1n858dExm+5nUBzhABgl1tz1jYeQ/fo4Uzq2AsfT3egV3FNKVVmWkSUqkD27NrOPQtPcPjMJcIb3sLBpoEEdxlKq6AQWlkdTlVJWkSUqgiMYeu8VwjbOo3uHo8y4r6H6Bg8yOpUSmkRUcrVZaaeI/H98USd/YkYn648PnES9evrECTKNWgRUcqFndoXR/rM0bTJTubngPvpNf5FPD3011a5Dv00KuWi1u9LYeHnC3ki5zwbe3/M9TcMszqSUr+jRUQpF2Oy0vl2ybc8ur46wX59ODfyXroFNrE6llIF0iKilAtJO7GfEx+N4qZLe7m95ec8e0cPavp4Wh1LqUJpEVHKRRyJ+QbfxZPxs2WzMvIlXr5tgE5NoFyeFhGlXMDuOc/SfPtb7JVALgz9iJs6dLY6klJ20SKilIWyc2z864dEam09SFiNG2gzcQYt6/tZHUspu2kRUcoiZ3f9wms/7OTT5EaM6fokEwdH4O2pv5KqYtFPrFLOZgyHvn+DRuv/wVDTnKgRXzOiU6DVqZQqFS0iSjmRST/PgU/uJeTYUta4d6beXR8yopkWEFVxaRFRyknSzxzl3Ds3EZhxiLl17+HGe1+iTg2dSVBVbFpElHKCgylpTP5sF+MuheLWbgq3Db8DNze9fVdVfFpElCpP2Rkc/OoZJu7swDH8qX/XB/QLa2B1KqUcRouIUuXEdvoAxz+6g6CL2xlSbTJDJj5HkF91q2Mp5VBaRJRyNGO4uGUefPsINXKy+TjwBe4Z+yDVvNytTqaUw2kRUcqBcmyGmK/+RdcdL7LDBLOr93Tuvr6XDl+iKi0tIko5QsZFtibs5M8rL3HoaCiPNniQ7n94kv9rqpNHqcpNi4hSZWEMZzbNhh/+glemD+d8XuflO3twS+QI7X2oKkGLiFKllHFkOyfmPELg2U3sMCH82u4vLLv1er32oaoULSJKlZAxho2rFtNx+VhqGW9m1X+EnqOe4g7/mlZHU8rptIgoZS9jSNqbwF9XXmDtbhvP176dVrc+zh1tW1udTCnLaBFRyg4XDv7K6TkP43thH3vlTf5yazSjrhuMp7ub1dGUspQWEaWKYEs7y545z9Is6QtqmeosD5jMN38YgF8t/dKgUqBFRKlCbU1IIHDOQFrYzvFj9ZsJGvEyw5uHWB1LKZeiRUSpa5w4fph/rjjJvC3J/K16X4L73cVNfW7UW3aVKkCFOKErIk+KiBER/3zLpojIHhFJFJEBVuZTlUPGxdPEvz8J37fb8+vWLTzQtwUj/vQhffvepAVEqUK4fE9ERAKBG4GD+Za1AUYBEUATYJmItDLG5FiTUlVoNhsJS9+h0cZptDHnWVlrCJ+MuonApk2tTqaUy3P5IgK8DjwNLMy3bCjwpTEmA9gvInuALsA6C/KpCmzv8bNkf3QL4Rnb2OYWxt4Bn3F91z5Wx1KqwnDpIiIiQ4DDxphfrzmd0BRYn+95ct6ywtqZBEwCCAoKKoekqqIxxvD6j7t4Z+VeHvMI42jESLoPewgvT/22uVIlYXkREZFlQKMCVj0LPAPcVNBuBSwzhR3DGDMDmAHQqVOnQrdTVcf6r98gdnMGg6IGMOKWt2hQU6epVao0LC8ixpj+BS0XkUggFLjcCwkANotIF3J7HoH5Ng8AjpRzVFVJ7P11FR3jX+Cx2tfRYeSfdJpapcrAZe/OMsbEG2MaGGNCjDEh5BaODsaYY8AiYJSIeItIKNAS2GhhXFVBXDx3mmoL7uW01CV0wkdaQJQqI8t7IqVhjNkuInOAHUA28KDemaWKY2w2dn8wgUjbCXbePJu29Qs6i6qUKokKU0TyeiP5n78IvGhNGlURrVoyiz4XlvNL6EP0uK6gS21KqZJy2dNZSjlS4rEL3Le+Lq/X+wvXjfm71XGUqjS0iKhKL+3iOf7++RJ8vb24a/wfcXfX23iVcpQKczpLqdLa/sF9vHVhFYkjV+mtvEo5mPZEVKW2acHbdD67hJ2Bo+ga0cLqOEpVOlpEVKV1IDGONlueZ4dXJJ3vnmZ1HKUqJS0iqlJKv5RKzuxxZIoX/uM+w93D0+pISlVKWkRUpfTS4u3EZAZzqM9rNGgaanUcpSotLSKq0lm4JZlPY06yv+crRPUbaXUcpSo1LSKqUjm0L4GgBcMY0vQiT9zYyuo4SlV6WkRUpZGefom0mWNpIck8OyQaD3f9eCtV3vS3TFUaMR8+RuucXezvPo2Gwa2tjqNUlaBFRFUKm36YRc+Ts4hpMJyom8ZaHUepKqNURUREvBwdRKnSOpiSRvov77HfvRlRE96yOo5SVUppeyJ/u/yDiPRwUBZVAf28aiXzlywhPcuakfgzsnN4aNZmHpGn8Br3NV4+1S3JoVRVVdqxs77Pm7c8FQgHfnFcJFVRpGdmEf7TeBpLCv+OuYemAx5lRMcAp17Q/mbmdA4kN+KVMX1oGqTzgyjlbCX+bReRD4EhwCggwhjzF4enUhXCltXf0VhS2OfXh8Q6PfnzvHjufm02y9fHYEz5T2W/afl8hu+byttBKxgQoQVEKSuUuCdijLlHRKoBHYDOIvKeMeY+x0dTri4rbjZp+BB87xe8512DH3ccx3f+GDou2cw3K26h8eBn6Nw2rFyOfST5ACErHyXZI4BOOi6WUpYpticiIuNE5JSInBaRT0WkpjHmkjHmF2PMG1pAqqZzF9Npe34Ve+v1wd3HFxHhpohGdH3wIw4HDWVQ+jeEz+3D/NceZPu+ZIceOys7mxP/G0dNk4rHyE/wrl7Loe0rpexnz+msvwI3AmHAAeClck2kKoSlCScYkPFPPPtffTbTvW4gze75iJz713GiYS+Gnf+cpR9O5cEvNrP/VKpDjr3yf38jOmsLuzr8hSatOzmkTaVU6dhzOuu8MWZL3s9/FZEN5RlIVQwL447g6x9A6/DIAtd7Nwqj2QNfcTFpE9473Fi+/gQXt//A0FAbPUY8TMM6vqU67s87jzNldxg5ofczYMgjZXkJSikHsKcn0lhEJolILxGpD+iY2lXc8ZOnmHjwT0xqloKIFLmtb0hnHhrUkZVP9eORhr8yPHkaaa93YuHM6ZxLyyjRcY+dPMmTszfj3ziIPhNegmKOrZQqf/YUkalAFPAPIBFoKyLficjLInJHuaZTLilh+Uyud4+jd4t6du9Tv6Y3Hf74BSdu+RhvHx+G7n6WI690ZdGC2XZ9xyQ7O4dDH4zmTdvLvHVHND6eOk+6Uq7AniJyCHjJGNPHGFMPaAZMB84Ct5RjNuWiau5ewHG3hjSN7FuyHUVo0Hk4Tf4US3Lf1/F3v8SaTZvo+68VfLnhANk5tkJ3XfnZC3TOWE/tqEE0a1CzbC9AKeUw9hSR24BFInJIRH4EHgHqAYuBceUZTrmepAP7aZcZx9GgwaU/neTmTkDfCdSfEs/w8U/RuI4POxe9yoaXB7B6zcrffcckdt3P9Ep6kx01exI5/E8OeBVKKUcptogYYyYaYzoB/wZ2AfuBfsBGcu/WUlXI/hWf4yE2Ano7YJBDDy+ua9GQeZO7M7JLCO2yt9Pjx6Gs+OcwYrdsBuDEyRM0/P5+zrrVJXTi//Q6iFIupiRfNhxvjGl3+YmIvA085fhIylUZY1hzxOBdrS/dm0U7rF0Roc3QJ8jpP55d816g+97PkQX9+XzlvWx0b8/jxg2GvUu12v4OO6ZSyjFKMuzJeRHpePmJMSYW0KnjqpD4w+f48FxHDl3/33Jp371GPcLGvI55eAt7A4ax6lwDFh2uSezgpQS0u75cjqmUKpuS9EQmAJ+LyA4gFogEssollXJJa9etoY57FgPbNi7X4/jUCyD83g95NT2LxGMX6Bxi/11gSinnsrsnYozZDXQHvgMaAgnAoHLKpVxMTo6NAduf5ota/6V2Ned8VaiWj6cWEKVcXIkGYDTG5ABz8x6qComPXU00h4kPm2R1FKWUC9HpcZVdzm38gizjTsu+d1kdRSnlQrSIqGKlZ2QSduoHdtXsik/t+lbHUUq5EC0iqlhx65bRkNO4tfuD1VGUUi5Gi4gq1qfJDRjl/iote42wOopSysVoEVFFupCexbKdJwlr1x0Pn9IN366UqrxcvoiIyB9FJFFEtovIK/mWTxGRPXnrBliZsTKL+3kuL8tbDA+rZnUUpZQLKvEc684kIv2AoUCUMSZDRBrkLW8DjAIigCbAMhFplXcLsnIgr/hZXO++lTrNAqyOopRyQa7eE5kM/NMYkwFgjDmRt3wo8KUxJsMYsx/YA3SxKGOldfLUSdqlrSOp0U2Ih5fVcZRSLsjVi0groJeIbBCRlSLSOW95U3LnObksOW9ZgfJmZowRkZiTJ0+WY9zKZefyL/CRLOp3H211FKWUi7L8dJaILAMaFbDqWXLz1QWuAzoDc0SkGVDQeOCmgGW5K4yZAcwA6NSpU6Hbqav57lrAUbeGBET2sTqKUspFWV5EjDH9C1snIpOBeSZ3lqKNImID/MnteQTm2zQAOFKuQauYA6cusuFSE6LC+tJY5/BQShXC1U9nLQCuBxCRVoAXcApYBIwSEW8RCQVakjtJlnKQRb8eZVrOnYTe+rTVUZRSLszynkgxPgI+EpFtQCYwLq9Xsl1E5gA7gGzgQb0zy3GMMSTELqdrcBsa19Zbe5VShXPpImKMyQQKvKprjHkReNG5iaqG3Tu38nbaU8QG/wnoaXUcpZQLc/XTWcoCR9d8js0IzXvfYXUUpZSL0yKirmLLsRF8eDG7q0VSp3Go1XGUUi5Oi4i6yrbNawjhMOlht1kdRSlVAWgRUVc5tWlu7uRT/XTyKaVU8bSIqCsys208eWIQb4a+S3WdfEopZQctIi4gO8dGepb1dyiv3HWS0+k2OnTra3UUpVQFoUXEBfztmx30+ddyjp1LtzRH9rIXeLraQnq28Lc0h1Kq4tAiYrFLmTnM25xMt4s/8e0HU8nMtlmS42JqKt1T5tG9zjk83fVjoZSyj/61sNj6mA0MyvmJyX6bufv8DL6Y/bklObat/Irakkr1TqMsOb5SqmLSImKx9I2f8rLnB7S4+13OVAtk8K5n+X5tjNNzuG+byxlq0aLrYKcfWylVcWkRsdC51Ewiz/xIUq0uuPuFUmf8HKq7ZdH4+0kkJp8ovgEHOZ1yiqjU9exrOAA3D0+nHVcpVfFpEbHQpjVLCZBTeLX/AwCeDcPIvPVtomQvsz57j/PpWU7JsfzXPSy1dca/+1inHE8pVXloEbGQbescMvAisNuIK8vqdBhO/JAlfHahI0/M+RWbrfzn0Poi0cY7flMIbte73I+llKpctIhY5MT5dLwvHCDJrzfiU+uqdZEduvPMoHCOJ6xl3rcLyjXH4SOHOXtwG0Oim5TrcZRSlZMWEYt8s/Uo4zL/jMdt7xW4fkK3QN6r+QE9Yx9jw9aEcsux76cP+cn7Kf4vKLPcjqGUqry0iFjk27hDRDSpRfMmBX+xT9w9qDP6U+pIKp7zxnM45Xy55Kif9A37PJrTpFmbcmlfKVW5aRGxwIFjp3jv5Bj+XH9tkdtVC4rm3I2v0YEEYt9/iIxsxw6Nsnfnr4Tl7OJMs6EObVcpVXVoEbHA9hVzaSBnCY9oX+y2DXuMIanFOIakL+SLmR87NMeRNZ9hM0Kz68c5tF2lVNWhRcTJjDH47lnIGbe6+Le9wa59Qu74N0tDnuaFhIbMjTnkkBw2m6HR4R9J9ImibqMQh7SplKp6tIg42c6kQ3TNiuF44CBwc7dvJ3dP+o+ZQtdm9fnPgtUk7D1Q5hybD55h+KW/kNzz5TK3pZSqurSIONneVbPxliya9BxTov083N3478g2zPGYyvmZ4zh78VKZciyIO0yWpy/du3QtUztKqapNi4gT2WyGz5Ib8XWdCdRqcV2J9/evU4us7o/S1baF1TMeL/UXEbOysrgx7hEeDdpHDW+PUrWhlFKgRcSpYg6cYcOFenj0fRJEStVG8E0PsjtgOLee/4LFcz8oVRvb1n5HH2LpHlS9VPsrpdRlWkScaNvqBfT33Er/sAZlaqfFuLc56BNG3x1/Zf2m9SXeP33zl6TiQ1jv28uUQymltIg4SVaOjU773ub56nOp4VO2kXLFsxoNJs5hvVc3nlx8lEOn0+zeNy3tIhFnV5BYty9e1XzLlEMppbSIOEnMls1EsZtLrYY5pD0f/2Ba3f8556jOQ5+tJz0z2679tq/4ilqSRvUOOvmUUqrstIg4Scr6WQAE9ynZXVlFCfarwX+Ht+BvKU/w8wdTMKb4C+2rDqaz0q0Lrbrp5FNKqbLTIuIElzJzaH1yKfurR+HlH+LQtvtGNsOnQQsGHH+fFd/NLnLbM6mZvHMwiLWd3tTJp5RSDqFFxAlWxe2gJqmYtiOK37ikRGg58WMOe4XQfuPjJOzYWuima9avpZbtnA77rpRyGC0iTjA3IYNhXu8R3P++cmnf3ceX2uNm4ybgPncsp8+eLXC74A3Ps6D6C7RpXKvA9UopVVJaRMrZudRMVu86xi3tAnH38im349QOaE3KgOmILYu/fbmKnGu+iHjscBIRGXEcDxiIlPI7KkopdS0tIuVs45qlrPJ4kFFNT5X7sUK7DWfLoG9YmOTO6z/uumrd3uWf4S6Gpr10HnWllONoESlntl/nUEfSaB4e7ZTjjezajDEd6xO0+ik2rfz2ynL//QvZ69GcJi2dk0MpVTVoESlHx89epGPqSg4UMI96eXp2YEu6e+0h9OcHOJi0h/17E2mds5vToTr5lFLKsbSIlKPNKxbgL+ep2eUOpx7Xp2Zd3O+YSXVJJ/WzO5m1/RK9Mv5D6I3lc2FfKVV1uXQREZFoEVkvInEiEiMiXfKtmyIie0QkUUQGWJmzMN4J87goNWjc8VanH7txy/bs7/kq4TmJBG38OyEtwvFv0MjpOZRSlZtLFxHgFeBvxpho4Lm854hIG2AUEAEMBN4WETtneHKO/adSef/CdcS2fhI8vC3JENF/LHFNRjHSfQWjW9o3LIpSSpWEqxcRA1y+mFAbOJL381DgS2NMhjFmP7AH6FLA/pZZFHeE9SaCVjdPtjRH1IS3OHTLZ9zUs+TzlyilVHFcfUaiR4HvReRVcgte97zlTYH8Y6An5y0rkIhMAiYBBAUFlUvQ/IwxpMV8xm1Nw2lcu1q5H68obh4eNO9yi6UZlFKVl+U9ERFZJiLbCngMBSYDjxljAoHHgA8v71ZAU4WOPmiMmWGM6WSM6VS/fn3Hv4hrJOw/xOOX3mJijTXlfiyllLKS5T0RY0z/wtaJyKfAI3lP5wKXp/JLBgLzbRrAb6e6LLdv9WzaSDZNejluxF6llHJFlvdEinEE6JP38/XA7ryfFwGjRMRbREKBlsBGC/L9js1maJC0iOMeTajVvKvVcZRSqlxZ3hMpxr3Af0TEA0gn77qGMWa7iMwBdgDZwIPGmBzrYv4mLiGRjrZ49rS8j4Y6RpVSqpJz6SJijFkDdCxk3YvAi85NVLytMasIw4ugvjpGlVKq8nP101kVSma2jTcOhPDXVgup1iTC6jhKKVXutIg40JpdxziblsWgDs2sjqKUUk7h0qezKpq0n19lic/PNA/5xeooSinlFNoTcZC0jCzCTn6PT/WaeFWrYXUcpZRyCi0iDrJh/WpaSDJElsM86kop5aK0iDhI+uYvycad4F53WR1FKaWcRouIA5y5mE7k2Z9Iqt0FN19/q+MopZTT6IV1B1gaf4SYrBE81KO31VGUUsqptIg4wML445zwu5mQzn2K31gppSoRPZ1VRkdPnyfswCxGhXsjOsyJUqqK0Z5IGcWtXMjznv/jmJ/2QpRSVY/2RMro8jzqjToOtjqKUko5nRaRMth39CRdMtZyuPGNls2jrpRSVtIiUgY7VszFV9Kp32201VGUUsoSWkRKyRhDyr44Utz8qRdxvdVxlFLKElpESmnb4fNMvTCUn/ovBjd3q+MopZQltIiU0qItB/F0F26K1mHflVJVl97iWwo5NkO/zX+kR72m1Kk+yOo4SillGe2JlMKW7TvpaoujcZNAq6MopZSltIiUwtF1X+AuhqA+Y6yOopRSltIiUkKZ2TaCjywh2buFzqOulKrytIiU0KbNsUSxm/SwYVZHUUopy2kRKaFvEi/yH7mL4D5jrY6ilFKW0yJSAqkZ2SzclcGJqMl41guyOo5SSllOi0gJrNu0gf45qxna1s/qKEop5RK0iJRA5sb/8brXO3Rq7GV1FKWUcglaROx05mI67c79RFLtrjqPulJK5dEiYqdNq5fQVE7h3X6k1VGUUsplaBGxk4n/inS8CLjuNqujKKWUy9AiYocTF9Kpe3EPB/z7ID61rI6jlFIuQwdgtEODmj5ceGgZNT0yrY6ilFIuRYuInZo3qGl1BKWUcjl6OksppVSpaRFRSilVai5RRETkdhHZLiI2Eel0zbopIrJHRBJFZEC+5R1FJD5v3ZsiIs5PrpRSVZtLFBFgGzAcWJV/oYi0AUYBEcBA4G0RuTyh+TvAJKBl3mOg09IqpZQCXKSIGGMSjDGJBawaCnxpjMkwxuwH9gBdRKQxUMsYs84YY4BPgf9zXmKllFLgIkWkCE2BQ/meJ+cta5r387XLlVJKOZHTbvEVkWVAowJWPWuMWVjYbgUsM0UsL+zYk8g99UVQkA7hrpRSjuK0ImKM6V+K3ZKBwHzPA4AjecsDClhe2LFnADMAOnXqVGixUUopVTKu/mXDRcAXIvIa0ITcC+gbjTE5InJBRK4DNgBjgf/a02BsbOwpETmQ99QfOFUOuZ2lIuevyNlB81upImeHipk/uLAVLlFERGQYuUWgPrBYROKMMQOMMdtFZA6wA8gGHjTG5OTtNhn4BKgGLMl7FMsYUz/fcWOMMZ2K2t6VVeT8FTk7aH4rVeTsUPHzX8sliogxZj4wv5B1LwIvFrA8BmhbztGUUkoVwdXvzlJKKeXCqnoRmWF1gDKqyPkrcnbQ/FaqyNmh4ue/iuR+V08ppZQquareE1FKKVUGWkSUUkqVWpUoIiIyMG8U4D0i8ucC1oeJyDoRyRCRJ63IWBg7st8lIlvzHmtFpJ0VOQtjR/6hednjRCRGRHpakbMwxeXPt11nEckRkRHOzFcUO977viJyLu+9jxOR56zIWRh73vu81xCXNwr4SmdnLIod7/9T+d77bXmfn3pWZC0TY0ylfgDuwF6gGeAF/Aq0uWabBkBncm8lftLqzCXM3h2om/fzzcAGq3OXML8vv12biwJ2Wp27JPnzbfcz8B0wwurcJXjv+wLfWp21DPnrkPsdsqC85w2szl3Sz06+7W8FfrY6d2keVaEn0gXYY4zZZ4zJBL4kd3TgK4wxJ4wxm4AsKwIWwZ7sa40xZ/Kerufq4WCsZk/+iybvtwioQRFjoFmg2Px5/gh8DZxwZrhi2JvdVdmT/05gnjHmIOT+Hjs5Y1FK+v7fAcxySjIHqwpFpLCRgCuCkma/Bzu/ue8kduUXkWEishNYDExwUjZ7FJtfRJoCw4B3nZjLHvZ+drqJyK8iskREIpwTzS725G8F1BWRFSISKyJjnZaueHb/7opIdXLnQ/raCbkcziW+sV7OSjTir4uxO7uI9CO3iLjSNQW78pu8EQtEpDfwAlCawTrLgz353wD+ZHLHcyv/RPazJ/tmINgYc1FEBgELyB2fzhXYk98D6AjcQO7wR+tEZL0xZld5h7NDSf7u3Ar8Yow5XY55yk1VKCKFjQRcEdiVXUSigA+Am40xKU7KZo8SvffGmFUi0lxE/I0xrjBAnT35OwFf5hUQf2CQiGQbYxY4JWHhis1ujDmf7+fvROTtCvbeJwOnjDGpQKqIrALaAa5QREry2R9FBT2VBVSJC+sewD4glN8ucEUUsu3zuNaF9WKzA0HkzvjY3eq8pczfgt8urHcADl9+bvWjJJ+dvO0/wXUurNvz3jfK9953AQ5WpPceCAd+ytu2OrnTbLe1OntJPjtAbeA0UMPqzKV9VPqeiDEmW0QeAr4n946Jj0zu6MD3561/V0QaATFALcAmIo+SeyfF+cLadQZ7sgPPAX7kzj8PkG1cZIRQO/PfBowVkSzgEvAHk/fbZTU787skO7OPACaLSDa57/2oivTeG2MSRGQpsBWwAR8YY7ZZl/o3JfjsDAN+MLm9qQpJhz1RSilValXh7iyllFLlRIuIUkqpUtMiopRSqtS0iCillCo1LSJKKaVKTYuIUmUgIiEiUuRtpflGy/0u3z5GRF7It42/iGSJyPRS5pgpIqddaRRhVTVoEVHKOVYbYwble74PGJzv+e3A9tI2boy5C1hU2v2VKi0tIko5iIg0E5EtItLZjs0vAQkicvmLoX8A5uRr6xMReVdEVovILhEZnLfcXUReFZH4vHlY/uj4V6KU/Sr9N9aVcgYRaU3ucN/jjTFxdu72JTBKRI4BOeSOrdQk3/oQoA/QHFguIi2A8eQOpdE+71vRFW8SI1WpaBFRquzqAwuB24wxJTkltZTcUYuPA7MLWD/HGGMDdovIPiCM3BGO3zXGZAOYCjryq6o89HSWUmV3jty5I3qUZCeTO1lRLPAEBc8lce2YRIbcIcZ1rCLlMrSIKFV2mcD/kTuQ5J0l3Pff5M5HUtAQ/reLiJuINCd3mtVE4AfgfhHxANDTWcpqejpLKQcwxqTmXfz+UURSjTEL7dxvO4XflZUIrAQaAvcbY9JF5ANyZ/Tbmjfy8ftAqW4LVsoRdBRfpcqZiPQld56awcVsmn+fT4BvjTFflec+SpWVns5SqvxlAm0vf9mwPIjITHLv5Eovr2MoVRDtiSillCo17YkopZQqNS0iSimlSk2LiFJKqVLTIqKUUqrUtIgopZQqtf8HS8KiJJGDL0QAAAAASUVORK5CYII=\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "plot(k, numerical[0], label=r'${d P_k}/{d \\Omega_c}$ Finite Difference')\n", + "plot(k, jacobian[:,0],'--',label=r'${d P_k}/{d \\Omega_c}$ Automatic Differentiation')\n", + "xlabel('k [Mpc]')\n", + "ylabel('$dP_k$')\n", + "legend()\n", + "savefig('psfinit.png',dpi=175)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "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.8.5" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +}