{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "97336ec7-efa5-4102-b3e6-1414f8f5ae54",
   "metadata": {},
   "source": [
    "# Zivich (2026): Pooled Logistic Regression for Survival Analysis\n",
    "\n",
    "In the corresponding pre-print, Zivich et al. propose a novel implementation of estimating equations for causal survival analysis. This example reviews those key results and demonstrates how to use ``ee_plogit``. See the corresponding paper for finer details on the full approach and how it works. This page primarily focuses on how to apply the code (and not the identification assumptions or underlying g-computation algorithm described in the paper).\n",
    "\n",
    "## Setup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "61e92251-fc93-40f6-8ec3-0e7a7e52c832",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Versions\n",
      "NumPy:        2.3.5\n",
      "SciPy:        1.16.3\n",
      "pandas:       2.3.3\n",
      "Matplotlib:   3.10.8\n",
      "Delicatessen: 4.1\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import scipy as sp\n",
    "import pandas as pd\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import delicatessen as deli\n",
    "from delicatessen import MEstimator\n",
    "from delicatessen.estimating_equations import ee_plogit\n",
    "from delicatessen.utilities import inverse_logit, spline, plogit_predict\n",
    "\n",
    "print(\"Versions\")\n",
    "print(\"NumPy:       \", np.__version__)\n",
    "print(\"SciPy:       \", sp.__version__)\n",
    "print(\"pandas:      \", pd.__version__)\n",
    "print(\"Matplotlib:  \", mpl.__version__)\n",
    "print(\"Delicatessen:\", deli.__version__)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "29ed5962-b056-477a-9594-da7fcf2408b7",
   "metadata": {},
   "source": [
    "## Example 1: Bladder Cancer Time-to-Recurrence\n",
    "\n",
    "The first example comes from a data set in Collett (2015). The data consists of 86 patients who received either placebo or the chemotherapeutic thiotepa following removal of superficial bladder tumors, with time measured in months. Baseline variables included the initial number of tumors and the diameter (in centimeters) of the largest initial tumor. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "6744141d-c1bc-4a28-9e35-1cd6c7f7d723",
   "metadata": {},
   "outputs": [],
   "source": [
    "d = pd.read_csv(\"data/collett.dat\", sep=r'\\s+',\n",
    "                names=['patient', 'time', 'delta', 'treat', 'init', 'size'])\n",
    "d['novel'] = d['treat'] - 1\n",
    "d['intercept'] = 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "6d244e7b-f572-409e-8686-9231448bd894",
   "metadata": {},
   "outputs": [],
   "source": [
    "t = np.asarray(d['time'])\n",
    "y = np.asarray(d['delta'])\n",
    "W = np.asarray(d[['novel', 'init', 'size', ]])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6a244b58-d2c5-4427-98ca-2de6742f83d4",
   "metadata": {},
   "source": [
    "Here, the parameter of interest was the effect comparing novel to standard treatment on disease-free survival at 59 months (end of follow-up). While treatment was randomized, both number of tumors and tumor diameter are adjusted for in this illustration. Here, we will review how pooled logistic regression can be used for these types of survival analyses.\n",
    "\n",
    "### Pooled Logistic Regression\n",
    "\n",
    "First, examine how a pooled logistic regression model can be fit to the data. Here, time is modeled using disjoint indicators which is handled automatically by ``ee_pooled_logistic`` by default."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "7fd92b94-779e-4b4b-8e41-3bc0a41b3969",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Number of unique event times (for specifying the initial values)\n",
    "unique_event_times = list(np.unique(d.loc[d['delta'] == 1, 'time']))\n",
    "params_plr = len(unique_event_times)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "1a4a1e8a-d12a-4fd5-a934-c594b59f4551",
   "metadata": {},
   "outputs": [],
   "source": [
    "def psi_plogit(theta):\n",
    "    ee_plog = ee_plogit(theta, t=t, delta=y, X=W)\n",
    "    return ee_plog"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "fcf6190a-d19e-4536-b761-2248ea389bcc",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "==============================================================\n",
      "              Estimation Method: M-estimator\n",
      "--------------------------------------------------------------\n",
      "No. Observations:          86 | No. Parameters:             24\n",
      "Solving algorithm:         lm | Max Iterations:           5000\n",
      "Solving tolerance:      1e-09 | Allow P-Inverse:             1\n",
      "Derivative Method:     approx | Deriv Approx:            1e-09\n",
      "Small N Correction:      None | Distribution:           Z-stat\n",
      "==============================================================\n",
      "   Theta   StdErr  Z-score      LCL      UCL  P-value  S-value \n",
      "--------------------------------------------------------------\n",
      "   -0.55     0.33    -1.66    -1.19     0.10     0.10     3.37 \n",
      "    0.26     0.08     3.15     0.10     0.42     0.00     9.28 \n",
      "    0.07     0.09     0.79    -0.11     0.26     0.43     1.22 \n",
      "==============================================================\n"
     ]
    }
   ],
   "source": [
    "inits = ([0., ]*W.shape[1] +          # Coefs for baseline vars\n",
    "         [-4., ] +                    # Coef for intercept\n",
    "         [0., ]*(params_plr - 1))  # Coefs for time terms\n",
    "estr = MEstimator(psi_plogit, init=inits)\n",
    "estr.estimate()\n",
    "estr.print_results(subset=[x for x in range(W.shape[1])])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "393865df-658c-4d57-977c-835ae1b99e2b",
   "metadata": {},
   "source": [
    "Here, we examine the 3 parameters corresponding to the baseline covariates. As detailed elsewhere, these coefficients can be interpreted as approximating the hazard ratio. Specifically, the hazard of recurrence among those on the novel drug is $\\exp(-0.55) = 0.58$ that of those on placebo, conditional on tumor size and count. \n",
    "\n",
    "Rather than model time non-parametrically, one can also use parametric specifications for how the discrete-time hazard of the event changes over time. As in the paper, we consider the use of splines to model time. Here, time was modeled restricted quadratic splines with knots at 10, 20, 30, and 40 months."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "da769605-cd19-47ad-bac7-0076fc4adf64",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Creating time design matrix\n",
    "t_steps = np.asarray(range(1, np.max(d['time']+1)))\n",
    "intercept = np.ones(t_steps.shape)[:, None]\n",
    "time_splines = spline(t_steps, knots=[10, 20, 30, 40],\n",
    "                      power=2, restricted=True, normalized=False)\n",
    "s_matrix = np.concatenate([intercept, t_steps[:, None], time_splines], axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "65e419cd-4e22-4845-97ce-583fc3b2da97",
   "metadata": {},
   "outputs": [],
   "source": [
    "def psi_plogit(theta):\n",
    "    ee_plog = ee_plogit(theta, t=t, delta=y, X=W, S=s_matrix)\n",
    "    return ee_plog"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "8b66c078-7589-4517-92a7-1f559cb129e9",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "==============================================================\n",
      "              Estimation Method: M-estimator\n",
      "--------------------------------------------------------------\n",
      "No. Observations:          86 | No. Parameters:              8\n",
      "Solving algorithm:         lm | Max Iterations:           5000\n",
      "Solving tolerance:      1e-09 | Allow P-Inverse:             1\n",
      "Derivative Method:     approx | Deriv Approx:            1e-09\n",
      "Small N Correction:      None | Distribution:           Z-stat\n",
      "==============================================================\n",
      "   Theta   StdErr  Z-score      LCL      UCL  P-value  S-value \n",
      "--------------------------------------------------------------\n",
      "   -0.54     0.32    -1.68    -1.16     0.09     0.09     3.41 \n",
      "    0.25     0.08     3.30     0.10     0.40     0.00     9.99 \n",
      "    0.07     0.09     0.78    -0.11     0.26     0.43     1.21 \n",
      "==============================================================\n"
     ]
    }
   ],
   "source": [
    "inits = ([0., ]*W.shape[1] +          # Coefs for baseline vars\n",
    "         [-4., ] +                    # Coef for intercept\n",
    "         [0., 0., 0., 0.])            # Coefs for time terms\n",
    "estr = MEstimator(psi_plogit, init=inits)\n",
    "estr.estimate()\n",
    "estr.print_results(subset=[x for x in range(W.shape[1])])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "141b0c9a-1508-4319-bbd4-94544692d56d",
   "metadata": {},
   "source": [
    "Here, we see the coefficients for the covariates are largely the same. While this is the case in this example, this may not always be true (especially when relying on more restrictive functional forms for time). \n",
    "\n",
    "### G-computation\n",
    "\n",
    "As also detailed elsewhere, the hazard ratio has some major challenges for interpretation, including non-collapsibility. Now we show how to transform these pooled logistic regression results into marginal risk functions, which are easier to interpret and avoid some of these challenges. We will use a g-computation algorithm here. Briefly, we will fit a pooled logistic regression model, then we will use those coefficients to project what would have happened if everyone as given thiotepa (placebo). Both tumor count and tumor size were modeled as linear relationships. For flexibility, we will fit pooled logistic regression models stratified by treatment. This approach to modeling does not rely on a proportional hazards assumption for the treatment (unlike the previous approach)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "00e28095-ca04-4f51-952b-78f1e37f0a81",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Data under the plans of interest (all-treat, all-placebo)\n",
    "d1 = d.copy()\n",
    "d1['novel'] = 1\n",
    "d0 = d.copy()\n",
    "d0['novel'] = 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "bbd5b25d-8e21-42e6-941d-a86557d4cb40",
   "metadata": {},
   "outputs": [],
   "source": [
    "a = np.asarray(d['novel'])\n",
    "t = np.asarray(d['time'])\n",
    "y = np.asarray(d['delta'])\n",
    "W = np.asarray(d[['init', 'size', ]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "e460502f-1604-4221-bfe4-ef41bcef74d9",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Counting up events to define number of parameters later\n",
    "event_times = [0, ] + list(np.unique(d.loc[d['delta'] == 1, 'time'])) + [59, ]\n",
    "event_times_a1 = list(np.unique(d.loc[(d['delta'] == 1) & (d['novel'] == 1), 'time']))\n",
    "event_times_a0 = list(np.unique(d.loc[(d['delta'] == 1) & (d['novel'] == 0), 'time']))\n",
    "params_rd = len(event_times)\n",
    "params_plr_a1 = len(event_times_a1)\n",
    "params_plr_a0 = len(event_times_a0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "af27e483-1dd0-4025-9d19-5dd1fdf91f1a",
   "metadata": {},
   "outputs": [],
   "source": [
    "def psi_rd(theta):\n",
    "    # Extracting parameters\n",
    "    rds = theta[:params_rd]\n",
    "    idPLR = params_rd + W.shape[1] + params_plr_a1\n",
    "    beta1 = theta[params_rd: idPLR]\n",
    "    beta0 = theta[idPLR:]\n",
    "\n",
    "    # Nuisance models\n",
    "    ee_plog1 = ee_plogit(beta1, t=t, delta=y, X=W, unique_times=event_times_a1)\n",
    "    ee_plog1 = ee_plog1 * (a == 1)[None, :]\n",
    "\n",
    "    ee_plog0 = ee_plogit(beta0, t=t, delta=y, X=W, unique_times=event_times_a0)\n",
    "    ee_plog0 = ee_plog0 * (a == 0)[None, :]\n",
    "\n",
    "    # Predictions to get risk differences\n",
    "    risk1 = plogit_predict(theta=beta1, delta=y, t=t, X=W,\n",
    "                           times_to_predict=event_times, \n",
    "                           measure='risk', \n",
    "                           unique_times=event_times_a1)\n",
    "    risk0 = plogit_predict(theta=beta0, delta=y, t=t, X=W,\n",
    "                           times_to_predict=event_times, \n",
    "                           measure='risk', \n",
    "                           unique_times=event_times_a0)\n",
    "    ee_rd = (risk1 - risk0) - np.asarray(rds)[:, None]\n",
    "\n",
    "    # Returning stacked estimating equations\n",
    "    return np.vstack([ee_rd, ee_plog1, ee_plog0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "53c39f08-e02f-48ed-935b-6abb66842365",
   "metadata": {},
   "outputs": [],
   "source": [
    "inits = ([0., ]*params_rd\n",
    "         + [0., ]*W.shape[1] + [-4., ] + [0., ]*(params_plr_a1 - 1)\n",
    "         + [0., ]*W.shape[1] + [-4., ] + [0., ]*(params_plr_a0 - 1))\n",
    "estr = MEstimator(psi_rd, init=inits)\n",
    "estr.estimate()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0982697e-b69c-4f9d-b03d-2f46ecf73642",
   "metadata": {},
   "source": [
    "We can then examine the table for the causal risk difference at 59 months."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "409409de-9e90-4a3b-ab3f-41d15c30b8d1",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "==============================================================\n",
      "              Estimation Method: M-estimator\n",
      "--------------------------------------------------------------\n",
      "No. Observations:          86 | No. Parameters:             54\n",
      "Solving algorithm:         lm | Max Iterations:           5000\n",
      "Solving tolerance:      1e-09 | Allow P-Inverse:             1\n",
      "Derivative Method:     approx | Deriv Approx:            1e-09\n",
      "Small N Correction:      None | Distribution:           Z-stat\n",
      "==============================================================\n",
      "   Theta   StdErr  Z-score      LCL      UCL  P-value  S-value \n",
      "--------------------------------------------------------------\n",
      "   -0.19     0.12    -1.63    -0.42     0.04     0.10     3.28 \n",
      "==============================================================\n"
     ]
    }
   ],
   "source": [
    "estr.print_results(subset=[params_rd-1])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e56a9319-846a-440c-86bc-d29163c1efe1",
   "metadata": {},
   "source": [
    "Perhaps more informative is to examine a plot of the results. In this plot, inference is implicitly for the underlying causal risk difference *function*, so rather than confidence *intervals* the confidence *bands* are presented. The confidence bands can be straightforwardly computed using `delicatessen`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "5fd68347-0d47-499d-b389-90e562053ff1",
   "metadata": {},
   "outputs": [],
   "source": [
    "cb = estr.confidence_bands(method='supt', subset=[x for x in range(1, params_rd)], seed=7777)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "d6ee347a-7378-4ac5-abab-99c9061f83cd",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Formatting RD results for figure\n",
    "rd = estr.theta[:params_rd]\n",
    "rd_ci = estr.confidence_intervals()[:params_rd, :]\n",
    "rd_results = pd.DataFrame()\n",
    "rd_results['time'] = event_times\n",
    "rd_results['rd'] = rd\n",
    "rd_results['lcb'] = [0.,] + list(cb[:, 0])\n",
    "rd_results['ucb'] = [0.,] + list(cb[:, 1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "490a4bde-3edb-41ab-ab41-575641f24037",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'Risk Difference')"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjUAAAFSCAYAAAAKDC7nAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjgsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvwVt1zgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAOxVJREFUeJzt3QucjGX7wPFrnc+LHJZyjBxyjGjlLW82x4/wShQ55PAJKVEOFUISlbeUHDo4vCnixVu9JaJSrEOESKKI5JhYh1qy8/9c9///zH9md2Z2ZmdmZ+eZ3/fzeezOM888M/PMY55rr/u67zvO4XA4BAAAIMrlivQLAAAACAWCGgAAYAsENQAAwBYIagAAgC0Q1AAAAFsgqAEAALZAUAMAAGyBoAYAANgCQQ0AALAFghoAAGALURXUrF+/Xjp06CDly5eXuLg4WblyZaaP+fzzz+Wmm26S/PnzS7Vq1WT+/PkZtpk5c6ZUrlxZChQoIE2bNpUtW7aE6R0AAIBwiaqg5uLFi1K/fn0ThPjj4MGD0r59e/n73/8uO3bskGHDhkn//v3lk08+cW6zZMkSGT58uIwfP162b99u9t+6dWs5efJkGN8JAAAItbhondBSMzUrVqyQTp06ed1m1KhR8t///ld2797tXNe9e3c5e/asrFq1ytzWzMzNN98sr776qrmdlpYmFSpUkKFDh8ro0aOz4Z0AAIBQyCM2lpycLElJSW7rNAujGRt1+fJl2bZtm4wZM8Z5f65cucxj9LHepKammsWigdCZM2fkmmuuMcEWAADwj+ZWzp8/b0pL9BocDFsHNcePH5eyZcu6rdPbKSkp8scff8jvv/8uV69e9bjN999/73W/U6ZMkQkTJoTtdQMAEGuOHDki1113XVD7sHVQEy6a2dE6HMu5c+ekYsWK5gMpVqxYRF8bAADRRBMNWvZRtGjRoPdl66AmISFBTpw44bZOb2vgUbBgQcmdO7dZPG2jj/VGe1Lpkp7ul6AGAIDAhaJ8I6p6PwUqMTFR1q5d67ZuzZo1Zr3Kly+fNGrUyG0brY/R29Y2AAAgOkRVUHPhwgXTNVsXq8u2/n748GFns1CvXr2c2z/44IPy008/yciRI02NzGuvvSbvvfeePProo85ttBnp9ddflwULFsjevXtl0KBBput43759I/AOAQBATDQ/ff3112bMGYtV19K7d28zqN6xY8ecAY6qUqWK6dKtQczLL79sCpDeeOMN0wPK0q1bNzl16pSMGzfOFBY3aNDAdPdOXzwMAABytqgdpyanFTnFx8ebgmFqagAAiMw1NKqanwAAALwhqAEAALZAUAMAAGyBoAYAANgCQQ0AALAFghoAAGALBDUAAMAWCGoAAIAtENQAAABbIKgBAAC2QFADAABsgaAGAADYAkENAACwBYIaAABgCwQ1AADAFghqAACALRDUAAAAWyCoAQAAtkBQAwAAbIGgBgAA2AJBDQAAsAWCGgAAYAsENQAAwBaiLqiZOXOmVK5cWQoUKCBNmzaVLVu2eN22RYsWEhcXl2Fp3769c5s+ffpkuL9NmzbZ9G4AAECo5JEosmTJEhk+fLjMnj3bBDQvvfSStG7dWvbt2ydlypTJsP3y5cvl8uXLztu//fab1K9fX7p27eq2nQYx8+bNc97Onz9/mN8JAACI6UzN9OnTZcCAAdK3b1+pXbu2CW4KFSokb731lsftS5YsKQkJCc5lzZo1Zvv0QY0GMa7blShRIpveEQAAiLmgRjMu27Ztk6SkJOe6XLlymdvJycl+7ePNN9+U7t27S+HChd3Wf/755ybTU6NGDRk0aJDJ6PiSmpoqKSkpbgsAAIisqAlqTp8+LVevXpWyZcu6rdfbx48fz/TxWnuze/du6d+/f4amp4ULF8ratWtl6tSp8sUXX0jbtm3Nc3kzZcoUiY+Pdy4VKlQI4p0BAICYq6kJhmZp6tatK02aNHFbr5kbi95fr149uf766032pmXLlh73NWbMGFPbY9FMDYENAACRFTWZmlKlSknu3LnlxIkTbuv1ttbB+HLx4kVZvHix9OvXL9PnqVq1qnmuAwcOeN1Ga3CKFSvmtgAAgMiKmqAmX7580qhRI9NMZElLSzO3ExMTfT526dKlpg6mZ8+emT7PL7/8YmpqypUrF5LXDQAAskfUBDVKm3xef/11WbBggezdu9cU9WoWRntDqV69epmmIU9NT506dZJrrrnGbf2FCxfk8ccfl02bNsmhQ4dMgNSxY0epVq2a6SoOAACiR1TV1HTr1k1OnTol48aNM8XBDRo0kFWrVjmLhw8fPmx6RLnSMWy++uorWb16dYb9aXPWrl27TJB09uxZKV++vLRq1UomTZrEWDUAAESZOIfD4Yj0i4h2WiisvaDOnTtHfQ0AABG6hkZV8xMAAIA3BDUAAMAWCGoAAIAtENQAAABbIKgBAAC2QFADAABsgaAGAADYAkENAACwBYIaAABgCwQ1AADAFghqAACALRDUAAAAWyCoAQAAtkBQAwAAbIGgBgAA2AJBDQAAsAWCGgAAYAsENQAAwBYIagAAgC0Q1AAAAFsgqAEAALZAUAMAAGwh6oKamTNnSuXKlaVAgQLStGlT2bJli9dt58+fL3FxcW6LPs6Vw+GQcePGSbly5aRgwYKSlJQk+/fvz4Z3AgAAYjaoWbJkiQwfPlzGjx8v27dvl/r160vr1q3l5MmTXh9TrFgxOXbsmHP5+eef3e6fNm2azJgxQ2bPni2bN2+WwoULm33++eef2fCOAABATAY106dPlwEDBkjfvn2ldu3aJhApVKiQvPXWW14fo9mZhIQE51K2bFm3LM1LL70kTz31lHTs2FHq1asnCxculF9//VVWrlyZTe8KAADEVFBz+fJl2bZtm2kesuTKlcvcTk5O9vq4CxcuSKVKlaRChQomcNmzZ4/zvoMHD8rx48fd9hkfH2+atXztEwAA5DxRE9ScPn1arl696pZpUXpbAxNPatSoYbI4//nPf+Ttt9+WtLQ0adasmfzyyy/mfutxgexTpaamSkpKitsCAAAiK2qCmqxITEyUXr16SYMGDeT222+X5cuXS+nSpWXOnDlB7XfKlCkmo2MtmgUCAACRFTVBTalSpSR37txy4sQJt/V6W2tl/JE3b15p2LChHDhwwNy2HhfoPseMGSPnzp1zLkeOHMnCOwIAADEZ1OTLl08aNWoka9euda7T5iS9rRkZf2jz1bfffmu6b6sqVaqY4MV1n9qUpL2gfO0zf/78pleV6wIAACIrj0QR7c7du3dvady4sTRp0sT0XLp48aLpDaW0qenaa681zUNq4sSJcsstt0i1atXk7Nmz8vzzz5su3f3793f2jBo2bJg888wzUr16dRPkjB07VsqXLy+dOnWK6HsFAAA2Dmq6desmp06dMoPlaSGv1sqsWrXKWeh7+PBh0yPK8vvvv5su4LptiRIlTKZn48aNpju4ZeTIkSYwGjhwoAl8mjdvbvaZfpA+AACQs8U5dLAWBEWbrLRgWOtraIoCACAy19CoqakBAADwhaAGAADYAkENAACwBYIaAABgCwQ1AADAFghqAACALRDUAAAAWyCoAQAAtkBQAwAAbIGgBgAAxG5Q89dff8mnn34qc+bMkfPnz5t1v/76q1y4cCHUrw8AACA8E1rqLNdt2rQxk0empqbKnXfeKUWLFpWpU6ea27Nnzw50lwAAANmfqXnkkUekcePGZgbsggULOtd37txZ1q5dG/wrAgAAyI5MzZdffikbN26UfPnyua2vXLmyHD16NCuvAQAAIPszNWlpaXL16tUM63/55RfTDAUAABAVQU2rVq3kpZdect6Oi4szBcLjx4+Xdu3ahfr1AQAA+CXO4XA4JACakWndurXow/bv32/qa/RnqVKlZP369VKmTBmJNSkpKRIfHy/nzp2TYsWKRfrlAAAQk9fQgIMaq0v3kiVLZOfOnSZLc9NNN0mPHj3cCodjCUENAABRGtQgfB8IAACxJCWE19CAa2qmTJkib731Vob1uk7HqgEAAIiEgIMaHUW4Zs2aGdbfeOONDLwHAACiJ6g5fvy4lCtXLsP60qVLy7Fjx0L1ugAAAMIb1FSoUEE2bNiQYb2uK1++fKC7AwAAiExQM2DAABk2bJjMmzfPzAOli9bTPProo+a+cJs5c6YZvbhAgQLStGlT2bJli9dtX3/9dfnb3/4mJUqUMEtSUlKG7fv06WPG2nFddG4rAABg82kSHn/8cfntt99k8ODBcvnyZbNOA4xRo0bJmDFjJJy0G/nw4cNN7Y4GNDoIoI6Zs2/fPo/j43z++edy7733SrNmzcxr1EJmHTxwz549cu211zq30yBGgzRL/vz5w/o+AABA6GW5S7eOT7N3714zNk316tWzJRDQQObmm2+WV1991TllgzaHDR06VEaPHp3p43V6B83Y6ON79erlzNScPXtWVq5cmeXXRZduAACisEu3pUiRIibAqFOnTrYENJoV2rZtm2lCsuTKlcvcTk5O9msfly5dkitXrkjJkiUzZHQ001OjRg0ZNGiQyUT5kpqaaj4E1wUAAERZ89PFixflueeek7Vr18rJkydNtsTVTz/9JOFw+vRpk2kpW7as23q9/f333/u1D20i02Jm18BIm57+8Y9/SJUqVeTHH3+UJ554Qtq2bWsCpdy5c3sdq2fChAlBviMAABDRoKZ///7yxRdfyP3332+6dmthbTTQQGzx4sUmK6P1NZbu3bs7f69bt67Uq1dPrr/+erNdy5YtPe5La4e0tseimRptBssqnXbCVyugHuM8eQL+qAAAiCkBXyk//vhj+e9//yu33nqrZCedMFMzJydOnHBbr7cTEhJ8PvaFF14wQc2nn35qghZfqlatap7rwIEDXoMabW4LVZObBjRHjhwxzWLe5M2b1wRNBDYAAISwpkYLbdPXpGSHfPnySaNGjUyzl0WbvvR2YmKi18dNmzZNJk2aJKtWrTIzivszC7nW1HgaYDAcNEOjAY0GbFaw5Lroer2fKboAAAhxUKMBwrhx40zRbXbTJh8de2bBggWm55UW9WqNT9++fc392qPJtVu5duEeO3asGUdHx7bR0ZB10Z5bSn9qF/VNmzbJoUOHTIDUsWNHqVatmukqnp00C6MZmfQL2RkAAMLU/PTiiy+aglot0NVAQS+8rrZv3y7h0q1bNzl16pQJqjQ4adCggcnAWMXDhw8fNj2iLLNmzTK9pu6++263/YwfP16efvppkwXZtWuXCZK0W7cWEes4Nhq4xcpYNdTzAABidpyazHr9aMAQa4LpY69NSwcPHjRBVPoA0bpfu5Br7yxP9weDeh4AgJ3GqQk4UxOLQUs0Z1t89ZxyrefxtI3ul3oeAEC0yFI/YW2qWbZsmWmG0poULRzWZidtBnKdfgDhl1m2xZ+eU1Y9jyc6NpA+R1bQFR3pz9VgCt45nwCEPKjRGhQdvE5TRVpcq5NYalCzfPlyU9OycOHCQHeJIPjKtliZFm/ZlsyCFb2IaE2SBk1ZQVd0+8lqVtCfps7McD4BCHlQoz2QdL4k7SpdtGhR5/p27drJfffdF+juECKesi16kfnjjz98BiUatLgORph+n9q+mZW/rmm6sp9gsoKZNXX689w0hQLITMDfLlu3bpU5c+ZkWK/NTtojCTmHP0GJBjSZNU1llTZdIbaygpkFwL6aOjPD+QQgMwFfsbSXjqcJHH/44QcpXbp0oLtDmDHODcJxTnkKTLIj6PDWZEq9DYAsBTV33XWXTJw4Ud577z3nl4nW0uhkkV26dOGowq+LEBei2Pq8s1ps7m99F/U2ALI8+J4OZlemTBlTr3H77bebZiedqmDy5Mkc1Si7WISLP0XGXIiiqwdSZgFqZp+3r/qtYJpSqbcBkOWgRns9rVmzRjZs2CA7d+40Uw3cdNNNpkcUAr+QhKIHUjAXi0jV83Ahyn6h6IHk7VwLRf1WZnw9lnobAOZ7IpDDoF+GBQsWlB07dphZurN7pm67XkiC7YEU7MUiXDJ7TcGMgaOoowhMsD2QMjvXIn0O5tSsZSTP42DHBgJCJbu+rwN6Bm0uqFixIn8VhfhCEs4eSDlVsGPg5OTmq5w+n1YwPZDsei5FUrjO41Bk5oBQya7v64D3/uSTT8oTTzwh//rXv8yge4jNC0mwghkDJyc3XzGfVvSdS5EUzvM4FJk5IBSy8/s64DP91VdflQMHDpgZrStVqiSFCxfOtlm6YS/BftHmxDoK5tOKjGi+aIf7POYPKuQE2fV9HfA3QadOncLzSgAbCdd8WrFYcwIA/mKWbkS1SF2os1oXE+n6j5zYUw4AQoVZuhGVIh0cZLXoLdL1Hzm1pxwAhAKzdCMqRTI4CLbojaACns6pUKO5EbGIWboRtSIZHOTEImVEn3BnHGluRKxhlm4ghH8J89cxclLGkeZGxBpm6QZC/Nc1fx0joC9hapyAkMmV1Vm6rVEqmaUbsfrXdf78+T0ueh8XKgCIgqBGZ+nWSSxdZ+muVq2aFC1alFm6s0DTzpcuXfK5RONIqbEyDo2nhYAGACKDWbqzKWjRwtL0g7HpfZ07d5Y9e/b43MfNN98sK1asMFkxAAAQRKZG53g6ffq0+f2BBx6Q8+fPmxm6Bw8eLCNHjpSkpCTJLjNnzpTKlSubArimTZvKli1bfG6/dOlSqVmzptm+bt268tFHH2UILMaNGyflypUzM5Dre9m/f39IXqvuu0WLFtKwYUOpXbu2VK9e3W254YYbMg1o1NatW01WDAAABBnUaOFjSkqK+X3BggXy559/SiQsWbJEhg8fLuPHjzdzTNWvX19at24tJ0+e9Lj9xo0b5d5775V+/frJN998Y6Z40GX37t3ObaZNmyYzZsyQ2bNny+bNm81cVrrPULxHzdAkJydnut2NN94oP/zwgwmmXJedO3cG/RoAAIgVcQ4/CjbuvPNOOXHihDRq1MgENd26dTNZDU/eeustCRfNzGhTjE6qqdLS0syorkOHDpXRo0dn2F5f58WLF+XDDz90rrvlllukQYMGJojRt64Tc44YMUIee+wxc/+5c+ekbNmyMn/+fOnevbtfr0sDvvj4eDPRpz7Waib6/fffzetTX3/9tdnGCnaUZo9y5cpljqU1oJvetoax1+00o6P27dsnRYoUMb/rdum3VZrN0fekxao6M6/S/WpQqq/J9TPLyrb58uVz1otoc1pqampQ22rgqJ+hVYsS6LZ62wo+CxUq5NxWH6/70efX1xHotvr6rcyYvgbr89Rjo8cokG312OoxtqT/7APdNpDP3s7nibfPM9jzxPXzDPY88fZ5ZvU88fZ5BnueRNtnz3dE9H1HXLlyxXy+CQkJznPd9bPXxIRei/X6qx0twp6pefvtt6Vdu3amQFjpE+sF29MSLnqAtm3b5tbUpR+E3vaWDdH16ZvGNAtjbX/w4EE5fvy42zYaeGjw5CvDoh+OBjKui7IKpjX40MUKaKwvSGvRwEqzTHq89LZ+8BpEaQBjBVfpacBkee+998y22vznSpu6dP23337rXPf++++bdX379nXbVj9PXa/ZKcunn35q1qUP5rp06WLWf/755851GzZsMOu0N5yr+++/36z/+OOPnes0q6br0n8WAwYMMOuXL1/uXLd3716zrnnz5m7bPvzww2b9okWLnOsOHTpk1mmw7WrUqFFm/Ztvvulcp0G5rqtVq5bbthMmTDDrX3nlFec6/TytJkLXcWemTp1q1ulPi95vbWudB0r3p+t0/670+XW9vh6Lvk5dp6/blb4vXa/v06LvX9fp8XClx0vX6/Gz6HHVdXqcXennoOv1c7Ho56Xr9PNzpZ+vrtfP26Lnga7T88KVnje6Xs8ji55fuk7PN1d6Pup6PT8tet7qOj2PXel5ruv1vLdoJlPX6f9VV/r/R9fr/yfL0aNHzbp69eq5bfvUU0+Z9bNmzXKuO3PmjPPzdDV58mSzbvr06W5f5ta2rs3Duo2u08e4srbV57Doc+s6fS2u9LXqen3tFm/fEXoMdL1rsznfEf+L74ic9R3Rpk0bc210TTR8+eWXJqDJ1kJhzT4899xz5vcqVarIv/71L7nmmmskO2lNj0Z0+lrSv7bvv//e42M0YPG0va637rfWedvGkylTpmS4WPmix8r66w4AAESw+UkLhbXmo1SpUqZQ+OWXXzYZiez066+/yrXXXmvqZBITE53rtVD5iy++cMs4WDSQ0OYyrauxvPbaayYg0b+UdV9a8Kz71kJhyz333GOyJ1rD4y1To4tF/0LXrMyuXbvcmp80u6TpU00RakrNSp/6mzKk+YnmJ3/OE0Xz0/+i+YnmJ5qor8R081OeQAqFNajRIEHT79kd1Ohz68FwTdsrva0HyhNd72t766eucw1q9LbW3XhjDbKWXqVKlTJ8IPphajOXK9e2eovrf0RPrIuar2091TnpieNp3JRgt9XPwtP7CGRb1/9EWdlWj4mnbT19NoFsq/8RPW2rQXL6jFsg26pgtw3ks7fzeeLt8wz2PPH2eQZ7noTzsw/2PIm2z57viOj+jsib7vn089QOOqHiV1CjmRHtNaRt/BqRaXt+dhcK6398ff61a9ea16K0gExvP/TQQ15ft94/bNgw57o1a9Y4Mz3alKaBjW5jBTEavGnWZ9CgQWF5HwAAIDzy+Fso/M9//lN+/PFH89eJpogi0a1bu3P37t1bGjduLE2aNJGXXnrJ9G6yimB79eplmqi05kU98sgjZsRjHQW5ffv2snjxYtMLae7cueZ+fS8a8DzzzDOmqEmDnLFjx5o0mBU4AQCA6BA1hcJWF+1Tp06ZwfK0kFezK6tWrXIW+h4+fNitmaZZs2byzjvvmJ4FTzzxhAlcVq5cKXXq1HGrydHAaODAgXL27FlTIa779JT2BAAAUV4oDP/GqfFU5GTV1GibvK+aGU9cC4W1u6an9lPXsTI8cR0TI9DHZvZ4AAAyYxUKa1LE03XQ1zU0LJkaHXFXMxmavdDffUk/fgbCR4MSbSbTJrVA543y57G+Hg8AQFRmajS60oufNjnp7153FhcnP/30k8SaSGVqXO/3Radb8PRYHQDQH96yRAAARF2mxrVLcvruycgZ0gcurkFLZsFLZkGPNa5OoGi6AgBkJ7+CGgTPdbh9TxkuT336A2FNweAaUGjTkc7w7Ytuoxk4X81L/mZ0PO2bpisAQHYJ6EqqvYR04D2dL0Ln1NALoaaT7r77bjMfCU0UGekx0nSbpt905ERvgxtqyi3YwCb982pAkdVCYH+DIl/0sfr8nBcAgOzg91VUL7w65svu3bulbdu20qFDB1NsqpNj6cRtOtnV+vXrA64bsTsNVHQKBW+lS5rBOXLkiNf7g+FtxNNQBkWeBFKvAwBAtgc1OpvsL7/8YuovatSo4XafTiipM+vOnj1bhg4dGrIXZxehyMB4qmvJaq1LdgRFAABkN7+vttrkpKPtpg9oVM2aNeXJJ5+UZcuWEdSECZkPRDN/xkTKDIXnAEIW1Hz33XcmG+PN3//+d5k4caK/u4MfAin29TYXV6T5yiZxkYoN/o6JlBkKzwGELKjRKQR8TY2g92kfc0iOKfbN6RkmLlLRJaujV+tjgg1oFIXnAEIW1OiM2DpFuDc655K33j2IrboWfzNMXKRiY/TqzMZEygyF5wBCHtTol1rLli29Fr36GocFsSWzDFMoBvbL6Rkqu/En2+JPkJp+PKVAeTtfOBcABBTUjB8/PtNtunTpwlFFQBmmYAqgab7KeaNXh5u35+FcABDyoAbwRygG9lM0X2VfDyTXDImvbEs4hh7w53zhXACgmCYB2S6Ygf0UNRaR64GUmXBkbHydL5wLAFwR1CAiorEAOlTZkEjUf4SqB5Kn4QOyY+gBf86XcA9GGS7UAwGhQ1ADZHM2JNL1H1npgeTrApxThh6I1gEqw3U+hGLAQyAUdO7D1NTUsEwHlB5BTYzM4m1Xkfrr3NMF2t9sSKTrP4LtgZSTMm+hqs+KpHCcD9nV3AgE4vfff5d8+fJJOAV8pdT5n6677jqP923atEluueWWULyumBGpWbztIlJ/nWf217WnbEiourJnRbQ2zYS7PiuSwlkPFKrmRiDaBHyVbNWqlXz11VdSsmRJt/UbNmyQ9u3bm5GHER2zeEernPDXeWZ/XWeWDYnWppKcyK71WTmhuREIZfNTdpyHAQc1monRwOazzz6TokWLmnXr16+XDh06yNNPPx2O12h7ZGCi569zX9mWzLIhOSEYy8nzhCF6mhuBQIManZEgO+oIAw5q3njjDbn77rtNEPPJJ5/Ixo0b5a677pJnnnlGHnnkkfC8SiAH/nUeaLYlJzSV0NMGgJ3lCvgBuXLJ4sWLTR3IHXfcYQKaKVOmhD2gOXPmjPTo0cPUlhQvXlz69esnFy5c8Ln90KFDpUaNGuaLvGLFivLwww9nmHRTLzTpF31/gK9sS1azIVYwFqmFaSUA2JlfmZpdu3ZlWKdNTffee6/07NlTbrvtNuc29erVC/2rFDEBzbFjx2TNmjUmldW3b18ZOHCgvPPOOx63//XXX83ywgsvSO3ateXnn3+WBx980KxbtmyZ27bz5s2TNm3aOG9r0ARkNdtCNgSRLuS2a2E4kJk4hx8VqJqd0S9z101db1u/689wzNS9d+9eE5hoLULjxo3NulWrVkm7du1Mb6zy5cv7tZ+lS5eaIOzixYvOOhbrIqXdH7MqJSVF4uPjTRZIM0mhpAHcwYMHJX/+/CY7BsAeNPCoXr162J9n//79EW+qRWy78n+FwlWqVPF4HQvlNdSvTI1eVCMpOTnZZE+sgEYlJSWZYGvz5s3SuXNnv/ZjHbD0hblDhgyR/v37S9WqVU02R7NApOkBhFN2FI5TGI5Y41dQU6lSJYmk48ePS5kyZdzWaWCi3cr1Pn+cPn1aJk2aZJqsXE2cONHUBulfMqtXr5bBgwebWh2tv/FGI05dXKNMAMhpheM0hSLWBNz7acGCBVKqVCkzJo0aOXKkzJ071zQPvfvuuwEFQKNHj5apU6dm2vQULA069PXqa0zf7Xzs2LHO3xs2bGiapp5//nmfQY0WRk+YMCHo1wUgtuWEXnxATPd+evbZZ509O7RZ6NVXX5Vp06aZQOfRRx8NaF8jRowwQYuvRZuEEhIS5OTJkxkGpdMeTnqfL+fPnzdFwDqmjv5VlFldStOmTU2djmsmJr0xY8aYpixr0cHxAABAlGVq9AJerVo18/vKlSvNmDXapHPrrbdKixYtAtpX6dKlzZKZxMREM1Lxtm3bpFGjRmbdunXrJC0tzQQhvjI0rVu3NkW277//vhQoUCDT59qxY4eUKFHCPMYbvc/X/QAAIAoyNUWKFJHffvvN/K41KHfeeaf5XQOGcLUN16pVy2RbBgwYIFu2bDFTMjz00EPSvXt3Z8+no0ePSs2aNc39VkCjIx9rc9Kbb75pbmv9jS5WD60PPvjADCa4e/duOXDggMyaNctkonR8GwAAYPNMjQYx2lNI609++OEH061a7dmzRypXrizhsmjRIhPItGzZ0vR66tKli8yYMcOty9i+ffuc4zNs377d9IxSVmbJtTeXvlZtipo5c6ZpNtMu6brd9OnTTfAEAABsOE6NK20Geuqpp0wz1KBBg5yD1o0fP95MKf7kk09KrGGcGgAAIj9OTcBBDTIiqAEAIEoG39MpEOrUqWOafTxNmeAqXNMkwDvtCeYtNtUuo8wCDgCIBX4FNQ0aNHAOgKe/e5syIVzTJMB3QKNRrjb9eXL58mWPoygDAGA3fk+TYHW9jvSUCXCnwaQGNBUqVMgQuGjAo7VPmvrLaisjmR4AgG2nSfA1YnA4h/uGbxrQpG+r1IBE12lQk9UMGpkeAEC0CEmbhBYA6cjCOr2Av3MxIXsCHc3gZDVLY2V6qCUHANgqqNHARedNWrNmjWnu0DmfOnXqJPPmzTPduHPnzh3wNAkIv1DU0mhwkxU0XQEAspPfV7xx48bJnDlzJCkpSTZu3Chdu3aVvn37yqZNm8yAdXpbAxuEh7fAIqsBhz+Cbb6i6QoAkCODmqVLl8rChQvlrrvuMtMKaNdtvaDu3LnTXPwQucBC7w/HZxBM8xVNVwCAHBvU6MzV1mSSOmaNTuiozU0ENOHlT2ARzmYeuoIDAKKF31dCzRK4joWiFzud3BLhR2ABAEAIgxrNFPTp08dkaNSff/4pDz74oBQuXNhtu+XLl/u7SwAAgOwPanr37u12u2fPnqF7FbAtX4XM9I4CAEQkqNGu20AoC5zpHQUACCUmBEJECpzpHQUACDWCGkS0wDmYcXZovgIAuCKoQUQwL1Xs0QA2mCk3CGIBZIagBhHBvFSxF9CkpKS4DQsRKGqwAGSGoAYRw/g79sq2+Mqk6GM0oNE6q6x87tRgAfAHQQ2AkGRb/Mmk6H3aKw4AwoGgBoBffGVbsiuT4q2wnHobAIqgBkBAIpFtyaywnHobAOb7icOAaBZMl/Bg+MoMZNbLh6xCaAvLqbcBYMklUeLMmTPSo0cP02ZfvHhx6devn1y4cMHnY1q0aGEuIK6Lzlfl6vDhw9K+fXspVKiQlClTRh5//PGIXSgR+F/u+ld7ampqti9aW+LpPLHqTrLy2Oygz6vZjqws/rxmb/sPxfu1MkTpFwrOAURdpkYDmmPHjsmaNWvMl2Tfvn1l4MCB8s477/h83IABA2TixInO2xq8WPSCqAFNQkKCbNy40ey/V69e5ovy2WefDev7QWS7hAfDV2Ygs14+kcwqhKJbtf7f0IAyK+MOeXtsqETrHyNk7oAYC2r27t0rq1atkq1bt0rjxo3NuldeeUXatWsnL7zwgpQvX97rYzWI0aDFk9WrV8t3330nn376qZQtW1YaNGggkyZNklGjRsnTTz8d1Jc/wi/Sf6F7y9Tk1F4+wXar9nUB9ifIDNfFOxQDOUYS9UBAjAU1ycnJpsnJCmhUUlKS5MqVSzZv3iydO3f2+thFixbJ22+/bQKbDh06yNixY53ZGt1v3bp1TUBjad26tQwaNEj27NkjDRs29LhPqxnBon/9InZkdhH1JyMRiaxCuAOuSAWZkczaBYt6ICAGg5rjx4+bepf0X2QlS5Y093lz3333SaVKlUwmZ9euXSYDs2/fPlm+fLlzv64BjbJu+9rvlClTZMKECUG+K0SrzC6ivjISkc4qhLsJKFazdnadmgIIhez8Iy6i3wSjR4+WqVOnZtr0lFVac2PRjEy5cuWkZcuW8uOPP8r111+f5f2OGTNGhg8f7pap0YscYkdWL6KRzipQvxE7QlFDBUTbH1QRDWpGjBghffr08blN1apVTdPRyZMnM/yH1R5R3uplPGnatKn5eeDAARPU6GO3bNnits2JEyfMT1/7zZ8/v1mArCCrgGipoQKi7Q+qiJ7ppUuXNktmEhMT5ezZs7Jt2zZp1KiRWbdu3TpJS0tzBir+2LFjh/mpGRtrv5MnTzYBk9W8pb2rtNt47dq1s/iuACDnyIlF60BMj1NTq1YtadOmjemerZmVDRs2yEMPPSTdu3d39nw6evSo1KxZ05l50SYm7cmkgdChQ4fk/fffN921b7vtNqlXr57ZplWrViZ4uf/++2Xnzp3yySefyFNPPSVDhgwhEwMAQJSJiqDG6sWkQYvWxGhX7ubNm8vcuXOd92vhpRYBX7p0ydzWtKt21dbARR+nTV1dunSRDz74wPmY3Llzy4cffmh+atamZ8+eJvBxHdcGAABEhzgHpfFB02K8+Ph4OXfunGm6AgB/6B9jBw8eNH9YhbreQOsOtYddlSpVaH5CzFxDqR4DgAgJdxd/u3bhB7whqAGACAl3F3+68CPWENQAQATR3RqIwUJhAAAAXwhqAACALRDUAAAAWyCoAQAAtkBQAwAAbIGgBgAA2AJBDQAAsAWCGgAAYAsENQAAwBYIagAAgC0Q1AAAAFsgqAEAALZAUAMAAGyBoAYAANgCQQ0AALAFghoAAGALBDUAAMAWCGoAAIAtENQAAABbIKgBAAC2EDVBzZkzZ6RHjx5SrFgxKV68uPTr108uXLjgdftDhw5JXFycx2Xp0qXO7Tzdv3jx4mx6VwAAIFTySJTQgObYsWOyZs0auXLlivTt21cGDhwo77zzjsftK1SoYLZ3NXfuXHn++eelbdu2buvnzZsnbdq0cd7WoAkAAESXqAhq9u7dK6tWrZKtW7dK48aNzbpXXnlF2rVrJy+88IKUL18+w2Ny584tCQkJbutWrFgh99xzjxQpUsRtvQYx6bcFAADRJSqan5KTk03gYQU0KikpSXLlyiWbN2/2ax/btm2THTt2mGar9IYMGSKlSpWSJk2ayFtvvSUOh8PnvlJTUyUlJcVtAQAAkRUVmZrjx49LmTJl3NblyZNHSpYsae7zx5tvvim1atWSZs2aua2fOHGi3HHHHVKoUCFZvXq1DB482NTqPPzww173NWXKFJkwYUIW3w0AALBdpmb06NFei3mt5fvvvw/6ef744w9Te+MpSzN27Fi59dZbpWHDhjJq1CgZOXKkqbvxZcyYMXLu3DnncuTIkaBfIwAAiOJMzYgRI6RPnz4+t6lataqpdzl58qTb+r/++sv0iPKnFmbZsmVy6dIl6dWrV6bbNm3aVCZNmmSamPLnz+9xG13v7T4AABCDQU3p0qXNkpnExEQ5e/asqYtp1KiRWbdu3TpJS0szQYg/TU933XWXX8+ldTclSpQgaAEAIMpERU2N1sJol+sBAwbI7NmzTZfuhx56SLp37+7s+XT06FFp2bKlLFy40BT8Wg4cOCDr16+Xjz76KMN+P/jgAzlx4oTccsstUqBAAdNd/Nlnn5XHHnssW98fAACIkaBGLVq0yAQyGrhor6cuXbrIjBkznPdroLNv3z7TzORKezNdd9110qpVqwz7zJs3r8ycOVMeffRR0+OpWrVqMn36dBM8AQCA6BLnyKz/MjKlXbrj4+NN0bCOeAwAALL/GhoV49QAAABkhqAGAADYAkENAACwBYIaAABgCwQ1AADAFghqAACALRDUAAAAWyCoAQAAtkBQAwAAbIGgBgAA2AJBDQAAsAWCGgAAYAsENQAAwBYIagAAgC0Q1AAAAFsgqAEAALZAUAMAAGyBoAYAANgCQQ0AALAFghoAAGALBDUAAMAWCGoAAIAtENQAAABbiJqgZvLkydKsWTMpVKiQFC9e3K/HOBwOGTdunJQrV04KFiwoSUlJsn//frdtzpw5Iz169JBixYqZ/fbr108uXLgQpncBAAAk1oOay5cvS9euXWXQoEF+P2batGkyY8YMmT17tmzevFkKFy4srVu3lj///NO5jQY0e/bskTVr1siHH34o69evl4EDB4bpXQAAgHCJc2g6I4rMnz9fhg0bJmfPnvW5nb6t8uXLy4gRI+Sxxx4z686dOydly5Y1++jevbvs3btXateuLVu3bpXGjRubbVatWiXt2rWTX375xTzeHykpKRIfH2/2rxkfAAAg2X4NzSM2dfDgQTl+/LhpcrLoQWvatKkkJyeboEZ/apOTFdAo3T5Xrlwms9O5c2eP+05NTTWLRT8I64MBAAD+s66docix2Dao0YBGaWbGld627tOfZcqUcbs/T548UrJkSec2nkyZMkUmTJiQYX2FChVC9OoBAIgtv/32m0k+RG1QM3r0aJk6darPbbSJqGbNmpKTjBkzRoYPH+68rU1hlSpVksOHDwf9gcRSZK5B4JEjR2iy47hxruVA/B/lmGUXbe2oWLGiSSgEK6JBjda79OnTx+c2VatWzdK+ExISzM8TJ06Y3k8Wvd2gQQPnNidPnnR73F9//WV6RFmP9yR//vxmSU8DGmpqAqPHi2MWOI4bxyy7cK5xzLKLln5EdVBTunRps4RDlSpVTGCydu1aZxCjf3lorYzVgyoxMdFkWbZt2yaNGjUy69atWydpaWmm9gYAAESPqOnSrU07O3bsMD+vXr1qftfFdUwZbaZasWKF+T0uLs70knrmmWfk/fffl2+//VZ69eplejR16tTJbFOrVi1p06aNDBgwQLZs2SIbNmyQhx56yBQR+9vzCQAA5AxRUyisg+gtWLDAebthw4bm52effSYtWrQwv+/bt8/ZE0mNHDlSLl68aMad0YxM8+bNTZftAgUKOLdZtGiRCWRatmxpUl9dunQxY9sEQpuixo8f77FJChyzUOJc45hlF841jlk0nmtRN04NAABAVDc/AQAA+EJQAwAAbIGgBgAA2AJBDQAAsAWCmiDNnDlTKleubHpU6dg22jUc/09nPe/QoYPpIq/d7FeuXOl2eLROXXu26QCJBQsWNHNv7d+/P6YPoU7DcfPNN0vRokXNNB46BIH27HOlM80PGTJErrnmGilSpIjptacDS8aqWbNmSb169ZwDxekYVB9//LHzfo6Xf5577jnncBgcO8+efvppc4xcF9dR7znXvDt69Kj07NnTfG/p933dunXl66+/Dun1gKAmCEuWLDHTJWhXtO3bt0v9+vWldevWGUYpjmXapV6PiwZ/nkybNs10oZ89e7YZGLFw4cLmGOoXQ6z64osvTMCyadMmWbNmjVy5ckVatWpljqXl0UcflQ8++ECWLl1qtv/111/lH//4h8Sq6667zlyQdSBN/ZK84447pGPHjrJnzx5zP8crc1u3bpU5c+aY4NAVxy6jG2+8UY4dO+ZcvvrqK45XJn7//Xe59dZbJW/evOYPju+++05efPFFKVGiRGivB9qlG1nTpEkTx5AhQ5y3r1696ihfvrxjypQpHFIP9HRbsWKF83ZaWpojISHB8fzzzzvXnT171pE/f37Hu+++yzH8PydPnjTH7osvvnAeo7x58zqWLl3qPEZ79+412yQnJ3Pc/k+JEiUcb7zxBsfLD+fPn3dUr17dsWbNGsftt9/ueOSRRzjXvBg/fryjfv36Hu/j/6Z3o0aNcjRv3tzr/aG6HpCpyaLLly+bvwo1PWbRwfv0dnJyclZ3G1MOHjxoZkN3PYY6f5Y243EM/581oKQ12Zued5q9cT1umv7WCeE4bmJGHF+8eLHJbGkzFMcrc5oZbN++vds5xbnmnTaJaJO6zk3Yo0cPM9I9x8s3Hdm/cePG0rVrV9OsrgPovv766yG/HhDUZNHp06fNl2fZsmXd1utt/WCQOes4cQy903nItL5B07Z16tRxHrd8+fJJ8eLFOfdc6FQoWl+ko5I++OCDZsqU2rVrc7wyoQGgNp9rLZen/6Oca+70Ijt//nwzOr3WcunF+G9/+5ucP3+e4+XDTz/9ZI5X9erV5ZNPPjFzMD788MPOmQJCdT2ImmkSgFj9C3r37t1ubfbwrEaNGmY+OM1sLVu2THr37m3qjeDdkSNH5JFHHjG1W67Tx8C7tm3bOn/X+iMNcipVqiTvvfeeKW6F9z/QNFPz7LPPmtuaqdHvNq2f0f+roUKmJotKlSoluXPnztDjRG/r7ODInHWcOIae6ZxkH374oZnfTAthXY+bNn/qfGace/9PMwrVqlWTRo0amayDFqi//PLLHC8ftGlOOzbcdNNNkidPHrNoIKjFmvq7/pXMueabZkxvuOEGOXDgAOeaD9qjSTOnrnRSaavpLlTXA4KaIL5A9ctz7dq1bpGo3tZ2fGSuSpUq5mR1PYYpKSmm6j2Wj6HWVGtAo80n69atM8fJlZ532oPA9bhpl2/9cojl45ae/n9MTU3lePmgE/lqs51muKxF/5rWOhHrd8413y5cuCA//vijuWjzf9M7bUJPPzTFDz/8YLJcIb0e+F1SjAwWL15sKrPnz5/v+O677xwDBw50FC9e3HH8+HGOlkuvim+++cYserpNnz7d/P7zzz+b+5977jlzzP7zn/84du3a5ejYsaOjSpUqjj/++CNmj+GgQYMc8fHxjs8//9xx7Ngx53Lp0iXnNg8++KCjYsWKjnXr1jm+/vprR2Jiolli1ejRo03vsIMHD5rzSG/HxcU5Vq9ebe7nePnPtfcTxy6jESNGmP+beq5t2LDBkZSU5ChVqpTppcjx8m7Lli2OPHnyOCZPnuzYv3+/Y9GiRY5ChQo53n77bec2obgeENQE6ZVXXjEXl3z58pku3ps2bQp2l7by2WefmWAm/dK7d29nN76xY8c6ypYtawLEli1bOvbt2+eIZZ6Oly7z5s1zbqP/yQcPHmy6LesXQ+fOnU3gE6seeOABR6VKlcz/w9KlS5vzyApoFMcr60ENx85dt27dHOXKlTPn2rXXXmtuHzhwgOPlhw8++MBRp04d811fs2ZNx9y5c93uD8X1IE7/8T+vAwAAkDNRUwMAAGyBoAYAANgCQQ0AALAFghoAAGALBDUAAMAWCGoAAIAtENQAAABbIKgBAAC2QFADIMv69OkjnTp1itgRvP/++52z/trx+HTv3l1efPHFkL4mwM4IagB4FBcX53N5+umnzSzY8+fPj8gR3Llzp3z00Ufy8MMPS6QdOnTIHBOdBDKUnnrqKZk8ebKcO3cupPsF7CpPpF8AgJzp2LFjzt+XLFki48aNc5tlt0iRImaJlFdeeUW6du0a0dcQbnXq1JHrr79e3n77bRkyZEikXw6Q45GpAeBRQkKCc4mPjzeZCNd1Gkykb15p0aKFDB06VIYNGyYlSpSQsmXLyuuvvy4XL16Uvn37StGiRaVatWry8ccfuz3X7t27pW3btmaf+hhtVjp9+rTXT+bq1auybNky6dChg9v6ypUryzPPPCO9evUy+6pUqZK8//77curUKenYsaNZV69ePfn666/dHvfvf/9bbrzxRsmfP7/ZR/omH12nzVwPPPCAeQ8VK1aUuXPnOu+vUqWK+dmwYUNznPQ4uHrhhRekXLlycs0115jg5MqVK877XnvtNalevboUKFDAvPe7777b7bH6HhcvXsxZCviBoAZASC1YsEBKlSolW7ZsMQHOoEGDTEalWbNmsn37dmnVqpUJWi5dumS2P3v2rNxxxx0mINBgY9WqVXLixAm55557vD7Hrl27TJNM48aNM9z3z3/+U2699Vb55ptvpH379ua5NMjp2bOneX7NfOhtay7fbdu2mefS+pVvv/3WNKuNHTs2Q7OaBjr6fLrfwYMHm/dlZa70vapPP/3UZLiWL1/ufNxnn30mP/74o/mpx0b3a+1b3682n02cONHsS9/7bbfd5va8TZo0MftPTU0N4lMBYkRAc3oDiEnz5s1zxMfHZ1jfu3dvR8eOHZ23b7/9dkfz5s2dt//66y9H4cKFHffff79z3bFjxzSacCQnJ5vbkyZNcrRq1cptv0eOHDHb7Nu3z+PrWbFihSN37tyOtLQ0t/WVKlVy9OzZM8NzjR071rlOn1fX6X3qvvvuc9x5551u+3n88ccdtWvX9rpffd4yZco4Zs2aZW4fPHjQ7PObb77JcHz0sXocLF27dnV069bN/P7vf//bUaxYMUdKSorDm507d5p9Hzp0yOs2AP4XmRoAIaXNO5bcuXObJpe6des612kTizp58qSz4FezGFaNji41a9Y092mGw5M//vjDNBVpU4+v57eey9fz792712R2XOnt/fv3m2YuT/u1muKsffiizVp6HCzaDGU97s477zRNZFWrVjUZpUWLFjkzWJaCBQuan+nXA8iIoAZASOXNm9fttgYAruusQCQtLc38vHDhgqkb0Z5DrosGFembYizavKUX+cuXL/t8fuu5fD1/MO/Ln334epzW52iT2LvvvmuCHS3Grl+/vmmSs5w5c8b8LF26dECvF4hFBDUAIuqmm26SPXv2mGJcLSJ2XQoXLuzxMQ0aNDA/v/vuu6Cfv1atWrJhwwa3dXr7hhtucMuw+JIvXz7z0zWz4688efJIUlKSTJs2zdQKaffwdevWuRVRX3fddSaQA+AbQQ2AiNLeQJqNuPfee2Xr1q2myemTTz4xvaW8BQmatdBg6Kuvvgr6+UeMGCFr166VSZMmyQ8//GCKeV999VV57LHH/N5HmTJlTDORVeTs77gyH374ocyYMcNkpn7++WdZuHChyeLUqFHDuc2XX35piqsBZI6gBkBElS9f3mRGNIDRi7fWv2iX8OLFi0uuXN6/ovr3729qUIKlwdF7771nuk3ruDDaBKS9kbS7eiDZFg1O5syZY96Pdh/3h75H7Smlvb80YzR79mzTFKV1OOrPP/+UlStXyoABA7L8/oBYEqfVwpF+EQAQKC0W1oyGDgyYmJhoywM4a9YsWbFihaxevTrSLwWICmRqAEQlbe7R5hpfg/RFOy0y1pGTAfiHTA0AALAFMjUAAMAWCGoAAIAtENQAAABbIKgBAAC2QFADAABsgaAGAADYAkENAACwBYIaAABgCwQ1AABA7OB/AMpGNZf1nwWbAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 600x350 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(figsize=[6, 3.5])\n",
    "ax.axhline(y=0, linestyle=':', color='k')\n",
    "ax.fill_between(rd_results['time'], rd_results['lcb'], rd_results['ucb'], color='k', alpha=0.1, step='post')\n",
    "ax.step(rd_results['time'], rd_results['rd'], color='k', where='post')\n",
    "ax.set_xlim([0, 60])\n",
    "ax.set_ylim([-1, 1])\n",
    "ax.set_xlabel(\"Time (months)\")\n",
    "ax.set_ylabel(\"Risk Difference\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "03404a84-77b4-4875-a3f4-128a60c4c6c4",
   "metadata": {},
   "source": [
    "The same process can also be done when modeling time using splines. The following code illustrates this process"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "d342d52c-9c5f-4fbe-aea4-e2c38b2a69a9",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Creating time design matrix\n",
    "t_steps = np.asarray(range(1, 60))\n",
    "intercept = np.ones(t_steps.shape)[:, None]\n",
    "time_splines = spline(t_steps, knots=[10, 20, 30, 40],\n",
    "                      power=2, restricted=True, normalized=False)\n",
    "s_matrix = np.concatenate([intercept, t_steps[:, None], time_splines], axis=1)\n",
    "tp_intervals = list(range(0, 60))\n",
    "params_rd = len(tp_intervals)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "694b6aaa-006c-4501-a12e-1b2c64131de1",
   "metadata": {},
   "outputs": [],
   "source": [
    "def psi_rd(theta):\n",
    "    # Extracting parameters\n",
    "    risks = theta[:params_rd]\n",
    "    idPLRM = params_rd + 7\n",
    "    beta1 = theta[params_rd:idPLRM]\n",
    "    beta0 = theta[idPLRM:]\n",
    "\n",
    "    # Nuisance models\n",
    "    ee_plog1 = ee_plogit(beta1, t=t, delta=y, X=W, S=s_matrix)\n",
    "    ee_plog1 = ee_plog1 * (a == 1)[None, :]\n",
    "    ee_plog0 = ee_plogit(beta0, t=t, delta=y, X=W, S=s_matrix)\n",
    "    ee_plog0 = ee_plog0 * (a == 0)[None, :]\n",
    "\n",
    "    # Predictions to get risk differences\n",
    "    risk1 = plogit_predict(theta=beta1, t=t, delta=y, X=W, S=s_matrix,\n",
    "                           times_to_predict=tp_intervals, measure='risk')\n",
    "    risk0 = plogit_predict(theta=beta0, t=t, delta=y, X=W, S=s_matrix,\n",
    "                           times_to_predict=tp_intervals, measure='risk')\n",
    "    ee_rd = (risk1 - risk0) - np.asarray(risks)[:, None]\n",
    "\n",
    "    # Returning stacked estimating equations\n",
    "    return np.vstack([ee_rd, ee_plog1, ee_plog0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "c70dfcff-0d5a-433c-831a-bf6af05af568",
   "metadata": {},
   "outputs": [],
   "source": [
    "inits = [0., ]*params_rd + [0., 0., -4., ] + [0., ]*4 + [0., 0., -4., ] + [0., ]*4\n",
    "estr = MEstimator(psi_rd, init=inits)\n",
    "estr.estimate()\n",
    "cb = estr.confidence_bands(method='supt', subset=[x for x in range(1, params_rd)], seed=7777)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "92f7f072-8630-4dda-92db-79b073147567",
   "metadata": {},
   "outputs": [],
   "source": [
    "rd_results = pd.DataFrame()\n",
    "rd_results['time'] = tp_intervals\n",
    "rd_results['rd'] = estr.theta[:params_rd]\n",
    "rd_results['lcl'] = [0.,] + list(cb[:, 0])\n",
    "rd_results['ucl'] = [0.,] + list(cb[:, 1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "c8311966-b997-4631-b04c-7ae1e7aa5c76",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjUAAAFSCAYAAAAKDC7nAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjgsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvwVt1zgAAAAlwSFlzAAAPYQAAD2EBqD+naQAARaFJREFUeJzt3QvcjHX+//HPjRyLlHM5RqFyiLLUVhvRYYu2WtrksGGTdFJhNySVztvJpiPa7FItttpSosMmIoVIoijKMSEKlfk/3t/9fed/zZi577nP91zzej4el5m55pq5Z665zPWe7zErEolEDAAAIM2VKu4XAAAAUBAINQAAIBQINQAAIBQINQAAIBQINQAAIBQINQAAIBQINQAAIBQINQAAIBQINQAAIBQINQAAIBTSKtS88847dt5551mdOnUsKyvLZsyYkeNj3nrrLTvhhBOsXLly1rhxY5s4ceIB24wbN84aNGhg5cuXt3bt2tmCBQsK6R0AAIDCklahZvfu3dayZUsXQlKxZs0aO/fcc+03v/mNLV682K699lrr16+fvfbaa9Ftpk6datdff72NGjXKPvzwQ/f8Xbp0sc2bNxfiOwEAAAUtK10ntFRJzfTp061bt25Jtxk6dKj95z//sWXLlkXX9ejRw7Zv324zZ850t1Uyc+KJJ9ojjzzibu/fv9/q1q1rgwcPtmHDhhXBOwEAAAWhjIXYvHnzrFOnTjHrVAqjEhvZt2+fLVq0yIYPHx69v1SpUu4xemwye/fudYunILRt2zY7/PDDXdgCAACpUdnK999/75qW6BycH6EONRs3brSaNWvGrNPtnTt32o8//mjfffed/fLLLwm3+fTTT5M+79ixY2306NGF9roBAMg069atsyOPPDJfzxHqUFNYVLKjdjjejh07rF69eu4DqVy5crG+NgAA0okKGtTs45BDDsn3c4U61NSqVcs2bdoUs063FTwqVKhgpUuXdkuibfTYZNSTSks8PS+hBgCA3CuI5htp1fspt9q3b2+zZ8+OWTdr1iy3XsqWLWtt2rSJ2UbtY3TbbwMAANJDWoWaXbt2ua7ZWnyXbV3/6quvotVCvXr1im5/xRVX2BdffGE33XSTayPzt7/9zZ577jm77rrrotuoGumJJ56wSZMm2YoVK2zgwIGu63jfvn2L4R0CAICMqH764IMP3Jgznm/X0rt3bzeo3oYNG6IBRxo2bOi6dCvEPPjgg64B0pNPPul6QHndu3e3LVu22MiRI13D4latWrnu3vGNhwEAQMmWtuPUlLRGTlWqVHENhmlTAwBA8ZxD06r6CQAAIBlCDQAACAVCDQAACAVCDQAACAVCDQAACAVCDQAACAVCDQAACAVCDQAACAVCDQAACAVCDQAACAVCDQAACAVCDQAACAVCDQAACAVCDQAACAVCDQAACAVCDQAACAVCDQAACAVCDQAACAVCDQAACAVCDQAACAVCDQAACAVCDQAACAVCDQAACIW0CzXjxo2zBg0aWPny5a1du3a2YMGCpNuefvrplpWVdcBy7rnnRrfp06fPAfefddZZRfRuAABAQSljaWTq1Kl2/fXX2/jx412geeCBB6xLly62cuVKq1GjxgHbT5s2zfbt2xe9/e2331rLli3t4osvjtlOIWbChAnR2+XKlSvkdwIAADK6pOb++++3/v37W9++fa158+Yu3FSsWNGefvrphNsfdthhVqtWregya9Yst318qFGICW5XtWrVInpHAAAg40KNSlwWLVpknTp1iq4rVaqUuz1v3ryUnuOpp56yHj16WKVKlWLWv/XWW66k55hjjrGBAwe6Ep3s7N2713bu3BmzAACA4pU2oWbr1q32yy+/WM2aNWPW6/bGjRtzfLza3ixbtsz69et3QNXTM888Y7Nnz7a77rrL3n77bTv77LPd30pm7NixVqVKlehSt27dfLwzAACQcW1q8kOlNMcff7yddNJJMetVcuPp/hYtWthRRx3lSm86duyY8LmGDx/u2vZ4Kqkh2AAAULzSpqSmWrVqVrp0adu0aVPMet1WO5js7N6926ZMmWKXX355jn+nUaNG7m+tXr066TZqg1O5cuWYBQAAFK+0CTVly5a1Nm3auGoib//+/e52+/bts33s888/79rB9OzZM8e/s379etempnbt2gXyugEAQNFIm1AjqvJ54oknbNKkSbZixQrXqFelMOoNJb169XJVQ4mqnrp162aHH354zPpdu3bZjTfeaPPnz7e1a9e6gNS1a1dr3Lix6yoOAADSR1q1qenevbtt2bLFRo4c6RoHt2rVymbOnBltPPzVV1+5HlFBGsPm3Xfftddff/2A51N11tKlS11I2r59u9WpU8c6d+5sY8aMYawaAADSTFYkEokU94tId2oorF5QO3bsoH0NAADFdA5Nq+onAACAZAg1AAAgFAg1AAAgFAg1AAAgFAg1AAAgFAg1AAAgFAg1AAAgFAg1AAAgFAg1AAAgFAg1AAAgFAg1AAAgFAg1AAAgFAg1AAAgFAg1AAAgFAg1AAAgFAg1AAAgFAg1AAAgFAg1AAAgFAg1AAAgFAg1AAAgFAg1AAAgFAg1AAAgFNIu1IwbN84aNGhg5cuXt3bt2tmCBQuSbjtx4kTLysqKWfS4oEgkYiNHjrTatWtbhQoVrFOnTrZq1aoieCcAACBjQ83UqVPt+uuvt1GjRtmHH35oLVu2tC5dutjmzZuTPqZy5cq2YcOG6PLll1/G3H/33XfbQw89ZOPHj7f333/fKlWq5J5zz549RfCOAABARoaa+++/3/r37299+/a15s2buyBSsWJFe/rpp5M+RqUztWrVii41a9aMKaV54IEH7Oabb7auXbtaixYt7JlnnrFvvvnGZsyYUUTvCgAAZFSo2bdvny1atMhVD3mlSpVyt+fNm5f0cbt27bL69etb3bp1XXBZvnx59L41a9bYxo0bY56zSpUqrloru+cEAAAlT9qEmq1bt9ovv/wSU9Iiuq1gksgxxxzjSnH+/e9/27PPPmv79++3Dh062Pr16939/nG5eU7Zu3ev7dy5M2YBAADFK21CTV60b9/eevXqZa1atbLTTjvNpk2bZtWrV7fHHnssX887duxYV6LjF5UCAQCA4pU2oaZatWpWunRp27RpU8x63VZbmVQcdNBB1rp1a1u9erW77R+X2+ccPny47dixI7qsW7cuD+8IAABkZKgpW7astWnTxmbPnh1dp+ok3VaJTCpUffXxxx+77tvSsGFDF16Cz6mqJPWCyu45y5Ur53pVBRcAAFC8ylgaUXfu3r17W9u2be2kk05yPZd2797tekOJqpqOOOIIVz0kt956q/3qV7+yxo0b2/bt2+2ee+5xXbr79esX7Rl17bXX2m233WZNmjRxIWfEiBFWp04d69atW7G+VwAAEOJQ0717d9uyZYsbLE8NedVWZubMmdGGvl999ZXrEeV99913rgu4tq1ataor6Xnvvfdcd3DvpptucsFowIABLviccsop7jnjB+kDAAAlW1ZEg7UgX1RlpQbDal9DVRQAAMVzDk2bNjUAAADZIdQAAIBQINQAAIBQINQAAIBQINQAAIBQINQAAIBQINQAAIBQINQAAIBQINQAAIBQINQAAIDMDTU///yzvfHGG/bYY4/Z999/79Z98803tmvXroJ+fQAAAIUzoaVmuT7rrLPc5JF79+61M8880w455BC766673O3x48fn9ikBAACKvqTmmmuusbZt27oZsCtUqBBdf8EFF9js2bPz/4oAAACKoqTmv//9r7333ntWtmzZmPUNGjSwr7/+Oi+vAQAAoOhLavbv32+//PLLAevXr1/vqqEAAADSItR07tzZHnjggejtrKws10B41KhRds455xT06wMAAEhJViQSiVguqESmS5cupoetWrXKta/RZbVq1eydd96xGjVqWKbZuXOnValSxXbs2GGVK1cu7pcDAEBGnkNzHWp8l+6pU6fakiVLXCnNCSecYJdeemlMw+FMQqgBACBNQw0K7wMBACCT7CzAc2iu29SMHTvWnn766QPWa53GqgEAACgOuQ41GkW4adOmB6w/9thjGXgPAACkT6jZuHGj1a5d+4D11atXtw0bNhTU6wIAACjcUFO3bl2bO3fuAeu1rk6dOrl9OgAAgOIJNf3797drr73WJkyY4OaB0qL2NNddd527r7CNGzfOjV5cvnx5a9eunS1YsCDptk888YT9+te/tqpVq7qlU6dOB2zfp08fN9ZOcNHcVgAAIOTTJNx444327bff2pVXXmn79u1z6xQwhg4dasOHD7fCpG7k119/vWu7o0CjQQA1Zs7KlSsTjo/z1ltv2SWXXGIdOnRwr1ENmTV44PLly+2II46IbqcQo5DmlStXrlDfBwAAKHh57tKt8WlWrFjhxqZp0qRJkQQBBZkTTzzRHnnkkeiUDaoOGzx4sA0bNizHx2t6B5XY6PG9evWKltRs377dZsyYkefXRZduAADSsEu3d/DBB7uAcdxxxxVJoFGp0KJFi1wVkleqVCl3e968eSk9xw8//GA//fSTHXbYYQeU6Kik55hjjrGBAwe6kqjs7N27130IwQUAAKRZ9dPu3bvtzjvvtNmzZ9vmzZtdaUnQF198YYVh69atrqSlZs2aMet1+9NPP03pOVRFpsbMwWCkqqff/e531rBhQ/v888/tz3/+s5199tkuKJUuXTrpWD2jR4/O5zsCAADFGmr69etnb7/9tl122WWua7ca1qYDBbEpU6a4Uhm1r/F69OgRvX788cdbixYt7KijjnLbdezYMeFzqe2Q2vZ4KqlRNRgAAEijUPPqq6/af/7zHzv55JOtKGnCTJWcbNq0KWa9bteqVSvbx957770u1LzxxhsutGSnUaNG7m+tXr06aahRdRuNiQEAKFly3aZGDW3j26QUhbJly1qbNm1ctZenqi/dbt++fdLH3X333TZmzBibOXOmm1E8lVnI1aYm0QCDAAAgRKFGAWHkyJGu0W1RU5WPxp6ZNGmS63mlRr1q49O3b193v3o0BbuVqwv3iBEj3Dg6GttGoyFrUc8t0aW6qM+fP9/Wrl3rAlLXrl2tcePGrqs4AAAIcfXTfffd5xrUqoGugsJBBx0Uc/+HH35ohaV79+62ZcsWF6oUTlq1auVKYHzj4a+++sr1iPIeffRR12vqoosuinmeUaNG2S233OKqs5YuXepCkrp1qxGxxrFRcKN6CQCAkI9Tk1OvHwWGTMM4NQAAFP85NM+D76FwPhAAADLJzuIefE9VNU8++aRrv7Jt27ZotdPXX3+drxcDAABQZG1q1AZFg9cpValxrSaxVG+oadOmuTYtzzzzTJ5fDAAAQJGV1KgHkuZLWrVqVcwgduecc4698847eX4hAAAARRpqFi5caH/6058OWK9Zr9UjCQAAIC1Cjbo6J5rA8bPPPrPq1asX1OsCAAAo3DY1559/vt1666323HPPudua+0ltaTRZ5IUXXpjbp0MK1EFNk3lqBOVEi+7X+Dz6LBJdatGYPOkyTxcAAEXSpVtdrjSY3QcffGDff/+9G7BO1U6aquCVV16xSpUqWaYpjC7dP//8s+3Zs8eNmKyRj+ODTJA+QgWW4GUw1CjQ6FJTTagdVJkyZdy64GVw0EIAANLxHJrrkhr94VmzZtncuXNtyZIl7oR7wgknuB5RyJ+ffvrJBRntU4WZvXv3umCiKj8fQHxYyanURcHGl+L4Uh49pw4eXffP4Z9Xf0OBRyNEa52/JOwAANJFmdyedCtUqGCLFy92s3QX9UzdYaTQoRCjUi+FDk3roCChkHHooYfmucpIj1NYEYWTZH9bJUJaNJeXUrKnx2hR6Y4+c10Gww5VWcgvhW2/6BgMVqf62/7+YAllMNgHFx3vOk59+Ndl8DoBHQi/XIUandDq1avnvmSQf/ri3rp1q5sVXF+4KimpWLFikQUG/R19pvHzdwXDzo8//uhKd3yVlg87eq2JSnaA+GM8uOiHkUogFd6DoSW+FjxR+zBfteqDjz9W/WUw+Piq2GCo0TGqY1bBxx+vvqQSCKtIJBJd4m+nui7Z7eD/ueCPkkT36Ye6mqsU9vkt12ehv/zlL/bnP//Z/v73v7tB95A3KhnZvHmzK6VRHWJJCgTZhR2dlHQSUtD57rvv3DrfLsefNHTwBsMOJ43w0zGhY0OLwouCi6pSfZDR/cESFb/4kpXC+KLzX6w+OPmqXf+l7INO/HFLaSQKQ7IAkdOS7DGJAsT+BMEi0Q+AnJ4/+P8xu9vBy2CziPjrwZLYwj4f5PpM+sgjj9jq1atd4qpfv/4BDYMLc5buMNABoWoeBRp90VatWjVtqnL0OnUSElVJeTpQ/clMQcf/R/K/hH3Jjm8bFFzS5b0jNtgGS11UmucDjf/sfWDw7bWKI9gGA1QivvRI78E3xveBXsemjnFfGukXqrDCKadgEQwK2W2TrHdq8PGp/D2/XU7fj4mqYSW726lsX9D0gyK+g0uJCTXdunUrnFeSAfSl76ub9IV58MEHWxj4NkDxfFsJH3Z8tWWw55Ue538h+/X+PgJP8QpWGfmqyGDpi6+S9J+djud0KpXzx1qy0kj9+PBz28W3MQsGHcJ50UglaGQXPnypXfAyGDwS/Y1EgcRXgya6ntPit0vWLiy+BAS5xyzdRdQdrSRXNxU1fUnoxBgMPZ5vt+OrJnz7h2CVBaGncNu9qL2LjlcfXvwvLF/6komlbPHtgjy/PwjnsZKFjFRuxwcOX3URHzwSBRh/v8Qfn/HhIb4naaqhBHkvqWnYsGHCHz7F2qXbz9L9wgsv2Oeff2433nija1ujaqeaNWu66RIQSx+UxvLRh5pO1U2FxY+Zk1OPGJ1YfRf0+OoEX0Lgfz0n6uniLzO9ykBf8vE9jXy1ixZ/25ekBatgiqvqKN1KdfRjRd+Lvr1BMIQH2+okOz79Cbao5aa6JdG6YFuJYE81HyyShZjg308UPrIryUjUgJzQAY9ZuguZiuw3bdrkvsAOOeSQwv5zaS+7NhASPDnrhKzg479Ug7+4gicMPx5PcIk/mcSPwBz8Ai2pkv2y9b9ufZsXX9oS7Gnk91GwGpDwkvc2ZvGC3dQVePx+D57Eg+HGP18wnPtSy/hjMFn4CYaFYK8Tf5msvYd/bLIQE//8wb+fKFz41+v/3yW7HygRocbP0n333XfHnKQ1S/cf/vCHgn59aU1fCmpDoy83jTmDwg89wZKJ4Be3b9AaXOe/nH3bkGS/AIMBKHjSyanRXVD8iSDRySJ4Pf4EE98WwJesBBsoBgOLBF9vsNQl00sKS8JxGh9I/eesqr9kbTn844LHrT+Wcjre8tPGg+MFoQ41mqX7scceO2A9s3QfSI1jVfVEoCla/hdibgRPIvHdIhUgsiuCD37pJzrBxP+dZCeinN5TfMjyv+x9tUZxVWEg94KlMgCKMdQwS3dqVCrgezlR3Jo+JxkAQPoqlddZun0PAJ0MmKU7cbWTipaD47kAAIASFGruu+8+1/CtRo0arjTitNNOs8aNG7v2NbfffnvhvMo0rXaiYTAAAEWHWboLqdpJczhR7QQAQAkrqdE4NKpOkT/+8Y9uRmnN0H3llVfaTTfdZJ06dbKiMm7cOGvQoIEbvrxdu3a2YMGCbLd//vnnrWnTpm77448/3l555ZWY+9Vgc+TIkVa7dm1XVaT3smrVqjy9NlU3+d5O+nsAAKCEhRp1M9QgaDJp0iQ3OmBxmDp1qutSPmrUKDfYX8uWLa1Lly5upN5E3nvvPbvkkkvs8ssvt48++shN8aBl2bJl0W3UNf2hhx6y8ePH2/vvv+/mstJz5uU9qsqpIEZEBAAAhTRNwplnnukGkGvTpo0LNd27d0/aAPbpp5+2wqKSmRNPPNFNqukb5NatW9cGDx5sw4YNO2B7vc7du3fbyy+/HF33q1/9ylq1auVCjN66JuYcMmSI3XDDDe5+hRKNjDxx4kTr0aNHSq/LD/H8zjvvuOfz+0ZhUKU2fnAzTwPGiUpzfBWVnxBQt4OlPLnZVlVfek/BgdT8sPdq0B38zPKyrQYa812l/eB3+dnWD50dnBE8N9vqtg+fqu7z9Hg9jx9tOLfb6vXrfYheg+8m7T/P3GybymdfEMdJos8zzMdJss8zv8dJ8PPM73GS7PPM63HCdwTfEen6HbHn//5f1qpVy63zkxv7//cqmNC5syAKBVIqqXn22Wfd4HpqICz6w2oMm2gpLNpBixYtiqnq0geh2/PmzUv4GK2PrxpTKYzffs2aNW76guA2CicKT8me03/BKcgEFzn11FPtvPPOs169ernpIzQYYZMmTdxghevXr48OotWiRQu3/uuvv44+p0KU1vlw5em1aH2wSuy5555z61T9F3T66ae79R9//HF03YsvvujW9e3bN2ZbfZ5ar9Ip74033nDr4sPchRde6Na/9dZb0XVz585169QbLuiyyy5z61999dXoOpWqaV38Z9G/f3+3ftq0adF1K1ascOtOOeWUmG2vvvpqt37y5MnRdWvXrnXrFLaDhg4d6tY/9dRT0XUK5VrXrFmzmG1Hjx7t1j/88MPRdfo8tU5LcF6qu+66y63Tpaf7/bb+OBA9n9bp+YP097Ver8fT69Q6ve4gvS+t1/v09P61TvsjSPtL67X/PO1XrdN+DtLnoPX6XDx9Xlqnzy9In6/W6/P2dBxonY6LIB03Wq/jyNPxpXU63oJ0PGq9jk9Px63W6TgO0nGu9TruPf1/0Dr9/wjS/x+t1/8nT//PtE7/74Juvvlmt/7RRx+NrtMElv7zDFInCK27//77Y77M/bY+3Ii20br4jhN+Wz9Jpuhva51eSxDfEf/Dd0S4viM6d+7sJr4NFjT897//dYGmSBsKq+TizjvvdNc1IdXf//53O/zww60o+S7Sei3xr+3TTz9N+BgFlkTba72/369Ltk0iY8eOPeBk5emACR40ohIcfflq3id9Wfnu8Bs2bLAjjzySAdMAACiq6ic1FP7ss8+sWrVqrqHwgw8+WOTdlb/55hs3arHaybRv3z66Xg2V33777ZgSB0/FxKouU7sa729/+5sLJPqlrOdSg2c9txoKe7///e9d0FAbnkT8RICefqGrGuzxxx93xWoqStPzKxjpUr8Utf+Cv/o9BZ3WrVu7arW2bdu60KMk61H99D9UP5XcouX4bal+ovpJqKKmijqrGKqfUgo1OskuXbrUGjVq5L7YdLKuXr26FSV9iaqOW7ODq7Gv17t3bzc77r///e8DHlOvXj3XsPjaa6+NrlMj4xkzZtiSJUvsiy++sKOOOso1IlY7G09j7+i2wltu2tQoWKmUx9flBykEqURJf1fF7Nqfuh0fdNRQWaU6Cm5a1GMrt0P+AwBQUvhQo5qeRCO3+3NoQYSalM6WOrkqSKj+ThlIdXVF3VBYv/7092fPnh0NNdpJun3VVVclfd26PxhqZs2aFS3p0Q5WctQ2PtRo5yqcDBw4MNevUXM8qbu7SrbiKZmqt5aWYNBRVdUHH3zg2vDMnz/fBbQ5c+a4xQdKNW5W0FKbHYUw5vcBACCPoUYNhf/617/a559/7k6oSlPF0a1bpS4qmVE1zUknnWQPPPCA693kG8Gqga6qqNTmRa655hoXBjQK8rnnnmtTpkxxAULVRKL3osBz2223uUZNCjkjRoxwxWDB0qBUqSpJoUTF8KlMj6CgozClpV+/fi6kffLJJy7g+JCjfa2Gl77xpd6fDzhq+KW/CQAAUqx+CtKJX8GgqBsKe+rOfc8997gqMIUBjTHje0Co14QG5gv2fNDge+pZ4HvKaFyaYE8MvX1VSSnoKJAoKKjdzdFHH53yawoWnakNgxoAK2zkt0RFdY0KOWodrnZDGmhQ1XCenl8lPwo5eu8nnHACVVUAgIytfsp1qIFl+4GoTcy6detc1VJBN6ZWCZBKbxRw1KNq5cqVMffrYFAoU8DRolIdAACKU4kLNSoNGTBggOtBoevZie8bnwniPxC1q9HM5Qo1hdnIVyVCCjcKOVpU0hSkiUZ9KY7a5QQHEwMAICNDTbDKSdeTPllWlutRlGniPxDtUgUODUZYVG1eVFWlHlUa9EiLBk3SQRRsaK12SAo5Wpo3b06DYwBA5oUaZC/RB6IPUaU1Kqkpjskt9VreffddF3BUihMcvVjUJf/Xv/61Czi6jB+AEACAgkCoSTPJUua3335bYI2G80O5VT3XVFWlkKNBB4PDussxxxzj2uMo4KjLe3AAQAAAQhdq1H1a895ovgj1JtKJWi/yoosucnOuZGqbjWQfiAbWU6Nh9YgqSSFBjZg1j5ZKcNSzStVWwcNAB51GOVbAUdDR9eAkagAApHWoUVfiDh062LJly+zss8+2pk2buhOhBo+bOXOm606skoBEo+mGXXYfiBoNK9go1JTUkYHV9kelNwo4WoKTo4mqzzQ2kD5/LepG7mcqBgAgrUYU9rPJaqZpDfOvqoogDfevHjbjx4+3wYMH5+sFhY3CjD4s32i4JI4GrNelwQm1iEKY2uMopGoQwC1btrjbWkQDC6rRsQKOelVpvipCDgCguKVcUqMGpZrocdCgQQnvf/jhh928TKrSyDQ5pUyVcmnSTE1OqakU0okOj9WrV7tp5VWao5Czbdu2A0pyVEWlQRC1aDoLjdcDAMCeklj9pN4yamR67LHHJrxf1VK/+c1v3K/6TJPKB6IPVT2Q1M6mqGc4L0g6MDXon8KNQo7myYoPOTpoNRGnSnNUbaWQozm2AACZZ09JDDVqK6NqiWQnJ/XyqV+/fsww/pki1Q9EDa0VbEqVKhWaRtW+JEfhRoumclA1ZbwjjzwyGnB02axZs4xsfwUAmWZPSQw1eiGab0klNols2rTJTQSpQeAyTW4+EG2rYKPeRMUxfk1R0PvzAUe9rNTmKjgQoOi9a+4uP3O5rterV69EtjkCAIQs1Kh04bjjjkvag0fVKsuXLyfUpPCBaPwaBUQ1Is6E0opdu3bZRx995Eal1kjHCjo6eBM1WA4GHTVA1qCABB0ASF97SmKoGT16dEpPqBmvM01uPxDtcrU92rx5s3tcog85zHRwazBABRz1ptOiQKzxfOJVq1bNteMKLo0aNcq4fRbWLzoFfLXJ2rp1qxv+QINCxi9qYK9tRZ+7Fv240g8tf10/DtQIX8H4sMMOi7lUGzaCMVB8SmSoQXJ5+UD0AavKTl/q+jLWF3Qm04CAGvNo8eLF0aCzatWqA6qtfNWV2uRorKQmTZq4IQaOPvpoq127NievEkDt6lQSqWpI9frzlzreFV58iFEJXlFQ6FE4VvVm3bp1Xdu/4KXaCWb6/z+gMBFq0kxeU6aq7NTAWo8rqWPYFCf9SldPK/WsU0mOlk8++eSAKR487XsfcnSpEh0tOnFlQjVfUVG7OR23X375pZvfzF9qUXhRCWSqv5V84FCpiv4PaQyk4KIG9br07c/0t+MX/T9SkNIs9QpMGhPKXyY7VoI0xpICjo6Zo446yl1qhnstDE0A5B+hJs3kp+jMj2GjX62U2ORMJzGNeKyA89lnn7nQo8s1a9Ykbc+lE6eCjQKOij/9pXpkHXHEEUwBkaSnnsKKX7TPfYBRyUuiqsIgNYRXyZn2rzoQ6FIlIgowPsToUv9fCjPMK9Qo3KiUyAcv9eLU+9Cleupl17lB7yEYcnRdi147P0KAWCpZ14+M+EXfJ/p/dvLJJ1P9lA7yWx+oYKNft/qlqfp/ShXyVn31xRdfuIDjFwUdLb49RiI6MakxskKPFgUdXepErPU6EStshu0EppO9wokWndh1gtd1H2JUPZQdHaO+CkfVOn7R/tO+O/zww9Nin+kLVz8qdOxoaAJVeaq9ly6z2wc6JnzAUdhRSNaifRLWXo0oPDrh6zygHwvZLb5UUpfx6+Ov/5xkm0SXydYlu63XG79eS6LmAp7OjSpBpU1NGiiIRk46SPQlqkVfiipyR/7pP5nad+ikpYDjL1XyoJN5KtUT+jwUcPziSxxU2qBFJ3DfMLW4S9v0ftVOSyUTCspagtdVbaT3ncogmXovDRo0iAaX4HXtg7A31lYJj4KODztadF0BMFn1moKcSqW0r7Qo6GifKej5sEf7naL/P6Eg4Bf9APIBwt/21/2i28nWx1/3ocFvH78+0br4+7MLA+kuKyvL/QhSb199B5W4UKMvRP0aS2T+/PluLqBMU1AfiD4KldboJKQvvpI0s3cYaX/rl4NOUr4qwl/XyV9hSCe23NDnpjCgEjd9fokWtRPxPXb8ottq26FLva5Ev5T8pXoJ6ZiLX3T86TLVsaLUXkT/l4OLShp0MlZw0TGNAykIqzTHhx1dKiQrLKuYPTv6jBUIFXBUtaVLjf0Vtl5bvuTBBwhdBq8nW6dS1dw8Lj6oJLpf/2fSjU78+l7w3wn+Mng9/vsj/nqZJNv450m0bfw2yW771+cvtT54PbhofYluU9O8eXM3saH+4wVpbiBNiKiTcqYpyA9E1L5GwUb/MfW86frFFgb6z+hLOxRytPheawpEwUXHQUmhkqQaNWq4kiVd+us6ofq2RDROL1j6KlVJq6/29EHHh+TcNKAWnRD0GSl8+kbT8YtKEYMnlfhLhWz9zWSLTviJqhJ8FUN8iUR8SUYwWPhAUtKDhA8LqSzx2wZvJ7ruPwu1KfNhIni/DyXB68HtfGgIW2nenpI4S7enkpjOnTvbm2++GZ3DSLM5n3feeXbLLbfk68Xgf/RrXge3Tp4qKcjEsWxKCp00fHuRnOjLXJ+XFgXT+EUlLPolr3FXktWNa1GIze4Xk/7f6T9+cNEx4q8r0NAuq+jpc1OpixbNexZPn63+TyvgqB2PX3xA9seOrqs0SMeGqgnTfT497Red5LXoxB28zGldom3it0t0XzCI+Ns+OIQtMCCfJTVKWxdddJH7j/faa6+5SQ3PP/98u+222+yaa66xTFTQJTXBL0F9oWlf+y6uAMLP99rSohCcbNEvYF+6Et/N3TfcVKhItvhSnWRVCzmVWCQKHPHrtS2lzZltT0mufhL9wlRVk/5TLV261MaOHWtXXXWVFSad2AcPHmwvvfSSS9oXXnihPfjgg0nbnWh7jW78+uuvu26c+vXUrVs3GzNmTExbgUT/2f75z39ajx49ij3UiA4E/+tN+13PT6kNACBd7Clp1U8KLvFU1XTJJZdYz5497dRTT41uo/l6CsOll17qim1nzZrlSjD69u1rAwYMsH/84x8Jt/dFu/fee69rB6RuqldccYVb98ILL8RsO2HCBDvrrLOit9XQs6RQgFOPCdWrq5habZb0C0j16fz6AQAglyU1OrHqBBrcNHjbX9dlYczSreHzFUwWLlxobdu2detmzpxp55xzjuuxoh4EqXj++eddCFO7Bj8xp17z9OnTXSlOXhVmSU2Q9rHaZahKSsXTKqVSwAEAoKQqcSU1asFfnObNm+dKT3ygkU6dOrmw9f7779sFF1yQ0vP4HRY/0/igQYOsX79+bqRZleaoFKgkloLoNen1q22N73Gjg0XhhsZvAIBMl1Ko0eBRxUndaNUlNUjBRN3KdV8q1NVS7WlUZRV066232hlnnOGqc9T+5sorr3Q9Va6++uqkz+XHQfCKuiuvGt6pe67CjN6XqqTUII8qKQBAJst137ZJkybZf/7zn+jtm266yZWidOjQwbVbyY1hw4Zl2zJfy6effmr5pdChhs2qworvdj5ixAg3H0Xr1q1t6NCh7v3cc8892T6fGkarqMwvGrCsOPjB0/T3FfLUoFhVa0y8DgDIRLnu/aQZkB999FFXuqFqoY4dO9oDDzxgL7/8sjuxTps2LeXnUtsQNX7NjqqEnn32WRsyZEjM6K7qrqgxRNROJrvqJ7VB6dKliyvF0GvMaV4WBbbf/va3rlpHpR+pltQoWBR2m5rsqC2TSphUJaVgQ8kNAKAkKHFtaoI0OqYmcJMZM2a4MWtUpaPSjtNPPz1Xz+UHqspJ+/btXRXLokWLrE2bNm7dnDlz3E5q165d0sdpRynQ6AT/4osvpjTR3OLFi90onskCjfiBn0oSHSg6KFQl5cONQiDhBgCQKXJd/aSTpi9dURuUM888011XYEhlcsC8aNasmety3b9/f1uwYIGbkkHj4mgsGd/zSTMMN23a1N3vA41GPlapxVNPPeVu+2HufQ8tjXnz5JNP2rJly9z8LSqBuuOOO9x4OOnKhxuNgBusltL7L8lDlwMAkF+5LqlRiFFPIbVB+eyzz1y3alm+fLmbCK+wTJ482QUZVXf5wfceeuih6P0au2blypVuQED58MMPXc8o8SVLwd5ceq1qcDtu3Di77rrrXDsUbXf//fe78JTugiU32icq6VLAU6BT76lUSq0AAAh1mxqdHG+++WZXDTVw4MDooHUavVdjpvzlL3+xTFNU49Tkhz5m1WuqakqvU9f9IH50BwcAZOw0CSi8D6QoqFRLpTYKqL5kSyU3an9TEsfnAQCkrz0lraGwpkA47rjj3C/6RFMmBBXWNAkoOKp2Uzd8HTxqB+VnkFbbG7XBUcBhpGIAQLpJKdS0atUqOgCeriebMqGwpklA4VBI1Vg3WjS/lEptFG4UcrQo2Kj9DRNoAgDSQcrTJPiu18U9ZQIKh0poVHKjRWPwqARH1VMKNwqqqppSCQ4BBwAQmmkSspsyobC6dKNo+XF4VMepulB9rqrzDAYcLfFzaAEAkFbj1CSiX/b33XefawSE8FB1oqqfNMeWwqy6wR9xxBGuTY6qqjTAn4IO498AAEqCMrkJLpo3adasWa6theZI6tatm02YMMF141a1hMZ7QbgDjhY1MlYJjhaV4KgkRz2qVHJTEkdbBgBkhpRDzciRI+2xxx6zTp062XvvvWcXX3yx9e3b1+bPn+8GrNNt2ltkZsDZt29fdAwcleBo0Ta+HQ7dxAEAJSrUaOLIZ555xs4//3w3rYC6bqvaYcmSJZy0MpgPL74NjkpsFHAUbBRy1NhYVLpHOxwAQIkINevXr49OJqkxa3SCUnUTv8IRpPY2Wg455BCrVq1atJpKXcUVdBSEVU3lQw7HDwCgyEONer0EB2TTiUnzCgHJqDrSj4OjxsbBairfZVx0XGlRGAIAoNBDjQbX69OnT7QRqE5OV1xxhTthBU2bNi3PLwaZU02lEhtfiqOQ4y81IKAvxWFOKgBAoYSa3r17x9zu2bNnrv4QEHPg/V9JnxaNZqxSHPWwUxWV5qVSryrNFaLSG1+SQ1UVAKBAQo26bgOFXYqjEY1V1amAo0UBR6U4vqrKt9kh5AAA4jEkLEpkW5yKFSu6pWrVqq6qSgFHpTm+qkolOqoS9W1xCDkAAEIN0qKqSovabynkqNu4r65SSY4uVZKjkKPtFHJokwMAmYdQg7Tjq6B8ryqV5CjkaPHVVepCrmos3/DYhx3a5QBAeBFqEJqSHFVXaYRjhRkfcnzjYx94gqU5WlTVRdABgHAg1CB0FFT8NA6iIONLc1R1pTFyFHR0qfW6X49RyKFEp/CpV5v2uV+Ct73g7fj7JFEQ1bpEi0rrgpcAwotQg9DTicyXzIhKc3QiVaBRyPFTO2jRdQUe3a/H+VIghR5dMnbO//gwElzi1wVDSTBMBANGMHj4+8TPI6f1wW2Dzxnkg49K6eJfS/C+4OvRor/j/4au+wVAeiLUICP5tjbBUbJ9iY6CjS/ZUfWVwo6u+x5X4k9+eh4fetI98ASDgAJAMCT4QBDkg4APHcEAGAyBPrQEw0miUOMDTX5LU5KVAgXfS3yo9ddVeqf7vGDQ8e+J0h6g5CLUAElKdOLDjl+Cva98d3NfjeUFQ09wKepSAH9iT1SC4YNLsCTFv2bf1kizrPuG1tm9p5JUteODkqS6v/1nrH2ixV8Pfs4Kt8F95kMOgQcoOQg1QB7DjvhwEDwh+l//vgTAr/PBwguGifhSjOxeS7AkIr5UIj5YBEtJ/N/wY/v4NkTBYBK8nknVMNl9xuJLdvzn6UvyFHQIPEDJkTahZtu2bTZ48GB76aWX3BfuhRdeaA8++GC2k2qefvrp9vbbb8es+9Of/mTjx4+P3v7qq69s4MCB9uabb7rn0nQQY8eOdV/2QE58EEh2MowvGQmWkASrRnRi9PclquoJNpaNbwsSH1oSrY/fBvmvrsxr4PFVdcHqS//ZAMiftDlzX3rppbZhwwabNWuW+9Xbt29fGzBggP3jH//I9nH9+/e3W2+9NXpb3X49fdGce+65VqtWLXvvvffc8/fq1cudoO64445CfT/IDD6A+N5VyMzAE6zWCo6rpOtqx+NDj+gyvvSMYBpOiRq853abvDwm1W3iJdsmp+fXcZ7o/0dhyIqk8k6K2YoVK6x58+a2cOFCa9u2rVs3c+ZMO+ecc2z9+vVWp06dpCU1rVq1sgceeCDh/a+++qr99re/tW+++cZq1qzp1qkUZ+jQobZly5aUPwRNvqiZp3fs2OHmLgKAVPiwEyzF84sPPsGSveDiBUvecupV5q8na5Sd11K8gjoZ57RdYd0XfN/xt+Mfl922ub0t+b1dkNvEy27ohNyu04+6GjVqJHzOgjyHpkVJzbx581w3XB9opFOnTu4/7Pvvv28XXHBB0sdOnjzZnn32WVcac95559mIESOipTV63uOPPz4aaKRLly6uOmr58uXWunXrhM/pJ1sMfiAAkFupdCFPFGjiu84He3cFqzuD1V7BkiAfirI78ft1qQadnE6k2T1PotAVDF7xwS1+20Tbpfq4RK8t1fsK43ZxPadXUEG3uKRFqNm4caNLeEGqh9YQ+bovmT/84Q9Wv359V5KzdOlSVwKzcuVKmzZtWvR5g4FG/O3snldtbkaPHp3PdwUAuavCzI/4YBO8jL+e6LZ/Ldm9zpzW5SY8AGkXaoYNG2Z33XVXjlVPeaU2N55KZGrXrm0dO3a0zz//3I466qg8P+/w4cPt+uuvjympqVu3bp6fDwAKW6IqAiBsijXUDBkyxPr06ZPtNo0aNXJVR5s3b45Zr4ZH6hGl+1LVrl07d7l69WoXavTYBQsWxGyzadMmd5nd82oGaC0AAKDkKNZQU716dbfkpH379rZ9+3ZbtGiRtWnTxq2bM2eOqxf2QSUVixcvdpcqsfHPe/vtt7vA5Ku31LtKDZXUMBkAAKSPtBgYoVmzZnbWWWe57tkqWZk7d65dddVV1qNHj2jPp6+//tqaNm0aLXlRFdOYMWNcEFq7dq29+OKLrrv2qaeeai1atHDbdO7c2YWXyy67zJYsWWKvvfaa3XzzzTZo0CBKYgAASDNpEWp8LyaFFrWJUVfuU045xR5//PHo/Rq7Ro2ANT+PqDv2G2+84YKLHqeqLg3Yp8H7PDW8e/nll92lSm169uzpgk9wXBsAAJAe0mKcmpKOcWoAACj+c2jalNQAAABkh1ADAABCgVADAABCgVADAABCgVADAABCgVADAABCgVADAABCgVADAABCgVADAABCgVADAABCgVADAABCgVADAABCgVADAABCgVADAABCgVADAABCgVADAABCgVADAABCgVADAABCgVADAABCgVADAABCgVADAABCgVADAABCgVADAABCIW1CzbZt2+zSSy+1ypUr26GHHmqXX3657dq1K+n2a9eutaysrITL888/H90u0f1TpkwponcFAAAKShlLEwo0GzZssFmzZtlPP/1kffv2tQEDBtg//vGPhNvXrVvXbR/0+OOP2z333GNnn312zPoJEybYWWedFb2t0AQAANJLWoSaFStW2MyZM23hwoXWtm1bt+7hhx+2c845x+69916rU6fOAY8pXbq01apVK2bd9OnT7fe//70dfPDBMesVYuK3BQAA6SUtqp/mzZvngocPNNKpUycrVaqUvf/++yk9x6JFi2zx4sWu2ireoEGDrFq1anbSSSfZ008/bZFIJNvn2rt3r+3cuTNmAQAAxSstSmo2btxoNWrUiFlXpkwZO+yww9x9qXjqqaesWbNm1qFDh5j1t956q51xxhlWsWJFe/311+3KK690bXWuvvrqpM81duxYGz16dB7fDQAACF1JzbBhw5I25vXLp59+mu+/8+OPP7q2N4lKaUaMGGEnn3yytW7d2oYOHWo33XSTa3eTneHDh9uOHTuiy7p16/L9GgEAQBqX1AwZMsT69OmT7TaNGjVy7V02b94cs/7nn392PaJSaQvzwgsv2A8//GC9evXKcdt27drZmDFjXBVTuXLlEm6j9cnuAwAAGRhqqlev7pactG/f3rZv3+7axbRp08atmzNnju3fv9+FkFSqns4///yU/pba3VStWpXQAgBAmkmLNjVqC6Mu1/3797fx48e7Lt1XXXWV9ejRI9rz6euvv7aOHTvaM8884xr8eqtXr7Z33nnHXnnllQOe96WXXrJNmzbZr371KytfvrzrLn7HHXfYDTfcUKTvDwAAZEiokcmTJ7sgo+CiXk8XXnihPfTQQ9H7FXRWrlzpqpmC1JvpyCOPtM6dOx/wnAcddJCNGzfOrrvuOtfjqXHjxnb//fe78AQAANJLViSn/svIkbp0V6lSxTUa1ojHAACg6M+haTFODQAAQE4INQAAIBQINQAAIBQINQAAIBQINQAAIBQINQAAIBQINQAAIBQINQAAIBQINQAAIBQINQAAIBQINQAAIBQINQAAIBQINQAAIBQINQAAIBQINQAAIBQINQAAIBQINQAAIBQINQAAIBQINQAAIBQINQAAIBQINQAAIBQINQAAIBQINQAAIBTSJtTcfvvt1qFDB6tYsaIdeuihKT0mEonYyJEjrXbt2lahQgXr1KmTrVq1Kmabbdu22aWXXmqVK1d2z3v55Zfbrl27CuldAAAAy/RQs2/fPrv44ott4MCBKT/m7rvvtoceesjGjx9v77//vlWqVMm6dOlie/bsiW6jQLN8+XKbNWuWvfzyy/bOO+/YgAEDCuldAACAwpIVUXFGGpk4caJde+21tn379my309uqU6eODRkyxG644Qa3bseOHVazZk33HD169LAVK1ZY8+bNbeHChda2bVu3zcyZM+2cc86x9evXu8enYufOnValShX3/CrxAQAAVuTn0DIWUmvWrLGNGze6KidPO61du3Y2b948F2p0qSonH2hE25cqVcqV7FxwwQUJn3vv3r1u8fRB+A8GAACkzp87C6KMJbShRoFGVDITpNv+Pl3WqFEj5v4yZcrYYYcdFt0mkbFjx9ro0aMPWF+3bt0CevUAAGSWb7/91hU+pG2oGTZsmN11113ZbqMqoqZNm1pJMnz4cLv++uujt1UVVr9+ffvqq6/y/YFkUjJXCFy3bh1Vduw3jrUSiP+j7LOiotqOevXquQKF/CrWUKP2Ln369Ml2m0aNGuXpuWvVquUuN23a5Ho/ebrdqlWr6DabN2+OedzPP//sekT5xydSrlw5t8RToKFNTe5of7HPco/9xj4rKhxr7LOioqYfaR1qqlev7pbC0LBhQxdMZs+eHQ0x+uWhtjK+B1X79u1dKcuiRYusTZs2bt2cOXNs//79ru0NAABIH2nTpVtVO4sXL3aXv/zyi7uuJTimjKqppk+f7q5nZWW5XlK33Xabvfjii/bxxx9br169XI+mbt26uW2aNWtmZ511lvXv398WLFhgc+fOtauuuso1Ik615xMAACgZ0qahsAbRmzRpUvR269at3eWbb75pp59+uru+cuXKaE8kuemmm2z37t1u3BmVyJxyyimuy3b58uWj20yePNkFmY4dO7qirwsvvNCNbZMbqooaNWpUwiopsM8KEsca+6yocKyxz9LxWEu7cWoAAADSuvoJAAAgO4QaAAAQCoQaAAAQCoQaAAAQCoSafBo3bpw1aNDA9ajS2DbqGo7/T7Oen3feea6LvLrZz5gxI2b3qJ26erZpgMQKFSq4ubdWrVqV0btQ03CceOKJdsghh7hpPDQEgXr2BWmm+UGDBtnhhx9uBx98sOu1p4ElM9Wjjz5qLVq0iA4UpzGoXn311ej97K/U3HnnndHhMNh3id1yyy1uHwWX4Kj3HGvJff3119azZ0/3vaXv++OPP94++OCDAj0fEGryYerUqW66BHVF+/DDD61ly5bWpUuXA0YpzmTqUq/9ovCXyN133+260I8fP94NjFipUiW3D/XFkKnefvttF1jmz59vs2bNsp9++sk6d+7s9qV33XXX2UsvvWTPP/+82/6bb76x3/3ud5apjjzySHdC1kCa+pI844wzrGvXrrZ8+XJ3P/srZwsXLrTHHnvMhcMg9t2Bjj32WNuwYUN0effdd9lfOfjuu+/s5JNPtoMOOsj94Pjkk0/svvvus6pVqxbs+UBdupE3J510UmTQoEHR27/88kukTp06kbFjx7JLE9DhNn369Ojt/fv3R2rVqhW55557ouu2b98eKVeuXOSf//wn+/D/bN682e27t99+O7qPDjrooMjzzz8f3UcrVqxw28ybN4/99n+qVq0aefLJJ9lfKfj+++8jTZo0icyaNSty2mmnRa655hqOtSRGjRoVadmyZcL7+L+Z3NChQyOnnHJK0vsL6nxASU0e7du3z/0qVPGYp8H7dHvevHl5fdqMsmbNGjcbenAfav4sVeOxD/8/P6Ckn+xNx51Kb4L7TcXfmhCO/WZuxPEpU6a4ki1VQ7G/cqaSwXPPPTfmmOJYS05VIqpS19yEl156qRvpnv2VPY3s37ZtW7v44otdtboG0H3iiScK/HxAqMmjrVu3ui/PmjVrxqzXbX0wyJnfT+zD5DQPmdo3qNj2uOOOi+63smXL2qGHHsqxF6CpUNS+SKOSXnHFFW7KlObNm7O/cqAAqOpzteVK9H+UYy2WTrITJ050o9OrLZdOxr/+9a/t+++/Z39l44svvnD7q0mTJvbaa6+5ORivvvrq6EwBBXU+SJtpEoBM/QW9bNmymDp7JHbMMce4+eBUsvXCCy9Y7969XXsjJLdu3Tq75pprXNut4PQxSO7ss8+OXlf7I4Wc+vXr23PPPecatyL5DzSV1Nxxxx3utkpq9N2m9jP6v1pQKKnJo2rVqlnp0qUP6HGi25odHDnz+4l9mJjmJHv55Zfd/GZqCBvcb6r+1HxmHHv/n0oUGjdubG3atHGlDmqg/uCDD7K/sqGqOXVsOOGEE6xMmTJuURBUY01d169kjrXsqcT06KOPttWrV3OsZUM9mlRyGqRJpX3VXUGdDwg1+fgC1Zfn7NmzY5KobqseHzlr2LChO1iD+3Dnzp2u1Xsm70O1qVagUfXJnDlz3H4K0nGnHgTB/aYu3/pyyOT9Fk//H/fu3cv+yoYm8lW1nUq4/KJf02on4q9zrGVv165d9vnnn7uTNv83k1MVevzQFJ999pkr5SrQ80HKTYpxgClTpriW2RMnTox88sknkQEDBkQOPfTQyMaNG9lbgV4VH330kVt0uN1///3u+pdffunuv/POO90++/e//x1ZunRppGvXrpGGDRtGfvzxx4zdhwMHDoxUqVIl8tZbb0U2bNgQXX744YfoNldccUWkXr16kTlz5kQ++OCDSPv27d2SqYYNG+Z6h61Zs8YdR7qdlZUVef3119397K/UBXs/se8ONGTIEPd/U8fa3LlzI506dYpUq1bN9VJkfyW3YMGCSJkyZSK33357ZNWqVZHJkydHKlasGHn22Wej2xTE+YBQk08PP/ywO7mULVvWdfGeP39+fp8yVN58800XZuKX3r17R7vxjRgxIlKzZk0XEDt27BhZuXJlJJMl2l9aJkyYEN1G/8mvvPJK121ZXwwXXHCBCz6Z6o9//GOkfv367v9h9erV3XHkA42wv/Ieath3sbp37x6pXbu2O9aOOOIId3v16tXsrxS89NJLkeOOO8591zdt2jTy+OOPx9xfEOeDLP2TerkOAABAyUSbGgAAEAqEGgAAEAqEGgAAEAqEGgAAEAqEGgAAEAqEGgAAEAqEGgAAEAqEGgAAEAqEGgB51qdPH+vWrVux7cHLLrssOutvGPdPjx497L777ivQ1wSEGaEGQEJZWVnZLrfccoubBXvixInFsgeXLFlir7zyil199dVW3NauXev2iSaBLEg333yz3X777bZjx44CfV4grMoU9wsAUDJt2LAhen3q1Kk2cuTImFl2Dz74YLcUl4cfftguvvjiYn0Nhe24446zo446yp599lkbNGhQcb8coMSjpAZAQrVq1YouVapUcSURwXUKE/HVK6effroNHjzYrr32WqtatarVrFnTnnjiCdu9e7f17dvXDjnkEGvcuLG9+uqrMX9r2bJldvbZZ7vn1GNUrbR169akn8wvv/xiL7zwgp133nkx6xs0aGC33Xab9erVyz1X/fr17cUXX7QtW7ZY165d3boWLVrYBx98EPO4f/3rX3bsscdauXLl3HPEV/lonaq5/vjHP7r3UK9ePXv88cej9zds2NBdtm7d2u0n7Yege++912rXrm2HH364Cyc//fRT9L6//e1v1qRJEytfvrx77xdddFHMY/Uep0yZwlEKpIBQA6BATZo0yapVq2YLFixwAWfgwIGuRKVDhw724YcfWufOnV1o+eGHH9z227dvtzPOOMMFAoWNmTNn2qZNm+z3v/990r+xdOlSVyXTtm3bA+7761//aieffLJ99NFHdu6557q/pZDTs2dP9/dV8qHbfi7fRYsWub+l9isff/yxq1YbMWLEAdVqCjr6e3reK6+80r0vX3Kl9ypvvPGGK+GaNm1a9HFvvvmmff755+5S+0bP659b71fVZ7feeqt7Lr33U089NebvnnTSSe759+7dm49PBcgQuZrTG0BGmjBhQqRKlSoHrO/du3eka9eu0dunnXZa5JRTTone/vnnnyOVKlWKXHbZZdF1GzZsUJqIzJs3z90eM2ZMpHPnzjHPu27dOrfNypUrE76e6dOnR0qXLh3Zv39/zPr69etHevbsecDfGjFiRHSd/q7W6T75wx/+EDnzzDNjnufGG2+MNG/ePOnz6u/WqFEj8uijj7rba9ascc/50UcfHbB/9FjtB+/iiy+OdO/e3V3/17/+FalcuXJk586dkWSWLFninnvt2rVJtwHwP5TUAChQqt7xSpcu7apcjj/++Og6VbHI5s2bow1+VYrh2+hoadq0qbtPJRyJ/Pjjj66qSFU92f19/7ey+/srVqxwJTtBur1q1SpXzZXoeX1VnH+O7KhaS/vBUzWUf9yZZ57pqsgaNWrkSpQmT54cLcHyKlSo4C7j1wM4EKEGQIE66KCDYm4rAATX+SCyf/9+d7lr1y7XbkQ9h4KLQkV8VYyn6i2d5Pft25ft3/d/K7u/n5/3lcpzZPc4tc9Rldg///lPF3bUGLtly5auSs7btm2bu6xevXquXi+QiQg1AIrVCSecYMuXL3eNcdWIOLhUqlQp4WNatWrlLj/55JN8//1mzZrZ3LlzY9bp9tFHHx1TwpKdsmXLustgyU6qypQpY506dbK7777btRVS9/A5c+bENKI+8sgjXZADkD1CDYBipd5AKo245JJLbOHCha7K6bXXXnO9pZKFBJVaKAy9++67+f77Q4YMsdmzZ9uYMWPss88+c415H3nkEbvhhhtSfo4aNWq4aiLfyDnVcWVefvlle+ihh1zJ1JdffmnPPPOMK8U55phjotv897//dY2rAeSMUAOgWNWpU8eVjCjA6OSt9i/qEn7ooYdaqVLJv6L69evn2qDkl8LRc88957pNa1wYVQGpN5K6q+emtEXh5LHHHnPvR93HU6H3qJ5S6v2lEqPx48e7qii1w5E9e/bYjBkzrH///nl+f0AmyVJr4eJ+EQCQW2osrBINDQzYvn37UO7ARx991KZPn26vv/56cb8UIC1QUgMgLam6R9U12Q3Sl+7UyFgjJwNIDSU1AAAgFCipAQAAoUCoAQAAoUCoAQAAoUCoAQAAoUCoAQAAoUCoAQAAoUCoAQAAoUCoAQAAoUCoAQAAFgb/DxmwuIuYnvAiAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 600x350 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(figsize=[6, 3.5])\n",
    "ax.axhline(y=0, linestyle=':', color='k')\n",
    "ax.fill_between(rd_results['time'], rd_results['lcl'], rd_results['ucl'], color='k', alpha=0.1)\n",
    "ax.plot(rd_results['time'], rd_results['rd'], color='k')\n",
    "ax.set_xlim([0, 60])\n",
    "ax.set_ylim([-1, 1])\n",
    "ax.set_xlabel(\"Time (months)\")\n",
    "ax.set_ylabel(\"Risk Difference\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f09786da-2d97-4fc4-8fe7-6a7e6789997e",
   "metadata": {},
   "source": [
    "The results between the two approaches are largely the same. These examples highlight how pooled logistic regression can be used to flexibily conduct causal survival analyses. Importantly, the estimating equation approach taken here can substantially reduce run-times and avoids reliance on the bootstrap for inference. Further details on this general approach and the benefits of estimating equations can be found in the following references.\n",
    "\n",
    "## References\n",
    "\n",
    "Abbott RD. (1985). Logistic regression in survival analysis. *American Journal of Epidemiology*, 121(3), 465-471.\n",
    "\n",
    "D'Agostino RB, Lee ML, Belanger AJ, Cupples LA, Anderson K, & Kannel WB. (1990). Relation of pooled logistic regression to time dependent Cox regression analysis: the Framingham Heart Study. *Statistics in Medicine*, 9(12), 1501-1515.\n",
    "\n",
    "Dumas E, & Stensrud MJ. (2025). How hazard ratios can mislead and why it matters in practice *European Journal of Epidemiology*, 40(6), 603-609.\n",
    "\n",
    "Hernán MA. (2010). The hazards of hazard ratios. *Epidemiology*, 21(1), 13-15.\n",
    "\n",
    "Zivich PN, Cole SR, Shook-Sa BE, DeMonte JB, & Edwards JK. (2025). Estimating equations for survival analysis with pooled logistic regression. *arXiv:2504.13291*\n",
    "\n",
    "Zivich PN, Klose M, DeMonte JB, Shook-Sa BE, Cole SR, & Edwards JK. (2026). An Improved Pooled Logistic Regression Implementation. *Epidemiology*, In-Press."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.13.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}