research.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "a445e139-a70f-42a7-b622-ce7b59ecce1c",
   "metadata": {},
   "outputs": [],
   "source": [
    "#! /usr/bin/python\n",
    "# -*- coding: utf-8 -*-\n",
    "# @author izhangxm\n",
    "# Copyright 2017 izhangxm@gmail.com. All Rights Reserved.\n",
    "#\n",
    "# Licensed under the Apache License, Version 2.0 (the \"License\");\n",
    "# you may not use this file except in compliance with the License.\n",
    "# You may obtain a copy of the License at\n",
    "#\n",
    "#     http://www.apache.org/licenses/LICENSE-2.0\n",
    "#\n",
    "# Unless required by applicable law or agreed to in writing, software\n",
    "# distributed under the License is distributed on an \"AS IS\" BASIS,\n",
    "# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.\n",
    "# See the License for the specific language governing permissions and\n",
    "# limitations under the License.\n",
    "# =============================================================================="
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "09fda08e-37e5-470c-aca7-49be59cc4cf9",
   "metadata": {},
   "outputs": [],
   "source": [
    "import arviz as az\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import pymc as pm\n",
    "from itertools import product\n",
    "import os.path as osp\n",
    "from scipy.optimize import leastsq\n",
    "import time\n",
    "%config InlineBackend.figure_format = 'retina'\n",
    "az.style.use(\"arviz-darkgrid\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "f0c5995b-88f0-4e33-abd7-f5c5873d83d9",
   "metadata": {},
   "outputs": [],
   "source": [
    "class MyDataset(object):\n",
    "\n",
    "    def __init__(self, dataset_path):\n",
    "\n",
    "        self.dataset_path = dataset_path\n",
    "\n",
    "        print(f\"loading dataset from {self.dataset_path}....\")\n",
    "        df = pd.read_csv(dataset_path)\n",
    "        # 数据初步处理\n",
    "        # 计算反应速率rate，初始速率固定设置为0\n",
    "\n",
    "        for c_i, (col_name, col_sr) in enumerate(df.items()):\n",
    "            if \"error\" in col_name or \"time\" in col_name or \"rate\" in col_name:\n",
    "                continue\n",
    "            rate_col_name = f\"{col_name}_rate\"\n",
    "            rates = []\n",
    "            pre_t = None\n",
    "            pre_v = None\n",
    "            for th, (index, value) in zip(df['time'],col_sr.items()):\n",
    "                if int(index) == 0:\n",
    "                    rates.append(0.0)\n",
    "                    pre_t = th\n",
    "                    pre_v = value\n",
    "                    continue\n",
    "\n",
    "                delta_t = th-pre_t\n",
    "                delta_value = value - pre_v\n",
    "                # print(col_name, index, pre_t, th, pre_v ,value)\n",
    "                rates.append(delta_value/delta_t)\n",
    "                pre_t = th\n",
    "                pre_v = value\n",
    "            df[rate_col_name] = rates\n",
    "\n",
    "        self.df = df\n",
    "         # 准备输出值 Y\n",
    "        self.cct_names = []\n",
    "        for x in self.df.columns:\n",
    "            if \"time\" in x or \"error\" in x or \"rate\" in x:\n",
    "                continue\n",
    "            self.cct_names.append(x)\n",
    "        self.rates_names = [f\"{x}_rate\" for x in self.cct_names]\n",
    "        self.error_names = [f\"{x}-error\" for x in self.cct_names]\n",
    "\n",
    "        self.cct = self.df[self.cct_names].values\n",
    "        self.rates = self.df[self.rates_names].values\n",
    "        self.errors = self.df[self.error_names].values\n",
    "\n",
    "        print(\"done!\")\n",
    "\n",
    "    def get_rates(self):\n",
    "        return self.rates\n",
    "\n",
    "    def get_df(self):\n",
    "        return self.df\n",
    "\n",
    "    def get_errors(self):\n",
    "        return self.errors\n",
    "\n",
    "    def get_cct(self):\n",
    "        return self.cct\n",
    "\n",
    "    def get_var_col_names(self):\n",
    "        return self.cct_names, self.rates_names, self.error_names\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "1c450d25-38a7-4cdb-9e14-483b5c628fc6",
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_target(ks, df, k_kinetics):\n",
    "\n",
    "    def _swt_k_kinetics(k_i, sr_vs):\n",
    "        kk = k_kinetics[k_i]\n",
    "        k = ks[k_i]\n",
    "        if kk == 0:\n",
    "            return kk\n",
    "        elif kk ==1:\n",
    "            return k * sr_v\n",
    "        elif kk == 2:\n",
    "            return k * sr_v ** 2\n",
    "        raise Exception(f\"unknown kinetics:{kk}\")\n",
    "\n",
    "    target = [] # 'xNH3', 'xNO3', 'xNO2', 'xNOrg', 'xN2', 'ANH3', 'ANO3', 'ANO2', 'ANOrg', 'AN2'\n",
    "    for i in range(0,len(df)):\n",
    "        sr = df.iloc[i]\n",
    "\n",
    "        r1 = ks[0] * sr['xN2'] if k_kinetics[0] == 1 else ks[0]\n",
    "        r2 = ks[1] * sr['xNH3'] if k_kinetics[1] == 1 else ks[1]\n",
    "        r3 = ks[2] * sr['xNO2'] if k_kinetics[2] == 1 else ks[2]\n",
    "        r4 = ks[3] * sr['xNO3'] if k_kinetics[3] == 1 else ks[3]\n",
    "        r5 = ks[4] * sr['xNO2'] if k_kinetics[4] == 1 else ks[4]\n",
    "        r6 = ks[5] * sr['xNO2'] * sr['xNO3'] if k_kinetics[5] == 1 else ks[5]\n",
    "        r7 = ks[6] * sr['xNO3'] if k_kinetics[6] == 1 else ks[6]\n",
    "        r8 = ks[7] * sr['xNO3'] if k_kinetics[7] == 1 else ks[7]\n",
    "        r9 = ks[8] * sr['xNH3'] if k_kinetics[8] == 1 else ks[8]\n",
    "        r10 = ks[9] * sr['xNOrg'] if k_kinetics[9] == 1 else ks[9]\n",
    "        r11 = ks[10] * sr['xNOrg'] if k_kinetics[10] == 1 else ks[10]\n",
    "\n",
    "        xNH3_rate =  2*r1 + r7 + r10 - r2 - r6 - r9\n",
    "        xNO3_rate = r3 - r7 - r4 - r8 + r11\n",
    "        xNO2_rate = r2 + r4 - r3 - r6 - 2*r5\n",
    "        xNOrg_rate = r8 + r9 - r10 -r11\n",
    "        xN2_rate = r5 + r6 - r1\n",
    "        ANH3_rate = (2*r1*(sr['AN2'] - sr['ANH3']) + (sr['ANO3']-sr['ANH3'])*r7 + (sr['ANOrg']-sr['ANH3'])*r10 )/sr['xNH3']\n",
    "        ANO3_rate = ( (sr['ANO2'] - sr['ANO3'])*r2 + (sr['ANOrg'] - sr['ANO3'])*r11 ) / sr['xNO3']\n",
    "        ANO2_rate = ( (sr['ANH3']-sr['ANO2'] )*r2 + (sr['ANO3']-sr['ANO2'])*r4 ) / sr['xNO2']\n",
    "        ANOrg_rate = ( (sr['ANO3']-sr['ANOrg'] )*r8 + (sr['ANH3']-sr['ANOrg'])*r9 ) / sr['xNOrg']\n",
    "        AN2_rate = ( (sr['ANO2']-sr['AN2'] )*r5 + (sr['ANO2']*sr['ANH3'] - sr['AN2'])*r6 ) / sr['xN2']\n",
    "\n",
    "        line_rate = [xNH3_rate, xNO3_rate,xNO2_rate, xNOrg_rate, xN2_rate, ANH3_rate,ANO3_rate, ANO2_rate,ANOrg_rate,AN2_rate]\n",
    "        target.append(line_rate)\n",
    "    target = np.array(target)\n",
    "    return target\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "a80b5701-2324-4e1e-bfb9-fbf50a15b016",
   "metadata": {},
   "outputs": [],
   "source": [
    "def r2_loss(pred, y):\n",
    "    r2_loss = 1 - np.square(pred - y).sum() / np.square(y - np.mean(y)).sum()\n",
    "    return r2_loss\n",
    "\n",
    "def get_model(dataset, k_kinetics, k_sigma_priors = 0.1, kf_type=0):\n",
    "    # 定义参数优化模型\n",
    "    mcmc_model = pm.Model()\n",
    "    ## 参数个数\n",
    "    params_n = 11\n",
    "\n",
    "    ks = []\n",
    "    with mcmc_model:\n",
    "        for ki in range(1, params_n+1):\n",
    "            if kf_type == 0:\n",
    "                p_dense = pm.HalfNormal(f\"k{ki}\", sigma=k_sigma_priors)\n",
    "            else:\n",
    "                p_dense = pm.Normal(f\"k{ki}\",mu=0, sigma=k_sigma_priors)\n",
    "            ks.append(p_dense)\n",
    "\n",
    "    df = dataset.get_df()\n",
    "    errors = dataset.get_errors()\n",
    "    rates = dataset.get_rates()\n",
    "\n",
    "\n",
    "    target= get_target(ks, df, k_kinetics)\n",
    "    target = np.array(target)[1:].reshape(-1).tolist()\n",
    "    sigma_Y = errors[1:].reshape(-1).tolist()\n",
    "    rata_Y = rates[1:].reshape(-1).tolist()\n",
    "\n",
    "    with mcmc_model:\n",
    "        sigma = pm.HalfCauchy('sigma', beta=10, initval=0.1)\n",
    "        y_obs = pm.Normal(f\"rates\", mu=target, sigma=sigma, observed=rata_Y, shape=len(rata_Y))\n",
    "\n",
    "    return mcmc_model\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "6fb28131-62e8-4431-80b6-d274fbde85bb",
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_predict_ks(idata):\n",
    "    parames_summary = az.summary(idata, round_to=10)\n",
    "    ks_names = [f\"k{x+1}\" for x in range(11)]\n",
    "\n",
    "    predict_ks = []\n",
    "    for k_name in ks_names:\n",
    "        k_v = parames_summary[\"mean\"][k_name]\n",
    "        predict_ks.append(k_v)\n",
    "    return np.array(predict_ks)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "8f1b8490-11a4-47be-9d5f-ae6781749648",
   "metadata": {},
   "outputs": [],
   "source": [
    "def opt_model(dataset, k_kinetics, k_sigma_priors=0.01,  kf_type=0, draws=10000, tune=2000, chains=4, cores=4):\n",
    "    mcmc_model = get_model(dataset, k_kinetics, k_sigma_priors=k_sigma_priors,kf_type=kf_type)\n",
    "    idata = pm.sample(draws=draws,model=mcmc_model, chains=chains, cores=cores, tune=tune)\n",
    "    return idata\n",
    "\n",
    "def eval_model(idata, dataset):\n",
    "    predict_ks = get_predict_ks(idata)\n",
    "    predict = get_target(predict_ks, df, k_kinetics)\n",
    "    rates_y  = dataset.get_rates()\n",
    "    r2 = r2_loss(predict[1:],rates_y[1:])\n",
    "    return r2\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "7b389085-39f5-49a8-b27a-f0022036688f",
   "metadata": {},
   "outputs": [],
   "source": [
    "def ltq_fit(dataset,k_kinetics):\n",
    "    def _error_loss(ks, dataset):\n",
    "        df = dataset.get_df()\n",
    "        rates_y = dataset.get_rates()\n",
    "        predict= get_target(ks, df, k_kinetics)\n",
    "\n",
    "        r2 = r2_loss(predict[1:],rates_y[1:])\n",
    "        res =  (rates_y[1:] - predict[1:]).reshape(-1)\n",
    "\n",
    "        is_nagative = False\n",
    "        for x in ks:\n",
    "            if x<=0:\n",
    "                is_nagative = True\n",
    "                break\n",
    "        # if is_nagative:\n",
    "        #     res = res + 1000\n",
    "        # print('step', ks,r2, res.shape)\n",
    "        return res\n",
    "\n",
    "    ks_o = np.repeat(1,11).tolist()\n",
    "    ks_res =leastsq(_error_loss, ks_o, args=(dataset,))[0]\n",
    "    return ks_res\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "1c81a624-0229-427d-a14b-60724b2498a7",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "loading dataset from dataset/data.csv....\n",
      "done!\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>time</th>\n",
       "      <th>xNH3_rate</th>\n",
       "      <th>xNO3_rate</th>\n",
       "      <th>xNO2_rate</th>\n",
       "      <th>xNOrg_rate</th>\n",
       "      <th>xN2_rate</th>\n",
       "      <th>ANH3_rate</th>\n",
       "      <th>ANO3_rate</th>\n",
       "      <th>ANO2_rate</th>\n",
       "      <th>ANOrg_rate</th>\n",
       "      <th>AN2_rate</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0.5</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000e+00</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>48.0</td>\n",
       "      <td>-0.011153</td>\n",
       "      <td>0.011779</td>\n",
       "      <td>0.000381</td>\n",
       "      <td>0.000674</td>\n",
       "      <td>2.547368e-07</td>\n",
       "      <td>-0.050772</td>\n",
       "      <td>0.032181</td>\n",
       "      <td>0.000227</td>\n",
       "      <td>0.000528</td>\n",
       "      <td>0.000022</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>96.0</td>\n",
       "      <td>-0.007440</td>\n",
       "      <td>0.005084</td>\n",
       "      <td>-0.000385</td>\n",
       "      <td>-0.021708</td>\n",
       "      <td>3.500000e-07</td>\n",
       "      <td>-0.061275</td>\n",
       "      <td>0.010557</td>\n",
       "      <td>-0.000206</td>\n",
       "      <td>0.000374</td>\n",
       "      <td>0.000065</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>144.0</td>\n",
       "      <td>-0.001736</td>\n",
       "      <td>0.004588</td>\n",
       "      <td>0.000058</td>\n",
       "      <td>0.021417</td>\n",
       "      <td>4.375000e-07</td>\n",
       "      <td>-0.027376</td>\n",
       "      <td>-0.001019</td>\n",
       "      <td>0.000081</td>\n",
       "      <td>0.000273</td>\n",
       "      <td>0.000380</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "    time  xNH3_rate  xNO3_rate  xNO2_rate  xNOrg_rate      xN2_rate  \\\n",
       "0    0.5   0.000000   0.000000   0.000000    0.000000  0.000000e+00   \n",
       "1   48.0  -0.011153   0.011779   0.000381    0.000674  2.547368e-07   \n",
       "2   96.0  -0.007440   0.005084  -0.000385   -0.021708  3.500000e-07   \n",
       "3  144.0  -0.001736   0.004588   0.000058    0.021417  4.375000e-07   \n",
       "\n",
       "   ANH3_rate  ANO3_rate  ANO2_rate  ANOrg_rate  AN2_rate  \n",
       "0   0.000000   0.000000   0.000000    0.000000  0.000000  \n",
       "1  -0.050772   0.032181   0.000227    0.000528  0.000022  \n",
       "2  -0.061275   0.010557  -0.000206    0.000374  0.000065  \n",
       "3  -0.027376  -0.001019   0.000081    0.000273  0.000380  "
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABbcAAAPXCAYAAAAYJXYaAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAB7CAAAewgFu0HU+AABvyUlEQVR4nOzdeXRV9b3w/08CBgnQhEGCAipiRcEJFRQnFL3VVh+rjw+1TnWmWG8drl1i22sd2uVQL7aVanvF61Br7bUtrVq9DlXEttiFXJwWIlgsM0EEAwYCRJLfH/485YTpHDgnyZe8Xmu5Vr4ne+/zybQNbzb7lDQ2NjYGAAAAAAAkpLSlBwAAAAAAgHyJ2wAAAAAAJEfcBgAAAAAgOeI2AAAAAADJEbcBAAAAAEiOuA0AAAAAQHLEbQAAAAAAkiNuAwAAAACQHHEbAAAAAIDkiNsAAAAAACRH3AYAAAAAIDniNgAAAAAAyRG3AQAAAABIjrgNAAAAAEByxG0AAAAAAJIjbgMAAAAAkBxxGwAAAACA5IjbAAAAAAAkR9wGAAAAACA57Vt6gLbmo48+aukRWlRJSUlUVlZGRERNTU00Nja27EBAq+V8AeTK+QLIh3MGkCvnCyi8rl27FvR4rtwGAAAAACA54jYAAAAAAMkRtwEAAAAASI64DQAAAABAcsRtAAAAAACSI24DAAAAAJAccRsAAAAAgOSI2wAAAAAAJEfcBgAAAAAgOeI2AAAAAADJEbcBAAAAAEiOuA0AAAAAQHLEbQAAAAAAkiNuAwAAAACQHHEbAAAAAIDkiNsAAAAAACRH3AYAAAAAIDniNgAAAAAAyRG3AQAAAABIjrgNAAAAAEByxG0AAAAAAJIjbgMAAAAAkBxxGwAAAACA5IjbAAAAAAAkR9wGAAAAACA54jYAAAAAAMlp35xP1tDQEPPmzYu5c+fGkiVLYuXKlbFu3booLy+PysrK2HfffePzn/98tGvXrjnHioiIVatWxdSpU2PJkiVRU1MT3bp1i969e8ehhx4aZWVlzT4PAAAAAACbV/S4vXz58viv//qvmDZtWsyYMSPq6uq2uH1FRUWcdtppcckll8Suu+5a7PFi6dKlcdddd8Wzzz4bq1ev3uj9lZWVcdppp8VVV10VnTt3Lvo8AAAAAABsXdFvS7Jw4cK4//77Y9q0aVsN2xERK1asiEceeSS+9KUvxYQJE4o62+TJk+O0006LCRMmbDJsR0TU1NTEL37xizjjjDPi3XffLeo8AAAAAADkpllvSxIR0aNHj9hnn31ijz32iIqKimjXrl3U1NTEjBkz4o033oiGhoaIiFi9enV8+9vfjvr6+jjrrLMKPseMGTPiiiuuyIraPXv2jGOPPTZ69OgRixcvjokTJ8bKlSsjImLevHlx2WWXxW9/+9uoqqoq+DwAAAAAAOSu6HG7Xbt2MWTIkDjppJPiqKOOir322muz2y5cuDBuueWWePnllzOP3XrrrTFs2LDYfffdCzbT2rVrNwrbF198cVxzzTVZ99eura2NG264IZ555pmIiPjggw/i6quvjscee6xgswAAAAAAkL+i35Zk4MCB8ctf/jLOP//8LYbtiIjevXvHvffeG0cddVTmsTVr1sSjjz5a0JkeffTRWLhwYWZ95plnxpgxYzZ64cjOnTvH2LFjY9iwYZnHpk2bFi+++GJB5wEAAAAAID9Fj9v5ateuXVx77bVZj/35z38u2PHr6+tj/PjxmXWXLl1izJgxm92+tLQ0br755igt/een6p577inYPAAAAAAA5K/Vxe2IiEGDBkV5eXlmvXjx4oIde8qUKbF8+fLM+tRTT42Kioot7rPHHnvEkUcemVlPnz495s+fX7CZAAAAAADIT6uM2xERnTp1yrzd2NhYsOO+9NJLWeuTTjopp/1OPvnkrLVbkwAAAAAAtJxWGbfXrFkTNTU1mXXfvn0LduypU6dm3m7Xrl0ceOCBOe03ePDgzR4HAAAAAIDm1Srj9rPPPhv19fWZ9fHHH1+Q4zY0NMScOXMy6z322CPrCvEt6d+/f+y8886Z9ezZswsyEwAAAAAA+Wt1cfu9996LO+64I7Pu2rVrXHDBBQU59sKFC2PNmjWZ9W677ZbzviUlJdGrV6/Mev78+fHJJ58UZC4AAAAAAPLTvqUHaGxsjNra2pg1a1Y8//zz8dhjj8XatWsjIqK8vDzGjRsX3bt3L8hzLVmyJGtdVVWV1/5VVVWZK7/r6+tj2bJleR8DAAAAAIDt1+xx+/33349TTz01s25oaNjkC0Yed9xx8e1vfzv23HPPgj33qlWrsta53pJkc9s3PV4uSkpK8t5nR7Lhx9/WPxfAljlfALlyvgDy4ZwB5Mr5Alq/Zo/bjY2NsX79+s2+v7S0NM4999y47LLLCn5VdF1dXda6Q4cOee3fdPvVq1fnPUNlZWXe++yoKioqWnoEIBHOF0CunC+AfDhnALlyvoDWqcVvS9JUQ0NDPPLII/HrX/86zj///LjmmmuirKysIMfe8H7bEZH3cZtu3/R4AAAAbOzQsbldWPS/167NabufHtExp+3+9W91W98oD32nXZHTdg+/mNt2tXNvyWm7037665y2aw71N/Ta+kYRsdP3q4s8CQC0QNzu379/zJw5M7Net25d1NTUxIwZM+LZZ5+Np556Kurr66O+vj4eeOCBmDVrVvzsZz8rSOBueuV1fX19XvuvW7dui8fLRU1NTd777EhKSkoyf9u5YsWKTd6SBiDC+QLInfMF7DgK/eelTR0vxXNGa/pzZK4392xNM8O2SvF8Aa1doe9q0eJXbpeVlUXPnj2jZ8+eMXz48Ljgggti9OjRsXjx4oiI+Mtf/hL33HNPXHPNNdv9XOXl5VnrfK+8/uyFLjd3vFw4Ef5TY2OjzweQE+cLIFfOF5C2Qv/8bu14qZwzUpixqRRnhi1J5XwBbU1pSw/Q1L777hvjx4+PnXbaKfPYQw89VJC/9W36gpD53jN7e1+QEgAAAACAwmh1cTsi4vOf/3x86UtfyqzXrFkTL7/88nYft+kLVFZX53cPsCVLlmTebt++fXTv3n27ZwIAAAAAIH+tMm5HRBx55JFZ6w3v072tevfunXWf7EWLFuW8b2NjY1YM79u3b9bV5QAAAAAANJ9WG7d79OiRta6trd3uY5aWlsaee+6ZWc+dOzfnW5PMnj076x7d/fv33+55AAAAAADYNq02bjeN2Z/73OcKctzDDjss8/b69evjzTffzGm/119/PWs9ZMiQgswDAAAAAED+Wm3cfuedd7LWu+66a0GOO2LEiKz1s88+m9N+Tbc74YQTCjIPAAAAAAD5a5Vxe82aNfHUU09lPdb0HtzbaujQodG1a9fM+umnn46VK1ducZ+5c+fG5MmTM+tBgwZF3759CzIPAAAAAAD5K2rcXrduXbz77rt57dPQ0BA33nhj1os9HnTQQbHXXnttcvsFCxbEgAEDMv81vTK7qbKysrj00ksz648//jhuv/32rc7T0NCQeewb3/hGrh8OAAAAAABFUNS4vWbNmjj99NPjyiuvjIkTJ8a6deu2uP2bb74ZX/va1+IPf/jDPwcsLY3vfve7BZ3rvPPOy7rNye9+97u44447NpqvtrY2rr322nj11Vczjw0ePDhOPPHEgs4DAAAAAEB+2hf7CRobG+O5556L5557Ljp27Bj77rtv7L333lFRUREdO3aMVatWRXV1dbz99tsxf/78rH1LSkriBz/4QRx00EEFnWnnnXeOe+65J84999yoq6uLiIgHHngg/vjHP8bw4cOje/fuUV1dHS+99FLWLUt22WWX+PGPf1zQWQAAAAAAyF/R4/aG6urq4vXXX4/XX399q9tWVVXFzTffHMcff3xRZhk0aFCMGzcuvvWtb0VNTU1ERHzwwQfxm9/8ZpPb9+nTJ8aNGxe9evUqyjwAAAAAAOSuqLcl6dSpU9xxxx1x6qmnRlVVVU77DBw4ML773e/GM888U7Sw/ZljjjkmnnrqqTj99NOjY8eOm9ymoqIizj///HjiiSdi4MCBRZ0HAAAAAIDcFPXK7Xbt2sXpp58ep59+ekR8emX07NmzY8GCBbFy5cpYs2ZNlJeXR+fOnaNPnz4xaNCg+NznPpfXc/Tp0ydmzpy5zTP27Nkz7rjjjvje974XU6dOjcWLF8eKFSuiW7du0bt37zjssMOirKxsm48PAAAAAEDhNettSXr27Bk9e/ZszqfMWadOnWL48OEtPQYAAAAAADko6m1JAAAAAACgGMRtAAAAAACSI24DAAAAAJAccRsAAAAAgOSI2wAAAAAAJEfcBgAAAAAgOeI2AAAAAADJEbcBAAAAAEiOuA0AAAAAQHLEbQAAAAAAkiNuAwAAAACQHHEbAAAAAIDkiNsAAAAAACRH3AYAAAAAIDniNgAAAAAAyRG3AQAAAABIjrgNAAAAAEByxG0AAAAAAJIjbgMAAAAAkBxxGwAAAACA5IjbAAAAAAAkR9wGAAAAACA54jYAAAAAAMkRtwEAAAAASI64DQAAAABAcsRtAAAAAACSI24DAAAAAJAccRsAAAAAgOSI2wAAAAAAJEfcBgAAAAAgOeI2AAAAAADJEbcBAAAAAEiOuA0AAAAAQHLEbQAAAAAAkiNuAwAAAACQHHEbAAAAAIDkiNsAAAAAACRH3AYAAAAAIDniNgAAAAAAyRG3AQAAAABIjrgNAAAAAEByxG0AAAAAAJIjbgMAAAAAkBxxGwAAAACA5IjbAAAAAAAkR9wGAAAAACA54jYAAAAAAMkRtwEAAAAASI64DQAAAABAcsRtAAAAAACSI24DAAAAAJAccRsAAAAAgOSI2wAAAAAAJEfcBgAAAAAgOeI2AAAAAADJEbcBAAAAAEiOuA0AAAAAQHLEbQAAAAAAkiNuAwAAAACQHHEbAAAAAIDkiNsAAAAAACRH3AYAAAAAIDniNgAAAAAAyRG3AQAAAABIjrgNAAAAAEByxG0AAAAAAJIjbgMAAAAAkBxxGwAAAACA5IjbAAAAAAAkR9wGAAAAACA54jYAAAAAAMkRtwEAAAAASI64DQAAAABAcsRtAAAAAACSI24DAAAAAJAccRsAAAAAgOSI2wAAAAAAJEfcBgAAAAAgOeI2AAAAAADJEbcBAAAAAEiOuA0AAAAAQHLEbQAAAAAAkiNuAwAAAACQHHEbAAAAAIDkiNsAAAAAACRH3AYAAAAAIDniNgAAAAAAyRG3AQAAAABIjrgNAAAAAEByxG0AAAAAAJIjbgMAAAAAkBxxGwAAAACA5IjbAAAAAAAkR9wGAAAAACA54jYAAAAAAMkRtwEAAAAASI64DQAAAABAcsRtAAAAAACSI24DAAAAAJAccRsAAAAAgOSI2wAAAAAAJEfcBgAAAAAgOeI2AAAAAADJEbcBAAAAAEiOuA0AAAAAQHLEbQAAAAAAkiNuAwAAAACQHHEbAAAAAIDkiNsAAAAAACRH3AYAAAAAIDntm/sJa2pqYtasWTF37tyoqamJxsbGqKioiN122y0OPvjg6NKlS3OPBAAAAABAYooetxsaGmLq1KnxwgsvxN/+9reYNWvWZrctKSmJYcOGxYUXXhjDhw8v2kznn39+TJkyZZv2ffDBB+PII48s8EQAAAAAAOSj6HH75JNPjrlz5+a0bWNjY0yePDkmT54cp5xyStxyyy3RuXPnIk8IAAAAAEBqih63ly9fvtFje+65Zxx44IHRo0eP6NChQ1RXV8err74a1dXVmW2efvrpWLp0adx///3RoUOHos7Yrl27nLctKSkp4iQAAAAAAOSi2e653bt37xg5cmScccYZ0atXr43ev379+nj88cfjtttui7Vr10ZExJQpU+LHP/5xjBkzpmhzDR06NB555JGiHR8AAAAAgMIretzebbfd4oILLojTTz99i1dIt2vXLs4+++zYbbfdYvTo0dHQ0BAREY888khceOGFUVVVVexRAQAAAABIRGmxn2DChAlx5pln5nzrj+HDh8cpp5ySWdfX18eLL75YrPEAAAAAAEhQ0eN2+/b5Xxy+YdyOiHj77bcLNQ4AAAAAADuAosftbbH77rtnrT/88MMWmgQAAAAAgNaoVcbtVatWZa235epvAAAAAAB2XK0ybs+cOTNr3atXrxaaBAAAAACA1qhVXhL95JNPZq2POOKIoj3XokWLYsyYMfH222/H0qVL45NPPomuXbtGVVVVHHbYYXHcccfFoYceWrTnBwAAAAAgf60ubk+ZMiWmTJmSWXfp0iWOPvrooj3fggULYsGCBVmPrV69OhYuXBjTpk2L++67LwYPHhzf/e5344ADDijaHAAAAAAA5K5Vxe3Vq1fHDTfckPXYRRddFJ06dWqhiT71+uuvx9lnnx033nhjjBw5cruOVVJSUqCp0rThx9/WPxfAljlfALlyvoAdR6F/hjd1vBTPGanMuaEUZ4amUjxfQFvTquL2TTfdFHPmzMms99prr7j00kuL8lyVlZUxYsSIOOaYY2LAgAFRVVUVO++8c6xcuTLee++9ePnll+Pxxx+P2traiIior6+P733ve9G1a9c48cQTt+t5+VRFRUVLjwAkwvkCyJXzBaSt0H9e2trxUjlntKY/R9bnuF1rmhkKIZXzBbQ1rSZuP/jgg/HEE09k1mVlZXHnnXdGhw4dCv5cV155ZRx44IGbPHa3bt3i8MMPj8MPPzwuueSSuOaaazK3SWloaIgxY8bEn/70p+jatWvB5wIAAAAAIDetIm4/88wz8cMf/jDrsVtuuSX233//ojzfkCFDctquR48ecd9998U555wT77zzTkRE1NbWxvjx4+O6667bpueuqanZpv12FCUlJZm/7VyxYkU0Nja28ERAa+V8AeTK+QJ2HIX+89KmjpfiOaM1/Tky15uGtqaZYVuleL6A1q7Q/7KnxeP25MmT47rrrouGhobMY9dee22cccYZLTjVP3Xs2DFuvPHGOOusszKPPfPMM9sct50I/6mxsdHnA8iJ8wWQK+cLSFuhf363drxUzhkpzNhUijPDlqRyvoC2prQln/zNN9+MK664Iurr/3nXrksuuSRGjRrVglNt7OCDD4699947s168eHHWvcEBAAAAAGheLRa3Z82aFaNGjYrVq1dnHhs5cuQ2XxFdbAcffHDWetGiRS0zCAAAAAAALRO3582bFxdffHHWPbi++MUvxi233NIS4+Ske/fuWeuPPvqohSYBAAAAAKDZ4/aSJUviwgsvjKVLl2YeGz58eNx5551RWtqid0nZorq6uqx1hw4dWmgSAAAAAACatSYvX748Lrzwwli4cGHmsaFDh8a4ceNip512as5R8jZ//vysdbdu3VpoEgAAAAAAmi1u19bWxqWXXhrvv/9+5rGDDjoofv7zn7f6q6Dr6uritddey6zbt28f++yzTwtOBAAAAADQtjVL3F6zZk2MHj06pk+fnnls3333jfHjx0enTp2aY4Tt8vDDD0dtbW1mfdBBB0Xnzp1bcCIAAAAAgLat6HH7k08+iauuuirryud+/frFAw88EBUVFdt9/AULFsSAAQMy/40YMWKL28+cOTOv4//lL3+Jn/70p1mPXXjhhfmOCQAAAABAARU1bjc2Nsb1118fL7/8cuaxPn36xMMPPxzdu3cv5lNv1nnnnRejRo2KF154IdauXbvZ7VauXBl33313fP3rX4/6+vrM40OHDo0vfOELzTEqAAAAAACb0b6YB1+0aFE89dRTGz12/PHH53Wc3r17xwsvvFCQmRobG2PSpEkxadKk6NixYwwYMCD23nvvqKysjJ133jlWrlwZs2fPjmnTpkVdXV3Wvv369Ytx48YVZA4AAAAAALZdUeN2Y2PjRo81NDTkfZz169cXYpyN1NXVxRtvvBFvvPHGVrc96aST4vvf/35BbqUCAAAAAMD2KWrcbo2+/vWvx1/+8pd46623YvXq1VvctqysLI455pg4//zzY9iwYc00IQAAAAAAW1PUuN2nT5+8X8Cx2M9x2WWXxWWXXRYNDQ0xd+7cmDt3blRXV8fHH38c9fX1UV5eHhUVFdGvX78YOHBglJWVFXF6AAAAAAC2RZu7cvszpaWl0a9fv+jXr19LjwIAAAAAQJ5KW3oAAAAAAADIl7gNAAAAAEByxG0AAAAAAJIjbgMAAAAAkBxxGwAAAACA5IjbAAAAAAAkR9wGAAAAACA54jYAAAAAAMkRtwEAAAAASI64DQAAAABAcsRtAAAAAACSI24DAAAAAJAccRsAAAAAgOSI2wAAAAAAJEfcBgAAAAAgOeI2AAAAAADJEbcBAAAAAEiOuA0AAAAAQHLEbQAAAAAAkiNuAwAAAACQHHEbAAAAAIDkiNsAAAAAACRH3AYAAAAAIDniNgAAAAAAyRG3AQAAAABIjrgNAAAAAEByxG0AAAAAAJIjbgMAAAAAkBxxGwAAAACA5IjbAAAAAAAkR9wGAAAAACA54jYAAAAAAMkRtwEAAAAASI64DQAAAABAcsRtAAAAAACSI24DAAAAAJAccRsAAAAAgOSI2wAAAAAAJEfcBgAAAAAgOeI2AAAAAADJEbcBAAAAAEiOuA0AAAAAQHLEbQAAAAAAkiNuAwAAAACQHHEbAAAAAIDkiNsAAAAAACRH3AYAAAAAIDniNgAAAAAAyRG3AQAAAABIjrgNAAAAAEByxG0AAAAAAJIjbgMAAAAAkBxxGwAAAACA5IjbAAAAAAAkR9wGAAAAACA54jYAAAAAAMkRtwEAAAAASI64DQAAAABAcsRtAAAAAACSI24DAAAAAJAccRsAAAAAgOSI2wAAAAAAJEfcBgAAAAAgOeI2AAAAAADJEbcBAAAAAEiOuA0AAAAAQHLEbQAAAAAAkiNuAwAAAACQHHEbAAAAAIDkiNsAAAAAACRH3AYAAAAAIDniNgAAAAAAyRG3AQAAAABIjrgNAAAAAEByxG0AAAAAAJIjbgMAAAAAkBxxGwAAAACA5IjbAAAAAAAkR9wGAAAAACA54jYAAAAAAMkRtwEAAAAASI64DQAAAABAcsRtAAAAAACSI24DAAAAAJAccRsAAAAAgOSI2wAAAAAAJEfcBgAAAAAgOeI2AAAAAADJEbcBAAAAAEiOuA0AAAAAQHLEbQAAAAAAkiNuAwAAAACQHHEbAAAAAIDkiNsAAAAAACRH3AYAAAAAIDniNgAAAAAAyRG3AQAAAABIjrgNAAAAAEByxG0AAAAAAJIjbgMAAAAAkBxxGwAAAACA5IjbAAAAAAAkR9wGAAAAACA54jYAAAAAAMkRtwEAAAAASI64DQAAAABAcsRtAAAAAACSI24DAAAAAJAccRsAAAAAgOSI2wAAAAAAJKd9cz9hTU1NzJo1K+bOnRs1NTXR2NgYFRUVsdtuu8XBBx8cXbp0ae6RIiJi1apVMXXq1FiyZEnU1NREt27donfv3nHooYdGWVlZi8wEAAAAAMCmFT1uNzQ0xNSpU+OFF16Iv/3tbzFr1qzNbltSUhLDhg2LCy+8MIYPH17s0SIiYunSpXHXXXfFs88+G6tXr97o/ZWVlXHaaafFVVddFZ07d26WmQAAAAAA2LKix+2TTz455s6dm9O2jY2NMXny5Jg8eXKccsopccsttxQ1KE+ePDmuvfbaWL58+Wa3qampiV/84hfx8ssvx7hx42Lfffct2jwAAAAAAOSm6HF7U+F4zz33jAMPPDB69OgRHTp0iOrq6nj11Vejuro6s83TTz8dS5cujfvvvz86dOhQ8LlmzJgRV1xxRdbV2j179oxjjz02evToEYsXL46JEyfGypUrIyJi3rx5cdlll8Vvf/vbqKqqKvg8AAAAAADkrtnuud27d+8YOXJknHHGGdGrV6+N3r9+/fp4/PHH47bbbou1a9dGRMSUKVPixz/+cYwZM6ags6xdu3ajsH3xxRfHNddck3V/7dra2rjhhhvimWeeiYiIDz74IK6++up47LHHCjoPAAAAAAD5KS32E+y2225x6623xgsvvBCXX375JsN2RES7du3i7LPPjnHjxkVp6T/HeuSRR2LJkiUFnenRRx+NhQsXZtZnnnlmjBkzZqMXjuzcuXOMHTs2hg0blnls2rRp8eKLLxZ0HgAAAAAA8lP0uD1hwoQ488wzo127djltP3z48DjllFMy6/r6+oLG5Pr6+hg/fnxm3aVLly1eGV5aWho333xzVnC/5557CjYPAAAAAAD5K3rcbt8+/zufbBi3IyLefvvtQo0TU6ZMyboP+KmnnhoVFRVb3GePPfaII488MrOePn16zJ8/v2AzAQAAAACQn6LH7W2x++67Z60//PDDgh37pZdeylqfdNJJOe138sknZ63dmgQAAAAAoOW0yri9atWqrPW2XP29OVOnTs283a5duzjwwANz2m/w4MGbPQ4AAAAAAM2rVcbtmTNnZq039yKU+WpoaIg5c+Zk1nvssUd06tQpp3379+8fO++8c2Y9e/bsgswEAAAAAED+WmXcfvLJJ7PWRxxxREGOu3DhwlizZk1mvdtuu+W8b0lJSVZknz9/fnzyyScFmQsAAAAAgPy0urg9ZcqUmDJlSmbdpUuXOProowty7CVLlmStq6qq8tp/w+3r6+tj2bJlBZkLAAAAAID8tKq4vXr16rjhhhuyHrvoootyvnXI1jS9l3e+x226fdPjAQAAAADQPAr3So0FcNNNN2XdE3uvvfaKSy+9tGDHr6ury1p36NAhr/2bbr969eq8ZygpKcl7nx3Jhh9/W/9cAFvmfAHkyvkCdhyF/hne1PFSPGekMueGUpwZmkrxfAFtTauJ2w8++GA88cQTmXVZWVnceeedeQfoLdnwftufPUc+mm7f9Hi5qKyszHufHVVFRUVLjwAkwvkCyJXzBaSt0H9e2trxUjlntKY/R9bnuF1rmhkKIZXzBbQ1rSJuP/PMM/HDH/4w67Fbbrkl9t9//4I+T9NQXl+f6/+WP7Vu3botHg8AAHY0gw7K7XVmpr/ZvciTAABAthaP25MnT47rrrsuGhoaMo9de+21ccYZZxT8ucrLy7PW+V55vXbt2i0eLxc1NTV577MjKSkpyfxt54oVK6KxsbGFJwJaK+cLIFfOF61DW/89l8Io9PfRpo6X4jmjNf185frKVa1pZthWKZ4voLUr9L/sadG4/eabb8YVV1yRdQX1JZdcEqNGjSrK8zV9Qch875m9vS9IGRFOhBtobGz0+QBy4nwB5Mr5ouX4vFMIhf4+2trxUjlnpDBjUynODFuSyvkC2prSlnriWbNmxahRo7IC88iRI+O6664r2nNWVVVlraurq/Paf8mSJZm327dvH927+6eXAAAAAAAtoUXi9rx58+Liiy/O+mdKX/ziF+OWW24p6vP27t076z7ZixYtynnfxsbGrBjet2/f2GmnnQo6HwAAAAAAuWn2uL1kyZK48MILY+nSpZnHhg8fHnfeeWeUlhZ3nNLS0thzzz0z67lz5+Z8a5LZs2dn3aO7f//+hR4PAAAAAIAcNWvcXr58eVx44YWxcOHCzGNDhw6NcePGNdtV0Icddljm7fXr18ebb76Z036vv/561nrIkCEFnQsAAAAAgNw1W9yura2NSy+9NN5///3MYwcddFD8/Oc/z7pVSLGNGDEia/3ss8/mtF/T7U444YSCzQQAAAAAQH6aJW6vWbMmRo8eHdOnT888tu+++8b48eOjU6dOzTFCxtChQ6Nr166Z9dNPPx0rV67c4j5z586NyZMnZ9aDBg2Kvn37Fm1GAAAAAAC2rOhx+5NPPomrrroqXnvttcxj/fr1iwceeCAqKiq2+/gLFiyIAQMGZP5remV2U2VlZXHppZdm1h9//HHcfvvtm92+oaEhbrzxxmhoaMg89o1vfGO75wYAAAAAYNsVNW43NjbG9ddfHy+//HLmsT59+sTDDz8c3bt3L+ZTb9F5550Xu+66a2b9u9/9Lu64445Yt25d1na1tbVx7bXXxquvvpp5bPDgwXHiiSc226wAAAAAAGysfTEPvmjRonjqqac2euz444/P6zi9e/eOF154oWBz7bzzznHPPffEueeeG3V1dRER8cADD8Qf//jHGD58eHTv3j2qq6vjpZdeyrplyS677BI//vGPCzYHAAAAAADbpqhxu7GxcaPHNry9R67Wr19fiHGyDBo0KMaNGxff+ta3oqamJiIiPvjgg/jNb36zye379OkT48aNi169ehV8FgAAAAAA8tMsLyjZWh1zzDHx1FNPxemnnx4dO3bc5DYVFRVx/vnnxxNPPBEDBw5s5gkBAAAAANiUol653adPn5g5c2Yxn2K7n6Nnz55xxx13xPe+972YOnVqLF68OFasWBHdunWL3r17x2GHHRZlZWUFnBgAAAAAgO1V1Lidkk6dOsXw4cNbegwAAAAAAHLQpm9LAgAAAABAmsRtAAAAAACSI24DAAAAAJAccRsAAAAAgOSI2wAAAAAAJEfcBgAAAAAgOeI2AAAAAADJEbcBAAAAAEiOuA0AAAAAQHLEbQAAAAAAkiNuAwAAAACQHHEbAAAAAIDkiNsAAAAAACRH3AYAAAAAIDniNgAAAAAAyRG3AQAAAABIjrgNAAAAAEByxG0AAAAAAJIjbgMAAAAAkBxxGwAAAACA5IjbAAAAAAAkR9wGAAAAACA54jYAAAAAAMkRtwEAAAAASI64DQAAAABAcsRtAAAAAACSI24DAAAAAJAccRsAAAAAgOSI2wAAAAAAJEfcBgAAAAAgOeI2AAAAAADJEbcBAAAAAEiOuA0AAAAAQHLEbQAAAAAAkiNuAwAAAACQHHEbAAAAAIDkiNsAAAAAACRH3AYAAAAAIDniNgAAAAAAyRG3AQAAAABIjrgNAAAAAEByxG0AAAAAAJIjbgMAAAAAkBxxGwAAAACA5IjbAAAAAAAkR9wGAAAAACA54jYAAAAAAMkRtwEAAAAASI64DQAAAABAcsRtAAAAAACSI24DAAAAAJAccRsAAAAAgOSI2wAAAAAAJEfcBgAAAAAgOeI2AAAAAADJEbcBAAAAAEiOuA0AAAAAQHLEbQAAAAAAkiNuAwAAAACQHHEbAAAAAIDkiNsAAAAAACRH3AYAAAAAIDniNgAAAAAAyRG3AQAAAABIjrgNAAAAAEByxG0AAAAAAJIjbgMAAAAAkBxxGwAAAACA5IjbAAAAAAAkR9wGAAAAACA54jYAAAAAAMkRtwEAAAAASI64DQAAAABAcsRtAAAAAACSI24DAAAAAJAccRsAAAAAgOSI2wAAAAAAJEfcBgAAAAAgOeI2AAAAAADJEbcBAAAAAEiOuA0AAAAAQHLEbQAAAAAAkiNuAwAAAACQHHEbAAAAAIDkiNsAAAAAACRH3AYAAAAAIDniNgAAAAAAyRG3AQAAAABIjrgNAAAAAEByxG0AAAAAAJIjbgMAAAAAkBxxGwAAAACA5IjbAAAAAAAkR9wGAAAAACA54jYAAAAAAMkRtwEAAAAASI64DQAAAABAcsRtAAAAAACSI24DAAAAAJAccRsAAAAAgOSI2wAAAAAAJEfcBgAAAAAgOeI2AAAAAADJEbcBAAAAAEiOuA0AAAAAQHLEbQAAAAAAkiNuAwAAAACQHHEbAAAAAIDkiNsAAAAAACRH3AYAAAAAIDniNgAAAAAAyRG3AQAAAABIjrgNAAAAAEBy2jfXE61atSreeeedeOutt+Ktt96Kt99+OxYuXJh5f+/eveOll14q+hzXX399/P73v9+mfX/wgx/EyJEjCzwRAAAAAAD5KnrcfvDBB2PChAnx97//PRoaGor9dAAAAAAAtAFFj9uvvfZazJo1q9hPs83atWuX87YlJSVFnAQAAAAAgFw1221JNlReXh6DBg2K6dOnx+rVq1tihIhovluhAAAAAABQWEWP2x06dIgDDzwwDjjggDjggANi//33j/79+0dpaWmMGDGiReM2AAAAAABpKnrc/tGPflTspwAAAAAAoI0pbekBAAAAAAAgX+I2AAAAAADJEbcBAAAAAEiOuA0AAAAAQHKK/oKSrVltbW1873vfizfeeCOWLFkSa9eujYqKiujZs2cMHjw4jjnmmDj66KOjpKSkpUcFAAAAAGADbTpur1ixIv77v/8767G6urqorq6Ot956Kx5++OHYZ599YsyYMXH00UcX5Dnbeijf8ONv658LYMucL4BcOV+0Dj73FEKhv482dbwUzxmpzLmhFGeGplI8X0Bb06bjdi5mzZoVl156afzrv/5r/Ou//ut2H6+ysnL7h9pBVFRUtPQIQCKcL4BcOV8Uw7KctvJ7LoVQ6O+jrR0vlXNGa/r5qs9xu9Y0MxRCKucLaGvaZNwuLy+P4447Lo499tjYb7/9onfv3tGxY8f4+OOPY+7cufGXv/wlfvWrX8WyZZ/+It/Y2Bjjxo2LLl26xAUXXNDC0wMAAAAA0Obi9le/+tX493//9+jcufNG7+vatWt07do1Dj744LjoooviO9/5Tjz77LOZ999xxx1x9NFHR//+/bf5+WtqarZ53x1BSUlJ5m87V6xYEY2NjS08EdBaOV8AuXK+aB3a+u+5FEahv482dbwUzxmt6eerU47btaaZYVuleL6A1q7Q/7KnzcXtgw8+OKftOnXqFD/60Y+irq4uJk2aFBER69evj7vvvjt+8pOfbPPzOxH+U2Njo88HkBPnCyBXzhctx+edQij099HWjpfKOSOFGZtKcWbYklTOF9DWlLb0AK1ZaWlp3HjjjbHTTjtlHnvppZeirq6uBacCAAAAAEDc3orevXvHkUcemVmvW7cupk2b1oITAQAAAAAgbueg6a1MFi1a1DKDAAAAAAAQEeJ2Trp37561/uijj1poEgAAAAAAIsTtnDS9x3ZZWVkLTQIAAAAAQIS4nZMFCxZkrZteyQ0AAAAAQPMSt7eisbEx/vznP2c9tt9++7XQNAAAAAAARIjbW/XHP/4x5syZk1nvuuuusffee7fcQAAAAAAARPuWHmB7DRgwIGs9c+bMzW47c+bMjbbfkhkzZsRNN92U9dgFF1yQ13wAAAAAABRem7py+7rrrouvfvWr8cQTT0Rtbe1mt1uzZk08/PDDcfbZZ2dt169fvzj33HObY1QAAAAAALag6FduL1y4MP7lX/5lk+9bv3591nYDBw7c5HYPPfRQDB06tCDzvP766/H666/HTjvtFPvss0/ss88+0bVr1ygvL49Vq1bFnDlzYurUqfHxxx9n7dejR48YP358lJWVFWQOAAAAAAC2XdHjdmNjY1bE3pLNbdfY2FjIkSIior6+PqZPnx7Tp0/f6rZDhw6NH/7wh7HrrrsWfA4AAAAAAPKX/D2383HeeefFCy+8EK+//nqsXLlyi9u2b98+hgwZEuedd16MGDEiSkvb1B1cAAAAAABataLH7T59+mzxRR63Vz7HHjlyZIwcOTIiIhYsWBD/+Mc/orq6OlasWBFr166N8vLy+NznPhd9+/aNAw44IDp27FissQEAAAAA2A5t6srtDfXp0yf69OnT0mMAAAAAALAN3GsDAAAAAIDkiNsAAAAAACRH3AYAAAAAIDniNgAAAAAAyRG3AQAAAABIjrgNAAAAAEByxG0AAAAAAJIjbgMAAAAAkBxxGwAAAACA5IjbAAAAAAAkR9wGAAAAACA54jYAAAAAAMkRtwEAAAAASI64DQAAAABAcsRtAAAAAACSI24DAAAAAJAccRsAAAAAgOSI2wAAAAAAJEfcBgAAAAAgOeI2AAAAAADJEbcBAAAAAEiOuA0AAAAAQHLEbQAAAAAAkiNuAwAAAACQHHEbAAAAAIDkiNsAAAAAACRH3AYAAAAAIDniNgAAAAAAyRG3AQAAAABIjrgNAAAAAEByxG0AAAAAAJIjbgMAAAAAkBxxGwAAAACA5IjbAAAAAAAkR9wGAAAAACA54jYAAAAAAMkRtwEAAAAASI64DQAAAABAcsRtAAAAAACSI24DAAAAAJAccRsAAAAAgOSI2wAAAAAAJEfcBgAAAAAgOeI2AAAAAADJEbcBAAAAAEiOuA0AAAAAQHLEbQAAAAAAkiNuAwAAAACQHHEbAAAAAIDkiNsAAAAAACRH3AYAAAAAIDniNgAAAAAAyRG3AQAAAABIjrgNAAAAAEByxG0AAAAAAJIjbgMAAAAAkBxxGwAAAACA5IjbAAAAAAAkR9wGAAAAACA54jYAAAAAAMkRtwEAAAAASI64DQAAAABAcsRtAAAAAACSI24DAAAAAJAccRsAAAAAgOSI2wAAAAAAJEfcBgAAAAAgOeI2AAAAAADJEbcBAAAAAEiOuA0AAAAAQHLEbQAAAAAAkiNuAwAAAACQHHEbAAAAAIDkiNsAAAAAACRH3AYAAAAAIDniNgAAAAAAyRG3AQAAAABIjrgNAAAAAEByxG0AAAAAAJIjbgMAAAAAkBxxGwAAAACA5IjbAAAAAAAkR9wGAAAAACA54jYAAAAAAMkRtwEAAAAASI64DQAAAABAcsRtAAAAAACSI24DAAAAAJAccRsAAAAAgOSI2wAAAAAAJEfcBgAAAAAgOeI2AAAAAADJEbcBAAAAAEiOuA0AAAAAQHLEbQAAAAAAkiNuAwAAAACQHHEbAAAAAIDkiNsAAAAAACRH3AYAAAAAIDniNgAAAAAAyRG3AQAAAABIjrgNAAAAAEByxG0AAAAAAJIjbgMAAAAAkBxxGwAAAACA5IjbAAAAAAAkR9wGAAAAACA54jYAAAAAAMkRtwEAAAAASE77lh6gNVi3bl387//+byxcuDCWL18elZWV0atXrzjssMOivLy8pccDAAAAAKCJZovbq1atinfeeSfeeuuteOutt+Ltt9+OhQsXZt7fu3fveOmll5prnIiIqK2tjZ/85Cfx5JNPRk1NzUbvLy8vj5NPPjmuvfba6NGjR7POBgAAAADA5hU9bj/44IMxYcKE+Pvf/x4NDQ3FfrqcvfPOO/HNb34zFixYsNltVq9eHRMmTIhJkybF2LFjY9iwYc04IQAAAAAAm1P0uP3aa6/FrFmziv00eVm8eHGMGjUqli5dmnmsoqIijj/++OjVq1csXbo0Xnnllcz7ly1bFt/4xjfi17/+dQwYMKClxgYAAAAA4P/XIvfcLi8vj0GDBsX06dNj9erVzfrcjY2NcdVVV2WF7VNPPTVuueWW6NSpU+axdevWxdixY+Ohhx6KiE+v4r7iiivimWeeibKysmadGQAAAACAbEWP2x06dIgDDzwwDjjggDjggANi//33j/79+0dpaWmMGDGi2eP2888/H2+++WZmfeSRR8Z//Md/RElJSdZ2ZWVl8e1vfztWrlwZEyZMiIiI+fPnx69//ev42te+1qwzAwAAAACQrehx+0c/+lGxnyIvP/vZzzJvl5aWxk033bRR2N7Q9ddfH88//3zU1tZGRMR9990X55xzTrRv3yIXvQMAAAAAEBGlLT1Ac5o7d27MmDEjsz766KNjjz322OI+FRUVccopp2TWS5cujalTpxZtRgAAAAAAtq5Nxe0XX3wxa33SSSfltN8Xv/jFLR4HAAAAAIDm1abidtMrrgcPHpzTfgcccEC0a9dus8cBAAAAAKB5tam4PXv27MzbO++8c+y111457de5c+fYfffdM+t//OMf0djYWPD5AAAAAADITZuJ2/X19bFgwYLMulevXlt8Icmmdtttt8zbdXV1sWjRooLOBwAAAABA7tpM3F62bFl88sknmXWvXr3y2r+qqiprXV1dXZC5AAAAAADIX/uWHqC5rFq1KmvdqVOnvPZvuv3q1au3aY58rhbfEW348bf1zwWwZc4XQK6cL1oHn3sKodDfR5s6XornjFTm3FCKM0NTKZ4voK1pM3G7aYzu0KFDXvs33X5b43ZlZeU27bcjqqioaOkRgEQ4XwC5cr4ohmU5beX3XAqh0N9HWzteKueM1vTzVZ/jdq1pZiiEVM4X0Na0mbi9du3arHVZWVle+zfdfs2aNds9E2kZtag2p+3u261zkSehOfm6w47Pz3kaXunyk5y2O/bjq4o8Cez4Bh2U219oTH+ze8GfO9ef9RhU2OddPG9ojlsOKewT02Jy/ZrvuvuUIk/S9vh/Os3Jz/qOr83E7aZxur4+179v/tS6deuy1vle+f2ZmpqabdpvR1FSUpL5284VK1ZEY2NjC0+Uj9x+XNr613jH4+veUtI+X5AWP+c7EueLluNnpO1J4Wu+qRlT/B2jNX2uc73BZ2uaeXvsKB9HilrD5z7F8wXbpjV8v7UVhf6XPW0mbje9Z3a+V143vfK7vLx8m+ZwIvynxsbGHfLzsSN+TGydr3tx7ajnC9LiezANzhctx+e97Unha761GVM5Z6QwY1MpzrwpO8rHkaLW9rlP5XzBtvG1TVdpSw/QXJrG6Hzvmd30BSm3NW4DAAAAALD92kzc7t69e7Rv/88L1RcvXpzX/kuWLMla9+rVqyBzAQAAAACQvzYTt8vKyqJPnz6ZdXV1dV7/5GDRokWZt3feeefo3bt3QecDAAAAACB3bSZuR0TstddembfXrFkT77//fk771dbWxrx58zLrfv36RUlJScHnAwAAAAAgN20qbh922GFZ69dffz2n/d56661Yv379Zo8DAAAAAEDzalNxe8SIEVnrZ599Nqf9mm53wgknFGwmAAAAAADy16bidr9+/WLAgAGZ9V//+teYO3fuFvdZsWJFPPPMM5l1jx49YsiQIUWbEQAAAACArUs+bg8YMCDrv625/PLLM283NDTETTfdtMUXlrz99tvj448/zqwvu+yyaN++/fYNDQAAAADAdkk+bufr5JNPjgMOOCCznjx5cnzrW9+KVatWZW23bt26uO2222LChAmZx3r37h1nn312s80KAAAAAMCmFf0S5IULF8a//Mu/bPJ9G75I48KFC2PgwIGb3O6hhx6KoUOHFmSekpKS+MlPfhIjR46MZcuWRUTEH//4x/jzn/8cI0aMiKqqqvjwww9j0qRJsXTp0sx+5eXlce+990aHDh0KMgcAAAAAANuu6HG7sbExK2Jvyea229JtQ7ZF796947777osrr7wyFi5cGBGf3lv797///Sa379q1a4wdOzb23Xffgs4BAAAAAMC2aXO3JfnM/vvvH0888UScf/75UVFRscltOnbsGKeffno89dRTcdRRRzXzhAAAAAAAbE7Rr9zu06dPzJw5s2jH355jd+nSJf793/89rrvuupg6dWosXLgwli9fHhUVFbHrrrvGYYcdFp06dSrgtAAAAAAAFELR43YKysrK4sgjj2zpMQAAAAAAyFGbvS0JAAAAAADpErcBAAAAAEiOuA0AAAAAQHLEbQAAAAAAkiNuAwAAAACQHHEbAAAAAIDkiNsAAAAAACRH3AYAAAAAIDniNgAAAAAAyRG3AQAAAABIjrgNAAAAAEByxG0AAAAAAJIjbgMAAAAAkBxxGwAAAACA5IjbAAAAAAAkR9wGAAAAACA54jYAAAAAAMkRtwEAAAAASI64DQAAAABAcsRtAAAAAACSI24DAAAAAJAccRsAAAAAgOSI2wAAAAAAJEfcBgAAAAAgOeI2AAAAAADJEbcBAAAAAEiOuA0AAAAAQHLEbQAAAAAAkiNuAwAAAACQHHEbAAAAAIDkiNsAAAAAACRH3AYAAAAAIDniNgAAAAAAyRG3AQAAAABIjrgNAAAAAEByxG0AAAAAAJIjbgMAAAAAkBxxGwAAAACA5IjbAAAAAAAkR9wGAAAAACA54jYAAAAAAMkRtwEAAAAASI64DQAAAABAcsRtAAAAAACSI24DAAAAAJAccRsAAAAAgOSI2wAAAAAAJEfcBgAAAAAgOeI2AAAAAADJEbcBAAAAAEiOuA0AAAAAQHLEbQAAAAAAkiNuAwAAAACQHHEbAAAAAIDkiNsAAAAAACRH3AYAAAAAIDniNgAAAAAAyRG3AQAAAABIjrgNAAAAAEByxG0AAAAAAJIjbgMAAAAAkBxxGwAAAACA5IjbAAAAAAAkR9wGAAAAACA54jYAAAAAAMkRtwEAAAAASI64DQAAAABAcsRtAAAAAACSI24DAAAAAJAccRsAAAAAgOSI2wAAAAAAJEfcBgAAAAAgOeI2AAAAAADJEbcBAAAAAEiOuA0AAAAAQHLEbQAAAAAAkiNuAwAAAACQHHEbAAAAAIDkiNsAAAAAACRH3AYAAAAAIDniNgAAAAAAyRG3AQAAAABIjrgNAAAAAEByxG0AAAAAAJIjbgMAAAAAkBxxGwAAAACA5IjbAAAAAAAkR9wGAAAAACA54jYAAAAAAMkRtwEAAAAASI64DQAAAABAcsRtAAAAAACSI24DAAAAAJAccRsAAAAAgOSI2wAAAAAAJEfcBgAAAAAgOeI2AAAAAADJEbcBAAAAAEiOuA0AAAAAQHLEbQAAAAAAkiNuAwAAAACQHHEbAAAAAIDkiNsAAAAAACRH3AYAAAAAIDniNgAAAAAAyRG3AQAAAABIjrgNAAAAAEByxG0AAAAAAJIjbgMAAAAAkBxxGwAAAACA5LRvqSeuqamJadOmRXV1ddTW1kbPnj2jT58+ccghh0RpqeYOAAAAAMDmNXvcnjNnTowdOzYmTpwY9fX1G72/Z8+ecdZZZ8WoUaOirKysKDOcf/75MWXKlG3a98EHH4wjjzyywBMBAAAAAJCPZr1E+sknn4wzzjgjnn/++U2G7YiIDz74IMaNGxdf/epXY+HChc05HgAAAAAAiWi2K7dfeeWVuP7662P9+vWZx/bcc884/PDDo7KyMubNmxcTJ06MNWvWRETE9OnTY/To0fHYY49F586dizpbu3btct62pKSkiJMAAAAAAJCLZonbS5cujX/7t3/LhO2SkpIYM2ZMXHDBBVn3116+fHlcddVVmVuGzJo1K2688cYYO3Zs0WYbOnRoPPLII0U7PgAAAAAAhdcstyX5+c9/Hh9//HFm/c1vfjMuuuiijV44slu3bnH//fdH//79M489/fTT8e677zbHmAAAAAAAJKLocXvZsmXx+OOPZ9a77757jBo1arPbd+jQIW644YbMurGxMe69996izggAAAAAQFqKHrdffPHFWLduXWb9la98JXbaaact7jNs2LDo169fZj1p0qSoq6sr2owAAAAAAKSl6HH7pZdeylqffPLJOe234XZr1qyJv/71rwWdCwAAAACAdBU9bk+dOjXzdo8ePaJv37457Td48OCs9WuvvVbQuQAAAAAASFdR4/YHH3yQ9UKS++23X877Dhw4MGs9e/bsgs0FAAAAAEDa2hfz4O+//37Werfddst53x49esROO+0U9fX1mzxWoSxatCjGjBkTb7/9dixdujQ++eST6Nq1a1RVVcVhhx0Wxx13XBx66KFFeW4AAAAAALZNUeP2kiVLstZVVVU571tSUhJVVVWxYMGCTR6rUBYsWJB5js+sXr06Fi5cGNOmTYv77rsvBg8eHN/97nfjgAMOKMoMAAAAAADkp6hxe9WqVVnrTp065bX/htt/8sknsW7duigrKyvIbPl4/fXX4+yzz44bb7wxRo4cuV3HKikpKdBUadrw499RPxc76sfFlvm6F15bOF+QFt+HafB1ajk+921PCl/zTc2Y4u8Yqcy5oRRn3pQd5eNIUWv43Kd4vmDb+Pqmq6hxu66uLmvdoUOHvPZvuv2qVasKFrcrKytjxIgRccwxx8SAAQOiqqoqdt5551i5cmW899578fLLL8fjjz8etbW1ERFRX18f3/ve96Jr165x4oknbtfz8qmKioqWHiE/q2tz2szXeAfj694qJHe+IC1+zncozhfFsCynrfyM7Eh2nK/51mb87JxRt7IZhtkOrelzXZ/jdq1p5k3J9Wve2j+OHVlr+9z7HSNNftZ3fEWN22vWrMla5xumm26/du3a7Z4pIuLKK6+MAw88cJOxvVu3bnH44YfH4YcfHpdccklcc801MWXKlIiIaGhoiDFjxsSf/vSn6Nq1a0FmAQAAAAAgf0WN203j8WcvDpmrdevWZa0LddX2kCFDctquR48ecd9998U555wT77zzTkRE1NbWxvjx4+O6667bpueuqanZpv12FCUlJZm/7VyxYkU0Nja28ET5yO3Hpa1/jXc8vu4tJe3zBWnxc74jcb5oOX5G2p4UvuabmjHF3zFa0+c615uNtqaZt8eO8nGkqDV87lM8X7BtWsP3W1tR6Kvkixq3y8vLs9ZNr+TemqZXaud7z+5C6NixY9x4441x1llnZR575plntjluOxH+U2Nj4w75+dgRPya2zte9uHbU8wVp8T2YBueLluPz3vak8DXf2oypnDNSmLGpFGfelB3l40hRa/vcp3K+YNv42qartJgHbxq3V69endf+G74gZfv27fO+Z3ehHHzwwbH33ntn1osXL445c+a0yCwAAAAAABQ5bldVVWWtq6urc963sbExlixZstljNbeDDz44a71o0aKWGQQAAAAAgOLG7b322itrnU8Q/vDDD7Pu0d2vX7+CzbUtunfvnrX+6KOPWmgSAAAAAACKfuV2ly5dMusZM2bkvO9nL+D4mf79+xdsrm1RV1eXtW6pW6QAAAAAAFDkuB0Rceihh2be/vDDD2P+/Pk57Tdt2rSs9ZAhQwo6V76azt2tW7cWmgQAAAAAgKLH7REjRmSt/+d//ien/Z577rnM2x06dIijjjqqoHPlo66uLl577bXMun379rHPPvu02DwAAAAAAG1d0eP2CSecEDvttFNm/Zvf/CbrXtqb8uqrr8Y//vGPzHr48OFRXl5etBm35uGHH47a2trM+qCDDorOnTu32DwAAAAAAG1d0eN2jx49YuTIkZn1vHnz4r777tvs9mvXro0f/OAHmXVJSUlcfvnlm91+wYIFMWDAgMx/Ta8Ub2rmzJl5TB/xl7/8JX76059mPXbhhRfmdQwAAAAAAAqr6HE7ImL06NHRqVOnzHrcuHHx0EMPRUNDQ9Z2y5cvj0svvTT+/ve/Zx770pe+FAMHDizYLOedd16MGjUqXnjhhVi7du1mt1u5cmXcfffd8fWvfz3rSvOhQ4fGF77whYLNAwAAAABA/to3x5NUVVXFXXfdFZdffnk0NDREY2Nj3HbbbfHYY4/FEUccEZWVlTF37tyYOHFirFmzJrPf3nvvHbfccktBZ2lsbIxJkybFpEmTomPHjjFgwIDYe++9o7KyMnbeeedYuXJlzJ49O6ZNmxZ1dXVZ+/br1y/GjRtX0HkAAAAAAMhfs8TtiIjjjjsubrvttrjpppsy0XjOnDkxZ86cTW6/3377xU9/+tOi3tu6rq4u3njjjXjjjTe2uu1JJ50U3//+96OioqJo8wAAAAAAkJtmuS3JZ04//fSYMGFCnHjiiVkvMrmhXXbZJa644op4/PHHo0+fPgWf4etf/3occcQROb1AZVlZWZxwwgnx0EMPxd133y1sAwAAAAC0Es125fZn9tprr7jnnnvio48+imnTpkV1dXWsWrUqevToEX379o1DDjkk2rVrl/Px+vTpk9eLRF522WVx2WWXRUNDQ8ydOzfmzp0b1dXV8fHHH0d9fX2Ul5dHRUVF9OvXLwYOHBhlZWXb8mECAAAAAFBEzR63P9O1a9c44YQTWurpo7S0NPr16xf9+vVrsRkAAAAAANg2zXpbEgAAAAAAKARxGwAAAACA5IjbAAAAAAAkR9wGAAAAACA54jYAAAAAAMkRtwEAAAAASI64DQAAAABAcsRtAAAAAACSI24DAAAAAJAccRsAAAAAgOSI2wAAAAAAJEfcBgAAAAAgOeI2AAAAAADJEbcBAAAAAEiOuA0AAAAAQHLEbQAAAAAAkiNuAwAAAACQHHEbAAAAAIDkiNsAAAAAACRH3AYAAAAAIDniNgAAAAAAyRG3AQAAAABIjrgNAAAAAEByxG0AAAAAAJIjbgMAAAAAkBxxGwAAAACA5IjbAAAAAAAkR9wGAAAAACA54jYAAAAAAMkRtwEAAAAASI64DQAAAABAcsRtAAAAAACSI24DAAAAAJAccRsAAAAAgOSI2wAAAAAAJEfcBgAAAAAgOeI2AAAAAADJEbcBAAAAAEiOuA0AAAAAQHLEbQAAAAAAkiNuAwAAAACQHHEbAAAAAIDkiNsAAAAAACRH3AYAAAAAIDniNgAAAAAAyRG3AQAAAABIjrgNAAAAAEByxG0AAAAAAJIjbgMAAAAAkBxxGwAAAACA5IjbAAAAAAAkR9wGAAAAACA54jYAAAAAAMkRtwEAAAAASI64DQAAAABAcsRtAAAAAACSI24DAAAAAJAccRsAAAAAgOSI2wAAAAAAJEfcBgAAAAAgOeI2AAAAAADJEbcBAAAAAEiOuA0AAAAAQHLEbQAAAAAAkiNuAwAAAACQHHEbAAAAAIDkiNsAAAAAACRH3AYAAAAAIDniNgAAAAAAyRG3AQAAAABIjrgNAAAAAEByxG0AAAAAAJIjbgMAAAAAkBxxGwAAAACA5IjbAAAAAAAkR9wGAAAAACA54jYAAAAAAMkRtwEAAAAASI64DQAAAABAcsRtAAAAAACSI24DAAAAAJAccRsAAAAAgOSI2wAAAAAAJEfcBgAAAAAgOeI2AAAAAADJEbcBAAAAAEiOuA0AAAAAQHLEbQAAAAAAkiNuAwAAAACQHHEbAAAAAIDkiNsAAAAAACRH3AYAAAAAIDniNgAAAAAAyRG3AQAAAABIjrgNAAAAAEByxG0AAAAAAJIjbgMAAAAAkBxxGwAAAACA5IjbAAAAAAAkR9wGAAAAACA54jYAAAAAAMkRtwEAAAAASI64DQAAAABAcsRtAAAAAACSI24DAAAAAJAccRsAAAAAgOSI2wAAAAAAJEfcBgAAAAAgOeI2AAAAAADJEbcBAAAAAEhO+5Z64pqampg2bVpUV1dHbW1t9OzZM/r06ROHHHJIlJY2f3NftWpVTJ06NZYsWRI1NTXRrVu36N27dxx66KFRVlbW7PMAAAAAALB5zR6358yZE2PHjo2JEydGfX39Ru/v2bNnnHXWWTFq1KhmicpLly6Nu+66K5599tlYvXr1Ru+vrKyM0047La666qro3Llz0ecBAAAAAGDrmvUS6SeffDLOOOOMeP755zcZtiMiPvjggxg3blx89atfjYULFxZ1nsmTJ8dpp50WEyZM2GTYjvj0CvNf/OIXccYZZ8S7775b1HkAAAAAAMhNs125/corr8T1118f69evzzy25557xuGHHx6VlZUxb968mDhxYqxZsyYiIqZPnx6jR4+Oxx57rChXTM+YMSOuuOKKrKjds2fPOPbYY6NHjx6xePHimDhxYqxcuTIiIubNmxeXXXZZ/Pa3v42qqqqCzwMAAAAAQO6aJW4vXbo0/u3f/i0TtktKSmLMmDFxwQUXZN1fe/ny5XHVVVfFlClTIiJi1qxZceONN8bYsWMLOs/atWs3CtsXX3xxXHPNNVm3QqmtrY0bbrghnnnmmYj49Kryq6++Oh577LGCzgMAAAAAQH6a5bYkP//5z+Pjjz/OrL/5zW/GRRddtNELR3br1i3uv//+6N+/f+axp59+uuC3A3n00Uezbnly5plnxpgxYza6x3fnzp1j7NixMWzYsMxj06ZNixdffLGg8wAAAAAAkJ+ix+1ly5bF448/nlnvvvvuMWrUqM1u36FDh7jhhhsy68bGxrj33nsLNk99fX2MHz8+s+7SpUuMGTNms9uXlpbGzTffnBXi77nnnoLNAwAAAABA/ooet1988cVYt25dZv2Vr3wldtpppy3uM2zYsOjXr19mPWnSpKirqyvIPFOmTInly5dn1qeeempUVFRscZ899tgjjjzyyMx6+vTpMX/+/ILMAwAAAABA/ooet1966aWs9cknn5zTfhtut2bNmvjrX/9alHlOOumkvOeJCLcmAQAAAABoQUWP21OnTs283aNHj+jbt29O+w0ePDhr/dprrxV8nnbt2sWBBx64TfNseBwAAAAAAJpXUeP2Bx98kPVCkvvtt1/O+w4cODBrPXv27O2ep6GhIebMmZNZ77HHHtGpU6ec9u3fv3/svPPOBZ0HAAAAAIBtU9S4/f7772etd9ttt5z37dGjR9a9uZsea1ssXLgw1qxZs03zlJSURK9evTLr+fPnxyeffLLdMwEAAAAAkL+ixu0lS5ZkrauqqnLet6SkJGv7psdq7nmabl9fXx/Lli3b7pkAAAAAAMhf+2IefNWqVVnrXG8BsqntP/nkk1i3bl2UlZW1ink2dbxclJSU5L3PjmTDj39H/VzsqB8XW+brXnht4XxBWnwfpsHXqeX43Lc9KXzNNzVjir9jpDLnhlKceVN2lI8jRa3hc5/i+YJt4+ubrqLG7bq6uqx1hw4d8tq/6farVq3arrhd6HlWr16d9wyVlZV577OjqqioaOkR8rO6NqfNfI13ML7urUJy5wvS4ud8h+J8UQy5/WtFPyM7kh3na761GT87Z9StbIZhtkNr+lzX57hda5p5U3L9mrf2j2NH1to+937HSJOf9R1fUeP2hve3joi8w3TT7deuXduq5ml6PHZs9+3WuaVHoAX4usOOz895Go79+KqWHqHNmv5m95YegWZWjK/5/16b45/lrs1ts2OjsOeEXXefktN283fP8YCH5PrMv851w1Zjp+9Xt/QIBZHr15zC8/90mpOf9R1fUe+53fRK5/r6XP+O91Pr1q3LWm/PVdvFmCffK78BAAAAACiMosbt8vLyrHW+Vzo3vVI733tkF3uepscDAAAAAKB5NGvczvce1Ru+YGP79u23+0rppnF8e+bZ1PEAAAAAAGgeRY3bVVVVWevq6tzvzdXY2BhLlizZ7LGae56IyJqnffv20b27+w8CAAAAALSEosbtvfbaK2u9aNGinPf98MMPs+6J3a9fv+2ep3fv3llXf+czT2NjY1YM79u3b+y0007bPRMAAAAAAPkr+pXbXbp0yaxnzJiR877vvPNO1rp///7bPU9paWnsueeemfXcuXNzvjXJ7Nmzs+7RXYh5AAAAAADYNkWN2xERhx56aObtDz/8MObPn5/TftOmTctaDxkypCDzHHbYYZm3169fH2+++WZO+73++utFmQcAAAAAgPwVPW6PGDEia/0///M/Oe333HPPZd7u0KFDHHXUUUWZ59lnn81pv6bbnXDCCQWZBwAAAACA/BU9bp9wwglZ96b+zW9+k3Uv7U159dVX4x//+EdmPXz48CgvLy/IPEOHDo2uXbtm1k8//XSsXLlyi/vMnTs3Jk+enFkPGjQo+vbtW5B5AAAAAADIX/tiP0GPHj1i5MiR8atf/SoiIubNmxf33XdfXHHFFZvcfu3atfGDH/wgsy4pKYnLL798s8dfsGBB1lXUvXv3jpdeemmz25eVlcWll14ad955Z0REfPzxx3H77bfHrbfeusntGxoa4sYbb4yGhobMY9/4xjc2e/wd2axZs+K9996LJUuWRGlpafTq1SsGDhwYu+++e4vNNG/evHjnnXeiuro6GhoaoqqqKj7/+c/HPvvsU9DnWbduXbz99tsxe/bsqKmpiYaGhujSpUv06dMn9tlnn9h1110L+nywI2jL5wwgP23pfLF+/fqYM2dOvPfee/HBBx/EqlWrory8PCorK2PgwIGx9957R0lJSYE+CmidampqYtq0aVFdXR21tbXRs2fP6NOnTxxyyCFRWlr06682smrVqpg6dWosWbIkampqolu3btG7d+849NBDo6ysrCDPsWTJknjrrbdiyZIlUVdXF1VVVdGvX7844IADCnJ82FG1pfPF0qVL47333ot58+bFypUro6SkJCorK6Nv375x4IEHFuyiT9jRFD1uR0SMHj06nnjiiVi1alVERIwbNy46deoUX/va17JORsuXL4+rrroq/v73v2ce+9KXvhQDBw4s6DznnXde/PKXv4zFixdHRMTvfve7qKioiGuuuSbrZFRbWxs33HBDvPrqq5nHBg8eHCeeeGJB52ntnnzyybj//vtj5syZm3z/4MGD44orrohjjjmm2WZ65ZVX4t57793oXuifGTBgQFx22WXxf/7P/9mu51mwYEH853/+Zzz99NOZ799NqaqqimOPPTZuvvnmaNeu3XY9J6SuLZ8zmlq0aFGccsopG7148Ysvvhh9+vQp6HNBitrK+eLDDz+M559/Pl555ZV47bXXora2drPbduvWLb7yla/EBRdcEN26dct7fmjN5syZE2PHjo2JEydu8l/z9uzZM84666wYNWpUwaLylixdujTuuuuuePbZZzf6f3VERGVlZZx22mlx1VVXRefOnbfpOd566634yU9+EpMnT866YOoze+yxR5x//vlx3nnn+Yst2EBbOF/U19fH5MmT48UXX4y//e1vMXfu3M1u2759+zj++OPj4osvjkMOOWSbPw7YEZU0NjY2NscTvfzyy3H55Zdn/Q99zz33jCOOOCIqKytj7ty5MXHixFizZk3m/XvvvXf893//9xZPDPleuf2Z6dOnx7nnnht1dXWZx3r27BnDhw+P7t27R3V1dbz00ktZtyzZZZdd4re//W306tUr5487ZWvXro1vf/vb8fTTT29125KSkrjoooviuuuuK+ovZY2NjXH77bfHww8/HLl865566qlx2223bdP/7H7xi1/EnXfeGevWrct5n7feeis6dOiQ93PBjqCtnzM2ZdSoUTFp0qSNHhe3aeva0vniD3/4Q3znO9+J9evX5zVP9+7d49Zbb43jjjsur/2gtXryySfjxhtv3GQUamrQoEExbty46N27d9HmmTx5clx77bWxfPnyrW67++67x7hx42LffffN6zn+67/+K+6666745JNPtrrtUUcdFT/60Y+ioqIir+eAHVFbOF+sWrUqjj/++FixYkVes5SUlMR5550XY8aMyboFMLRlzRa3Iz795f6mm27KCsqbs99++8VPf/rTrf7hf1vjdkTEn//85/jWt74VNTU1W922T58+MW7cuIJfRd6afetb34qnnnoqsy4pKYmjjjoqBgwYEPX19fH2229vdFXTFVdcEVdeeWXRZvrRj34UP//5z7MeO+SQQ+KAAw6Idu3axbvvvhuvvvpq1h9KzzjjjLj99tvzep7vf//78ctf/jLrse7du8cRRxwRVVVV0aFDh/joo49i5syZMX369EwAF7dpy9ryOWNT/vjHP8a11167yfeJ27R1bel88dBDD8Vtt92W9Vh5eXkMHjw4dt999+jatWvU1dXFjBkz4rXXXsuK4O3bt4+7777bC5mTvFdeeSVGjx6d9f295557xuGHHx6VlZUxb968jS502meffeKxxx7b5iumt2TGjBlxzjnnZIWznj17xrHHHhs9evSIxYsXx8SJE7MudOrZs2f89re/jaqqqpye4/HHH48bbrgh67H99tsvDj300OjYsWPMnj07XnnllazwPWzYsLj//vujfftm+QfW0Cq1lfPFypUrY8iQIVmPlZSUxOc///nYf//9o3v37tGuXbtYsGBB/PWvf42PPvooa9tTTjkl/uM//qNFbs0CrU2zxu2IiPfffz/Gjh0bkyZN2uQ/Ldlll13iK1/5SowePTqnK+e2J25HRHzwwQcxduzYeO655zYZ3SsqKuK0006Lq6++uignytbqV7/6Vdx8882Z9a677ho/+9nPYr/99sva7q9//WtcddVV8fHHH0fEpyfj+++/P44++uiCzzRp0qQYNWpUZv25z30u7r777hg2bFjWdu+8805cfvnlUV1dnXns+9//fnzlK1/J6Xn+8z//M+66667MuqqqKr7zne/EF77whU3+j6O2tjZefPHFePTRR+ORRx4Rt2mT2vI5Y1NqamriS1/6UixbtixKS0ujW7du8eGHH2beL27TlrW188Vncfuzf0581llnxRFHHLHJq63mz58f3/nOd2LKlCmZx8rLy+O5556Lnj17bs+HCC1m6dKl8cUvfjHrZ3nMmDFxwQUXbPIWlRt+/5966qkxduzYgs6zdu3a+OIXvxgLFy7MPHbxxRdv9haVzzzzTOaxQw45JB577LGtPse7774b/+///b/Mn3fLysri1ltv3eh2RvPnz49vfOMbMWvWrMxjo0ePjmuuuWabPz5IWVs6X2wYt/fZZ58YOXJknHrqqZu8JdnatWtj/Pjxcc8992TdDeGGG26I8847b7s+RtgRNHvc/sxHH32UeVGAVatWRY8ePaJv375xyCGHtMg9iz97UYDFixfHihUrMi8KcNhhhzXL/Ztak7q6ujjxxBMzIaZDhw7xxBNPRL9+/Ta5/ZQpU+JrX/ta5kqmQYMGxYQJEwo6U2NjY3z5y1/O3JOzpKQkHnnkkY3+pvMz//jHP+LLX/5yrF27NiI+/ZvTP/3pT1sNz3//+9/j9NNPz/wi2q9fv3j00Ueje/fuBfxoYMfSls8Zm3P99dfH73//+4iIOPvss2P27NlZv3yL27RVbfF88cgjj8Qbb7wRV199dfTt23er86xbty4uvvjieO211zKPnXPOOXHjjTfm+6FBq9D0X0ReeeWVccUVV2xy27Vr18YZZ5wRs2fPjohPfx7/8Ic/5H07kC154IEH4o477siszzzzzLj11ls3uW1DQ0NcfPHFWa/BdO+99271X1OMHj06Jk6cmFnfcccdcfrpp29y2+XLl8epp54ay5Yti4iIjh07xgsvvBC77LJLrh8S7DDa0vli5cqVcc4558TVV1+d8+u6/frXv876faCysjL+/Oc/t7lmBU212L9f6Nq1a5xwwglx7rnnxqhRo+L//t//G0OGDGmxF+Pr1KlTDB8+PL761a/G17/+9Rg5cmQceeSRbfIk8fjjj2ddYXjppZdu9g+dERFDhw6N0047LbOePn161i9zhfDiiy9mvdjUl7/85c3+oTPi0yh9ySWXZNYffPBB/OY3v9nq89xwww2ZsN2hQ4cYP368sA1b0ZbPGZsyefLkTNjeZZddNntrEmiL2uL54uyzz46xY8fmFLYjPr3C8/vf/37W/cVfeOGFnO4DDq3NsmXL4vHHH8+sd99996x/JdFUhw4dsm7l0djYGPfee2/B5qmvr4/x48dn1l26dIkxY8ZsdvvS0tK4+eabs64Yveeee7b4HDNmzMg6Tw0ZMmSzYTvi0xeR/bd/+7fMuq6uLh544IEtPgf8f+3de1BU5/3H8Q8LCjFeqAiieK0Go+QC3uqlrbdMdaJpnWRqkxJttNQqJpOiSZxoLiZNjZJOIsW2M+TSOLZJI5gaY5pkWgxpibZqsCaCQ9Uo8QJ4v0SFCOzvD3+c8iwLexaWhN3zfs0ww/f4nOd8d8bz5exzznmeUOS0etG5c2dt3rzZ9sC2JN19993GYpLnzp0zBtMBp2JyHjTy7rvvWr+Hh4frRz/6kc99fvzjHxvxe++912Y5SVJqaqrPfe6++27jZomvnHbt2qWioiIrvu+++2x/EQWczKk1w5srV67oiSeesOJHH31UXbp08bsfIFQ5sV60ZO7cgQMHGuu8nDx5UsePH/e7H+Drlp+fbyzOPmvWLJ8LoI0dO9a46fXhhx/aWrPJjh07dhgLws2YMcPnAo79+/fXuHHjrLi4uFhHjhxpsr1nPfCsYd7ccccd6tq1qxW///77PvcBQo3T6oXL5WrRfNnTp0834k8//dTvPoBQw+A2DGfOnNGePXusOCUlxdaiKcnJyYqPj7figoICYwGI1qipqdE///lPK+7Vq5duueUWn/v17NlTycnJVlxUVNRoEYaGNm7caMR2vnADTufkmuHNb37zG+sC9tvf/naji0/AyagX/unfv78R109ZAAQTz7WQpk2bZmu/hu2qqqr00UcftUk+U6dO9Tsf6dognJ1jdOjQwdaCsJGRkZo4caIVHzt2TPv27bOVGxAqnFgvWqJfv35G3PCNOMCpGNyGoaioyFigICUlxfa+Db/knTt3Tvv37w9ITvv379f58+dbnVNtba3xZHZDNTU1xlMWSUlJSkhI8D9ZwGGcWjO8KS4u1rp16yRJUVFRWrFihe19ASegXvjn0qVLRtySJ8CBr9uuXbus3+vXWLLD81xsOAd9oPIJDw+3dTPLWz4N+2nowoULxuKQw4YNs71+R8OaIgXuMwPBwmn1oqW4PgAaY3AbhvrFGOo1fCXWl6SkJCP+7LPP2iSnoUOHtjgnz77q/fe//9Xly5et+KabbvIjQ8C5nFozPNXW1uqxxx6zniZduHAh0xoBHqgX9rndbmOATJLx9DoQDE6cOKGLFy9asT/nl2d9CMT5VVdXp8OHD1tx//79df3119vad9CgQYqKivKZT3utKUB758R60VIN1wmRuD4AJIlbPDB4flns1auX7X092wbqi6dnP717925xTocOHfLabu/evUY8ePBgSde+XBYWFmrTpk0qKSlRRUWFIiMjFRMTo5SUFE2ePFmTJk0yFn0CnMSpNcPTK6+8opKSEknX6kfDxeYAXEO9sG/nzp0qLy+34sTERHXv3j1g/QNfhdacXz169FCHDh2shd4Dcc4fO3ZMVVVVLconLCxM8fHx1mDXkSNHVFNT0+iJyUDWlEDVOSAYOLFetITb7daWLVuMbWPGjGl1v0Cw48ltGCorK43Yn7uAnm0rKiqCJqcDBw4YcY8ePVRRUaG5c+cqLS1NW7Zs0WeffabLly/r7NmzOnDggHJzc7Vw4ULdeeed1qAW4DROrRkNff7551q7dq2kaxezTz31lM/FbwAnol7YV19T6tmddxRoTzzPLztz7NcLCwsz2nv29VXn49n+6tWrXufBb80xYmNjjYVqA/GZgWDhxHrREps2bTIWqOzXr1+jtz4AJ2JwGwbP+ZvsvnrjrW3DaT7ae04N59usP+Y999yj7du3+zxGSUmJUlNTVVhYaDsvIFQ4tWY09MQTT1hPdtx1110aOXKk7eMBTkK9sCc3N1f//ve/rTg6Olr33ntvQPoGvkqtOb8829fU1OjLL79sN/l468/bts6dO9vu3+Vy6brrrrPiQNUUIBg4sV746+TJk1q9erWxbdGiRbxFDojBbXi4cuWKEdtdAMVb20BdkHnm1LFjx4Dn1HB+L0nKzMzU8ePHJV2b/y4rK0vbtm3Tp59+qg8++EC//OUvjVeTLl++rIyMDB09etR2bkAocGrNqLdx40brJlj37t318MMP2z4W4DROrxd2lJaW6le/+pWxbenSperWrVur+wa+aq055721b+3gUKDz8XbOt6ameB4jEINhQLBwYr3wR01NjZYsWaKzZ89a2771rW/pBz/4Qav6BUIFg9swNJxXSvLvgsyzrWdf7Tknzz82Fy5ckCTdcccdys3N1bRp0xQTE6OOHTuqd+/emjVrljZv3qxbb73V2GflypW2cwNCgVNrhiSdPn1amZmZVrx06VJFR0fbPhbgNE6uF3acOXNGixYtMr5QT506VXfeeWer+gW+Lq05v7y1r66ublf5eDvnA3mM1n5eIJg4sV74Y9WqVY3e6lq1ahVPbQP/jwUl26nJkyfr2LFjbdb/I4884nXBM887jFevXrVdyD1f/fH37mZTvOVkl92cvH3GgQMHauXKlcbcdw116dJFWVlZmj59unVneOvWrTp06JAGDhxoO0cgEKgZzedklz85PfPMMzp37pwkafTo0Zo5c6bt4wBfJ+pF8znZFcicLl++rAULFhjzaN54443cNEdQa835JTU+x/wdXGrrfLyd84E8RqDqHBAMnFgv7Hr55Ze1fv16K+7QoYNeeOEFvxa5BEIdT27D0KlTJyP25w6j591Rz74ClZM/d2Ht5uRte1pams8/ir169dJdd91lxW63WwUFBbbzA4KdU2tGQUGB/vrXv0q6dvH81FNP+Zkl4DxOrRe+XL16VQ888ID27NljbevTp49efPFFv+brBdqb1pzzUuNzzN85b9s6H2/nfGtqimf7QNU5IBg4sV7Y8Ze//EXPPfecFYeFhWnVqlUaN25ci/oDQhWD2zB4Fl1/5obynNeqrb54+jN/lt2cvP3xmzx5sq1jeLbbvXu3zeyA4OfEmnHp0iWtWLHCiufPn69vfvOb/iUJOJAT64UvdXV1euSRR4xFqWNjY/WHP/xBcXFxfvcHtCetOecl8xyLiIho9ZPMrV0E1s4Cc62pKXV1dcYAGoPbcBIn1gtf/v73v2v58uVyu93Wtscff1wzZszwuy8g1DEtSTuVlZXVpvOs9e3b1+v2nj17GnF5ebni4+Nt9VlRUWHEdvfzxTMnz+MEIifP7TExMerevbutYyQmJhpxZWWl7fyAQKFmNJ1ToGvGunXrVF5eLkkaMGCAfv7zn7cgS+DrQ71oOqe2uMZozpNPPmm9BSJJ3bp100svvaR+/fr53RfQ3rTm/HK73cY1tWdfX3U+knmNHxERoZiYmIAe4+TJk6qpqbHiQNU5IBg4sV40Z/v27crIyFBtba217cEHH1Rqaqpf/QBOweB2O3XzzTd/Lcf1fPqwvLxcKSkptvatH+xpqq9A5XT8+HHb+9rNadCgQUbctWtX28fwbFs/By/wVaJmNN1PoGtGw1XKy8rKlJyc7LPfhhemkvS9733PiP/2t78pISHBdp5Aa1Avmu6nLa4xmpKZmakNGzZYcadOnZSTk6Mbb7zRr36A9qo159epU6eMOW4DsZ5NQkKCIiMjrZt7/uTjdruNwa2+ffuqQ4cOjdoFsqawhg+cxIn1oimffPKJ0tPTjXm7582bp/T0dNt9AE7DtCQweA7ylpSU2N63uLi42b4CldO+fftanFNTXzxvuOEGI/ZcAKI5bbXIFRAMnFoz6rndbtXW1vr88eT57w1fNwRCldPrRUO///3v9fLLL1txZGSkfve739m6WQYEi549e6pLly5W7M/55VkfAnHOu1wuDRgwwIrLyspsTzVw8OBBY8qQpvIJZE0JVJ0DgoET64U3paWlSktLM441a9YsLV261HYfgBMxuA3D8OHD5XL977+FP/NH/+c//7F+j46O1uDBgwOSU2JiovF0tD85NWwbHh6u4cOHe203aNAg41Wh06dPG68FNsfzFSW705kAocCpNQOA/6gX1/zxj3/UmjVrrDgiIkJr1qzR2LFjbR8bCBYjRoywfj916pSOHDlia7+ioiIjHjVqVEDyGTlypPV7bW2tsZBrczxrQ1P5dOvWzXhopqSkxPY0UHaPAYQqp9ULT2VlZfrpT3+q8+fPW9tmzJjBwvWADQxuwxATE6Nbb73Vinfv3m1rDundu3cbg7wTJkxQRERgZr2JiIjQd7/7XSsuLy+39YelsrLS+DKckpLS5MCzy+UyFoasqqqy/USZ5x+voUOH2toPCAVOrBnLly9XaWmpXz+jR482+sjPzzf+vU+fPi3/wECQcGK98LRp0yY988wzVuxyuZSZmWl7EWsg2Hj+33733Xdt7ff+++9bv0dGRmr8+PFtks97771naz/PdlOmTLF1jKtXryo/P99n/9XV1frwww+tuHfv3ho2bJit3IBQ4cR6Ua+yslJz587VyZMnjf1Wr15tPBgAwDvOEjQybdo06/fa2lq98cYbPvd5/fXXm+wj0DlJ0muvveZznz//+c+qq6uzndP06dONOC8vz1Zunu0C9ccUCBZOrRkA/OfkepGfn6/ly5db0xCFhYXp6aefbnT9AYSSKVOmGHPN5ubmGnPjerN9+3YdOnTIiidMmKBOnToFJJ/Ro0frG9/4hhW/8847unDhQrP7lJWVadu2bVaclJTU5MK5Ustqyttvv23kMXXqVJ/7AKHGifVCuraez9y5c3Xs2DFr2/jx47VmzZqA3cwHQh2D22hk1qxZxhQdL730kvEHw9OOHTu0efNmKx42bJgmTZrU7DGys7M1ZMgQ6yc7O7vZ9lOmTFFiYqIVv/XWW9q5c2eT7Q8dOmTMZRkbG6sf/vCHzR5j7NixxqtHeXl5Pl9P3rBhg9Fm0KBBGjduXLP7AKHGqTUDgP+cWi/+9a9/KSMjw5jybNmyZdQZhLwePXoY/88///xz5eTkNNm+urraeLshLCxMCxcubLL90aNHjfPd11sQHTt2VFpamhVfvHhRq1atarJ9XV2dnnzySeNmlq9F3YYNG6YJEyZY8c6dO7Vp06Ym2585c0bPP/+8FUdFRWnevHnNHgMIRU6sF1988YXS0tJ08OBBa9uIESP029/+Vh07dmx2XwD/w+A2GunUqZNRhKurqzV37lyvizp89NFHSk9PNxZDy8jIUFhYWEBzcrlcWrx4sRW73W6lp6dr+/btjdqWlJTovvvuM+a3u//++xUVFeXzOA8//LDCw8MlXXuibP78+SooKGjUrq6uTuvXr9eKFSuM7UuWLOG1ITiOk2sGAP84sV7s3btX6enpxj6LFy/WnDlzAvQJgPZtwYIFuv766604Oztbr776qjEAJF0b5E1LS9OBAwesbbfffnvAp+e499571atXLyveuHGjVq9e3WiR+C+++EJLliwxakFKSopuu+02n8fIyMgwnrh8/PHHtWXLlkbtjhw5op/85Cc6ffq0tW3OnDmKi4vz6zMBocJJ9eLLL79Uenq69u7da2276aablJOTo+uuuy6AnwIIfWHuht8YgAYWL16sd955x4rDwsI0fvx4DRkyRDU1Nfrkk08aPdm8YMECZWRk+Ow7Oztba9euteL7779fDzzwgM/9fv3rX+vFF180to0YMUI333yzXC6XSktLtW3bNuOL8Pe//30999xzPvuu96c//UlPP/20sW3o0KEaOXKkunbtqlOnTqmwsNB4bUiSfvazn+mhhx6yfRwg1Di1Ztgxe/Zs7dixw4rz8/OZZxuO5qR68eijj+rNN980ttXfSPfHypUrNXPmTL/3A9qDgoICLVy40BigGjBggMaMGaPo6GiVlZXpgw8+UFVVlfXvgwcP1htvvKHOnTs32e/Ro0eN+WwTEhK0detWn/kUFxcrNTVVV65csbbFxcVpwoQJiomJUUVFhbZu3WpMQRAbG6u8vDzFx8fb+syvv/56owdh6r9TREVF6eDBg/rHP/5hvNExevRovfLKK8bUDIDTOKVe7NixQ7Nnzza2uVwuv2/ijxo1SuvWrfNrHyDUMIEPmvTss8+qtrbWWhDB7XarsLBQhYWFjdqGhYVpzpw5+sUvftGmOS1evFhVVVVav369te3jjz/Wxx9/7LX97bffbryqZEdqaqqqq6v1/PPPW3N87du3z+tTZdK1z/7ggw82+woU4AROrRkA/OekeuHtOZLa2lq/8/N8ag0IJhMnTtSzzz6rFStWWANEhw8f1uHDh722Hzp0qNauXdvsQFVrJCUlKTs7Ww899JDOnTsnSTpx4oRyc3O9tu/Tp4+ys7NtD2xL0j333KOLFy8qKyvLGsBu7jvFmDFjlJWVxcA2HM8p9cLb9UFL/tZzfQAwLQmaERkZqaysLK1evdqYi9JTcnKycnJytGzZsoC/KuzJ5XLpscceU05OjpKTk5tsl5iYqMzMTL3wwguKjIz0+zjz5s1TXl6eJk2a1OQFZkREhCZOnKi8vDwGtgE5u2YA8A/1AnCemTNn6s0339Rtt93W5PV1bGysFi1apA0bNrT5G07f+c539Pbbb2vmzJlNTgHQrVs3zZ49W2+99VaLpjuYP3++XnvtNY0bN67JqQv79u2rZcuW6dVXX1V0dLTfxwBCkRPrBYCWY1oS2FZaWqr9+/ersrJS4eHhiouLU1JSkvr37/+15VRWVqbi4mKdOHFCtbW16tmzp2644QYNGTIkYMc4c+aMdu/ercrKSl24cEFdu3ZVfHy8Ro0apS5dugTsOECocWrNAOA/6gXgLGfPnlVRUZEqKip06dIl9ejRQ3379tXw4cNbNG1Pa126dEm7du1SeXm5zp8/r+7duyshIUEjR44M2KJulZWV2rNnjyorK1VVVaW4uDgNHDhQt9xyS0D6B0KVE+sFAP8wuA0AAAAAAAAACDpMSwIAAAAAAAAACDoMbgMAAAAAAAAAgg6D2wAAAAAAAACAoMPgNgAAAAAAAAAg6DC4DQAAAAAAAAAIOgxuAwAAAAAAAACCDoPbAAAAAAAAAICgw+A2AAAAAAAAACDoMLgNAAAAAAAAAAg6DG4DAAAAAAAAAIIOg9sAAAAAAAAAgKDD4DYAAAAAAAAAIOgwuA0AAAAAAAAACDoMbgMAAAAAAAAAgg6D2wAAAAAAAACAoMPgNgAAAAAAAAAg6DC4DQAAAAAAAAAIOgxuAwAAAAAAAACCDoPbAAAAAAAAAICgw+A2AAAAAAAAACDoMLgNAAAAAAAAAAg6DG4DAAAAAAAAAIIOg9sAAAAAAAAAgKDD4DYAAAAAAAAAIOgwuA0AAAAAAAAACDoMbgMAAAAAAAAAgs7/AbELdr73z/BCAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 720x480 with 1 Axes>"
      ]
     },
     "metadata": {
      "image/png": {
       "height": 491,
       "width": 731
      }
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "dataset = MyDataset(\"dataset/data.csv\")\n",
    "df = dataset.get_df()\n",
    "cct_names,rates_names,error_names = dataset.get_var_col_names()\n",
    "\n",
    "rates = dataset.get_rates()\n",
    "_ = plt.hist(rates[1:])\n",
    "\n",
    "df[['time'] +  rates_names]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "c53d82b2-7aaf-4aad-9566-14825fda8006",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Auto-assigning NUTS sampler...\n",
      "Initializing NUTS using jitter+adapt_diag...\n",
      "Multiprocess sampling (4 chains in 4 jobs)\n",
      "NUTS: [k1, k2, k3, k4, k5, k6, k7, k8, k9, k10, k11, sigma]\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "\n",
       "<style>\n",
       "    /* Turns off some styling */\n",
       "    progress {\n",
       "        /* gets rid of default border in Firefox and Opera. */\n",
       "        border: none;\n",
       "        /* Needs to be in here for Safari polyfill so background images work as expected. */\n",
       "        background-size: auto;\n",
       "    }\n",
       "    progress:not([value]), progress:not([value])::-webkit-progress-bar {\n",
       "        background: repeating-linear-gradient(45deg, #7e7e7e, #7e7e7e 10px, #5c5c5c 10px, #5c5c5c 20px);\n",
       "    }\n",
       "    .progress-bar-interrupted, .progress-bar-interrupted::-webkit-progress-bar {\n",
       "        background: #F44336;\n",
       "    }\n",
       "</style>\n"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "\n",
       "    <div>\n",
       "      <progress value='48000' class='' max='48000' style='width:300px; height:20px; vertical-align: middle;'></progress>\n",
       "      100.00% [48000/48000 00:37&lt;00:00 Sampling 4 chains, 11 divergences]\n",
       "    </div>\n",
       "    "
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Sampling 4 chains for 2_000 tune and 10_000 draw iterations (8_000 + 40_000 draws total) took 37 seconds.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.7755015446026299 [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]\n"
     ]
    }
   ],
   "source": [
    "# 全部是1级动力学，不限制K的范围\n",
    "k_kinetics = np.repeat(1,11).astype(np.uint8).tolist()\n",
    "idata = opt_model(dataset, k_kinetics, k_sigma_priors=0.1, kf_type=1, draws=10000, tune=2000, chains=4, cores=4)\n",
    "r2 = eval_model(idata, dataset)\n",
    "print(r2, k_kinetics)\n",
    "# 0.7754844473588883 [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1] k_sigma_priors=0.1\n",
    "# 0.7535096123714303 [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1] k_sigma_priors=0.01\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "34d08305-990a-49ec-9bcf-f516a527daaf",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Auto-assigning NUTS sampler...\n",
      "Initializing NUTS using jitter+adapt_diag...\n",
      "Multiprocess sampling (4 chains in 4 jobs)\n",
      "NUTS: [k1, k2, k3, k4, k5, k6, k7, k8, k9, k10, k11, sigma]\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "\n",
       "<style>\n",
       "    /* Turns off some styling */\n",
       "    progress {\n",
       "        /* gets rid of default border in Firefox and Opera. */\n",
       "        border: none;\n",
       "        /* Needs to be in here for Safari polyfill so background images work as expected. */\n",
       "        background-size: auto;\n",
       "    }\n",
       "    progress:not([value]), progress:not([value])::-webkit-progress-bar {\n",
       "        background: repeating-linear-gradient(45deg, #7e7e7e, #7e7e7e 10px, #5c5c5c 10px, #5c5c5c 20px);\n",
       "    }\n",
       "    .progress-bar-interrupted, .progress-bar-interrupted::-webkit-progress-bar {\n",
       "        background: #F44336;\n",
       "    }\n",
       "</style>\n"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "\n",
       "    <div>\n",
       "      <progress value='48000' class='' max='48000' style='width:300px; height:20px; vertical-align: middle;'></progress>\n",
       "      100.00% [48000/48000 00:25&lt;00:00 Sampling 4 chains, 302 divergences]\n",
       "    </div>\n",
       "    "
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Sampling 4 chains for 2_000 tune and 10_000 draw iterations (8_000 + 40_000 draws total) took 26 seconds.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.6191166714084839 [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]\n"
     ]
    }
   ],
   "source": [
    "# 全部是1级动力学，限制 k>0\n",
    "k_kinetics = np.repeat(1,11).astype(np.uint8).tolist()\n",
    "idata = opt_model(dataset, k_kinetics, k_sigma_priors=0.01, kf_type=0, draws=10000, tune=2000, chains=4, cores=4)\n",
    "r2 = eval_model(idata, dataset)\n",
    "print(r2, k_kinetics)\n",
    "# 0.6412686764018853 [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1] k_sigma_priors=0.1\n",
    "# 0.6191166714084839 [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1] k_sigma_priors=0.01"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "117a7be0-84b8-424e-b741-1a389126880b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>mean</th>\n",
       "      <th>sd</th>\n",
       "      <th>hdi_3%</th>\n",
       "      <th>hdi_97%</th>\n",
       "      <th>mcse_mean</th>\n",
       "      <th>mcse_sd</th>\n",
       "      <th>ess_bulk</th>\n",
       "      <th>ess_tail</th>\n",
       "      <th>r_hat</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>k1</th>\n",
       "      <td>-2.382871e-04</td>\n",
       "      <td>0.099421</td>\n",
       "      <td>-0.186561</td>\n",
       "      <td>0.186902</td>\n",
       "      <td>4.471629e-04</td>\n",
       "      <td>5.296000e-04</td>\n",
       "      <td>49411.168032</td>\n",
       "      <td>29111.388085</td>\n",
       "      <td>0.999984</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>k2</th>\n",
       "      <td>3.163900e-06</td>\n",
       "      <td>0.000176</td>\n",
       "      <td>-0.000336</td>\n",
       "      <td>0.000327</td>\n",
       "      <td>9.605000e-07</td>\n",
       "      <td>8.908000e-07</td>\n",
       "      <td>33588.569267</td>\n",
       "      <td>28117.350692</td>\n",
       "      <td>1.000049</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>k3</th>\n",
       "      <td>2.320253e-02</td>\n",
       "      <td>0.079675</td>\n",
       "      <td>-0.123965</td>\n",
       "      <td>0.175754</td>\n",
       "      <td>3.809907e-04</td>\n",
       "      <td>3.953232e-04</td>\n",
       "      <td>43739.496054</td>\n",
       "      <td>29806.225075</td>\n",
       "      <td>1.000027</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>k4</th>\n",
       "      <td>-5.134000e-07</td>\n",
       "      <td>0.000043</td>\n",
       "      <td>-0.000082</td>\n",
       "      <td>0.000081</td>\n",
       "      <td>2.370000e-07</td>\n",
       "      <td>2.128000e-07</td>\n",
       "      <td>32771.879805</td>\n",
       "      <td>27241.594908</td>\n",
       "      <td>1.000106</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>k5</th>\n",
       "      <td>3.312190e-05</td>\n",
       "      <td>0.000865</td>\n",
       "      <td>-0.001536</td>\n",
       "      <td>0.001745</td>\n",
       "      <td>5.403600e-06</td>\n",
       "      <td>4.600300e-06</td>\n",
       "      <td>25767.048597</td>\n",
       "      <td>24346.971425</td>\n",
       "      <td>0.999993</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>k6</th>\n",
       "      <td>-1.897000e-06</td>\n",
       "      <td>0.000053</td>\n",
       "      <td>-0.000102</td>\n",
       "      <td>0.000097</td>\n",
       "      <td>3.269000e-07</td>\n",
       "      <td>2.789000e-07</td>\n",
       "      <td>25884.402209</td>\n",
       "      <td>24044.757112</td>\n",
       "      <td>1.000027</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>k7</th>\n",
       "      <td>2.247107e-04</td>\n",
       "      <td>0.000495</td>\n",
       "      <td>-0.000688</td>\n",
       "      <td>0.001185</td>\n",
       "      <td>2.604100e-06</td>\n",
       "      <td>2.291300e-06</td>\n",
       "      <td>36384.256434</td>\n",
       "      <td>27356.852293</td>\n",
       "      <td>0.999989</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>k8</th>\n",
       "      <td>-7.593300e-03</td>\n",
       "      <td>0.003649</td>\n",
       "      <td>-0.014491</td>\n",
       "      <td>-0.000795</td>\n",
       "      <td>2.326070e-05</td>\n",
       "      <td>1.681750e-05</td>\n",
       "      <td>24721.299855</td>\n",
       "      <td>25935.097742</td>\n",
       "      <td>1.000131</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>k9</th>\n",
       "      <td>2.448837e-02</td>\n",
       "      <td>0.011986</td>\n",
       "      <td>0.002189</td>\n",
       "      <td>0.047620</td>\n",
       "      <td>6.272020e-05</td>\n",
       "      <td>4.604150e-05</td>\n",
       "      <td>36677.550249</td>\n",
       "      <td>28688.205008</td>\n",
       "      <td>1.000027</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>k10</th>\n",
       "      <td>9.241510e-05</td>\n",
       "      <td>0.000013</td>\n",
       "      <td>0.000066</td>\n",
       "      <td>0.000117</td>\n",
       "      <td>7.320000e-08</td>\n",
       "      <td>5.190000e-08</td>\n",
       "      <td>34193.245750</td>\n",
       "      <td>28445.322695</td>\n",
       "      <td>1.000224</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>k11</th>\n",
       "      <td>-3.179210e-04</td>\n",
       "      <td>0.000170</td>\n",
       "      <td>-0.000632</td>\n",
       "      <td>0.000005</td>\n",
       "      <td>1.059500e-06</td>\n",
       "      <td>7.733000e-07</td>\n",
       "      <td>25776.283183</td>\n",
       "      <td>26651.829741</td>\n",
       "      <td>1.000070</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>sigma</th>\n",
       "      <td>1.067338e-02</td>\n",
       "      <td>0.001795</td>\n",
       "      <td>0.007594</td>\n",
       "      <td>0.014040</td>\n",
       "      <td>1.204340e-05</td>\n",
       "      <td>8.589600e-06</td>\n",
       "      <td>22448.982823</td>\n",
       "      <td>20505.893074</td>\n",
       "      <td>0.999995</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "               mean        sd    hdi_3%   hdi_97%     mcse_mean       mcse_sd  \\\n",
       "k1    -2.382871e-04  0.099421 -0.186561  0.186902  4.471629e-04  5.296000e-04   \n",
       "k2     3.163900e-06  0.000176 -0.000336  0.000327  9.605000e-07  8.908000e-07   \n",
       "k3     2.320253e-02  0.079675 -0.123965  0.175754  3.809907e-04  3.953232e-04   \n",
       "k4    -5.134000e-07  0.000043 -0.000082  0.000081  2.370000e-07  2.128000e-07   \n",
       "k5     3.312190e-05  0.000865 -0.001536  0.001745  5.403600e-06  4.600300e-06   \n",
       "k6    -1.897000e-06  0.000053 -0.000102  0.000097  3.269000e-07  2.789000e-07   \n",
       "k7     2.247107e-04  0.000495 -0.000688  0.001185  2.604100e-06  2.291300e-06   \n",
       "k8    -7.593300e-03  0.003649 -0.014491 -0.000795  2.326070e-05  1.681750e-05   \n",
       "k9     2.448837e-02  0.011986  0.002189  0.047620  6.272020e-05  4.604150e-05   \n",
       "k10    9.241510e-05  0.000013  0.000066  0.000117  7.320000e-08  5.190000e-08   \n",
       "k11   -3.179210e-04  0.000170 -0.000632  0.000005  1.059500e-06  7.733000e-07   \n",
       "sigma  1.067338e-02  0.001795  0.007594  0.014040  1.204340e-05  8.589600e-06   \n",
       "\n",
       "           ess_bulk      ess_tail     r_hat  \n",
       "k1     49411.168032  29111.388085  0.999984  \n",
       "k2     33588.569267  28117.350692  1.000049  \n",
       "k3     43739.496054  29806.225075  1.000027  \n",
       "k4     32771.879805  27241.594908  1.000106  \n",
       "k5     25767.048597  24346.971425  0.999993  \n",
       "k6     25884.402209  24044.757112  1.000027  \n",
       "k7     36384.256434  27356.852293  0.999989  \n",
       "k8     24721.299855  25935.097742  1.000131  \n",
       "k9     36677.550249  28688.205008  1.000027  \n",
       "k10    34193.245750  28445.322695  1.000224  \n",
       "k11    25776.283183  26651.829741  1.000070  \n",
       "sigma  22448.982823  20505.893074  0.999995  "
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# idata.posterior[\"k1\"].sel(draw=slice(990, 1000))\n",
    "az.summary(idata, round_to=10)\n",
    "# _ = az.plot_trace(idata, combined=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "7a3483de-dbda-4ce1-9155-a61466f14f29",
   "metadata": {},
   "outputs": [],
   "source": [
    "def solve_01(dataset, kk_list, save_path, k_sigma_priors=0.01):\n",
    "    if osp.exists(save_path):\n",
    "        f  = open(save_path, 'a')\n",
    "    else:\n",
    "        f  = open(save_path, 'w+')\n",
    "        f.write(f\"k_kinetics,k_sigma_priors, ks_res,r2\\n\")\n",
    "\n",
    "    for k_kinetics in kk_list:\n",
    "        idata = opt_model(dataset, k_kinetics, k_sigma_priors=k_sigma_priors, draws=10000, tune=2000, chains=4, cores=4)\n",
    "        r2 = eval_model(idata, dataset)\n",
    "        line_str = f'\"{k_kinetics}\",\"{k_sigma_priors}\",\"{ks_res}\",\"{r2}\"\\n'\n",
    "        f.write(line_str)\n",
    "        f.flush()\n",
    "        print(f\"k_kinetics={k_kinetics}, k_sigma_priors={k_sigma_priors}, r2={r2}\")\n",
    "    f.close()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "37f7c723-3a8f-4984-aa16-4ee7c1c40840",
   "metadata": {},
   "outputs": [],
   "source": [
    "def solve_02(dataset, kk_list, save_path):\n",
    "\n",
    "    if osp.exists(save_path):\n",
    "        f  = open(save_path, 'a')\n",
    "    else:\n",
    "        f  = open(save_path, 'w+')\n",
    "        f.write(f\"k_kinetics,ks_res,r2\\n\")\n",
    "\n",
    "    for k_kinetics in kk_list:\n",
    "        ks_res = ltq_fit(dataset,k_kinetics )\n",
    "        predict= get_target(ks_res, dataset.get_df(), k_kinetics)\n",
    "        rates_y = dataset.get_rates()\n",
    "        r2 = r2_loss(predict[1:],rates_y[1:])\n",
    "        line_str = f'\"{k_kinetics}\",\"{ks_res}\",\"{r2}\"\\n'\n",
    "        f.write(line_str)\n",
    "        f.flush()\n",
    "        print(f\"k_kinetics={list(k_kinetics)}, r2={r2}\")\n",
    "    f.close()\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "31726371-b0f3-4759-8eaa-a9cf743f0a65",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "loading dataset from dataset/data.csv....\n",
      "done!\n",
      "k_kinetics=[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], r2=0.8124185944967925\n",
      "k_kinetics=[1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], r2=0.8132599345385618\n",
      "k_kinetics=[0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0], r2=0.8124188475175166\n",
      "k_kinetics=[0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0], r2=0.8125377132530056\n",
      "k_kinetics=[0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0], r2=0.8124180715538909\n",
      "k_kinetics=[0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0], r2=0.8124229761423063\n",
      "k_kinetics=[0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0], r2=0.8124170382487262\n",
      "k_kinetics=[0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0], r2=0.8134512961052908\n",
      "k_kinetics=[0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0], r2=0.7983958940970844\n",
      "k_kinetics=[0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0], r2=0.793370355698932\n",
      "k_kinetics=[0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0], r2=0.8108134886555125\n",
      "k_kinetics=[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1], r2=0.8142171969928428\n",
      "k_kinetics=[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1], r2=0.7767488693235725\n"
     ]
    }
   ],
   "source": [
    "kk_list_all = list(product([0,1], repeat=11))\n",
    "kk_list = np.vstack([np.repeat(0,11).astype(np.uint8), np.eye(11,11, dtype=np.uint8), np.repeat(1,11).astype(np.uint8)]).tolist()\n",
    "dataset = MyDataset(\"dataset/data.csv\")\n",
    "\n",
    "solve_02(dataset,kk_list,save_path=\"lsq_temp.csv\" )\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a8f72c53-decf-45f4-91fd-9c2129f1435a",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "pymc4",
   "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.9.16"
  },
  "vscode": {
   "interpreter": {
    "hash": "1cc05d71c091818f3a20105c117b9e65164e6ad677e8fb3b94fb2c20e86bdea0"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}