{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Gradient vs finite difference"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[CpuDevice(id=0)]\n"
     ]
    }
   ],
   "source": [
    "# NBVAL_SKIP\n",
    "from jax import config\n",
    "import os\n",
    "import jax\n",
    "\n",
    "print(jax.devices())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# NBVAL_SKIP\n",
    "import os\n",
    "os.environ['SPS_HOME'] = '/home/annalena/sps_fsps'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Load ssp template from FSPS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2025-11-10 17:11:57,608 - rubix - INFO - \n",
      "   ___  __  _____  _____  __\n",
      "  / _ \\/ / / / _ )/  _/ |/_/\n",
      " / , _/ /_/ / _  |/ /_>  <\n",
      "/_/|_|\\____/____/___/_/|_|\n",
      "\n",
      "\n",
      "2025-11-10 17:11:57,608 - rubix - INFO - Rubix version: 0.0.post626+g42b4b7505.d20251110\n",
      "2025-11-10 17:11:57,609 - rubix - INFO - JAX version: 0.7.2\n",
      "2025-11-10 17:11:57,609 - rubix - INFO - Running on [CpuDevice(id=0)] devices\n",
      "2025-11-10 17:11:57,609 - rubix - WARNING - python-fsps is not installed. Please install it to use this function. Install using pip install fsps and check the installation page: https://dfm.io/python-fsps/current/installation/ for more details. Especially, make sure to set all necessary environment variables.\n"
     ]
    }
   ],
   "source": [
    "# NBVAL_SKIP\n",
    "from rubix.spectra.ssp.factory import get_ssp_template\n",
    "ssp_fsps = get_ssp_template(\"FSPS\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(107,)\n",
      "(12,)\n"
     ]
    }
   ],
   "source": [
    "# NBVAL_SKIP\n",
    "age_values = ssp_fsps.age\n",
    "print(age_values.shape)\n",
    "\n",
    "metallicity_values = ssp_fsps.metallicity\n",
    "print(metallicity_values.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "start age: 15.848933219909668, start metallicity: 0.025251565501093864\n",
      "target age: 3.1622776985168457, target metallicity: 0.014199999161064625\n"
     ]
    }
   ],
   "source": [
    "# NBVAL_SKIP\n",
    "index_age = 90\n",
    "index_metallicity = 9\n",
    "\n",
    "#initial_metallicity_index = 5\n",
    "#initial_age_index = 70\n",
    "initial_metallicity_index = 10\n",
    "initial_age_index = 104\n",
    "\n",
    "learning_all = 1e-2\n",
    "tol = 1e-10\n",
    "\n",
    "print(f\"start age: {age_values[initial_age_index]}, start metallicity: {metallicity_values[initial_metallicity_index]}\")\n",
    "print(f\"target age: {age_values[index_age]}, target metallicity: {metallicity_values[index_metallicity]}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Configure pipeline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# NBVAL_SKIP\n",
    "from rubix.core.pipeline import RubixPipeline\n",
    "import os\n",
    "config = {\n",
    "    \"pipeline\":{\"name\": \"calc_gradient\",},\n",
    "    \n",
    "    \"logger\": {\n",
    "        \"log_level\": \"DEBUG\",\n",
    "        \"log_file_path\": None,\n",
    "        \"format\": \"%(asctime)s - %(name)s - %(levelname)s - %(message)s\",\n",
    "    },\n",
    "    \"data\": {\n",
    "        \"name\": \"IllustrisAPI\",\n",
    "        \"args\": {\n",
    "            \"api_key\": os.environ.get(\"ILLUSTRIS_API_KEY\"),\n",
    "            \"particle_type\": [\"stars\"],\n",
    "            \"simulation\": \"TNG50-1\",\n",
    "            \"snapshot\": 99,\n",
    "            \"save_data_path\": \"data\",\n",
    "        },\n",
    "        \n",
    "        \"load_galaxy_args\": {\n",
    "        \"id\": 14,\n",
    "        \"reuse\": True,\n",
    "        },\n",
    "        \n",
    "        \"subset\": {\n",
    "            \"use_subset\": True,\n",
    "            \"subset_size\": 2,\n",
    "        },\n",
    "    },\n",
    "    \"simulation\": {\n",
    "        \"name\": \"IllustrisTNG\",\n",
    "        \"args\": {\n",
    "            \"path\": \"data/galaxy-id-14.hdf5\",\n",
    "        },\n",
    "    \n",
    "    },\n",
    "    \"output_path\": \"output\",\n",
    "\n",
    "    \"telescope\":\n",
    "        {\"name\": \"TESTGRADIENT\",\n",
    "         \"psf\": {\"name\": \"gaussian\", \"size\": 5, \"sigma\": 0.6},\n",
    "         \"lsf\": {\"sigma\": 1.2},\n",
    "         \"noise\": {\"signal_to_noise\": 100,\"noise_distribution\": \"normal\"},\n",
    "         },\n",
    "    \"cosmology\":\n",
    "        {\"name\": \"PLANCK15\"},\n",
    "        \n",
    "    \"galaxy\":\n",
    "        {\"dist_z\": 0.1,\n",
    "         \"rotation\": {\"type\": \"edge-on\"},\n",
    "        },\n",
    "        \n",
    "    \"ssp\": {\n",
    "        \"template\": {\n",
    "            \"name\": \"FSPS\"\n",
    "        },\n",
    "        \"dust\": {\n",
    "                \"extinction_model\": \"Cardelli89\",\n",
    "                \"dust_to_gas_ratio\": 0.01,\n",
    "                \"dust_to_metals_ratio\": 0.4,\n",
    "                \"dust_grain_density\": 3.5,\n",
    "                \"Rv\": 3.1,\n",
    "            },\n",
    "    },        \n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/annalena/.conda/envs/rubix/lib/python3.12/site-packages/rubix/telescope/factory.py:26: UserWarning: No telescope config provided, using default stored in /home/annalena/.conda/envs/rubix/lib/python3.12/site-packages/rubix/telescope/telescopes.yaml\n",
      "  warnings.warn(\n",
      "2025-11-10 17:11:58,256 - rubix - INFO - Getting rubix data...\n",
      "2025-11-10 17:11:58,257 - rubix - INFO - Rubix galaxy file already exists, skipping conversion\n",
      "/home/annalena/.conda/envs/rubix/lib/python3.12/site-packages/jax/_src/numpy/scalar_types.py:50: UserWarning: Explicitly requested dtype float64 requested in asarray is not available, and will be truncated to dtype float32. To enable more dtypes, set the jax_enable_x64 configuration option or the JAX_ENABLE_X64 shell environment variable. See https://github.com/jax-ml/jax#current-gotchas for more.\n",
      "  return asarray(x, dtype=self.dtype)\n",
      "/home/annalena/.conda/envs/rubix/lib/python3.12/site-packages/rubix/core/data.py:491: UserWarning: Explicitly requested dtype <class 'jax.numpy.float64'> requested in array is not available, and will be truncated to dtype float32. To enable more dtypes, set the jax_enable_x64 configuration option or the JAX_ENABLE_X64 shell environment variable. See https://github.com/jax-ml/jax#current-gotchas for more.\n",
      "  rubixdata.galaxy.center = jnp.array(data[\"subhalo_center\"], dtype=jnp.float64)\n",
      "2025-11-10 17:11:58,318 - rubix - INFO - Centering stars particles\n",
      "2025-11-10 17:11:59,305 - rubix - WARNING - The Subset value is set in config. Using only subset of size 2 for stars\n",
      "2025-11-10 17:11:59,305 - rubix - INFO - Data loaded with 2 star particles and 0 gas particles.\n",
      "2025-11-10 17:11:59,306 - rubix - INFO - Data preparation completed in 1.05 seconds.\n",
      "2025-11-10 17:11:59,306 - rubix - INFO - Setting up the pipeline...\n",
      "2025-11-10 17:11:59,307 - rubix - DEBUG - Pipeline Configuration: {'Transformers': {'rotate_galaxy': {'name': 'rotate_galaxy', 'depends_on': None, 'args': [], 'kwargs': {}}, 'spaxel_assignment': {'name': 'spaxel_assignment', 'depends_on': 'rotate_galaxy', 'args': [], 'kwargs': {}}, 'calculate_datacube_particlewise': {'name': 'calculate_datacube_particlewise', 'depends_on': 'spaxel_assignment', 'args': [], 'kwargs': {}}, 'convolve_psf': {'name': 'convolve_psf', 'depends_on': 'calculate_datacube_particlewise', 'args': [], 'kwargs': {}}, 'convolve_lsf': {'name': 'convolve_lsf', 'depends_on': 'convolve_psf', 'args': [], 'kwargs': {}}, 'apply_noise': {'name': 'apply_noise', 'depends_on': 'convolve_lsf', 'args': [], 'kwargs': {}}}}\n",
      "2025-11-10 17:11:59,307 - rubix - DEBUG - Rotation Type found: edge-on\n",
      "2025-11-10 17:11:59,310 - rubix - INFO - Calculating spatial bin edges...\n",
      "/home/annalena/.conda/envs/rubix/lib/python3.12/site-packages/rubix/telescope/factory.py:26: UserWarning: No telescope config provided, using default stored in /home/annalena/.conda/envs/rubix/lib/python3.12/site-packages/rubix/telescope/telescopes.yaml\n",
      "  warnings.warn(\n",
      "2025-11-10 17:11:59,337 - rubix - INFO - Getting cosmology...\n",
      "2025-11-10 17:11:59,547 - rubix - INFO - Calculating spatial bin edges...\n",
      "2025-11-10 17:11:59,556 - rubix - INFO - Getting cosmology...\n",
      "2025-11-10 17:11:59,567 - rubix - INFO - Getting cosmology...\n",
      "2025-11-10 17:11:59,652 - rubix - DEBUG - Method not defined, using default method: cubic\n",
      "/home/annalena/.conda/envs/rubix/lib/python3.12/site-packages/rubix/telescope/factory.py:26: UserWarning: No telescope config provided, using default stored in /home/annalena/.conda/envs/rubix/lib/python3.12/site-packages/rubix/telescope/telescopes.yaml\n",
      "  warnings.warn(\n",
      "2025-11-10 17:11:59,807 - rubix - DEBUG - Method not defined, using default method: cubic\n",
      "/home/annalena/.conda/envs/rubix/lib/python3.12/site-packages/rubix/telescope/factory.py:26: UserWarning: No telescope config provided, using default stored in /home/annalena/.conda/envs/rubix/lib/python3.12/site-packages/rubix/telescope/telescopes.yaml\n",
      "  warnings.warn(\n",
      "2025-11-10 17:12:00,094 - rubix - INFO - Assembling the pipeline...\n",
      "2025-11-10 17:12:00,095 - rubix - INFO - Compiling the expressions...\n",
      "2025-11-10 17:12:00,096 - rubix - INFO - Number of devices: 1\n",
      "2025-11-10 17:12:00,180 - rubix - INFO - Rotating galaxy with alpha=90.0, beta=0.0, gamma=0.0\n",
      "2025-11-10 17:12:00,181 - rubix - INFO - Rotating galaxy for simulation: IllustrisTNG\n",
      "2025-11-10 17:12:00,181 - rubix - WARNING - Gas not found in particle_type, only rotating stellar component.\n",
      "2025-11-10 17:12:00,286 - rubix - INFO - Assigning particles to spaxels...\n",
      "2025-11-10 17:12:00,318 - rubix - INFO - Calculating Data Cube (combined per‐particle)…\n",
      "2025-11-10 17:12:00,481 - rubix - DEBUG - Datacube shape: (1, 1, 466)\n",
      "2025-11-10 17:12:00,481 - rubix - INFO - Convolving with PSF...\n",
      "2025-11-10 17:12:00,486 - rubix - INFO - Convolving with LSF...\n",
      "2025-11-10 17:12:00,493 - rubix - INFO - Applying noise to datacube with signal to noise ratio: 100 and noise distribution: normal\n",
      "2025-11-10 17:12:06,414 - rubix - INFO - Total time for sharded pipeline run: 7.11 seconds.\n"
     ]
    }
   ],
   "source": [
    "# NBVAL_SKIP\n",
    "pipe = RubixPipeline(config)\n",
    "inputdata = pipe.prepare_data()\n",
    "output = pipe.run_sharded(inputdata)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Set target values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# NBVAL_SKIP\n",
    "import jax.numpy as jnp\n",
    "\n",
    "inputdata.stars.age = jnp.array([age_values[index_age], age_values[index_age]])\n",
    "inputdata.stars.metallicity = jnp.array([metallicity_values[index_metallicity], metallicity_values[index_metallicity]])\n",
    "inputdata.stars.mass = jnp.array([[1.0], [1.0]])\n",
    "inputdata.stars.velocity = jnp.array([[0.0, 0.0, 0.0], [0.0, 0.0, 0.0]])\n",
    "inputdata.stars.coords = jnp.array([[0.0, 0.0, 0.0], [0.0, 0.0, 0.0]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2025-11-10 17:12:06,474 - rubix - INFO - Setting up the pipeline...\n",
      "2025-11-10 17:12:06,475 - rubix - DEBUG - Pipeline Configuration: {'Transformers': {'rotate_galaxy': {'name': 'rotate_galaxy', 'depends_on': None, 'args': [], 'kwargs': {}}, 'spaxel_assignment': {'name': 'spaxel_assignment', 'depends_on': 'rotate_galaxy', 'args': [], 'kwargs': {}}, 'calculate_datacube_particlewise': {'name': 'calculate_datacube_particlewise', 'depends_on': 'spaxel_assignment', 'args': [], 'kwargs': {}}, 'convolve_psf': {'name': 'convolve_psf', 'depends_on': 'calculate_datacube_particlewise', 'args': [], 'kwargs': {}}, 'convolve_lsf': {'name': 'convolve_lsf', 'depends_on': 'convolve_psf', 'args': [], 'kwargs': {}}, 'apply_noise': {'name': 'apply_noise', 'depends_on': 'convolve_lsf', 'args': [], 'kwargs': {}}}}\n",
      "2025-11-10 17:12:06,476 - rubix - DEBUG - Rotation Type found: edge-on\n",
      "2025-11-10 17:12:06,479 - rubix - INFO - Calculating spatial bin edges...\n",
      "/home/annalena/.conda/envs/rubix/lib/python3.12/site-packages/rubix/telescope/factory.py:26: UserWarning: No telescope config provided, using default stored in /home/annalena/.conda/envs/rubix/lib/python3.12/site-packages/rubix/telescope/telescopes.yaml\n",
      "  warnings.warn(\n",
      "2025-11-10 17:12:06,491 - rubix - INFO - Getting cosmology...\n",
      "2025-11-10 17:12:06,502 - rubix - INFO - Calculating spatial bin edges...\n",
      "2025-11-10 17:12:06,624 - rubix - INFO - Getting cosmology...\n",
      "2025-11-10 17:12:06,635 - rubix - INFO - Getting cosmology...\n",
      "2025-11-10 17:12:06,681 - rubix - DEBUG - Method not defined, using default method: cubic\n",
      "/home/annalena/.conda/envs/rubix/lib/python3.12/site-packages/rubix/telescope/factory.py:26: UserWarning: No telescope config provided, using default stored in /home/annalena/.conda/envs/rubix/lib/python3.12/site-packages/rubix/telescope/telescopes.yaml\n",
      "  warnings.warn(\n",
      "2025-11-10 17:12:06,757 - rubix - DEBUG - Method not defined, using default method: cubic\n",
      "/home/annalena/.conda/envs/rubix/lib/python3.12/site-packages/rubix/telescope/factory.py:26: UserWarning: No telescope config provided, using default stored in /home/annalena/.conda/envs/rubix/lib/python3.12/site-packages/rubix/telescope/telescopes.yaml\n",
      "  warnings.warn(\n",
      "2025-11-10 17:12:06,804 - rubix - INFO - Assembling the pipeline...\n",
      "2025-11-10 17:12:06,804 - rubix - INFO - Compiling the expressions...\n",
      "2025-11-10 17:12:06,806 - rubix - INFO - Number of devices: 1\n",
      "2025-11-10 17:12:06,888 - rubix - INFO - Rotating galaxy with alpha=90.0, beta=0.0, gamma=0.0\n",
      "2025-11-10 17:12:06,889 - rubix - INFO - Rotating galaxy for simulation: IllustrisTNG\n",
      "2025-11-10 17:12:06,889 - rubix - WARNING - Gas not found in particle_type, only rotating stellar component.\n",
      "2025-11-10 17:12:06,961 - rubix - INFO - Assigning particles to spaxels...\n",
      "2025-11-10 17:12:06,963 - rubix - INFO - Calculating Data Cube (combined per‐particle)…\n",
      "2025-11-10 17:12:06,971 - rubix - DEBUG - Datacube shape: (1, 1, 466)\n",
      "2025-11-10 17:12:06,972 - rubix - INFO - Convolving with PSF...\n",
      "2025-11-10 17:12:06,974 - rubix - INFO - Convolving with LSF...\n",
      "2025-11-10 17:12:06,977 - rubix - INFO - Applying noise to datacube with signal to noise ratio: 100 and noise distribution: normal\n",
      "2025-11-10 17:12:12,833 - rubix - INFO - Total time for sharded pipeline run: 6.36 seconds.\n"
     ]
    }
   ],
   "source": [
    "# NBVAL_SKIP\n",
    "targetdata = pipe.run_sharded(inputdata)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(466,)\n"
     ]
    }
   ],
   "source": [
    "# NBVAL_SKIP\n",
    "print(targetdata[0,0,:].shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Set initial datracube"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "# NBVAL_SKIP\n",
    "inputdata.stars.age = jnp.array([age_values[initial_age_index], age_values[initial_age_index]])\n",
    "inputdata.stars.metallicity = jnp.array([metallicity_values[initial_metallicity_index], metallicity_values[initial_metallicity_index]])\n",
    "inputdata.stars.mass = jnp.array([[1.0], [1.0]])\n",
    "inputdata.stars.velocity = jnp.array([[0.0, 0.0, 0.0], [0.0, 0.0, 0.0]])\n",
    "inputdata.stars.coords = jnp.array([[0.0, 0.0, 0.0], [0.0, 0.0, 0.0]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2025-11-10 17:12:12,913 - rubix - INFO - Setting up the pipeline...\n",
      "2025-11-10 17:12:12,914 - rubix - DEBUG - Pipeline Configuration: {'Transformers': {'rotate_galaxy': {'name': 'rotate_galaxy', 'depends_on': None, 'args': [], 'kwargs': {}}, 'spaxel_assignment': {'name': 'spaxel_assignment', 'depends_on': 'rotate_galaxy', 'args': [], 'kwargs': {}}, 'calculate_datacube_particlewise': {'name': 'calculate_datacube_particlewise', 'depends_on': 'spaxel_assignment', 'args': [], 'kwargs': {}}, 'convolve_psf': {'name': 'convolve_psf', 'depends_on': 'calculate_datacube_particlewise', 'args': [], 'kwargs': {}}, 'convolve_lsf': {'name': 'convolve_lsf', 'depends_on': 'convolve_psf', 'args': [], 'kwargs': {}}, 'apply_noise': {'name': 'apply_noise', 'depends_on': 'convolve_lsf', 'args': [], 'kwargs': {}}}}\n",
      "2025-11-10 17:12:12,914 - rubix - DEBUG - Rotation Type found: edge-on\n",
      "2025-11-10 17:12:12,916 - rubix - INFO - Calculating spatial bin edges...\n",
      "/home/annalena/.conda/envs/rubix/lib/python3.12/site-packages/rubix/telescope/factory.py:26: UserWarning: No telescope config provided, using default stored in /home/annalena/.conda/envs/rubix/lib/python3.12/site-packages/rubix/telescope/telescopes.yaml\n",
      "  warnings.warn(\n",
      "2025-11-10 17:12:12,929 - rubix - INFO - Getting cosmology...\n",
      "2025-11-10 17:12:12,941 - rubix - INFO - Calculating spatial bin edges...\n",
      "2025-11-10 17:12:12,951 - rubix - INFO - Getting cosmology...\n",
      "2025-11-10 17:12:12,961 - rubix - INFO - Getting cosmology...\n",
      "2025-11-10 17:12:13,008 - rubix - DEBUG - Method not defined, using default method: cubic\n",
      "/home/annalena/.conda/envs/rubix/lib/python3.12/site-packages/rubix/telescope/factory.py:26: UserWarning: No telescope config provided, using default stored in /home/annalena/.conda/envs/rubix/lib/python3.12/site-packages/rubix/telescope/telescopes.yaml\n",
      "  warnings.warn(\n",
      "2025-11-10 17:12:13,099 - rubix - DEBUG - Method not defined, using default method: cubic\n",
      "/home/annalena/.conda/envs/rubix/lib/python3.12/site-packages/rubix/telescope/factory.py:26: UserWarning: No telescope config provided, using default stored in /home/annalena/.conda/envs/rubix/lib/python3.12/site-packages/rubix/telescope/telescopes.yaml\n",
      "  warnings.warn(\n",
      "2025-11-10 17:12:13,159 - rubix - INFO - Assembling the pipeline...\n",
      "2025-11-10 17:12:13,160 - rubix - INFO - Compiling the expressions...\n",
      "2025-11-10 17:12:13,164 - rubix - INFO - Number of devices: 1\n",
      "2025-11-10 17:12:13,243 - rubix - INFO - Rotating galaxy with alpha=90.0, beta=0.0, gamma=0.0\n",
      "2025-11-10 17:12:13,244 - rubix - INFO - Rotating galaxy for simulation: IllustrisTNG\n",
      "2025-11-10 17:12:13,244 - rubix - WARNING - Gas not found in particle_type, only rotating stellar component.\n",
      "2025-11-10 17:12:13,315 - rubix - INFO - Assigning particles to spaxels...\n",
      "2025-11-10 17:12:13,317 - rubix - INFO - Calculating Data Cube (combined per‐particle)…\n",
      "2025-11-10 17:12:13,328 - rubix - DEBUG - Datacube shape: (1, 1, 466)\n",
      "2025-11-10 17:12:13,328 - rubix - INFO - Convolving with PSF...\n",
      "2025-11-10 17:12:13,331 - rubix - INFO - Convolving with LSF...\n",
      "2025-11-10 17:12:13,334 - rubix - INFO - Applying noise to datacube with signal to noise ratio: 100 and noise distribution: normal\n",
      "2025-11-10 17:12:19,690 - rubix - INFO - Total time for sharded pipeline run: 6.78 seconds.\n"
     ]
    }
   ],
   "source": [
    "# NBVAL_SKIP\n",
    "initialdata = pipe.run_sharded(inputdata)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Adam optimizer"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/annalena/.conda/envs/rubix/lib/python3.12/site-packages/rubix/telescope/factory.py:26: UserWarning: No telescope config provided, using default stored in /home/annalena/.conda/envs/rubix/lib/python3.12/site-packages/rubix/telescope/telescopes.yaml\n",
      "  warnings.warn(\n",
      "2025-11-10 17:12:19,786 - rubix - INFO - Setting up the pipeline...\n",
      "2025-11-10 17:12:19,787 - rubix - DEBUG - Pipeline Configuration: {'Transformers': {'rotate_galaxy': {'name': 'rotate_galaxy', 'depends_on': None, 'args': [], 'kwargs': {}}, 'spaxel_assignment': {'name': 'spaxel_assignment', 'depends_on': 'rotate_galaxy', 'args': [], 'kwargs': {}}, 'calculate_datacube_particlewise': {'name': 'calculate_datacube_particlewise', 'depends_on': 'spaxel_assignment', 'args': [], 'kwargs': {}}, 'convolve_psf': {'name': 'convolve_psf', 'depends_on': 'calculate_datacube_particlewise', 'args': [], 'kwargs': {}}, 'convolve_lsf': {'name': 'convolve_lsf', 'depends_on': 'convolve_psf', 'args': [], 'kwargs': {}}, 'apply_noise': {'name': 'apply_noise', 'depends_on': 'convolve_lsf', 'args': [], 'kwargs': {}}}}\n",
      "2025-11-10 17:12:19,787 - rubix - DEBUG - Rotation Type found: edge-on\n",
      "2025-11-10 17:12:19,789 - rubix - INFO - Calculating spatial bin edges...\n",
      "/home/annalena/.conda/envs/rubix/lib/python3.12/site-packages/rubix/telescope/factory.py:26: UserWarning: No telescope config provided, using default stored in /home/annalena/.conda/envs/rubix/lib/python3.12/site-packages/rubix/telescope/telescopes.yaml\n",
      "  warnings.warn(\n",
      "2025-11-10 17:12:19,798 - rubix - INFO - Getting cosmology...\n",
      "2025-11-10 17:12:19,808 - rubix - INFO - Calculating spatial bin edges...\n",
      "2025-11-10 17:12:19,818 - rubix - INFO - Getting cosmology...\n",
      "2025-11-10 17:12:19,828 - rubix - INFO - Getting cosmology...\n",
      "2025-11-10 17:12:19,871 - rubix - DEBUG - Method not defined, using default method: cubic\n",
      "/home/annalena/.conda/envs/rubix/lib/python3.12/site-packages/rubix/telescope/factory.py:26: UserWarning: No telescope config provided, using default stored in /home/annalena/.conda/envs/rubix/lib/python3.12/site-packages/rubix/telescope/telescopes.yaml\n",
      "  warnings.warn(\n",
      "2025-11-10 17:12:19,938 - rubix - DEBUG - Method not defined, using default method: cubic\n",
      "/home/annalena/.conda/envs/rubix/lib/python3.12/site-packages/rubix/telescope/factory.py:26: UserWarning: No telescope config provided, using default stored in /home/annalena/.conda/envs/rubix/lib/python3.12/site-packages/rubix/telescope/telescopes.yaml\n",
      "  warnings.warn(\n"
     ]
    }
   ],
   "source": [
    "# NBVAL_SKIP\n",
    "from rubix.pipeline import linear_pipeline as pipeline\n",
    "\n",
    "pipeline_instance = RubixPipeline(config)\n",
    "\n",
    "pipeline_instance._pipeline = pipeline.LinearTransformerPipeline(\n",
    "    pipeline_instance.pipeline_config, \n",
    "    pipeline_instance._get_pipeline_functions()\n",
    ")\n",
    "pipeline_instance._pipeline.assemble()\n",
    "pipeline_instance.func = pipeline_instance._pipeline.compile_expression()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "# NBVAL_SKIP\n",
    "import optax\n",
    "\n",
    "def loss_only_wrt_age_metallicity(age, metallicity, base_data, target):\n",
    "        \n",
    "        base_data.stars.age = age*20\n",
    "        base_data.stars.metallicity = metallicity*0.05\n",
    "\n",
    "        output = pipeline_instance.func(base_data)\n",
    "        #loss = jnp.sum((output.stars.datacube - target) ** 2)\n",
    "        #loss = jnp.sum(optax.l2_loss(output.stars.datacube, target.stars.datacube))\n",
    "        #loss = jnp.sum(optax.huber_loss(output.stars.datacube, target.stars.datacube))\n",
    "        loss = jnp.sum(optax.cosine_distance(output.stars.datacube, target))\n",
    "        \n",
    "        return jnp.log10(loss) #loss#/0.03 #jnp.log10(loss #/5e-5)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "#NBVAL_SKIP\n",
    "import jax\n",
    "\n",
    "def compute_gradient(age, metallicity, base_data, target):\n",
    "    loss, grad_fn = jax.value_and_grad(loss_only_wrt_age_metallicity, argnums=(0,1))\n",
    "    grads = grad_fn(age, metallicity, base_data, target)\n",
    "    return grads, loss"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2025-11-10 17:12:20,133 - rubix - INFO - Rotating galaxy with alpha=90.0, beta=0.0, gamma=0.0\n",
      "2025-11-10 17:12:20,134 - rubix - INFO - Rotating galaxy for simulation: IllustrisTNG\n",
      "2025-11-10 17:12:20,134 - rubix - WARNING - Gas not found in particle_type, only rotating stellar component.\n",
      "2025-11-10 17:12:20,223 - rubix - INFO - Assigning particles to spaxels...\n",
      "2025-11-10 17:12:20,239 - rubix - INFO - Calculating Data Cube (combined per‐particle)…\n",
      "2025-11-10 17:12:20,375 - rubix - DEBUG - Datacube shape: (1, 1, 466)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Initial parameters: {'age': Array([0.7924467, 0.7924467], dtype=float32), 'metallicity': Array([0.5050313, 0.5050313], dtype=float32)}\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2025-11-10 17:12:20,375 - rubix - INFO - Convolving with PSF...\n",
      "2025-11-10 17:12:20,378 - rubix - INFO - Convolving with LSF...\n",
      "2025-11-10 17:12:20,383 - rubix - INFO - Applying noise to datacube with signal to noise ratio: 100 and noise distribution: normal\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "grads: {'age': Array([5.885172, 5.885172], dtype=float32), 'metallicity': Array([0.27147812, 0.27147812], dtype=float32)}\n",
      "loss: -2.057193\n"
     ]
    }
   ],
   "source": [
    "#NBVAL_SKIP\n",
    "#calculate gradient with jax\n",
    "age_init = jnp.array([age_values[initial_age_index]/20, age_values[initial_age_index]/20])\n",
    "metallicity_init = jnp.array([metallicity_values[initial_metallicity_index]/0.05, metallicity_values[initial_metallicity_index]/0.05])\n",
    "\n",
    "\n",
    "# Pack both initial parameters into a dictionary.\n",
    "params_init = {'age': age_init, 'metallicity': metallicity_init}\n",
    "print(f\"Initial parameters: {params_init}\")\n",
    "\n",
    "data = inputdata\n",
    "target_value = targetdata\n",
    "\n",
    "loss, grads = jax.value_and_grad(lambda p: loss_only_wrt_age_metallicity(p['age'], p['metallicity'], data, target_value))(params_init)\n",
    "\n",
    "print(\"grads:\", grads)\n",
    "print(\"loss:\", loss)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "grads_fd: {'age': Array([0.35352707, 0.35352707], dtype=float32), 'metallicity': Array([0.25906563, 0.25891066], dtype=float32)}\n"
     ]
    }
   ],
   "source": [
    "#NBVAL_SKIP\n",
    "#calculate finite differnce\n",
    "import jax\n",
    "import jax.numpy as jnp\n",
    "from jax.flatten_util import ravel_pytree\n",
    "\n",
    "# 1) Skalares Loss über das ganze Param-PyTree\n",
    "f = lambda p: loss_only_wrt_age_metallicity(p['age'], p['metallicity'], data, target_value)\n",
    "\n",
    "# 2) Finite-Difference-Gradient (zentral) für beliebiges PyTree\n",
    "def finite_diff_grad(f, params, eps=1e-5):\n",
    "    flat, unravel = ravel_pytree(params)\n",
    "    def f_flat(x): return f(unravel(x))\n",
    "\n",
    "    def fd_i(i):\n",
    "        e_i = jnp.zeros_like(flat).at[i].set(1.0)\n",
    "        return (f_flat(flat + eps*e_i) - f_flat(flat - eps*e_i)) / (2*eps)\n",
    "\n",
    "    g_flat = jax.vmap(fd_i)(jnp.arange(flat.size))\n",
    "    return unravel(g_flat)\n",
    "\n",
    "# 3) Anwenden: JAX-Grad + FD-Grad berechnen und vergleichen\n",
    "grads_fd = finite_diff_grad(f, params_init, eps=1e-2)\n",
    "print(\"grads_fd:\", grads_fd)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlgAAAHECAYAAADlKlR5AAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjYsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvq6yFwwAAAAlwSFlzAAAPYQAAD2EBqD+naQAAfihJREFUeJzt3Xd4U9X/B/D3TZruRXcLpS0ypJS2rLKnFARBRVzgAgQVBQVEFPmK4E8QXKigKMgUEUFQEdl7D4ECZRZoKaOlFOheaXJ/f4SGhiRN2t40aft+PU+fJveenPtJTpJ+es655wqiKIogIiIiIsnIrB0AERERUU3DBIuIiIhIYkywiIiIiCTGBIuIiIhIYkywiIiIiCTGBIuIiIhIYkywiIiIiCRmZ+0AKkOtVuPGjRtwc3ODIAjWDoeIiIhqOFEUkZ2djaCgIMhkxvupqnWCdePGDQQHB1s7DCIiIqplrl69inr16hndX60TLDc3NwCaJ+nu7m6RYyiVSmzevBm9evWCQqGwyDHINLaD7WBb2A62hW1gO9iOqmiLrKwsBAcHa3MQY6p1glUyLOju7m7RBMvZ2Rnu7u784FgR28F2sC1sB9vCNrAdbEdVtoWpqUmc5E5EREQkMSZYRERERBJjgkVEREQkMSZYRERERBJjgkVEREQkMasnWNevX8eLL74Ib29vODk5oXnz5vjvv/+sHRYRERFRhVl1mYa7d++iY8eO6N69OzZs2ABfX18kJCSgTp061gyLiIiIqFKsmmDNnDkTwcHBWLRokXZbWFiY0fKFhYUoLCzU3s/KygKgWfdCqVRaJMaSei1VP5mH7WA72Ba2g21hG9gOtqMq2sLcugVRFEWLRWFCeHg4evfujWvXrmHXrl2oW7cu3nzzTYwYMcJg+SlTpmDq1Kl625cvXw5nZ2dLh0tERES1XF5eHgYPHozMzMwyFzm3aoLl6OgIABg3bhyeeeYZHDlyBO+88w5+/PFHvPLKK3rlDfVgBQcHIz093aIruW/ZsgWxsbFcodeK2A62g21hO9gWtoHtYDuqoi2ysrLg4+NjMsGy6hChWq1G69atMX36dABAixYtEB8fbzTBcnBwgIODg952hUJh8Td1VRyDTGM72A62he1gW9gGtoPtsGRbmFuvVc8iDAwMRHh4uM62pk2bIjk52UoREREREVWeVROsjh074vz58zrbLly4gJCQECtFRERERFR5Vh0iHDt2LDp06IDp06fj2WefxeHDhzFv3jzMmzevXPXkFRXDrqhYb7tMEOCokOuUM8ZYWaWyGIUqzX2FKBgsm1+kggjDU9kECHCyr1jZAqUK6jKmyDnb21m9rJNCrr2ieGGxCiq1NGUd7eSQyTRli4rVyC/SbwdjZYvVaqP1OtjJIa9AWaVKDaXKeFl7uQx2clm5yxar1Cgqo6xCLoOiAmVVahGFxSqjZe1kMtjblb+sWi0ir4y2eLBsQRn1ymUCHOw073dRFJGvlKasFJ97qcta8jvCWFsAtes7oqzPclV8R6jUxtsBqF3fEWV97qvmO6Lstqjs576s8qVZdZI7AKxbtw4TJ05EQkICwsLCMG7cOKNnET4oKysLHh4eCB6zEjIH/bMIuzfxxaKhMdr7TT/aaLSB2oZ54ffX22vvt/y/LbiTW2SwbGQ9D6wd1Ul7v+OM7biekW+wbCM/V2wZ11V7P/brXUhIyzFYtq6nE/Z90EN7//E5e3HyWqbBsl4u9jj2Uaz2/nM/HcChxDsGyzop5Dj7f49q7w9ddBg7zt8yWBYAkmY8pr395q9Hsf5UqtGyZz7prf2yfXflCaw+ds1o2aP/6wlvV80cuo/+iscvB68YLbtnQncEe2nadPr6s5i3+7LRspvHdkFjfzcAwKwtF/DttgSjZf9+qyOigj0BAD/tuoTPNpwzWva3Ee3Q/iFvAMDSA0mY/Pdpo2UXDmmNHg/7AwBW/XcV7/1x0mjZ7we3xGORgQCAf0+m4K3lx4yW/eLpSDzTOhgAsP3cTQxbbHwR3k+eaIaX24cCAA5cuo1B8w8aLTuxz8N4vetDAIATVzPwxPf7jJZ955FGGBvbGABw4WY2es3abbTsa10a4MO+TQEAV+/kofPnO4yWfaldCP7vyQgAwO2cQrT6dKvRsgNb1sNXz0YB0HzBhU/eZLRs3+YB+OGFVtr7oR/8a7Rsdf+O6D97D05dzzJYlt8R91n6O0KpVGLiwg34I1FutCy/IzQs/R2hVCrx5z/rMeGw8f6jyn5H5OZk4+o3z9r2JHcA6NevH/r162ftMIiIiIgkY/UerMoo6cFKuXXbYBYpzRChEps2bUbv3r20Zw5wiPC+KhsiLCzUawdjZTlEaNnu/+z8AqNtwSFCw2Ut9R2RnVeADRs3GWwLoHZ9R1hziFCpVOKfdevxSC/D7QDUru8Iaw4RKpVK/PvvenSPNd4Wlf3cZ2VlIdDX2/Z7sKTgbG+n84Evq1x56gQApSDCQa65r1AYfnzpLzxTylO2dKNWh7IlHwSpy9rbySCIdibboaSsvZnnbpSnbOkvJinL2pX6IpWyrFwmmP1+L09Z2b2y5rSFrBz1CoJlygIV+9xLXdaS3xHmtEVJ2fLUa4mylvyOsMTnvjxl5TLz2gGoHd8RUpct33eE+W0BlP9zX2zu8zO7ViIiIiIyCxMsIiIiIokxwSIiIiKSGBMsIiIiIokxwSIiIiKSGBMsIiIiIokxwSIiIiKSGBMsIiIiIokxwSIiIiKSGBMsIiIiIokxwSIiIiKSGBMsIiIiIokxwSIiIiKSGBMsIiIiIokxwSIiIiKSGBMsIiIiIokxwSIiIiKSGBMsIiIiIokxwSIiIiKSGBMsIiIiIokxwSIiIiKSGBMsIiIiIokxwSIiIiKSGBMsIiIiIokxwSIiIiKSGBMsIiIiIokxwSIiIiKSGBMsIiIiIokxwSIiIiKSGBMsIiIiIokxwSIiIiKSGBMsIiIiIokxwSIiIiKSGBMsIiIiIokxwSIiIiKSGBMsIiIiIokxwSIiIiKSGBMsIiIiIokxwSIiIiKSGBMsIiIiIolZNcGaMmUKBEHQ+Xn44YetGRIRERFRpdlZO4BmzZph69at2vt2dlYPiYiIiKhSrJ7N2NnZISAgwNphEBEREUnG6glWQkICgoKC4OjoiPbt2+Ozzz5D/fr1DZYtLCxEYWGh9n5WVhYAQKlUQqlUWiS+knotVT+Zh+1gO9gWtoNtYRvYDrajKtrC3LoFURRFi0VhwoYNG5CTk4MmTZogJSUFU6dOxfXr1xEfHw83Nze98lOmTMHUqVP1ti9fvhzOzs5VETIRERHVYnl5eRg8eDAyMzPh7u5utJxVE6wHZWRkICQkBF9//TVeffVVvf2GerCCg4ORnp5e5pOsDKVSiS1btiA2NhYKhcIixyDT2A62g21hO9gWtoHtYDuqoi2ysrLg4+NjMsGy+hBhaZ6enmjcuDEuXrxocL+DgwMcHBz0tisUCou/qaviGGQa28F2sC1sB9vCNrAdbIcl28Lcem1qHaycnBxcunQJgYGB1g6FiIiIqMKsmmCNHz8eu3btQlJSEvbv348BAwZALpdj0KBB1gyLiIiIqFKsOkR47do1DBo0CLdv34avry86deqEgwcPwtfX15phEREREVWKVROsFStWWPPwRERERBZhU3OwiIiIiGoCJlhEREREEmOCRURERCQxJlhEREREEmOCRURERCQxJlhEREREEmOCRURERCQxJlhEREREEmOCRURERCQxJlhEREREEmOCRURERCQxJlhEREREEmOCRURERCQxJlhEREREEmOCRURERCQxJlhEREREEmOCRURERCQxJlhEREREEmOCRURERCQxJlhEREREEmOCRURERCQxJlhEREREEmOCRURERCQxJlhEREREEmOCRURERCQxJlhEREREEmOCRURERCQxJlhEREREEmOCRURERCQxJlhEREREEmOCRURERCQxJlhEREREEmOCRURERCQxJlhEREREEmOCRURERCQxJlhEREREEmOCRURERCQxJlhEREREEmOCRURERCQxJlhEREREEmOCRURERCQxJlhEREREEmOCRURERCQxm0mwZsyYAUEQMGbMGGuHQkRERFQpNpFgHTlyBD/99BMiIyOtHQoRERFRpdlZO4CcnBy88MILmD9/Pj799NMyyxYWFqKwsFB7PysrCwCgVCqhVCotEl9JvZaqn8zDdrAdbAvbwbawDWwH21EVbWFu3YIoiqLFojDDK6+8Ai8vL8yaNQvdunVDdHQ0vvnmG4Nlp0yZgqlTp+ptX758OZydnS0cKREREdV2eXl5GDx4MDIzM+Hu7m60nFV7sFasWIFjx47hyJEjZpWfOHEixo0bp72flZWF4OBg9OrVq8wnWRlKpRJbtmxBbGwsFAqFRY5BprEdbAfbwnawLWwD28F2VEVblIyemWK1BOvq1at45513sGXLFjg6Opr1GAcHBzg4OOhtVygUFn9TV8UxyDS2g+1gW9gOtoVtYDvYDku2hbn1Wi3BOnr0KNLS0tCyZUvtNpVKhd27d2POnDkoLCyEXC63VnhEREREFWa1BOuRRx7BqVOndLYNHToUDz/8MN5//30mV0RERFRtWS3BcnNzQ0REhM42FxcXeHt7620nIiIiqk5sYh0sIiIioprE6utglbZz505rh0BERERUaezBIiIiIpIYEywiIiIiiTHBIiIiIpIYEywiIiIiiTHBIiIiIpIYEywiIiIiiTHBIiIiIpIYEywiIiIiiTHBIiIiIpIYEywiIiIiiTHBIiIiIpIYEywiIiIiiTHBIiIiIpIYEywiIiIiiTHBIiIiIpIYEywiIiIiiTHBIiIiIpIYEywiIiIiiTHBIiIiIpIYEywiIiIiiTHBIiIiIpIYEywiIiIiidlZOwAiIrI9oiiiuLgYKpXK2qHYPKVSCTs7OxQUFPD1sjIp2kIul8POzg6CIFQqFiZYRESko6ioCCkpKcjLy7N2KNWCKIoICAjA1atXK/1HmSpHqrZwdnZGYGAg7O3tK1wHEywiItJSq9VITEyEXC5HUFAQ7O3tmTSYoFarkZOTA1dXV8hknHljTZVtC1EUUVRUhFu3biExMRGNGjWqcJsywSIiIq2ioiKo1WoEBwfD2dnZ2uFUC2q1GkVFRXB0dGSCZWVStIWTkxMUCgWuXLmirasi+E4gIiI9TBSoNpPi/c9PEBEREZHEmGARERERSYwJFhERWYRKLeLApdv4O+46Dly6DZVatHZINUJSUhIEQUBcXFyZ5c6fP4+AgABkZ2dLHsMHH3yA0aNHS15vTcIEi4iIJLcxPgWdZm7HoPkH8c6KOAyafxCdZm7HxvgUa4dWa0ycOBGjR4+Gm5sbAGDnzp0QBEHv53//+5/efplMBg8PD7Ro0QITJkxASopuu40fPx5LlizB5cuXq/x5VRdMsIiISFIb41MwctkxpGQW6GxPzSzAyGXHmGRBc7amJSUnJ2PdunUYMmSI3r7z588jJSVF+/PBBx/o7b9x4waOHDmC999/H1u3bkVERAROnTqlLePj44PevXtj7ty5Fn0e1RkTLCIiKpMoisgrKjbrJ7tAiY/XnoahwcCSbVPWnkF2gdJkXaJYviHFjRs3olOnTvD09IS3tzf69euHS5cu6ZTZv38/oqOj4ejoiNatW+Ovv/7SG26Lj49Hnz594OrqCn9/f7z00ktIT08v89hLlixBSEgInJ2dMWDAAHz99dfw9PTU7p8yZQqio6Px888/IywsTHvqvzkxHz58GC1atNDGfPz4cZOvxcqVKxEVFYW6devq7fPz80NAQID2x9XV1eD+xo0b4/nnn8e+ffvg6+uLkSNH6pTr378/VqxYYTKW2orrYBERUZnylSqET94kSV0igNSsAjSfstlk2TOf9Iazvfl/pnJzczFu3DhERkYiJycHkydPxoABAxAXFweZTIasrCz0798fffv2xfLly3HlyhWMGTNGp46MjAz06NEDw4cPx6xZs5Cfn4/3338fzz77LLZv327wuPv27cO4ceMwY8YMPPHEE9i6dSs++ugjvXIXL17E6tWrsWbNGsjlcrNizsnJQb9+/RAbG4tly5YhMTER77zzjsnXYs+ePWjdurXZr11ZnJyc8MYbb2Ds2LFIS0uDn58fACAmJgbXrl1DUlISQkNDJTlWTcIEi4iIaoSBAwfq3F+4cCF8fX1x5swZREREYPny5RAEAfPnz4ejoyPCw8Nx/fp1jBgxQvuYOXPmoEWLFpg+fbpOPcHBwbhw4QIaN26sd9w5c+agZ8+eePfddyGTydC4cWPs378f69at0ylXVFSEpUuXwtfXt1wxq9VqLFiwAI6OjmjWrBmuXbum15v0oCtXrhhNsOrVq6dX1tvbu8z6Hn74YQCaCfYlCVZQUJD28Uyw9FU4wVKr1bh48SLS0tKgVqt19nXp0qXSgRERkW1wUshx5pPeZpU9nHgHQxYdMVlu8dA2iAnzMnnc8khISMDkyZNx6NAhpKena/82JScnIyIiAufPn0dkZKTOytwxMTE6dZw4cQI7duzQGzYDgEuXLhlMsC5cuIA+ffrobIuJidFLsEJCQnSSK3NiPnv2rF7M7du3N/la5OfnG12BfM+ePdqJ7wBQp04dk/WVDNeWvmySk5MTAPCalUZUKME6ePAgBg8ejCtXruiNkQuCwKuJExHVIIIgmD1U17mRLwI9HJGaWWBwHpYAIMDDEZ0b+UIuk/Yah/3790dISAjmz5+PoKAgqNVqRERElGtCeU5ODvr374+ZM2fq7QsMDKxUfC4uLhaJ2RAfHx/cvXvX4L6wsDCd+WHmOHv2LADo9FTduXMHAPSSRtKoUIL1xhtvoHXr1vj3338RGBjIC4ESEREAQC4T8HH/cIxcdgwCoJNklfyl+Lh/uOTJ1e3bt3H+/HnMnz8fnTt3BgDs3btXp0yTJk2wbNkyFBYWwsHBAQBw5Ihub1vLli2xevVqhIaGws7OvD+RjRs3xrFjx3S2PVhvRWNu2rQpfvnlFxQUFGh7pA4ePGiy7hYtWuDMmTNmxW9Kfn4+5s2bhy5duugkU/Hx8VAoFGjWrJkkx6lpKnQWYUJCAqZPn46mTZvC09MTHh4eOj9ERFR7PRoRiLkvtkSAh+4QVYCHI+a+2BKPRlSuJ8iQOnXqwNvbG/PmzcPFixexfft2jBs3TqfM4MGDoVar8dprr+Hs2bPYtGkTvvzySwD3h77eeust3LlzB4MGDcKRI0dw6dIlbNq0CUOHDjU6OjNq1Chs2bIFs2bNQkJCAn766Sds2LDBZOeDuTELgoARI0bgzJkzWL9+vTbmsvTu3RsHDhyo0IhSWloaUlNTkZCQgBUrVqBjx45IT0/XW5Jhz5496Ny5s3aokHRVKMFq27YtLl68KHUsRERUQzwaEYi97/fAbyPa4dvno/HbiHbY+34PiyRXgObivCtWrMDRo0cRERGBsWPH4osvvtAp4+7ujn/++QdxcXGIjo7GpEmTMHnyZADQ9g4FBQVh3759UKlU6NWrF5o3b44xY8bA09PT6AWAO3bsiK+//hqzZs1CVFQUNm7ciLFjxxqdA1WemF1dXfHPP//g1KlTaNGiBSZNmmRw+PJBffr0gZ2dHbZu3Wqy7IOaNGmCoKAgtGrVCjNmzEDPnj0RHx+P8PBwnXIrVqzQOUGAdFVoiHD06NF49913kZqaiubNm0OhUOjsj4yMlCQ4IiKqvuQyAe0fKvvsNCn17NlTb1jswXnCHTp0wIkTJ7T3f/31VygUCtSvX1+7rVGjRlizZk25jv3KK69g9OjR2iRsxIgRaNiwoXb/lClTMGXKlArF3K5dO73L4phaI8zOzg4ffvghvv76a/TurTlBoVu3bmU+ztT+0jZs2ACZTIann37arPK1UYUSrJLTSocNG6bdJggCRFHkJHciIrJZS5cuRYMGDVC3bl2cOHFCu8ZVZYe5Zs+ejX79+sHNzQ0bNmzAkiVL8MMPP0gUdcW8/vrryMjIQHZ2ts5Zg1LIzc3FokWLzJ6nVhtV6JVJTEyU5OBz587F3LlzkZSUBABo1qwZJk+erHe6KxERkRRSU1MxefJkpKamIjAwEM888wymTZtW6XqPHTuG2bNnIzs7Gw0aNMB3332H4cOHSxBxxdnZ2WHSpEkWqZs9V6ZVKMEKCQmR5OD16tXDjBkz0KhRI4iiiCVLluCJJ57A8ePHeVYCERFJbsKECZgwYYLk9S5atAju7u5G52lR7VPhvr1Lly7hm2++0a6NER4ejnfeeQcPPfSQ2XX0799f5/60adMwd+5cHDx4kAkWERERVVsVSrA2bdqExx9/HNHR0ejYsSMAzbWYmjVrhn/++QexsbHlrlOlUmHVqlXIzc01ukptYWEhCgsLtfezsrIAAEqlEkqlsgLPxLSSei1VP5mH7WA72Ba2wxJtoVQqIYoi1Gq13lU6yLCSieElrxtZj1RtoVarIYoilEql9rqRJcz9vAlieS9XDs0CZr1798aMGTN0tn/wwQfYvHmz3oJrZTl16hTat2+PgoICuLq6Yvny5ejbt6/BslOmTMHUqVP1ti9fvhzOzs7lexJERKTHzs4OAQEBCA4Ohr29vbXDIbKKoqIiXL16FampqSguLtbZl5eXh8GDByMzMxPu7u5G66hQguXo6IhTp06hUaNGOtsvXLiAyMhIFBQUmF1XUVERkpOTkZmZiT/++AM///wzdu3apbfeBmC4Bys4OBjp6ellPsnKUCqV2LJlC2JjY/WWo6Cqw3awHWwL22GJtigoKMDVq1cRGhpqch0n0hBFUXumHq9sYl1StUVBQQGSkpIQHBys9znIysqCj4+PyQSrQkOEvr6+iIuL00uw4uLitFfZNpe9vb12rZBWrVrhyJEj+Pbbb/HTTz/plXVwcNBe3qA0hUIh7Rd9xlUg77bmdnExPPKSoEg/A0XJ6ajO3oBnsHTHI7NJ3tZUYWwL2yFlW6hUKgiCAJlMxgnbZioZiip53ch6pGoLmUwGQRAMfrbM/axVKMEaMWIEXnvtNVy+fBkdOnQAoJmDNXPmTL1l/stLrVbr9FJVuYyrwJxWQLEmBgWAbgBwvlQZOwdg1FEmWURERGRQhdK7jz76CJMnT8bs2bPRtWtXdO3aFXPmzMGUKVPwv//9z+x6Jk6ciN27dyMpKQmnTp3CxIkTsXPnTrzwwgsVCUsaebe1yZVRxYX3e7iIiEhXxlXgRpzxn4yrVgyu6gwZMgRPPvmk9n63bt0wZswY7f3Q0FB88803ZtW1ePFieHp6ShqfNSQlJUEQBL2V6R90/vx5BAQEIDs7W/IYPvjgA4wePVryeh9UoR4sQRAwduxYjB07VvvkK7JKbFpaGl5++WWkpKTAw8MDkZGR2LRpU4XOQiQiIhvwwCiAQdVsFKBbt26Ijo42Oxky15EjR+Di4mJW2eeee07nBLApU6bgr7/+MpmoVFcTJ07E6NGjtbnFzp070b17d71ykyZNwqeffqqzXxAEuLm5oUGDBoiNjcXYsWMRGHj/Gpjjx49HgwYNMHbsWDRo0MBiz6HSa9xXZvn9BQsWVPbwRERkS8ozClBNEixL8fX1Nbusk5NTpS/nI5WioiKLnmGanJyMdevWYfbs2Xr7zp8/rzOx3NXVVWf/2bNntZfui4uLw+eff44FCxZg586daN68OQDAx8cHvXv3xty5c/Uuri0ls4cIW7Zsibt37wLQLNPQsmVLoz9ERFSDiCJQlGveT3G+eXUW55uuq5wnuXfr1g2jR4/GmDFjUKdOHfj7+2P+/PnIzc3F0KFD4ebmhoYNG2LDhg06j4uPj0efPn3g6uoKf39/vPTSS0hPTwegGebbtWsXvv32WwiCAEEQkJSUBJVKhVdffRVhYWFwcXFBmzZt8N1335Ur3geHCDMyMvD666/D398fjo6OiIiIwLp16wDoDhEuXrwYU6dOxYkTJ7QxLV68GMOGDUO/fv10jqFUKuHn51dmh8b8+fMRHBwMZ2dnDBgwAF9//bXOcOSUKVMQHR2Nn3/+GWFhYdqz6jZu3IhOnTrB09MT3t7e6NevHy5duqRT9+HDh9GiRQs4OjqidevWOH78uMnXZeXKlYiKikLdunX19vn5+SEgIED782CC5efnB39/fzRu3BjPP/889u3bB19fX4wcOVKnXP/+/bFixQqTsVSG2T1YTzzxhPYMvieeeIKnohIR1RbKPGB6kLR1LnzUdJkPbwD25g2hlViyZAkmTJiAw4cP4/fff8fIkSPx559/YsCAAfjwww8xa9YsvPTSS0hOToazszMyMjLQo0cPDB8+HLNmzUJ+fr72AtDbt2/Ht99+iwsXLiAiIgKffPIJAE3Pk1qtRr169bBq1SrUqVMH27Ztw9ixYxEUFIRnn3223C+HWq1Gnz59kJ2djWXLluGhhx7CmTNn9Ba5BDTDhfHx8di4cSO2bt0KAPDw8EDjxo3RpUsXpKSkaIfE1q1bh7y8PDz33HMGj7tv3z688cYbmDlzJh5//HFs3boVH330kV65ixcvYvXq1VizZo02ptzcXIwbNw6RkZHIycnB5MmTMWDAAMTFxUEmkyEnJwf9+vVDbGwsli1bhsTERLzzzjsmX4s9e/agdevWZr92ZXFycsIbb7yBsWPHIi0tTbvSQUxMDK5du4akpCSEhoZKcqwHmZ1gffzxx9rbU6ZMsUQsRERElRIVFaU92WrixImYMWMGfHx8MGLECADA5MmTMXfuXJw8eRLt2rXDnDlz0KJFC0yfPl1bx8KFCxEcHIwLFy6gcePGsLe3h7OzMwICArRl5HK5duFrtVqNZ599FidOnMDKlSsrlGBt3boVhw8fxtmzZ9G4cWMAMDo/yMnJCa6urtpFYUt06NABTZo0wS+//KK93uKiRYvwzDPP6PX0lJg9ezb69OmD8ePHAwAaN26M/fv3a3vOShQVFWHp0qU6w5oDBw7UKbNw4UL4+vrizJkziIiIwPLly6FWq7FgwQI4OjqiWbNmuHbtml5v0oOuXLliNMGqV6+eXllvb+8y63v44YcBaCbYlyRYQUFB2sdbPcEqrUGDBjhy5Ijek8rIyEDLli1x+fJlSYIjIiIboHDW9CaZI/Wkeb1TwzYCAZGmj1tOkZH365TL5fD29tbOvQEAf39/AJqTrADgxIkT2LFjh8EE5NKlS9pkx5Dvv/8eCxcuRHJyMvLz81FUVITo6Ohyxwxo1pGsV69emcczx/DhwzFv3jxMmDABN2/exIYNG7B9+3aj5c+fP48BAwbobIuJidFLsEJCQvTmjCUkJGDy5Mk4dOgQ0tPTtWtQJScnIyIiAmfPnkVkZKTOQp3GLoVXWn5+vtFFbvfs2aMz97tOnTom6ytZT730yFvJfLa8vDyTj6+oCiVYJePPDyosLMS1a9cqHZRVOXtrznAxdQaMc9kZMxFRjSEI5g/V2Zk5EdvOqdzDf+Z4cBHIksUiS98H7i9ImZOTg/79+2PmzJl6dZU+8+xBK1aswPjx4/HVV1+hbdu2EAQBP/74Iw4fPlyhuKWawP7yyy/jgw8+wIEDB7B//36EhYWhc+fOla7X0NmO/fv3R0hICObPn4+goCCo1WpERESgqKioUsfy8fHRzvl+UFhYWLmXqzh79iwA6PRU3blzB0D5TjQor3IlWGvXrtXe3rRpEzw8PLT3VSoVtm3bhrCwMOmiswbPYM3pw/fWuSq+egR2G8ZDdA2AMPh3TRmu5E5EVCO0bNkSq1evRmhoKOzsDP9JtLe31+tU2LdvHzp06IA333wTarUaWVlZlRq9iYyMxLVr17TDkqYYigkAvL298eSTT2LRokU4cOAAhg4dWmY9TZo0wZEjR3S2PXjfkNu3b+P8+fOYP3++NoHbu3evTpmmTZvil19+QUFBgbZH6uDBgybrbtGiBc6cOWOynDny8/Mxb948dOnSRSeZio+Ph0KhQLNmzSQ5jiHlSrBKFkwTBAGvvPKKzj6FQoHQ0FB89dVXkgVnNZ7B2gRKdKsHbBgPIScVqBMCOJnujiQiqrWq2SjAW2+9hfnz52PQoEGYMGECvLy8cPHiRaxYsQI///wz5HI5QkNDcejQISQlJcHV1RVeXl5o1KgRli5dik2bNiEkJAQLFizAkSNHKtzJ0LVrV3Tp0gUDBw7E119/jYYNG+LcuXMQBAGPPqo/5BoaGorExETt0KKbm5v2RLThw4ejX79+UKlUen+rHzR69Gh06dIFX3/9Nfr374/t27djw4YNJk9kq1OnDry9vTFv3jwEBgYiOTkZH3zwgU6ZwYMHY9KkSRgxYgQmTpyIpKQkfPnllyZfi969e2P48OFQqVQGJ/mXJS0tDYIg4ObNmzh+/Dg+//xzpKenY82aNTrl9uzZg86dO1t06YtyreSuVquhVqtRv359pKWlae+XXN7m/PnzeqeIVnuOHsi1v3d9xZQT1o2FiMjWlYwCvLbL+I8NLTIaFBSEffv2QaVSoVevXmjevDnGjBkDT09P7bXsxo8fD7lcjvDwcPj6+iI5ORmvv/46nnrqKTz33HNo37497ty5Y3LytimrV69GmzZtMGjQIISHh2PChAkGe6kAzQTzRx99FN27d4evry9+++037b6ePXsiMDAQvXv31k7mNqZjx4748ccf8fXXXyMqKgobN27E2LFjTV7oWyaTYcWKFTh69CgiIiIwduxYvTWlXF1d8c8//+DUqVNo0aIFJk2aZHAo9kF9+vSBnZ2d9gzJ8mjatCmaNm2KNm3aYMaMGejZsyfi4+MRHh6uU27FihXaEx8sRRDFci40YkOysrLg4eFh8orWlaFUKpH2fR/UzTgCxH4CdDR9iilJT6lUYv369ejbty8vMGxlbAvbYYm2KCgoQGJios56R1S2kiFCd3d3m7jYc05ODurWrYtFixbhqaeeKvfjR4wYgXPnzmHPnj0WiM4833//PdauXYtNmzaV63HmtMWGDRvw7rvv4uTJk0aHhcv6HJibe1R4Jffc3Fzs2rULycnJehPa3n777YpWa5MynUI1CRZ7sIiIyEap1Wqkp6fjq6++gqenJx5//HGzHvfll18iNjYWLi4u2LBhA5YsWYIffvjBwtGW7fXXX0dGRgays7MrdcUYQ3Jzc7Fo0SKjyZVUKlT78ePH0bdvX+Tl5SE3NxdeXl5IT0+Hs7Mz/Pz8alyCleEcqrnBBIuIiGxUcnIywsLCUK9ePSxevNjsBOLw4cP4/PPPkZ2djQYNGuC7777D8OHDLRxt2ezs7DBp0iSL1P30009bpN4HVSjBGjt2LPr3748ff/wRHh4eOHjwIBQKBV588UWzVmmtbjKdQjQ3bl8ECrIAR8sMRxIREVVUaGgoKjLrZ+XKlRaIhio0WBwXF4d3330XMpkMcrkchYWFCA4Oxueff44PP/xQ6hitrkjhDtH93jWRUk9ZNxgiIiKyeRVKsBQKhXbymJ+fH5KTkwForoV09epV6aKzIWJAlOYGhwmJiIjIhAoNEbZo0QJHjhxBo0aN0LVrV0yePBnp6en45ZdfEBERIXWMNkEMiAQurAdS4qwdChEREdm4CvVgTZ8+XXsJgWnTpqFOnToYOXIkbt26hXnz5kkaoK0QS66ZxR4sIiIiMqFCPVilr3Lt5+eHjRs3ShaQrRID7w0Rpl8AinItcg0tIiIiqhkq1IP16aefIjExUepYbJurP+AaAIhq4OZpa0dDRERENqxCCdaqVavQsGFDdOjQAT/88APS09Oljss2lfRi3YizahhERFR7JSUlQRAExMXFlVnu/PnzCAgIQHZ2tsVj+uCDDzB69GiLH6c6qVCCdeLECZw8eRLdunXDl19+iaCgIDz22GNYvnw58vLypI7RdgTyTEIiIqoeJk6ciNGjR2tXQt+5cycEQUBGRoZOuWvXrsHe3t7gSWonTpyAvb091q5dq7N99erVcHR0RHx8PADN9RqXLFmCy5cvW+bJVEMVvmhSs2bNMH36dFy+fBk7duxAaGgoxowZg4CAACnjsy1B0ZrfTLCIiKgSHrzEnNSSk5Oxbt06DBkyxGTZxYsX49lnn0VWVhYOHTqksy8qKgqTJ0/Ga6+9htu3bwMA0tLS8MYbb2Dq1KnapMzHxwe9e/fG3LlzJX8u1ZUkV6V0cXGBk5MT7O3toVQqpajSNpX0YN06CygLrBsLEVEVyysqNvpToFRJWrYiNm7ciE6dOsHT0xPe3t7o168fLl26pFNm//79iI6OhqOjI1q3bo2//vpLb7gtPj4effr0gaurK/z9/fHSSy+ZnAqzZMkShISEwNnZGQMGDMDXX38NT09P7f4pU6YgOjoaP//8s84FhM2J+fDhw2jRooU25uPHj5t8LVauXImoqCjUrVu3zHKiKGLRokV46aWXMHjwYCxYsECvzMSJE1G/fn289dZbADTXCWzUqBHGjx+vU65///5YsWKFydhqiwpf6TAxMRHLly/H8uXLcf78eXTt2hVTp06tsmv8WIV7XcDZG8i7DaSdBuq2snZERERVJnzyJqP7ujfxxaKhMdr7rf5vK/IfSKRKtA3zwu+vt9fe7zRzB+7k6vboJM14rNzx5ebmYty4cYiMjEROTg4mT56MAQMGIC4uDjKZDFlZWejfvz/69u2L5cuX48qVKxgzZoxOHRkZGejRoweGDx+OWbNmIT8/H++//z6effZZbN++3eBx9+3bh3HjxmHGjBl44oknsHXrVnz00Ud65S5evIjVq1djzZo1kMvlZsWck5ODfv36ITY2FsuWLUNiYqJZl6Tbs2ePzhn/xuzYsQN5eXno2bMn6tatiw4dOmDWrFlwcbl/prxcLseSJUvQsmVLDB48GJs2bUJcXJz2OZSIiYnBtWvXkJSUhNDQUJPHrukqlGC1a9cOR44cQWRkJIYOHYpBgwaZzJJrBEHQ9GJd2q4ZJmSCRURkMwYOHKhzf+HChfD19cWZM2cQERGB5cuXQxAEzJ8/H46OjggPD8f169cxYsQI7WPmzJmDFi1aYPr06Tr1BAcH48KFC2jcuLHecefMmYOePXtqLyHXuHFj7N+/H+vWrdMpV1RUhKVLl8LX17dcMavVaixYsACOjo5o1qwZrl27hpEjR5b5Wly5csWsBGvBggV4/vnnIZfLERERgQYNGmDVqlV6Q4tNmzbFmDFjMGPGDMycOdPg6xAUFKQ9NhOsCiZYjzzyCBYuXIjw8HCp47F9gdH3EywiolrkzCe9je6TCYLO/aMf9TS77N73u1cusHsSEhIwefJkHDp0COnp6VCr1QA085EiIiJw/vx5REZGaofnAE2vS2knTpzAjh074Orqqlf/pUuXDCYWFy5cQJ8+fXS2xcTE6CVYISEhOsmVOTGfPXtWL+b27dvDlPz8fJ3HGJKRkYE1a9Zg79692m0vvvgiFixYoJdg5eTk4Pfff4ezszP27NmDCRMm6NXn5OQEADX7ZLdyqFCCNW3aNKnjqD54JiER1VLO9ub/ybBU2bL0798fISEhmD9/PoKCgqBWqxEREVGuCeU5OTno378/Zs6cqbev5AomFVV62E3KmA3x8fHB3bt3yyyzfPlyFBQUoG3bttptoihCrVbr9da99957cHR0xP79+9GuXTssXboUL7/8sk59d+7cAQC9JLK2MvtdPW7cOPzf//0fXFxcMG7cuDLLfv3115UOzGaVJFg3TwPFRYCdvXXjISIi3L59G+fPn8f8+fPRuXNnANDpmQGAJk2aYNmyZSgsLISDgwMA4MiRIzplWrZsidWrVyM0NBR2dub9iWzcuDGOHTums+3Beisac9OmTfHLL7+goKBA2yN18OBBk3W3aNECZ86cKbPMggUL8O677+r1Vr355ptYuHAhZsyYAQDYsmULfv75Z+zfvx9RUVH49NNPMWbMGMTGxuoknfHx8VAoFGjWrJnJ+GoDs88iPH78uPYMwWPHjuH48eMGf0wtfFbt1QkFHD0AVRFw65y1oyEiIgB16tSBt7c35s2bh4sXL2L79u16nQGDBw+GWq3Ga6+9hrNnz2LTpk348ssvAQDCvWHLt956C3fu3MGgQYNw5MgRXLp0CZs2bcLQoUOhUhmetD9q1Chs2bIFs2bNQkJCAn766Sds2LBBW2dlYxYEASNGjMCZM2ewfv16bcxl6d27Nw4cOGA05ri4OBw7dgzDhw9HRESEzs+gQYOwZMkSFBcXIysrC6+++iree+89tGnTBgAwduxYhIeH47XXXtOpc8+ePejcubN2qLDWE6uxzMxMEYCYmZlpsWMUFRWJf/31l1hUVHR/4+J+ovixuygeXWqx45Iug+1AVsG2sB2WaIv8/HzxzJkzYn5+vmR1VpUtW7aITZs2FR0cHMTIyEhx586dIgDxzz//1JbZt2+fGBkZKdrb24utWrUSly9fLgIQz507py1z4cIFccCAAaKnp6fo5OQkPvzww+KYMWNEtVpt8LgqlUr85ptvxLp164pOTk7ik08+KX766adiQECAtszHH38sRkVFVSjmAwcOiFFRUaK9vb0YHR0trl69WgQgHj9+3OhroVQqxaCgIHHjxo3abdu2bRMBiNnZ2eKoUaPE8PBwg49NSUkRZTKZ+Pfff4tDhw4VIyIixMLCQp0yFy5cEJ2dncUlS5ZotzVp0kT87bffjMZUFVQqlXj37l1RpVJVqp6yPgfm5h7lTrCKiopEuVwunjp1qrwPlZzVEqxNkzQJ1rp3LXZc0sU/6raDbWE7mGBV3rJly0SFQiHm5eVVuA5Df9SHDx8udurUSYoQK2zOnDlir169tPd/++030dXV1SLHWr9+vdi0aVNRqVRapH5z2VKCVe6ZhQqFAvXr1zfa7VgrBEZrfqfEWTMKIiIqp6VLl6JBgwaoW7cuTpw4oV3jqrLDWrNnz0a/fv3g5uaGDRs2YMmSJfjhhx8kirpiXn/9dWRkZCAzMxPXr1/HnDlz8Mgjj1jkWLm5uVi0aJHZ89Zqgwqt5D5p0iR8+OGH2jMGap2Sie6p8YCqYisOExFR1UtNTcWLL76Ipk2bYuzYsXjmmWcwb968Std77Ngx9O7dG82bN8ePP/6I7777DsOHD5cg4oqzs7PDpEmTsGPHDrRt2xYuLi747rvvLHKsp59+WudsRKrgMg1z5szBxYsXERQUhJCQEL1TTx88m6LG8XoIsHcFinKA2wmAX1NrR0RERGaYMGGCwTWcKmvRokVwd3eHTCbJFegk9eSTTyI7O9vaYdQ6FUqwnnzySYnDqGZkMiAgEkjeD9yIY4JFREREOiqUYH388cdSx1H9BEZpEqyUE0D0IGtHQ0RERDakwn2ZGRkZ+PnnnzFx4kTtXKxjx47h+vXrkgVn04KiNb+5ojsRERE9oEI9WCdPnkTPnj3h4eGBpKQkjBgxAl5eXlizZg2Sk5OxdOlSqeO0PdqJ7icBtVozbEhERESECvZgjRs3DkOGDEFCQoLOxST79u2L3bt3SxacTfNuBNg5aSa637lk7WiIiIjIhlQowTpy5Ahef/11ve1169ZFampqpYOqFuR2QECE5jaHCYmIiKiUCiVYDg4OyMrK0tt+4cKF2nUVbS44SkRERgwZMkTnrPtu3bphzJgxZj12586dEAQBGRkZFomtKgmCgL/++qvMMrdv34afnx+SkpIsHs+PP/6I/v37W/w4FUqwHn/8cXzyySfaiz8LgoDk5GS8//77GDhwoKQB2rSSeVjswSIiqrHKkxiVZc2aNfi///s/s8p26NABKSkp8PDwAAAsXrwYnp6elY7BVk2bNg1PPPEEQkNDAQBJSUkQBAFxcXE65fLz8+Hl5QUfHx8UFhbq7Ltx4wa8vb3x008/6Ww/dOgQFAoFNm/eDAAYNmwYjh07hj179ljs+QAVTLC++uor5OTkwNfXF/n5+ejatSsaNmwINzc3TJs2TeoYbVfpBEsUrRsLERHZNC8vL7i5uZlV1t7eHgEBARAEwcJRmaZSqaBWqy1Wf15eHhYsWIBXX33VZNnVq1ejWbNmePjhh/V6xYKCgvDtt9/ik08+QUJCAgBNQvbKK69g+PDh6NWrFwDNazt48GCLrWpfokIJloeHB7Zs2YJ///0X3333HUaNGoX169dj165dequ612i+DwNye6AgE7ibZO1oiIgsqyjX+I+yoBxl802XrYBu3bph9OjRGDNmDOrUqQN/f3/Mnz8fubm5GDp0KNzc3NCwYUNs2LBB53Hx8fHo06cPXF1d4e/vj5deegnp6ekANMN8u3btwrfffgtBECAIApKSkqBSqfDqq68iLCwMLi4uaNOmjck/2A/2hBUWFuL9999HcHAwHBwc0LBhQyxYsACA7hDhzp07MXToUGRmZmpjmDJlCj755BNEREToHSc6OhofffSR0TjWrl2LRo0awdHREd27d8eSJUt0hiNLesvWrl2L8PBwODg4IDk5GUeOHEFsbCx8fHzg4eGBrl276l25JSEhAV26dIGjoyPCw8OxZcuWMl8TAFi/fj0cHBzQrl07k2UXLFiAF198ES+++KL2tSrtxRdfRI8ePTBs2DCo1WpMnDgRSqUSX3zxhU65/v37Y+3atcjPz9erQyrlXqZBrVZj8eLFWLNmjbYLLywsDAEBARBFsVzZ9meffYY1a9bg3LlzcHJyQocOHTBz5kw0adKkvGFZh5094N8MuHFc04vlFWbtiIiILGd6kPF9jXoBL6y6f/+LhoAyz3DZkE7A0H/v3/+mOZB3W7fMlMwKhbhkyRJMmDABhw8fxu+//46RI0fizz//xIABA/Dhhx9i1qxZeOmll5CcnAxnZ2dkZGSgR48eGD58OGbNmoX8/HztBaC3b9+Ob7/9FhcuXEBERAQ++eQTAICvry/UajXq1auHVatWoU6dOti2bRvGjh2LoKAgPPvss2bF+vLLL+PAgQP47rvvEBUVhcTERG1iV1qHDh3wzTffYPLkyTh//jwAwNXVFRkZGZg6dSqOHDmCNm3aAACOHz+OkydPYs2aNQaPmZiYiKeffhrvvPMOhg8fjuPHj2P8+PF65fLy8jBz5kz8/PPP8Pb2hp+fHy5fvoxXXnkFs2fPhiiK+Oqrr9C3b18kJCTAzc0NarUaTz31FPz9/XHo0CFkZmaaNbS6Z88etGrVymS5S5cu4cCBA1izZg1EUcTYsWNx5coVhISE6JT7+uuv0bFjR7zwwgtYtWoVtm/fDldXV50yrVu3RnFxMQ4dOoRu3bqZPHZFlCvBEkURjz/+ONavX4+oqCg0b94coiji7NmzGDJkCNasWWNyIltpu3btwltvvYU2bdqguLgYH374IXr16oUzZ85Un56wwKj7CVazJ60dDRFRrRYVFYX//e9/AICJEydixowZ8PHxwYgRIwAAkydPxty5c3Hy5Em0a9cOc+bMQYsWLTB9+nRtHQsXLkRwcDAuXLiAxo0bw97eHs7OzggICNCWkcvlmDp1KgBNx8Ozzz6LEydOYOXKlWYlWBcuXMDKlSuxZcsW9OzZEwDQoEEDg2Xt7e3h4eEBQRB0YnB1dUXv3r2xaNEibYK1aNEidO3a1WhdP/30E5o0aaLt0WnSpAni4+P1pvcolUr88MMPiIqK0m7r0aOHTpl58+bB09MTu3btQr9+/bB161acO3cOmzZtQlCQJhmfPn06+vTpU+ZrceXKFW35sixcuBB9+vRBnTp1AED73KdMmaJTztfXF1OnTsWbb76JkSNHokuXLnp1OTs7w8PDA1euXDF53IoqV4K1ePFi7N69G9u2bUP37t119m3fvh1PPvkkli5dipdfftms+jZu3KhXv5+fH44ePWrwBSksLNSZ1FZyJqNSqdROuJdaSb3G6pf5RUAOQH3jOFQWioFMtwNVHbaF7bBEWyiVSoiiCLVarT/v5oNrxh8ok2sWXS7x7gXjZQWZbtm3DZwoVME5P82bN9fGLQgCvL29ERERod1WcqZ7amoq1Go14uLisGPHDr0eDkAz3NWwYUMA0L4mpf3www9YtGgRkpOTkZ+fj6KiIkRHR2vLiaKo97iS+8eOHYNcLkfnzp0Nzm8q2VbSDqXvl/bqq69i+PDh+PLLLyGTybB8+XJ89dVXRudMnTt3Dq1bt9bZ37p1a71j2dvb67xuAHDz5k189NFH2LVrF9LS0qBSqZCXl4crV65ArVbjzJkzCA4ORkBAgPZxbdu21anbkLy8PAQFBensf/D5q1QqLFmyBLNmzdLuGzx4MCZMmID//e9/2otsi6KoLevs7IyDBw+iqKgIdnb66Y6TkxNycnKMvv6iKEKpVEIul+vsM/fzVq4E67fffsOHH36ol1wBmsz2gw8+wK+//mp2gvWgzExNl7CXl5fB/Z999pn2P4bSNm/eDGdn5wod01zGxpE9c3PQFYAy+T9s/PdfwAYmJNZk5oznU9VgW9gOKdvCzs4OAQEByMnJQVFRUTkeqQJgbnkVkG/ij1SB/lJAphQXF0MURZ1lhEr+4D64tFBubi6ysrKQkZGBRx99VK8XBAD8/f2RlZWF4uJiFBUV6dSxevVqvPfee/i///s/xMTEwNXVFd999x2OHj2q889/cXGx9n7pesR7J0ZlZWVBoVDoHTsvTzO8mp2dDZlMhoKCAr3nBgBdu3aFvb09li9fDnt7exQVFaFXr14Gl1IqiUGpVOrsN3QsR0dHZGdn6zz2pZdewp07dzBt2jTtvLGSY2VlZaGgoABqtVqn7pLb+fn5RmPy8PBAWlqazv6cnBwA99tp8+bNuH79OgYN0r32r0qlwj///KOTl8yePRuXL1/G9u3b0a9fP0yZMgUTJkzQO+6dO3fg6upqMK6ioiLk5+dj9+7dKC4u1tlX8nqZUq4E6+TJk/j888+N7u/Tp0+FZ+Wr1WqMGTMGHTt2NDhpD9B0944bN057PysrC8HBwejVqxfc3d0rdFxTlEoltmzZgtjYWIMfAhQXQPziUzgUZ6Nv52jAva5F4qjtTLYDVRm2he2wRFsUFBTg6tWrcHV11blSR3VgZ2cHe3t7nb8HMpkMjo6Oen8jnJyc4O7ujpiYGKxZswYREREGezlKysrlcp06jh8/jg4dOmDcuHEQRRHZ2dm4evWqTjmFQgE7Ozvt/dLxtW3bFmq1GsePH9cOEZZW0mng5uYGd3d3uLu7Q61WG/xb98orr+D333+Hvb09nn/+efj7+xt9jZo1a4YNGzbo1HPmzBmdYzk6OkIQBL1jHTp0CHPmzMHTTz8NALh69Spu376tfX2jo6Nx/fp15ObmIjAwEABw4MABndfbkJiYGPz66686+0t6FF1cXODu7o4VK1bgueeew4cffqjz2OnTp2PFihV44oknAGhOWJgxYwaWLVuGNm3a4IcffsALL7yAZ599FpGRkdrHXbp0CQUFBejQoYPBuAoKCuDk5KSdsF+asUTxQeVKsO7cuVNmw/n7++Pu3bvlqVLrrbfeQnx8PPbu3Wu0jIODAxwcHPS2KxQKi3/RGz2GQgH4NgVunoLi1mnAO9SicdR2VdHWZB62he2Qsi1UKhUEQYBMJtMOu1QnJbGb2lby/EaNGoWff/4ZL7zwAiZMmAAvLy9cvHgRK1aswM8//wy5XI6wsDAcPnwYycnJcHV1hZeXFxo3boxffvkFW7ZsQUhICBYsWID//vsPYWFh2mOVnPFX+tgl9xs0aKBdPqBkkvuVK1eQlpaGZ599VvuYkjgbNGiAnJwc7NixA1FRUXB2dtYmYSNGjEDTpk0BAPv27Suz3d544w3MmjULEydOxKuvvoq4uDgsWbIEgGZeWel2f7CeRo0a4ddff0VMTAyysrLw3nvvwcnJSfucevXqhcaNG2Po0KH44osvkJWVpT2bsaz306OPPooPP/wQmZmZ2vlVpWO4ffs21q1bh7Vr1+okSYAmuRwwYAAyMjLg7u6OYcOGoV+/fnjqqacgk8nwzDPPYM2aNRg2bBgOHz6sTaL37duHBg0aoFGjRgZjkslkEATB4GfL3M9auT49KpXKaIYPaBrnwa40c4waNQrr1q3Djh07UK9evXI/3upK1sO6EWfVMIiIqHyCgoKwb98+qFQq9OrVC82bN8eYMWPg6emp/SM/fvx4yOVyhIeHw9fXF8nJyXj99dfx1FNP4bnnnkP79u1x584djBw5slzHnjt3Lp5++mm8+eabePjhhzFixAjk5hpeoqJDhw5444038Nxzz8HX11dnNKlRo0bo0KEDHn74Ye2cJ2PCwsLwxx9/YM2aNYiMjMTcuXMxadIkADDYgVHaggULcPfuXbRs2RIvvfQS3n77bfj5+Wn3y2Qy/Pnnn8jPz0dMTAyGDx9u1tqYzZs3R8uWLbFy5UrttpJ5UXZ2dli6dClcXFzwyCOP6D32kUcegZOTE5YtW4bp06fj+vXreksyfP/990hJSdE5keG3337TnvhgKYIomr9CpkwmQ58+fYw2QmFhITZu3AiVSmVWfaIoYvTo0fjzzz+xc+dOo5mkMVlZWfDw8EBmZqZFhwjXr1+Pvn37Gs9aD80DNrwHNOoNvLDScBmqFLPagaoE28J2WKItCgoKkJiYiLCwsGo3RGgtJfOO3N3drdLrJ4oiGjVqhDfffFNnGo25pk2bhh9//BFXr161QHTm+ffff/Hee+8hPj4eMpkMBw8eRPv27XHr1i34+PiYXY85bXH69Gn06NEDFy5c0K6U/6CyPgfm5h7lGiJ85ZVXTJYpzwT3t956C8uXL8fff/8NNzc37YWiPTw84OTkVJ7QrCsoWvObl8whIqIqdOvWLaxYsQKpqakYOnSoWY/54Ycf0KZNG3h7e2Pfvn344osvMGrUKAtHWrbHHnsMCQkJuH79OgoLC/HFF18gKiqqXMmVuVJSUrB06VKjyZVUypVgLVq0SNKDz507FwD0FvlatGgRhgwZIumxLMq/mea045xUIDsVcAsw/RgiIqJK8vPzg4+PD+bNm6edv2RKQkICPv30U9y5cwf169fHu+++i4kTJ1o4UtPGjBmDuLg4dOjQAdHR0Vi6dKlFjmPopAJLKPdK7lIqx+ikbbN3AXwaA7fOaXqxmGAREVEVqMjf0VmzZmHWrFkWiKbyoqOjzV4GwdZVv1NEbFVgtOY3hwmJiIhqPSZYUik5k5AJFhHVADVmhIGoAqR4/zPBkgqXaiCiGqDkbMSaMkxDVBEl7//KnJ1r1TlYNUpAc83vrGtAbjrgIv2ZD0REliaXy+Hp6Ym0tDQAmhXFBV4CrExqtRpFRUUoKCiolouz1iSVbQtRFJGXl4e0tDR4enrqXYewPJhgScXRHfBuCNy+qBkmbKi/IBoRUXUQEKA5UackyaKyiaKI/Px87armZD1StYWnp6f2c1BRTLCkFBjFBIuIqj1BEBAYGAg/Pz8olSYuykxQKpXYvXs3unTpwsV3rUyKtlAoFJXquSrBBEtKgVFA/GogJc7akRARVZpcLpfkD01NV3KZOEdHRyZYVmZLbcHBYinxTEIiIiICEyxplSRYd5OA/LtWDYWIiIishwmWlJzqAJ4hmtupp6wbCxEREVkNEyypcT0sIiKiWo8JltQ4D4uIiKjWY4IltaBozW8mWERERLUWEyypBdzrwbp9ESjMtm4sREREZBVMsKTm6gu41wUgcqI7ERFRLcUEyxI4D4uIiKhWY4JlCYHRmt9MsIiIiGolJliWwKUaiIiIajUmWJZQkmClnweK8qwbCxEREVU5JliW4B4IuPoDohq4edra0RAREVEVY4JlKdqJ7nFWDYOIiIiqHhMsS2GCRUREVGsxwbIULtVARERUazHBspSSpRrSzgLFhVYNhYiIiKoWEyxL8agHOHkB6mIg7Yy1oyEiIqIqxATLUgSB62ERERHVUkywLInzsIiIiGolJliWFBSt+c0Ei4iIqFZhgmVJJT1YN08DKqV1YyEiIqIqwwTLkuqEAQ4egKoQuHXO2tEQERFRFWGCZUmCAARGam5zmJCIiKjWYIJlaZzoTkREVOswwbK0kgVHuVQDERFRrcEEy9JKerBSTwFqlXVjISIioirBBMvSvB8CFC5AcT6QnmDtaIiIiKgKMMGyNJmcE92JiIhqGSZYVUE70T3OqmEQERFR1WCCVRV4JiEREVGtwgSrKpScSZhyElCrrRoKERERWR4TrKrg0xiwcwSKsoG7idaOhoiIiCyMCVZVkNsB/hGa2zeOWzcWIiIisjgmWFWF87CIiIhqDasmWLt370b//v0RFBQEQRDw119/WTMcywqK1vxmgkVERFTjWTXBys3NRVRUFL7//ntrhlE1SvdgiaJ1YyEiIiKLsrPmwfv06YM+ffpYM4Sq49sUkCmAggwg4wpQJ9TaEREREZGFWDXBKq/CwkIUFhZq72dlZQEAlEollEqlRY5ZUm/l6xdg59cUQupJFF87BtG1buWDq0WkaweqLLaF7WBb2Aa2g+2oirYwt25BFG1jvEoQBPz555948sknjZaZMmUKpk6dqrd9+fLlcHZ2tmB00ohKXojQ2ztxwb8/zgY9Y+1wiIiIqJzy8vIwePBgZGZmwt3d3Wi5apVgGerBCg4ORnp6eplPsjKUSiW2bNmC2NhYKBSKStUlO7oI8o3vQd3gEagG/S5RhLWDlO1AlcO2sB1sC9vAdrAdVdEWWVlZ8PHxMZlgVashQgcHBzg4OOhtVygUFn9TS3KMeq0AALLUE5DZ2QGCIEFktUtVtDWZh21hO9gWtoHtYDss2Rbm1st1sKqSfzggyIG8dCDrhrWjISIiIguxag9WTk4OLl68qL2fmJiIuLg4eHl5oX79+laMzEIUToBfU+BmvGa5Bg9OdCciIqqJrNqD9d9//6FFixZo0aIFAGDcuHFo0aIFJk+ebM2wLEu7HlacdpNKLeLApdv4O+46Dly6DZXaJqbFERERUQVZtQerW7dusJE59lUnMAqI+1W7ovvG+BRM/ecMUjIL7hfxcMTH/cPxaESgtaIkIiKiSuAcrKpWakX3jfEpGLnsmE5yBQCpmQUYuewYNsanWCFAIiIiqiwmWFUtoDkAAchOwey1+2Co/65k29R/znC4kIiIqBpiglXV7F0An8YAAN+cc0aLiQBSMgtwOPFOFQVGREREUqlW62DVGIFRQPp5RAhJ2IkWZRZNyy4oc79NybgK5N2GShRx+noW7uQVwcvZHs3qukMuCICzN+AZbO0oiYiILI4JljUERQOnViJClgSoyi7q5+ZYFRFVXsZVYE4roLgQcgCRhsrYOQCjjjLJIiKiGo9DhNZwb6J7pDzJaBEBmrMJY8K8qiamysq7DRQXll2muFBTjoiIqIZjgmUNAc0BAEG4BU9kGywiAvi4fzjksupxOR2VmcttmFuOiIioOmOCZQ2OHijyCAUARMiSUMdZ/7pGLet7Vqt1sE5fz5K0HBERUXXGOVhWIIoijhWFoB2S8IT/LSwZHYvDiXeQll2A/CIVJq45hWPJGdibkI5OjXysHa5Z7uQVSVqOiIioOmMPlhVsjE/FjqwgAEBf75uQywS0f8gbT0TXxfMx9fFKh1AAwMdr41FUrLZipObzcrY3q9z6UynYfu4milXV43kRERFVBBOsKpZfpML/rTuDeDEUAOByO16vzNjYxvB2scelW7lYvD+xiiOsmGZ13c0qd/pGFoYt/g8dZ27Hl5vOI/l2noUjIyIiqnpMsKrY9zsu4kZmAe64PazZcDcRyM/QKePhpMD7fTT7v92agJtZtr8WltrMDqmBzb1Qx1mBm1mFmLPjIrp8sQOD5h3E33HXUaA0sWYFERFRNcEEqwolpudi3u7LAIB3+rcDPOtrdqSe0iv7dMt6aFHfE7lFKkz792xVhlkhi09ko0DUn6z/oGFZc3FwbCt8P7glujT2hSAABy7fxjsr4hAzbSsm/x2P0zcyqyBiIiIiy2GCVUVEUcTUf06jSKVGl8a+6N3MX+fCzw+SyQT83xMREARg7YkbOHDJdtePOnUtEzP252CycggAQHTyxsk+f2Fntz9wsu9aqEbsBJ6cCzi4A6kn4fBLfzwWJmDpsBjsmdAdY3o2Ql1PJ2QVFGPpgSt47Lu96Dd7D345kITMfKXR46rUIg5cuo2/467jwKXbvG4jERHZDJ5FWEW2nk3DzvO3oJALmNI/HIIgaBKss/8AKXEGHxNR1wMvtK2PZQeTMWXtaax7uxMUctvKiQuUKoxbGQeVWsRzPglADiBEPY/Itt11C9ZtAQRGA78MANJOAwsfBV7+C/XqhGJMz8YY3aMR9l1Mx+//XcWW0zcRfz0L8ddP49N/z6Jv80AMbBGI0ktobYxPwdR/ziAl8/7waaCHIz7uH16tlrcgIqKaybb+WtdQBUoVpv5zGgAwvHMDNPB11ewIjNb8NtCDVWJ8ryao46zA+ZvZWHrgioUjLb9ZWy8gIS0H9V3UaFlwULOx+dOGC/uHA8M2AnVCNXPPFvQG0jTDn3KZgC6NffH94JY4+OEj+KhfOJr4u6GwWI0/j1/Hiwv/w6dxcvy46zJWHE7GyGXHdJIrAEjNLMDIZcewMT7Fgs+YiIjINCZYVWDuzku4djcfgR6OGN2j4f0dJUOE6QlAYY7Bx3o622PCo5oJ799suWBTF38+euWOdk7ZD62uQyguALwbAkFlXMDaKwwYtgnwCwdyUoFFfYBrR3WLuNjj1U5h2DimM/58swMGxQTDxV6O9AIBX229iA/WnIKhwcCSbVP/OcPhQiIisiomWBaWfDsPc3ddAgD877FwONuXGpV19QPcggCIBie6l3iudTCi6nkgu7AYM9afs3DE5skvUmH8qpMQReCplnURcXuzZkfzZwDBxOV93AKAIf8CdVsD+XeBpY8Dl3fpFRMEAS3q18FnT0Vi//tdMfghFRr7uZRZtQggJbMAhxPvVPCZERERVR4TLAv7ZN1pFBWr0bGhN/o2D9AvUMZE9xIymYBP7k14X3P8Oo4kWT95mLnxHBLTcxHg7ogpPXyByzs0O5o/Y14Fzl7Ay38DYV2Bohzg12eAc/8aL25vh7Z+It7o2sCs6m2pp4+IiGofJlgWtP3cTWw9mwY7mYCpjzfTTGx/kBkJFgBEBXviudbBAIDJf5+26kro+y+lY/H+JADAzKcj4X5pHSCqgaCWgPdD5lfk4Aq8sAp4uB+gKgR+fwmI+63Mh/i5OZhVtbO93Pw4iIiIJMYEy0I0E9vPAACGdQpDQz83wwWDojW/TSRYADDh0Yfh4aTA2ZQs/HooWaJIyyensBjvrToJABjctj66NvYFTq7U7Ix8tvwV2jkAzywBol8ARBXw1xvAwR+NFm8dUgeBHo4wMQiJMSvi8MWmc7iTy2sfEhFR1WOCZSHzd1/Gldt58Hd3wNuPNDJcKOMqILs3JyvtLJB8CLgRd/8n46pOcS8Xe4zv3QQA8OXm80jPKbRY/MZM+/cMrmfko14dJ3zYtylw5zJw/T9AkAHNnqpYpXI74PE5QLu3NPc3vg/snAGddRlKisoEfNw/HAD0kqyS+3U9HZFbpML3Oy6h08ztmLHhHG5b4bWyVVw/jIjI8rgOlgVcu5uH73deBAB82LcpXB0MvMwZV4E5rYDikj/8amBhL90ydg7AqKOAZ7B20+CY+vj9SDLir2fh843n8PnTURZ6Fvp2nk/Db4c1Sd8XT0dpntfBPzQ7w7oCbv4Vr1wmA3pPA5zqADs+BXZ+ppkA3/szzb5SHo0IxNwXW+qtgxVwbx2sXuEB2HL2Jr7bloDTN7Lw465LWLI/CS+1D8GIzg3ga+YwY03E9cOIiKoGEywL+HTdWRQo1Wgb5oXHo4IMF8q7XSq5MqK4UFOuVIIllwmY+ngEBs7dj5X/XcPzMfXRsn4dCaM3LDNPiQ9Wa850HNoxFO0f8tb0MJUMD5o7ub0sggB0fQ9w9AA2vAcc+hEoyNT0bj3g0YhAxIYH4HDiHaRlF8DPzRExYV6QyzT9WL2bBaBXuD+2nU3Dt9sScOp6JubtvoylB5LwQtsQvN61AfzcHCsfczWyMT4FI5cd01viomT9sLkvtqx2SZZKLRp9D/A4PA6PU/uOcyjxDo6mC/BOvIP2Df0schxzMcGS2O4Lt7DxdCrk2jP/pG/cViF18HSrevjj6DVM/jsef7/VyeJvoqn/nEZqVgHCfFwwofe9C1WnnABuJwByB6Bpf+kO1vY1TZL110jgxG9AQRbw5E96xeQyQZPoGSEIAnqG++ORpn7YcT4N325NwIlrmViwNxHLDl7B4Lb18UbXh+DvXvMTLZVaxI9rdyFcSDe4XwDw49psxIY/J8l7qSq+UKuqN47H4XFMHeeRJj5Vcpya9rpZ7jhyLE34z+q984IoGpjoUk1kZWXBw8MDmZmZcHd3t8gxlEol1q9fj759+0KhKPtixoXFKvT5Zg8up+diWMcwTL43V8igG3HAvK6mA3ht1/2J8KWk5xSi+5c7kV1QjE+fjMCL7UJM11VBm06n4vVfjkImAKve6IBWIfd6zDZNAg7MAcKfAJ5dKv2Bz28AVr4CqAqhDu2MDe4volf/gSbbwRhRFLHrwi18uy0Bx5MzAAD2djIMahOMN7o9hEAPJwmDL4eMq5qeSmOcvXV6McvrVnYhNu8/goH7n4SjYPzajgWiAj9ErEDjJuHwc3OEr5sD/Nwc4PLAELepz0RVfKEa640rSeGk6o2z9eOU5/upMscpr9p2nNnPR0F15ajZ7VDR49S01626HaeEubkHe7AktGBvIi6n58LH1QFjYo1MbJeIj6sD3o1tjCn/nMEXm86jb/NAeLnYS36cO7lFmPSnZmjwtS4P3U+u1CogfrXmdvMKnD1ojiZ9gBdXA789D1nSHnRwvg7kdQc8KjbXSxAEdGvih66NfbH3Yjq+3ZqA/67cxZIDV/Db4at4tk09jOzWEHU9dRMti/bG6M3FM8DAXDxDVGoRiek5OH0jC2dTsnEmJQtnU7JwK7sQzYREvOBgPLkCAEdBiW3HzuG7o7qxONvL4efmcC/hcoSXiwJ3bwjIP3YdAZ7O2mTMy8UeW86kVnoYUhRFFBarUahUI1+pQoFShYJiFfKLVChQqpFXWGxyNf8Jf5xE8p08yAQBgiBAgGYEWvNb0N6Gzj4BMuH+bTVEfLb+XJnHeX/1KWTkKSEImhFzESW/xVL3790WxVL7NfcBoFitxvc7LpV5nPf+OIlLt3J13ncCAJVahfM3BNzYmwS5XAbh3p8UQx3nalHE7O0XTb5u1zPyNa+bgXLaYxvpmRcEzXG+3nzB5HFSMgt0jvFg+dL/+os62zX3zH0+Je+DijLnOBP/Oo0efgJS9iVBLpPf26f7COPP5/5xftxp+n2QkJYDmSDgwb4RY/U/uE8tiliwN7HM44xfdRKnb2SVPQJjpG+mZKtaFLFkf5LJ45xPzYZcJpT6XOp+DgUB2s+wTCj1+S31Hv1i03mjxxGgubpHbHhAlQ8XsgfLBHP/Q7yRkY9HvtqFfKUKXz8bhada1iu74kr2YAFAsUqNfrP34lxqNgbF1MdnTzU3XV85vfXrMfx7KgWN/V3xz+hOcLC7t75U4m5gSX/NUN74BE0SYCnXj0FcNhBC/h2IPk0gPPnD/bMvDTGzx0cUNWfTfbMtQbvyu0Iu4JnWwXiz20OoV8fZ8r0xFXwfZBcocS41G2dTsnDmhiaROpeajcJi/fXRBAHo4nIdS4rfM3mYD3xm47JdI6RlFyAtuxB5RSqzn4pMAIKQDg9kGy2TZ+eBsAZNUKhSaxOmAm0SdW9bscrYd7dWENJRRzB+nLuiG26g8sM2PA6PQzXDbyPalTmlpDzYg1XFpq0/i3ylCm1C62BAi7pVckw7uQyfPBGBZ386gBVHkjEoJhiR9Twlq/+fEzfw76kU2MkEfPVM9P3kCrg/uT38CcsmVwBQtyWKX/oHxYv6wSn9PPDzI2WXN7PHRxAEdGjogw4NfXDw8m18uzUBBy7fxvJDyVh55Crahnlh3yX9oTtrTAo/e2QLEgqO4sadbKTczcbdnHzYQQW5oIYdVGgONVpABSd7EYFudgh0s4O/qx18Xezg4yyDPCcVMH41Jq1pTzaHvN79a0nmFhYjLbsQt7ILkZZdgFvZhUjNyMfxc5fg4OGLWzlFSM8pxO3cIgSI6djq8K7JYcge578y+4+eTACcFHI42cvhYCeHo0IGz6Kb+LXA9HHG+M6Hg0+obk8SAIia/671e5o0O0vuq+5ewU8Zpo/zuudPsKtT/16vkVCqpwz3etA0/4mjdA9aqf0Ft5Iw65bp44z1XwAX31Bt8ilChKhW49r166gbVBcymUz7X3zp/5tLbuXfSsJ36aaPM9p7Phx8QnQeW6YHChXevoI5t00fZ1Sp4+gtufJAz4mgs+9eHelJ+NaM122M73w4+oSaeBLGXbubhxtXLppM5BzcvBDZoC7kpc96fuCJle4XLP0UBQDJd/JwNfGCyeOENGiMEG8XvToePOCD+0ruJqXnGvxee1Cnhj4I83Ep88pnxnYJgoDLt3KwO8HwnM/S2j/kjRAvZ4gln0tofuPe57Dks6o28PlViyJuZOYj/nqWyeNY4+oeTLAksO9iOv49mQKZAEx93DIT242JCfPCgBZ18efx6/jo79P4c2QHyCToBk3LLsBHf8cDAN7q3hDN63nc31lcCJxZq7ktxdmD5vBtgj2N/ofYq19ByL5RdlkDZ1+a0q6BN9q95o3DiXfw3bYE7L2YbvRLqHS3c8+mmuFK9b0Pu+bn3heAGlBpt2m+EFTqe7eLiyC/fQ4OF/6FOf9TNT3+KZqW3lDWaHDuvZ8KkC+MBTzqaq6R6R4IF/cghLkFIcw9EKhTF6gfCKVjKNYXJ6Bv31baXl2lSo1NWzfB8YDpYcjnI5wREhGtTZicFHI43kuiHO9tc7x3WyEX9D5PJw/vguN608d5M8YLkTFlXHjchJOH88w6zrsdfRAZ06YSx8k06zgjW3siMkZ3WRZND/tV9O3b3OTcn5OHs8w6zuh2XoiMaWle8AaPk23Wcd6u9HHMez6VfR8cPXkSzVJNJ3JLAj/HsKcfrfAcrKMnT6LZDdPHOd1uO1pFRlboGABw4NJtJF46bzKRe6t75Xp8Dly6bVaC9XaPRpU+zqD5B02Ws8ZZ4zUjwSrKBYoMXBpFkAMKR91yxggyQOGkX1aphFxVqLkvKvTKFhWrMf3vY3BCAV5oUx/hPvIHjiMA9s6l6s0DIGoeL7cHVGWsNC6315QzFre95r+YiX0exu4zV3HhaipWHzyPZ1obSCzsS10kWVmgWTXdCFHhjIn35pZEBzpiVKdA3RjOrwcKMwHXAKB+h/vbiwsBdbHx56Nwvv9vlamydk73178qLgKK8lBk547iR2dCseol44/TPsd8w6+bnSNwb44EiosAte6XWUxdByx7OQK/HryCuesPwUPQ1CGHCnbQfc0yMl3RfFIGAKAICqigqdcOxVDg/nOzgwoPCTcQLktCcyEJTWVJaCpcg0MZX6QPSkJdCI7ucLS3g6ODA5wdHDRf5DK5ZrhUJte8X+wcNPcFAYBwb58dkH/n/py5sqiVwN0kzY8RdgB6y90hv17/XjIWAIVbAKKyc8x6Lt0buKN5uKcmrpLeT1EElPc+G1BpftTQ/AA6ZZsFGbkqwgOa1b3XdV/m5974d0QzP/P+UDbzU+gfQ+/7pOS5GXi8v3lzJ5v52gGF2ZrPBkTNa5afCYUyG8hNB+zs7nXDifd/K+6XbeZh3n/wzXztNP+giKLmeej83Pv8lvF9UqnXrRzfEc2CXM0/jlqt+32iLuOz98B3RLR7DuQmPquOghKNHTJ0/048SO6gWVQZAFRKve9+c48T7Z6j+7rJ7QH5vWOqijWXHTNGbo8Yr1zscHwXDjB+rELYwc6jE1Dy759aBRSX8f6RKQC7e+9jtRoozkdMXQe0dUtHcc5tnXe+CjIUww4CAIWrF2LqOhj/jBr8jtAVU9cBYe7AjSwlClHy2otwguZ1EAD4uzvqHqdcuYGBsmWVL/3QGjEH6wM3uDsY6LVp1EtzrbsS0wINNhAAIKQTMLTUxYY/b2D8zK6gFsBrOwEA83ZfQt9tvVDPyOnv8H0YeOvQ/fvftwVunTNc1jUA8GuquXByw1gg+wZw87Thss7ewITL2rsp3/RAYMZRw2UVzsCklPv3f30GSNhsuCyAVf3i8d4fJ2Evl+G/Jr/A/bLxizDjwxv3v2z/HAmcWG687HuXAJd7Q0P/vgsc+dl42XdOAnXunRm5+X/A/tnGy5bHiO1A3Vaa2/u+BbZMNlq0SJTDXjBvDtLQovewQ90Ccqjwhnwt3lOsMvmYTNEFd+GKUOGmybI7u/2BbkFqYHkZJxT0/RKIGaG5nbgHWNLPrNhN8n0YcPTUvB8zbwBiGYlxeTh7Ax71AAiaP6Q3442XdfEFvBtp/vjmZwFpZox3BkRp3pvJ+42XcXAHvBpAm5TcPGV08m65yOw0z09U30uE7mhuGyTAzME42yFT3E+6VEVl/sNWLt4N79UtANkpmsWGjfF6CLhzybx6Qztr5osKMiD9gvHvYAAIfxJwC9QkWTeOA1f2mXWIVLfmCMgu4335yMdASAfNeyNhC7Brhnmxm/LED0CLFzS3L2wy/R1Rr4158z7bvgH0mam5ff0oML+H8bJdPwC6T9TcTjsL/NDOrNBVggJysYykss1w4LGvNLdz04EvjF/rdquqBb4pfhoiAAcUYY3DVOP1Pnjm+xQP42UN5BFZObnwmJHNOViWdDOrAN9uTUBfqUYE5fZAzynAvB3Ape2A38NmP9TfwxHIkCaMT+5dQ3FMbCO4p9Wwt0j8GiDzmuYLND+jzKLmJlcAMC/8FGT5uyBLi4dgLIkv0W4kEPMaPOqEIWPtdOD45ybr93K2B1D1cwgAaP7glHx5pp4GfuxQZnGz5d0ue3mK0nJvaX7KI/WE6TKFWUBKXPnqNYe6GMgxnThrVLPkCii7B6gybl80v6y5yRUAJO0xv+yZv8wvW0qZyRUAbCvjD35l7PkKOLVS00Nm4jsNibs033/mKC7U9LTJzeuRrIgyk6ty6ik/jp7y45LVJ4Wa0YN164bhLFKCIUKlUolNmzajd+9e98fW75V9Z8Vx/B13AzH1HLFieFsjc5+MDBEaDkJTdukTwOWdmuw99hPjMet00+fj0KV0DFl8GDIBWPl6BzQLcjdS1vAQoSiKGLH0P2y9lIsW9T2x6vX2sFM/8N/pid+Bf8dq/tN8bZemXnOH/SoxRKgszNO0Q1QgFL9I1DNTwtlXc5kftwBNL6J7INTFSsj2f1Ox+hQuQEBzIDAKCIjU/K4Tev/5lBoqUF09AvmCniarVI3YCXlAc5Pd/9ovwwe79DOvAT92Mj0k/cbee71KpTzQ/a/Mz9L/TABAyklg0aMmnwtipwHeoQDkmtdBFDW9O8UF0PQkqR8Y6gIA4X7ZO4nA9jI+FyUe+RjwCtMMCQnC/WHT0r8FuWYIQpBptqkKcG8quuY4mz80fZx+3wA+jUrVKbv33eN0/1jFRfeqlZUqd29f+gVg1SumjzNoJRAQoVOvMj8X27ZvxSM9HoFCYV/qud07lr3r/W3XjpjXPsM2Az4NNcNNovp+e2hvqzTDaCX3lfmaz3JJmVtngT9fN32cvl9peqFK2llU33+fQbw3lFcMQA3t7ObS74vMa8DmSaaP03EsUKe+5rZapfkMFCs1z0NU3X9O6nu/S6ZuiypNDDmpwOk/TR4m3bkRvPyCIBNVmuOolZr4VcX36ldpen/VxZrjq4t1f4qLAEjUEyglmR1g56z5W2rvrPket3PSvA8VTve3Objf2+Z8f1rCAf2rcOh5+R+gbkvDa4uYMURo9vfO0I1A4L25a5UcIszKyoKHb1At6cGyd9FNIMoqV546AUBQQiV30Nwv9cfk4OXb+DvuBgQB+OjJ1pA5mjcfQCfZMqbjO5oEK2450H0S4Oxl+jEKJ7R9OBg9o9Lxz4kb+Gj9ZfzxhpEJ76XfLKX8evAKtl7KhYOdDF8+EwU7uQyQP1D27L3J7ZHPAw4PPGc7BwBmnlFYrrL2gCho2qF0ElyWYRsBFz8g6waQnaoZ3spKKfU7VTMEoVYCebc0P6WGqMy+Crqdo2bIODBa8zuoxb1hDvNqkLsFQCWzh1xtPPFRyewhd/HRJBhyMz+yMrnu+923CTD6WOUXNJXJAHsXg58Js9smrJPRpUfMciPOvATroR6VP445glpU7jhlJb2luflr5ryVZueCQoWnZp+pydXmto+dg+a9UGFm/s9er3XVtE+zJyt/HDMSrPjgl9Dx6Tchq+hCo+Yu2fLkj5rpE8UFmqRMVaj5h7W4ULNNVXT/ful9JbezbwJJu82PS10MFGVpfip48kyZlvbX/DNQOmmzc3zg9r0Ez1CZ/DvmHUfhZDwHKG9uYG9eIlwzEqwqplSp8fHfmrlRg2Pq655hJ4UG3TU9H6kngcPzgW7vm/3QSX2bYvvZmziWnIE1x6/j6VYm1uO6J/l2HqavPwsAmPDow3jI10DCmH1T08UMAM0Hmh2TVdg5Ad4PaX6MUas1CUfp5Cs7VZOU3ToPXDts+jhDN2j++6ooz2DI3z6G/afO46fdl5Gec/+PrY+rPV7v0gAdmjep1ErupY8lST1EZD1+TSufMJqTyL26VdOTWZSn6aVU5t3/0W7LvX9CkTJfd1vWDfOHZkX1vcdZIoOzHiZYFbD0wBWcv5mNOs4KvNe7ifQHEASg0xjgj2HA4Z+ADqPN6/kCEODhiLcfaYTPNpzDjA1nERvuDw+nsv+jUqtFjP/jBPKKVGgb5oWhHUINFzy9RvNBqNfm3sRgK3D21vyHbWrlc3P+A5fJAFdfzU+g7unvZn8JCWb3dRnnGYwOnYPRtmPVXBDVYqRsGx6Hx7GV49RWcgXgVEfzUxHmfoe+sk7T61+cr5m+oswvdTtP0ytXkrwZ2p51A7i0rWIxWhgTrHJKyy7AN1suAND09Hg6S395GgBA0ycAzxAg4woQ9+v9s8PMMLRjGFb+dxWXbuVi1pYLmPJ4szLLL9qfhMOJd+BsL8cXT0cZX0erZHHRqlr7yhCPeppFRC147T5rMXXxapvnGVw1bcPj8DhVeRwzEjlR7oAiOzOniVTiODUyYXRwA9wrsWDzjTgmWDXFjA3nkF1YjMh6HnjW0HpTUpHbaXqu1o8H9n8HtBpq9vwbezsZpj4egRcXHMLSA0l4rk0wmgYanoh36VYOPt+oOWV50mNNUd/bSE/Z7UvAjWOaCX/NBlToKUmGQ122q6rahsfhcarqOGYkcsX2Hsjfd9Lix6mqhLFGJnJWwASrHP5LuoM1x65DEIBPnoiw/PBNixeBnTOAjGTNqcPNnzb7oZ0a+aBv8wCsP5WKj/8+jd9fb6e3InaxSo13V55AYbEanRv5YHBMfeMVnrq3DkiDboCrX/mfS3XDLyEiKmEqkVMqAVQywTLnOFJgIldlmGCZqVilxkf3JrY/1zoY0cGelj+owglo+zqwYxqw7xsgYqDhU1mNmPRYOHacu4XDSXfwd9wNPPnANRLn7bmMuKsZcHOww8yBkcYv8SOK9xOsyDIWsatJqupLiIioqjGRqxJMsMy04r9rOJuSBQ8nBSY8av4CoJXWZjiw9xsg9ZRm8dGGJi50XEpdTyeM6tEQX2w6j2nrz+KRpn5wc9RMeD+XmoVZ9+aSffx4MwR5lnEK943jmgUA7ZyAhx+rzLOpXjgUSURUcVZI5JTFxdi3bx86duwIhd29FMdK/wwzwSqDSi3iUOId7EsVsO4/TTIyvncTeLlYaGK7Ic5eQKtXgIM/aHqxypFgAcDwzmH44+g1JKZrJrzHhgcgJSMf321PgFIlomdTPwxsWbfsSk79ofndpI9mQiIREZGtKJ3IKZXIdL6uOTO8omuSSUSCc8wr7/vvv0doaCgcHR3Rtm1bHD5sxvpDFrYxPgWdZm7Hiwv/w8pEOfKUatjJBHhb6qzBsrR7U7OibeJu4Pqxcj3UwU6Oj/uHAwAW7kvCoPkHMW7VCSTdzoMgALHh/saHBgHN6sMlFwm25tmDRERE1YjVE6zff/8d48aNw8cff4xjx44hKioKvXv3RlpamtVi2hifgpHLjiElU/fab8VqEW8tP4aN8SlGHmkhnsFAxL0J7vu+LffDC5SGV50VReCD1afKfj5JezSXinD0BBqavqQLERER2UCC9fXXX2PEiBEYOnQowsPD8eOPP8LZ2RkLFy60SjwqtYip/5wp82IPU/85A5W6ii/h2PFtze+zazVLJpip5PmUpcznc/Le5PZmT5a6ThgRERGVxapzsIqKinD06FFMnDhRu00mk6Fnz544cOCAXvnCwkIUFt4/UyArKwuA5oLMSqU0V+U+lHhHr+eqNBFASmYBDlxMQ9swM64RKBWvxpA/1BOyS1uh2vcd1H2+NOthlXo+xQWwO/s3BADF4U9BlOg1roiS9pWqnani2Ba2g21hG9gOtqMq2sLcuq2aYKWnp0OlUsHf319nu7+/P86dO6dX/rPPPsPUqVP1tm/evBnOzuZdSsaUo+kCALnJcpv3HMLts1Xbi+Uti0EnbAWO/4ptRa1QqDB9DcTKPJ/AjCOIKcxGnsILW07dAeLXVzR0yWzZssXaIdA9bAvbwbawDWwH22HJtsjLyzOrXLU6i3DixIkYN26c9n5WVhaCg4PRq1cvuLsbXqm8vLwT72Bpwn8my/Xq3LZqe7AAQOwD9ZLNkF//D7Hul6HuPsnkQyrzfOR/aC6N49D6BfTt0a9iMUtEqVRiy5YtiI2NhcLKZ4bUdmwL28G2sA1sB9tRFW1RMnpmilUTLB8fH8jlcty8eVNn+82bNxEQEKBX3sHBAQ4ODnrbFQqFZC9k+4Z+CPRwRGpmgcF5WAI0F1Ru39DPOhfi7TQG+P1FyI8thLzruyaXTajw88nPAC5uBgDIo56D3Ea+NKRsa6octoXtYFvYBraD7bBkW5hbr1Unudvb26NVq1bYtu3+hRrVajW2bduG9u3bWyUmuUzQLmvwYPpUcv/j/uHWSa4AoMljgHcjoCATOLrEZPEKP5+z/wCqIsC3KeAfUfm4iYiIahGrn0U4btw4zJ8/H0uWLMHZs2cxcuRI5ObmYujQoVaL6dGIQMx9sSUCPBx1tgd4OGLuiy3xaEQlrvxdWTLZ/TMKD3wPFBeZfEiFns8pzfAgmj9drsvzEBERkQ3MwXruuedw69YtTJ48GampqYiOjsbGjRv1Jr5XtUcjAhEbHoADF9Owec8h9Orc1nrDgg+KfA7YPg3IvgHE/wFEDzb5kJLnczjxDtKyC+Dn5oiYMC/DzycrBUjco7ldjgtMExERkYbVEywAGDVqFEaNGmXtMPTIZQLahnnh9lkRbY0lI9Zg5wC0Gwls/Viz8Gjk85qeLRPkMgHtHzLjgpfxqwGIQHBboE5opcMlIiKqbaw+REgV1Hoo4OAO3DoHJGyStu5T9xYX5aVxiIiIKoQJVnXl6KFJsoAKXT7HqPQEICUOEORAswHS1UtERFSLMMGqztq9CcjtgeQDQPIhaeos6b1q+Ajg4iNNnURERLUME6zqzC0AiHpec3vfN5WvTxSBkyVnD3J4kIiIqKKYYFV3Hd4GIADn1wO3zleuruvHgLuJgMIZaNJXkvCIiIhqIyZY1Z1PI+DhxzS3931XubpK1r5q0hdwcK1cXURERLUYE6yaoNNYze+TvwOZ1ytWh6oYiF+juR35rDRxERER1VJMsGqCeq2BkI6AWgkcmluxOhJ3AblpgJMX8FAPaeMjIiKqZZhg1RQdx2h+/7dYc6Hm8jr1h+Z3swGAnBcrJSIiqgwmWDVFo1jArxlQlA38t6B8j1Xmay7uDPDsQSIiIgkwwaopBAHo+I7m9sEfAWWB+Y+9sFGTmHnU11weh4iIiCqFCVZNEvEU4BGsmUt14jfzH3ey5NI4A826piERERGVjX9NaxK5Amj/lub2/u8Atcr0Y/LvAgmbNbc5PEhERCQJJlg1TcuXAac6wJ3LwLl1psuf+Vtz9qFfM8C/meXjIyIiqgWYYNU09i5AmxGa23u/0Vz+piwlZw82f9qiYREREdUmTLBqoravA3ZOwI1jQNIe4+UyrwNJezW3mWARERFJhglWTeTiA7R4UXN737fGy8WvBiAC9dsDnvWrJDQiIqLagAlWTdX+LUCQARe3AqmnDJc5VXL2ICe3ExERSYkJVk3lFaZZlR0w3It16zyQehKQ2QHhT1ZpaERERDUdE6yarGTh0fg1wN0ruvtKeq8a9gRcvKs2LiIiohqOCVZNFhgFNOgOiCrgwPf3t4sihweJiIgsiAlWTddpjOb3saVA7m3N7Wv/AXeTAIUL0KSPtSIjIiKqsZhg1XRhXYHAaKA4HzgyX7Pt1ErN74cf06ybRURERJJiglXTCcL9JRsOfA8k7QNO3kuw6rUGbsQBGVetFh4REVFNZGftAMjCMq4CmydpbhdmAYv73t+3YYLmt50DMOoo4Blc9fERERHVQOzBqunybgPFhWWXKS7UlCMiIiJJMMEiIiIikhgTLCIiIiKJMcEiIiIikhgTLCIiIiKJMcEiIiIikhgTLCIiIiKJMcGq6Zy9NetclcXOQVOOiIiIJMGFRms6z2DNIqJlrXPl7M1FRomIiCTEBKs28AxmAkVERFSFOERIREREJDEmWEREREQSY4JFREREJDEmWEREREQSY4JFREREJDEmWEREREQSY4JFREREJDEmWEREREQSY4JFREREJLFqvZK7KIoAgKysLIsdQ6lUIi8vD1lZWVAoFBY7DpWN7WA72Ba2g21hG9gOtqMq2qIk5yjJQYyp1glWdnY2ACA4mJeBISIioqqTnZ0NDw8Po/sF0VQKZsPUajVu3LgBNzc3CIKANm3a4MiRIwbLGtpnzrasrCwEBwfj6tWrcHd3l/5JlKGs52OpOswtb6qcsf3l2V56W3Vvh4rUY075irZDWftqelvUlM8EYL22sOXPhKly/ExIU4c1vp9s5TMhiiKys7MRFBQEmcz4TKtq3YMlk8lQr1497X25XG70BTW0z9xtAODu7l7lH5yyno+l6jC3vKlyxvaXZ7uhbdW1HSpSjznlK9oOZe2r6W1R0z4TQNW3hS1/JkyV42dCmjqs8f1kS5+JsnquStSoSe5vvfVWufaZu81apIilvHWYW95UOWP7y7PdVtpCqjgs0RYVbYey9tX0tuBnovJs+TNhqhw/E9LUYY3vJ1tuB0Oq9RBhVcjKyoKHhwcyMzOr/D8Tuo/tYDvYFraDbWEb2A62w5baokb1YFmCg4MDPv74Yzg4OFg7lFqN7WA72Ba2g21hG9gOtsOW2oI9WEREREQSYw8WERERkcSYYBERERFJjAkWERERkcSYYBERERFJjAkWERERkcSYYEkoMTER3bt3R3h4OJo3b47c3Fxrh1RrhYaGIjIyEtHR0ejevbu1w6nV8vLyEBISgvHjx1s7lForIyMDrVu3RnR0NCIiIjB//nxrh1RrXb16Fd26dUN4eDgiIyOxatUqa4dUaw0YMAB16tTB008/bZH6uUyDhLp27YpPP/0UnTt3xp07d+Du7g47u2p9NaJqKzQ0FPHx8XB1dbV2KLXepEmTcPHiRQQHB+PLL7+0dji1kkqlQmFhIZydnZGbm4uIiAj8999/8Pb2tnZotU5KSgpu3ryJ6OhopKamolWrVrhw4QJcXFysHVqts3PnTmRnZ2PJkiX4448/JK+fPVgSOX36NBQKBTp37gwA8PLyYnJFtV5CQgLOnTuHPn36WDuUWk0ul8PZ2RkAUFhYCFEUwf+trSMwMBDR0dEAgICAAPj4+ODOnTvWDaqW6tatG9zc3CxWf61JsHbv3o3+/fsjKCgIgiDgr7/+0ivz/fffIzQ0FI6Ojmjbti0OHz5sdv0JCQlwdXVF//790bJlS0yfPl3C6GsWS7cFAAiCgK5du6JNmzb49ddfJYq8ZqmKdhg/fjw+++wziSKuuaqiLTIyMhAVFYV69erhvffeg4+Pj0TR1yxV0RYljh49CpVKheDg4EpGXfNUZTtYSq1JsHJzcxEVFYXvv//e4P7ff/8d48aNw8cff4xjx44hKioKvXv3RlpamrZMyfyFB39u3LiB4uJi7NmzBz/88AMOHDiALVu2YMuWLVX19KoVS7cFAOzduxdHjx7F2rVrMX36dJw8ebJKnlt1Yul2+Pvvv9G4cWM0bty4qp5StVUVnwlPT0+cOHECiYmJWL58OW7evFklz626qYq2AIA7d+7g5Zdfxrx58yz+nKqjqmoHixJrIQDin3/+qbMtJiZGfOutt7T3VSqVGBQUJH722Wdm1bl//36xV69e2vuff/65+Pnnn0sSb01mibZ40Pjx48VFixZVIsqazxLt8MEHH4j16tUTQ0JCRG9vb9Hd3V2cOnWqlGHXSFXxmRg5cqS4atWqyoRZK1iqLQoKCsTOnTuLS5culSrUGs2Sn4kdO3aIAwcOlCJMPbWmB6ssRUVFOHr0KHr27KndJpPJ0LNnTxw4cMCsOtq0aYO0tDTcvXsXarUau3fvRtOmTS0Vco0lRVvk5uYiOzsbAJCTk4Pt27ejWbNmFom3ppKiHT777DNcvXoVSUlJ+PLLLzFixAhMnjzZUiHXWFK0xc2bN7WficzMTOzevRtNmjSxSLw1mRRtIYoihgwZgh49euCll16yVKg1mhTtUBU4CxtAeno6VCoV/P39dbb7+/vj3LlzZtVhZ2eH6dOno0uXLhBFEb169UK/fv0sEW6NJkVb3Lx5EwMGDACgOXtqxIgRaNOmjeSx1mRStANJQ4q2uHLlCl577TXt5PbRo0ejefPmlgi3RpOiLfbt24fff/8dkZGR2nlFv/zyC9ujHKT6furZsydOnDiB3Nxc1KtXD6tWrUL79u0li5MJloT69OnDs6VsQIMGDXDixAlrh0GlDBkyxNoh1GoxMTGIi4uzdhgEoFOnTlCr1dYOgwBs3brVovVziBCAj48P5HK53qTPmzdvIiAgwEpR1U5sC9vAdrAdbAvbwbawDdWlHZhgAbC3t0erVq2wbds27Ta1Wo1t27ZJ2l1IprEtbAPbwXawLWwH28I2VJd2qDVDhDk5Obh48aL2fmJiIuLi4uDl5YX69etj3LhxeOWVV9C6dWvExMTgm2++QW5uLoYOHWrFqGsmtoVtYDvYDraF7WBb2IYa0Q4WOTfRBu3YsUMEoPfzyiuvaMvMnj1brF+/vmhvby/GxMSIBw8etF7ANRjbwjawHWwH28J2sC1sQ01oB16LkIiIiEhinINFREREJDEmWEREREQSY4JFREREJDEmWEREREQSY4JFREREJDEmWEREREQSY4JFREREJDEmWEREREQSY4JFREREJDEmWEREZhoyZAiefPJJa4dBRNUAEywisqpbt25h5MiRqF+/PhwcHBAQEIDevXtj37592jKCIOCvv/6yXpD3fPvtt1i8eLG1wyCiasDO2gEQUe02cOBAFBUVYcmSJWjQoAFu3ryJbdu24fbt29YOTY+Hh4e1QyCiaoI9WERkNRkZGdizZw9mzpyJ7t27IyQkBDExMZg4cSIef/xxAEBoaCgAYMCAARAEQXsfAP7++2+0bNkSjo6OaNCgAaZOnYri4mLtfkEQMHfuXPTp0wdOTk5o0KAB/vjjjzJj+uOPP9C8eXM4OTnB29sbPXv2RG5uLgDdIcKkpCQIgqD3061bN21de/fuRefOneHk5ITg4GC8/fbb2rqIqGZjgkVEVuPq6gpXV1f89ddfKCwsNFjmyJEjAIBFixYhJSVFe3/Pnj14+eWX8c477+DMmTP46aefsHjxYkybNk3n8R999BEGDhyIEydO4IUXXsDzzz+Ps2fPGjxWSkoKBg0ahGHDhuHs2bPYuXMnnnrqKYiiqFc2ODgYKSkp2p/jx4/D29sbXbp0AQBcunQJjz76KAYOHIiTJ0/i999/x969ezFq1KgKv15EVH0IoqFvDiKiKrJ69WqMGDEC+fn5aNmyJbp27Yrnn38ekZGR2jKCIODPP//UmWDes2dPPPLII5g4caJ227JlyzBhwgTcuHFD+7g33ngDc+fO1ZZp164dWrZsiR9++EEvlmPHjqFVq1ZISkpCSEiI3v4hQ4YgIyNDbz5YQUEBunXrBl9fX/z999+QyWQYPnw45HI5fvrpJ225vXv3omvXrsjNzYWjo2O5Xysiqj7Yg0VEVjVw4EDcuHEDa9euxaOPPoqdO3eiZcuWJieTnzhxAp988om2F8zV1RUjRoxASkoK8vLytOXat2+v87j27dsb7cGKiorCI488gubNm+OZZ57B/PnzcffuXZPPYdiwYcjOzsby5cshk8m08S1evFgnvt69e0OtViMxMdFknURUvXGSOxFZnaOjI2JjYxEbG4uPPvoIw4cPx8cff4whQ4YYfUxOTg6mTp2Kp556ymB9FSGXy7Flyxbs378fmzdvxuzZszFp0iQcOnQIYWFhBh/z6aefYtOmTTh8+DDc3Nx04nv99dfx9ttv6z2mfv36FYqPiKoP9mARkc0JDw/XmQyuUCigUql0yrRs2RLnz59Hw4YN9X5KepEA4ODBgzqPO3jwIJo2bWr02IIgoGPHjpg6dSqOHz8Oe3t7/PnnnwbLrl69Gp988glWrlyJhx56SC++M2fOGIzP3t7e7NeCiKon9mARkdXcvn0bzzzzDIYNG4bIyEi4ubnhv//+w+eff44nnnhCWy40NBTbtm1Dx44d4eDggDp16mDy5Mno168f6tevj6effhoymQwnTpxAfHw8Pv30U+1jV61ahdatW6NTp0749ddfcfjwYSxYsMBgPIcOHcK2bdvQq1cv+Pn54dChQ7h165bBhCw+Ph4vv/wy3n//fTRr1gypqakAAHt7e3h5eeH9999Hu3btMGrUKAwfPhwuLi44c+YMtmzZgjlz5kj8ShKRzRGJiKykoKBA/OCDD8SWLVuKHh4eorOzs9ikSRPxf//7n5iXl6ctt3btWrFhw4ainZ2dGBISot2+ceNGsUOHDqKTk5Po7u4uxsTEiPPmzdPuByB+//33YmxsrOjg4CCGhoaKv//+u9F4zpw5I/bu3Vv09fUVHRwcxMaNG4uzZ8/W7n/llVfEJ554QhRFUVy0aJEIQO+na9eu2vKHDx8WY2NjRVdXV9HFxUWMjIwUp02bVvkXjohsHs8iJKIay9DZh0REVYFzsIiIiIgkxgSLiIiISGKc5E5ENRZnQBCRtbAHi4iIiEhiTLCIiIiIJMYEi4iIiEhiTLCIiIiIJMYEi4iIiEhiTLCIiIiIJMYEi4iIiEhiTLCIiIiIJPb/KmcDqhcvb54AAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 700x500 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# NBVAL_SKIP\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# eps-Werte, über die wir scannen\n",
    "eps_values = jnp.logspace(-6, -1, 20)  # von 1e-6 bis 1e-1\n",
    "\n",
    "age_fd_values = []\n",
    "metal_fd_values = []\n",
    "\n",
    "for eps in eps_values:\n",
    "    g_fd = finite_diff_grad(f, params_init, eps=float(eps))\n",
    "    # g_fd hat die gleiche Struktur wie params_init:\n",
    "    # {'age': array([..,..]), 'metallicity': array([..,..])}\n",
    "    # Beispiel: nimm hier den Mittelwert pro Array\n",
    "    age_fd_values.append(float(jnp.mean(g_fd['age'])))\n",
    "    metal_fd_values.append(float(jnp.mean(g_fd['metallicity'])))\n",
    "\n",
    "plt.figure(figsize=(7,5))\n",
    "plt.semilogx(eps_values, age_fd_values, 'o-', label=\"age grad (FD)\")\n",
    "plt.semilogx(eps_values, metal_fd_values, 's-', label=\"metallicity grad (FD)\")\n",
    "\n",
    "# horizontale Linien = JAX-Gradient\n",
    "plt.axhline(float(grads['age'][0]), color='C0', linestyle='--', label=\"age grad (JAX)\")\n",
    "plt.axhline(float(grads['metallicity'][0]), color='C1', linestyle='--', label=\"metalicity grad (JAX)\")\n",
    "\n",
    "plt.xlabel(\"Step size\")\n",
    "plt.ylabel(\"Derivation\")\n",
    "# plt.title(\"Gradient vs finite difference step size\")\n",
    "plt.legend()\n",
    "plt.grid(True)\n",
    "plt.savefig(\"output/optimisation_finite_diff.jpg\", dpi=1000)\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "rubix",
   "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.12.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}