{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "cf4a3599",
   "metadata": {},
   "source": [
    "# Cluster-randomised field trial through savings groups: randomisation, balance, analysis\n",
    "\n",
    "**Worked example on synthetic data, by Ahmed Beddi.** The 80 savings groups, their members and the\n",
    "endline outcomes are generated below from a seeded random-number generator. No real programme,\n",
    "group, institution or dataset is involved.\n",
    "\n",
    "The intervention (illustrative): eight short financial-capability sessions delivered by a trained\n",
    "facilitator during regular savings-group meetings, with a pictorial rule-of-thumb card (keep business\n",
    "and household money apart; name one savings goal; deposit on meeting day). The group is the unit of\n",
    "randomisation because sessions are delivered to the whole group.\n",
    "\n",
    "Contents: (1) synthetic cluster list, identical to the website demo; (2) stratified randomisation\n",
    "and balance; (3) a pre-specified rerandomisation rule; (4) power and minimum detectable effect for the\n",
    "cluster design, checked by simulation; (5) a simulated endline and the pre-specified ANCOVA with\n",
    "cluster-robust standard errors; (6) randomisation inference; (7) heterogeneity; (8) implementation\n",
    "fidelity."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "dcc95da7",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-11-03T14:58:31.859566Z",
     "iopub.status.busy": "2026-11-03T14:58:31.859361Z",
     "iopub.status.idle": "2026-11-03T14:58:33.138707Z",
     "shell.execute_reply": "2026-11-03T14:58:33.137716Z"
    }
   },
   "outputs": [],
   "source": [
    "import math\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import statsmodels.formula.api as smf\n",
    "from scipy import stats\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "pd.set_option('display.width', 120)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cdd344fd",
   "metadata": {},
   "source": [
    "## 1. Synthetic cluster list\n",
    "\n",
    "The generator is a line-by-line port of the TypeScript used by the website demo (mulberry32 PRNG,\n",
    "Box–Muller normals, same draw order and rounding), so the list below is the list the demo shows."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "6ff29e51",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-11-03T14:58:33.141019Z",
     "iopub.status.busy": "2026-11-03T14:58:33.140646Z",
     "iopub.status.idle": "2026-11-03T14:58:33.161220Z",
     "shell.execute_reply": "2026-11-03T14:58:33.159895Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " id  region  rural  members  age  female   mm  savings  agentKm\n",
      "  1       1      1       21  4.6    0.66 0.16     1740      5.7\n",
      "  2       0      1       14  1.5    0.39 0.23     1650      7.7\n",
      "  3       2      1       15  1.6    0.64 0.14     3410      2.2\n",
      "  4       0      0       22  5.0    0.54 0.49     3360      2.4\n",
      "  5       2      1       18  4.2    0.84 0.02     1000      6.7\n",
      "  6       2      1       21  0.5    0.71 0.29     2380      3.0\n",
      "  7       1      1       16  6.9    0.73 0.39     1230      5.3\n",
      "  8       0      0       24  4.3    0.73 0.02     3120      1.0\n",
      "\n",
      "clusters: 80  members: 1474  mean size: 18.42  CV of size: 0.232\n"
     ]
    }
   ],
   "source": [
    "M32 = 0xFFFFFFFF\n",
    "\n",
    "def mulberry32(seed):\n",
    "    a = seed & M32\n",
    "    def rng():\n",
    "        nonlocal a\n",
    "        a = (a + 0x6D2B79F5) & M32\n",
    "        t = ((a ^ (a >> 15)) * (1 | a)) & M32\n",
    "        t = ((t + (((t ^ (t >> 7)) * (61 | t)) & M32)) & M32) ^ t\n",
    "        return ((t ^ (t >> 14)) & M32) / 4294967296\n",
    "    return rng\n",
    "\n",
    "def normal(rng):\n",
    "    u1 = 1 - rng(); u2 = rng()\n",
    "    return math.sqrt(-2 * math.log(u1)) * math.cos(2 * math.pi * u2)\n",
    "\n",
    "r1 = lambda x: math.floor(x * 10 + 0.5) / 10\n",
    "r2 = lambda x: math.floor(x * 100 + 0.5) / 100\n",
    "clip = lambda x, lo, hi: min(hi, max(lo, x))\n",
    "\n",
    "DATA_SEED, J = 2026, 80\n",
    "\n",
    "def make_clusters(seed=DATA_SEED, n=J):\n",
    "    rng = mulberry32(seed); rows = []\n",
    "    for i in range(n):\n",
    "        u = rng()\n",
    "        region = 0 if u < 0.4 else (1 if u < 0.75 else 2)\n",
    "        rural = 1 if rng() < [0.5, 0.7, 0.85][region] else 0\n",
    "        members = 12 + math.floor(rng() * 14)\n",
    "        age = r1(max(0.5, 3 + 2 * normal(rng)))\n",
    "        female = r2(clip(0.7 + 0.15 * normal(rng), 0.2, 1))\n",
    "        mm = r2(clip(0.35 - 0.12 * rural + 0.15 * normal(rng), 0.02, 0.95))\n",
    "        savings = math.floor(math.exp(math.log(2000) - 0.25 * rural + 0.4 * normal(rng)) / 10 + 0.5) * 10\n",
    "        agent_km = r1(2 + 6 * rng() if rural else 0.5 + 2 * rng())\n",
    "        rows.append(dict(id=i + 1, region=region, rural=rural, members=members, age=age, female=female,\n",
    "                         mm=mm, savings=savings, agentKm=agent_km))\n",
    "    return pd.DataFrame(rows)\n",
    "\n",
    "cl = make_clusters()\n",
    "cl.to_csv('synthetic_savings_groups.csv', index=False)\n",
    "print(cl.head(8).to_string(index=False))\n",
    "print('\\nclusters:', len(cl), ' members:', cl.members.sum(), ' mean size:', round(cl.members.mean(), 2),\n",
    "      ' CV of size:', round(cl.members.std(ddof=0) / cl.members.mean(), 3))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "261a4551",
   "metadata": {},
   "source": [
    "## 2. Stratified randomisation and balance\n",
    "\n",
    "Strata are formed from the selected variables (region, rural, large group = members at or above the\n",
    "median, high savings = savings at or above the median). Within each stratum, in sorted key order, the\n",
    "clusters are shuffled (Fisher–Yates) and the first half assigned to treatment; in an odd-sized stratum the\n",
    "remaining cluster is treated with probability ½. Balance is reported as the standardised difference\n",
    "\n",
    "$$d = \\frac{\\bar x_T - \\bar x_C}{\\sqrt{(s_T^2 + s_C^2)/2}} .$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "12d6e44a",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-11-03T14:58:33.163773Z",
     "iopub.status.busy": "2026-11-03T14:58:33.163442Z",
     "iopub.status.idle": "2026-11-03T14:58:33.190484Z",
     "shell.execute_reply": "2026-11-03T14:58:33.189120Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "strata: 6   treated clusters: 40   control: 40\n",
      "covariate   mean_T   mean_C  std_diff\n",
      "    rural    0.700    0.650     0.106\n",
      "  members   17.750   19.100    -0.316\n",
      "      age    2.997    3.622    -0.390\n",
      "   female    0.690    0.695    -0.038\n",
      "       mm    0.237    0.257    -0.134\n",
      "  savings 1697.000 1754.750    -0.083\n",
      "  agentKm    3.610    4.283    -0.300\n",
      "max |d| = 0.39\n"
     ]
    }
   ],
   "source": [
    "COVARS = ['rural', 'members', 'age', 'female', 'mm', 'savings', 'agentKm']\n",
    "\n",
    "def strata_keys(cl, vars_):\n",
    "    m_mem, m_sav = cl.members.median(), cl.savings.median()\n",
    "    cols = {'region': cl.region, 'rural': cl.rural,\n",
    "            'large': (cl.members >= m_mem).astype(int), 'highSav': (cl.savings >= m_sav).astype(int)}\n",
    "    return [tuple(int(cols[v].iloc[i]) for v in vars_) for i in range(len(cl))]\n",
    "\n",
    "def assign_once(keys, rng):\n",
    "    groups = {}\n",
    "    for i, k in enumerate(keys):\n",
    "        groups.setdefault(k, []).append(i)\n",
    "    z = np.zeros(len(keys), int)\n",
    "    for k in sorted(groups):\n",
    "        idx = list(groups[k])\n",
    "        for i in range(len(idx) - 1, 0, -1):\n",
    "            j = math.floor(rng() * (i + 1)); idx[i], idx[j] = idx[j], idx[i]\n",
    "        h = len(idx) // 2\n",
    "        for i in range(h):\n",
    "            z[idx[i]] = 1\n",
    "        if len(idx) % 2 == 1 and rng() < 0.5:\n",
    "            z[idx[-1]] = 1\n",
    "    return z\n",
    "\n",
    "def balance(cl, z):\n",
    "    out = []\n",
    "    for c in COVARS:\n",
    "        xt, xc = cl.loc[z == 1, c], cl.loc[z == 0, c]\n",
    "        sd = math.sqrt((xt.var(ddof=1) + xc.var(ddof=1)) / 2)\n",
    "        out.append(dict(covariate=c, mean_T=xt.mean(), mean_C=xc.mean(), std_diff=(xt.mean() - xc.mean()) / sd))\n",
    "    return pd.DataFrame(out)\n",
    "\n",
    "def randomise(cl, vars_, seed, threshold=None, cap=5000):\n",
    "    rng = mulberry32(seed); keys = strata_keys(cl, vars_); draws = 0\n",
    "    while True:\n",
    "        z = assign_once(keys, rng); b = balance(cl, z); draws += 1\n",
    "        if threshold is None or b.std_diff.abs().max() < threshold or draws >= cap:\n",
    "            return z, b, draws, len(set(keys))\n",
    "\n",
    "STRATA, SEED = ['region', 'rural'], 1\n",
    "z, bal, draws, ns = randomise(cl, STRATA, SEED)\n",
    "print(f'strata: {ns}   treated clusters: {z.sum()}   control: {len(z) - z.sum()}')\n",
    "print(bal.round(3).to_string(index=False))\n",
    "print('max |d| =', round(bal.std_diff.abs().max(), 3))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cc489015",
   "metadata": {},
   "source": [
    "Region and rural status are balanced by construction (every stratum is split in half). The other\n",
    "covariates are balanced only in expectation; with 40 clusters per arm the standard error of a\n",
    "standardised difference is about $\\sqrt{1/40 + 1/40} \\approx 0.22$, so differences of 0.2–0.3 on some\n",
    "covariate are expected in any single draw. That is not a failure of randomisation, and it is not a\n",
    "reason to redraw *unless a rule for redrawing was fixed in advance*.\n",
    "\n",
    "How often does a single stratified draw show a large imbalance? Over 2,000 draws:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "a96bda80",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-11-03T14:58:33.193573Z",
     "iopub.status.busy": "2026-11-03T14:58:33.193365Z",
     "iopub.status.idle": "2026-11-03T14:58:39.538866Z",
     "shell.execute_reply": "2026-11-03T14:58:39.538163Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "P(max |d| < 0.1) = 0.002   -> expected draws under the rule: 500.0\n",
      "P(max |d| < 0.2) = 0.073   -> expected draws under the rule: 13.6\n",
      "P(max |d| < 0.25) = 0.191   -> expected draws under the rule: 5.2\n",
      "P(max |d| < 0.3) = 0.349   -> expected draws under the rule: 2.9\n"
     ]
    }
   ],
   "source": [
    "rng = mulberry32(99); keys = strata_keys(cl, STRATA)\n",
    "mx = np.array([balance(cl, assign_once(keys, rng)).std_diff.abs().max() for _ in range(2000)])\n",
    "for th in [0.1, 0.2, 0.25, 0.3]:\n",
    "    print(f'P(max |d| < {th}) = {np.mean(mx < th):.3f}   -> expected draws under the rule: {1 / max(np.mean(mx < th), 1e-9):.1f}')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1bcc4c29",
   "metadata": {},
   "source": [
    "## 3. Pre-specified rerandomisation rule\n",
    "\n",
    "Rule written into the pre-analysis plan before the list is seen: *draw stratified assignments from the\n",
    "seeded stream; accept the first whose maximum absolute standardised difference across the seven\n",
    "covariates is below 0.25.* Because the rule is fixed and the seed logged, anyone can reproduce the\n",
    "accepted assignment, and inference can use the same rule (section 6). Morgan and Rubin (2012) show\n",
    "that rerandomisation under a pre-specified rule is valid if the analysis respects it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "c0d927c4",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-11-03T14:58:39.541093Z",
     "iopub.status.busy": "2026-11-03T14:58:39.540904Z",
     "iopub.status.idle": "2026-11-03T14:58:39.594781Z",
     "shell.execute_reply": "2026-11-03T14:58:39.593549Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "accepted after 11 draws; max |d| = 0.234\n",
      "covariate   mean_T   mean_C  std_diff\n",
      "    rural    0.675    0.675     0.000\n",
      "  members   18.550   18.300     0.058\n",
      "      age    3.448    3.173     0.169\n",
      "   female    0.685    0.700    -0.111\n",
      "       mm    0.265    0.229     0.234\n",
      "  savings 1651.750 1800.000    -0.215\n",
      "  agentKm    4.162    3.730     0.192\n"
     ]
    }
   ],
   "source": [
    "z, bal, draws, ns = randomise(cl, STRATA, SEED, threshold=0.25)\n",
    "print(f'accepted after {draws} draws; max |d| = {bal.std_diff.abs().max():.3f}')\n",
    "print(bal.round(3).to_string(index=False))\n",
    "cl['treat'] = z\n",
    "cl.to_csv('synthetic_savings_groups.csv', index=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "009beeab",
   "metadata": {},
   "source": [
    "## 4. Power and minimum detectable effect for the cluster design\n",
    "\n",
    "For $J$ clusters split in proportion $P$ to treatment, mean cluster size $\\bar m$, intra-cluster\n",
    "correlation $\\rho$ of the outcome and share $R^2$ of variance explained by the baseline covariate\n",
    "(assumed equal at both levels for simplicity), the minimum detectable effect in standard-deviation\n",
    "units is\n",
    "\n",
    "$$\\text{MDE} = (t_{1-\\alpha/2,\\,J-2} + t_{1-\\beta,\\,J-2})\\sqrt{\\frac{(1-R^2)\\,\\big[\\rho + (1-\\rho)/\\bar m\\big]}{P(1-P)\\,J}} .$$\n",
    "\n",
    "Unequal cluster sizes cost power; a standard adjustment (Eldridge, Ashby and Kerry, 2006) inflates the\n",
    "design effect to $1 + \\big[(\\text{CV}^2 + 1)\\bar m - 1\\big]\\rho$, with CV the coefficient of variation of cluster size."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "25cf694e",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-11-03T14:58:39.596841Z",
     "iopub.status.busy": "2026-11-03T14:58:39.596565Z",
     "iopub.status.idle": "2026-11-03T14:58:39.615760Z",
     "shell.execute_reply": "2026-11-03T14:58:39.613943Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "rho = 0.05:  MDE no covariate = 0.205 SD   with R2 = 0.4: 0.159 SD\n",
      "rho = 0.10:  MDE no covariate = 0.249 SD   with R2 = 0.4: 0.193 SD\n",
      "rho = 0.20:  MDE no covariate = 0.320 SD   with R2 = 0.4: 0.248 SD\n",
      "\n",
      "clusters needed for MDE 0.20 SD (rho 0.10, R2 0.4): 76\n",
      "adding members per group instead (J = 80, rho 0.10, R2 0.4): {10: np.float64(0.214), 20: np.float64(0.187), 40: np.float64(0.172), 1000: np.float64(0.156)}\n"
     ]
    }
   ],
   "source": [
    "def mde_cluster(J, m, rho, R2=0.0, P=0.5, alpha=0.05, power=0.8, cv=0.0):\n",
    "    df = J - 2\n",
    "    M = stats.t.ppf(1 - alpha / 2, df) + stats.t.ppf(power, df)\n",
    "    de_per_m = (1 + ((cv ** 2 + 1) * m - 1) * rho) / m   # = rho + (1-rho)/m when cv = 0\n",
    "    return M * math.sqrt((1 - R2) * de_per_m / (P * (1 - P) * J))\n",
    "\n",
    "mbar = cl.members.mean(); cv = cl.members.std(ddof=0) / mbar\n",
    "for rho in [0.05, 0.10, 0.20]:\n",
    "    print(f'rho = {rho:.2f}:  MDE no covariate = {mde_cluster(J, mbar, rho, cv=cv):.3f} SD   '\n",
    "          f'with R2 = 0.4: {mde_cluster(J, mbar, rho, 0.4, cv=cv):.3f} SD')\n",
    "print('\\nclusters needed for MDE 0.20 SD (rho 0.10, R2 0.4):',\n",
    "      next(j for j in range(10, 400, 2) if mde_cluster(j, mbar, 0.10, 0.4, cv=cv) <= 0.20))\n",
    "print('adding members per group instead (J = 80, rho 0.10, R2 0.4):',\n",
    "      {m: round(mde_cluster(J, m, 0.10, 0.4), 3) for m in [10, 20, 40, 1000]})"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c0700551",
   "metadata": {},
   "source": [
    "The last line is the practical lesson of cluster trials: beyond a point, more members per group buys\n",
    "almost nothing. With $\\rho = 0.10$ the MDE can never fall below $\\approx (t+t)\\sqrt{(1-R^2)\\rho/(P(1-P)J)}$\n",
    "however large the groups. More *clusters* is the only lever that keeps working.\n",
    "\n",
    "## 5. Simulated endline and the pre-specified analysis\n",
    "\n",
    "Data-generating process (synthetic; parameters chosen for illustration, not estimated from any\n",
    "programme). For member $i$ in group $j$, baseline monthly savings $x_{ij}$ (local currency units) varies\n",
    "around the group mean; endline savings (average monthly contribution recorded in the group ledger,\n",
    "months 4–9) is\n",
    "\n",
    "$$y_{ij} = 0.6\\,x_{ij} + 0.4\\,\\bar s_j + \\tau_{j} T_j\\cdot \\tfrac{\\text{sessions}_j}{8} + u_j + e_{ij},\\quad\n",
    "u_j\\sim N(0, 250^2),\\ e_{ij}\\sim N(0, 750^2),$$\n",
    "\n",
    "with a true effect of 150 units in peri-urban groups and 250 in rural groups (full dose), and\n",
    "sessions delivered drawn per treated group. Control groups receive the same ledger-photography visits\n",
    "as treated groups, so measurement is identical across arms."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "02f8c33b",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-11-03T14:58:39.618168Z",
     "iopub.status.busy": "2026-11-03T14:58:39.617927Z",
     "iopub.status.idle": "2026-11-03T14:58:39.682637Z",
     "shell.execute_reply": "2026-11-03T14:58:39.681206Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "members: 1474   control endline mean: 1713.1   control SD: 1200.5\n"
     ]
    }
   ],
   "source": [
    "def simulate_members(cl, z, seed=7, tau_urban=150.0, tau_rural=250.0, sd_u=250.0, sd_e=750.0, full_dose=False):\n",
    "    g = np.random.default_rng(seed); rows = []\n",
    "    sessions = np.where(z == 1, g.binomial(8, 0.85, len(cl)), 0)\n",
    "    if full_dose:\n",
    "        sessions = np.where(z == 1, 8, 0)\n",
    "    attendance = np.where(z == 1, g.beta(8, 2, len(cl)), 0.0)\n",
    "    u = g.normal(0, sd_u, len(cl))\n",
    "    for j, c in cl.reset_index(drop=True).iterrows():\n",
    "        m = int(c.members)\n",
    "        fem = (g.random(m) < c.female).astype(int)\n",
    "        x = c.savings * np.exp(g.normal(-0.08, 0.4, m))\n",
    "        tau = tau_rural if c.rural else tau_urban\n",
    "        y = 0.6 * x + 0.4 * c.savings + tau * z[j] * sessions[j] / 8 + u[j] + g.normal(0, sd_e, m)\n",
    "        for k in range(m):\n",
    "            rows.append(dict(group=int(c.id), treat=int(z[j]), region=int(c.region), rural=int(c.rural),\n",
    "                             female=fem[k], x0=x[k], y=y[k], sessions=int(sessions[j]),\n",
    "                             attendance=attendance[j]))\n",
    "    return pd.DataFrame(rows)\n",
    "\n",
    "df = simulate_members(cl, z)\n",
    "df['stratum'] = df.region.astype(str) + '_' + df.rural.astype(str)\n",
    "df.to_csv('synthetic_members_endline.csv', index=False, float_format='%.1f')\n",
    "sd_c = df.loc[df.treat == 0, 'y'].std()\n",
    "print('members:', len(df), '  control endline mean:', round(df.loc[df.treat == 0, 'y'].mean(), 1), '  control SD:', round(sd_c, 1))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d9a6559f",
   "metadata": {},
   "source": [
    "**Pre-specified primary model (ANCOVA).** Intention-to-treat at the group level:\n",
    "\n",
    "$$y_{ij} = \\alpha + \\tau T_j + \\beta x_{ij} + \\sum_s \\gamma_s S_{sj} + \\varepsilon_{ij},$$\n",
    "\n",
    "estimated by OLS with standard errors clustered by savings group (CR1) and t critical values with\n",
    "$G-1$ degrees of freedom ($G$ = number of groups). The ANCOVA uses the baseline\n",
    "outcome, which McKenzie (2012) shows is more efficient than difference-in-differences when outcomes\n",
    "are not highly autocorrelated."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "75d1a395",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-11-03T14:58:39.684750Z",
     "iopub.status.busy": "2026-11-03T14:58:39.684535Z",
     "iopub.status.idle": "2026-11-03T14:58:39.726237Z",
     "shell.execute_reply": "2026-11-03T14:58:39.725062Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ANCOVA, cluster-robust:  tau = 141.2  SE = 76.5  p = 0.0689  95% CI [-11.2, 293.6]\n",
      "same model, ignoring clustering: SE = 43.2  (too small by a factor 1.77)\n",
      "effect in control-SD units: 0.118\n",
      "true average effect among treated members (given sessions delivered): 177.7\n"
     ]
    }
   ],
   "source": [
    "fit = smf.ols('y ~ treat + x0 + C(stratum)', data=df).fit(cov_type='cluster', cov_kwds={'groups': df.group}, use_t=True)\n",
    "naive = smf.ols('y ~ treat + x0 + C(stratum)', data=df).fit(cov_type='HC1')\n",
    "tau_hat, se_hat = fit.params['treat'], fit.bse['treat']\n",
    "print(f'ANCOVA, cluster-robust:  tau = {tau_hat:.1f}  SE = {se_hat:.1f}  p = {fit.pvalues[\"treat\"]:.4f}  '\n",
    "      f'95% CI [{fit.conf_int().loc[\"treat\", 0]:.1f}, {fit.conf_int().loc[\"treat\", 1]:.1f}]')\n",
    "print(f'same model, ignoring clustering: SE = {naive.bse[\"treat\"]:.1f}  (too small by a factor {se_hat / naive.bse[\"treat\"]:.2f})')\n",
    "print(f'effect in control-SD units: {tau_hat / sd_c:.3f}')\n",
    "true_avg = (df.treat * np.where(df.rural == 1, 250, 150) * df.sessions / 8).sum() / df.treat.sum()\n",
    "print(f'true average effect among treated members (given sessions delivered): {true_avg:.1f}')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "08c59fed",
   "metadata": {},
   "source": [
    "Estimated residual ICC and $R^2$ from the control arm (what a real trial would use to update the\n",
    "power calculation for the next round):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "9427510f",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-11-03T14:58:39.728825Z",
     "iopub.status.busy": "2026-11-03T14:58:39.728594Z",
     "iopub.status.idle": "2026-11-03T14:58:41.354071Z",
     "shell.execute_reply": "2026-11-03T14:58:41.352547Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "one control arm (40 groups): residual ICC = 0.114   R2 of baseline + strata = 0.489\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "pooled over 30 endlines:     residual ICC = 0.122   R2 = 0.499   SD of y = 1147\n"
     ]
    }
   ],
   "source": [
    "ctrl = df[df.treat == 0].copy()\n",
    "r = smf.ols('y ~ x0 + C(stratum)', data=ctrl).fit()\n",
    "ctrl['res'] = r.resid\n",
    "import statsmodels.api as sm\n",
    "mix = sm.MixedLM.from_formula('res ~ 1', groups='group', data=ctrl).fit()\n",
    "var_u, var_e = float(mix.cov_re.iloc[0, 0]), float(mix.scale)\n",
    "icc_res = var_u / (var_u + var_e)\n",
    "r2 = 1 - r.resid.var() / ctrl.y.var()\n",
    "print(f'one control arm (40 groups): residual ICC = {icc_res:.3f}   R2 of baseline + strata = {r2:.3f}')\n",
    "\n",
    "# The same quantities averaged over 30 synthetic endlines without treatment (the \"true\" planning values\n",
    "# of this synthetic world; a single trial only sees one noisy draw, as above). ANOVA estimator of the ICC.\n",
    "pool = []\n",
    "for s_ in range(30):\n",
    "    d0 = simulate_members(cl, np.zeros(J, int), seed=5000 + s_)\n",
    "    d0['stratum'] = d0.region.astype(str) + '_' + d0.rural.astype(str)\n",
    "    f0 = smf.ols('y ~ x0 + C(stratum)', data=d0).fit(); d0['res'] = f0.resid\n",
    "    gm = d0.groupby('group').res.agg(['mean', 'count', 'var'])\n",
    "    s_u = gm['mean'].var() - gm['var'].mean() / gm['count'].mean()\n",
    "    pool.append((s_u / (s_u + gm['var'].mean()), 1 - f0.resid.var() / d0.y.var(), d0.y.std()))\n",
    "icc_true, r2_true, sd_true = np.mean(pool, axis=0)\n",
    "print(f'pooled over 30 endlines:     residual ICC = {icc_true:.3f}   R2 = {r2_true:.3f}   SD of y = {sd_true:.0f}')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b1a9e57d",
   "metadata": {},
   "source": [
    "## 6. Randomisation inference under the same rule\n",
    "\n",
    "The p-value from the sampling distribution of $\\hat\\tau$ across *re-draws of the assignment using the\n",
    "exact procedure used in the field* (same strata, same rerandomisation rule), holding outcomes fixed under\n",
    "the sharp null of no effect for any group. 1,000 re-draws."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "eabbc9dc",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-11-03T14:58:41.356191Z",
     "iopub.status.busy": "2026-11-03T14:58:41.355906Z",
     "iopub.status.idle": "2026-11-03T14:59:00.914302Z",
     "shell.execute_reply": "2026-11-03T14:59:00.913062Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "observed tau = 141.2;  randomisation-inference p = 0.082\n"
     ]
    }
   ],
   "source": [
    "def tau_ols(df, zc):\n",
    "    t = zc[df.group.values - 1]\n",
    "    X = np.column_stack([np.ones(len(df)), t, df.x0.values, pd.get_dummies(df.stratum, drop_first=True).values.astype(float)])\n",
    "    b = np.linalg.lstsq(X, df.y.values, rcond=None)[0]\n",
    "    return b[1]\n",
    "\n",
    "ri_rng = mulberry32(424242); keys = strata_keys(cl, STRATA); draws_ri = []\n",
    "obs = tau_ols(df, z)\n",
    "while len(draws_ri) < 1000:\n",
    "    zz = assign_once(keys, ri_rng)\n",
    "    if balance(cl, zz).std_diff.abs().max() < 0.25:\n",
    "        draws_ri.append(tau_ols(df, zz))\n",
    "draws_ri = np.array(draws_ri)\n",
    "print(f'observed tau = {obs:.1f};  randomisation-inference p = {np.mean(np.abs(draws_ri) >= abs(obs)):.3f}')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a3d6920b",
   "metadata": {},
   "source": [
    "## 7. Power by simulation of the full pipeline\n",
    "\n",
    "Repeat 300 times: draw a new assignment under the rule, a new endline, fit the ANCOVA with\n",
    "cluster-robust errors (t with $G-1$ degrees of freedom), record rejection at 5 %. The analytic MDE from\n",
    "section 4, using the pooled residual ICC and $R^2$ above, should give roughly 80 % power at an effect\n",
    "equal to the MDE."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "d2b2acfd",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-11-03T14:59:00.916651Z",
     "iopub.status.busy": "2026-11-03T14:59:00.916355Z",
     "iopub.status.idle": "2026-11-03T14:59:18.402929Z",
     "shell.execute_reply": "2026-11-03T14:59:18.401544Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "analytic MDE with pooled ICC and R2: 0.189 SD = 216 units\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "simulated power at the analytic MDE: 0.810  (Monte Carlo SE ~ 0.023)\n"
     ]
    }
   ],
   "source": [
    "mde_sd = mde_cluster(J, mbar, icc_true, r2_true, cv=cv)\n",
    "eff = mde_sd * sd_true\n",
    "print(f'analytic MDE with pooled ICC and R2: {mde_sd:.3f} SD = {eff:.0f} units')\n",
    "rej = []; sim_rng = mulberry32(777)\n",
    "for rep in range(300):\n",
    "    while True:\n",
    "        zz = assign_once(keys, sim_rng)\n",
    "        if balance(cl, zz).std_diff.abs().max() < 0.25:\n",
    "            break\n",
    "    d = simulate_members(cl, zz, seed=1000 + rep, tau_urban=eff, tau_rural=eff, full_dose=True)\n",
    "    d['stratum'] = d.region.astype(str) + '_' + d.rural.astype(str)\n",
    "    f = smf.ols('y ~ treat + x0 + C(stratum)', data=d).fit(cov_type='cluster', cov_kwds={'groups': d.group}, use_t=True)\n",
    "    rej.append(f.pvalues['treat'] < 0.05)\n",
    "print(f'simulated power at the analytic MDE: {np.mean(rej):.3f}  (Monte Carlo SE ~ {math.sqrt(0.8 * 0.2 / 300):.3f})')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3f502e01",
   "metadata": {},
   "source": [
    "Note: in this check every treated group receives the full dose (fidelity is not part of the power\n",
    "assumption). Realistic under-delivery dilutes the intention-to-treat effect, which is why the plan\n",
    "sets the MDE on the *delivered* programme and monitors fidelity.\n",
    "\n",
    "## 8. Heterogeneity (pre-specified, exploratory)\n",
    "\n",
    "Two pre-specified moderators: rural groups (cluster level, also a stratification variable, so the\n",
    "interaction is estimated within balanced strata) and member sex (individual level)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "b90ae436",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-11-03T14:59:18.405288Z",
     "iopub.status.busy": "2026-11-03T14:59:18.405085Z",
     "iopub.status.idle": "2026-11-03T14:59:18.436623Z",
     "shell.execute_reply": "2026-11-03T14:59:18.435034Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "treat:rural    coef =   120.3  SE =  164.0  p = 0.465\n",
      "treat:female   coef =  -186.0  SE =  108.7  p = 0.091\n"
     ]
    }
   ],
   "source": [
    "# rural is a stratification variable, so its main effect is absorbed by the stratum fixed effects\n",
    "for mod, rhs in [('rural', 'treat + treat:rural'), ('female', 'treat * female')]:\n",
    "    h = smf.ols(f'y ~ {rhs} + x0 + C(stratum)', data=df).fit(cov_type='cluster', cov_kwds={'groups': df.group}, use_t=True)\n",
    "    term = f'treat:{mod}'\n",
    "    print(f'{term:<14} coef = {h.params[term]:7.1f}  SE = {h.bse[term]:6.1f}  p = {h.pvalues[term]:.3f}')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cb827d6f",
   "metadata": {},
   "source": [
    "The true rural–peri-urban difference in this synthetic world is 100 units at full dose. With 80\n",
    "clusters an interaction of that size is estimated imprecisely: heterogeneity analyses in cluster\n",
    "trials need far more clusters than the main effect, and the plan labels them exploratory.\n",
    "\n",
    "## 9. Implementation fidelity\n",
    "\n",
    "Fidelity is measured, not assumed: sessions delivered (facilitator log, verified by spot checks),\n",
    "attendance (group register), and a five-item facilitator checklist per session. Below, the\n",
    "distribution of sessions delivered among treated groups and the descriptive dose–response. The\n",
    "dose–response is **not causal** (groups that received more sessions may differ), which is why the\n",
    "primary estimate is intention-to-treat."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "f378af88",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-11-03T14:59:18.438943Z",
     "iopub.status.busy": "2026-11-03T14:59:18.438713Z",
     "iopub.status.idle": "2026-11-03T14:59:18.587472Z",
     "shell.execute_reply": "2026-11-03T14:59:18.586035Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "sessions\n",
      "4     3\n",
      "5     3\n",
      "6    10\n",
      "7    16\n",
      "8     8\n",
      "Name: groups by sessions delivered (of 8), dtype: int64\n",
      "mean attendance share: 0.811\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAk4AAAFKCAYAAADrINB+AAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjExLjIsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvgI3uAAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAbCFJREFUeF7tnQeYE9X39w9t6b0jvSNdaQpSpEmVIiBIbwKiIIqAYEFpCggoCiJFBQWR3ntVQCkC0nsXXKr0pcz7fM/vP3kn2WQ32WyZhO/nefLsJpnM3Db3njntxjMMwxBCCCGEEBIp8SM9ghBCCCGEKBScCCGEEEK8hIITIYQQQoiXUHAihBBCCPESCk6EEEIIIV5CwYkQQgghxEsoOBFCCCGEeAkFJ0IIIYQQL6HgRAghhBDiJRScCCGEEEK8hIITIYQQQoiXUHAihBBCCPESCk4k2vn3339l9erV0Xbev/76Sw4cOBBt5yPESnSPVxL8/PPPP1KzZk3ZtGlT8Fc2lujTp48MHDgwlq7mHxScbEz27NklUaJEjle2bNmkdu3a8vnnn8uDBw9sW/J+/fpJrVq1ZM+ePT6X8fTp03Ljxg3H7wzDkOrVq0u9evV8Pld0gnKhL2bPnh2dpyU2wJ/x6srly5d1UY0Lzpw543Tv+Mv9+/d1zH/55Zf+nirowBy8ZcsWKV26dIzU7fbt2zrnPH78OEbObwe++eYbHV+oK8iVK5d89tlncuTIETsUL0IoONmYR48eSdmyZWXy5Mny3XffSc+ePSV16tTy/vvvS/ny5eXs2bO2LP27774rY8eOlWLFivlcvvz588v48eMdv4sXL57W/euvv/b5XNHNw4cPVZAjwYU/49WVN998U+rWrRsnDVSoUCG976ITjvnwhIWFybRp06RJkyaSMmXK6GxuB3PmzJHcuXNLaGhojJzfDmAuxfgyadmypWO+tzsJ7V7AJ518+fJJu3btnJph27Zt0qBBA+ndu7fMnTvXbRPdu3dPkiRJEmHz4YkyYcKEkiBBAo/HYXDjOG/OhZsgefLkUrRoUX15AudLnDix113btGlTr491B8oFASyieroDdYdmLyQkJNLre1MnPD1i0o2sLb3pO2/Akxz617WtvTn/3bt3tZxot4jw5lzR3Y537tyRZMmSRTq20W+RlR9ENl7RZ3gy9uZc3gCNKhbc+PGdn1u9GUPoU/Qn+tVXvDk/8PY4V3zpZ4wbHJ80aVKn00T12p7A/Yt2jmyOQ3lcy+IJuCGgD9u2bevxGIyZyOYNzAd4OMbYii68uR9v3rwpKVKkiJbxjHsR7RbRvXHr1i1dFyK7fzJlyiQvvfSSCo0jR46MriaJGQxiW7JkyWK0bt3abfm++uorqD6MP/74w+n78ePHGwULFjTixYtnpEuXzujTp49x8+ZNp2OWL19ulCxZ0ogfP76RMGFC48UXXzSOHj3qdMy5c+eMl19+2UidOrWe6+mnnzZ++eUXxzEPHz40UqRIYXz77bdG+/btjWTJkhmdO3fW70ePHq3fmXz99dd6nhMnThgVKlTQ86VPn954/fXXjdu3b+ths2bN0t+gTiEhIfo/PgPVqlUzXn31Vf3/u+++0+/27t0brl2KFi2qZTbZvXu3UblyZSNJkiRG4sSJjerVq7v9nSvXrl0zWrVqpWVOkCCBtg/aGWUzy2Qybtw4o0CBAlqnbNmyGW+++aZx584dp1MePnzYeOmll7QMOK5IkSLaB1bu3r1rvPXWW0bWrFn1Orly5TImTJhgPH782OvxuWXLFm2btWvXGlWrVtWyT548WX8fFhZm9OvXz8iRI4ee/6mnntJ+Qj9aQfvmyZNHj0F5W7Rooe1h0qhRI+OVV14xFi9e7DhX3rx5jW+++SZcOf/55x+jadOmRpo0aXSclSpVSn9n5ZNPPtH2Q/s+//zzehyOx+dWHj16ZHz88cdGpkyZtGwpU6bU9kK9rCxbtkyvg/OgLTBuzp8/H2EbuhuvuHf27dun7Yh2wPXeeeedcO1lcu/ePb0erov7yrw2wBjBeP/zzz+N4sWLa/lxbnD58mWjZcuWej38DvfZggULnMqLMYAyYkzgt+jX2rVrG2fOnNHj5syZE+7emTFjhuMc8+bNM0qUKKG/Qz3atm1rXLp0yekaGLO4H1FOjNHy5csbBw4c0OuNHTs2wva7evWq1sG8X3Cf4bdoA5QN4HrmPY0xgTbt37+/47yR3Ufm7631AqtXr9bz7tq1S9+b98DWrVuNJk2aGIkSJTKSJ09u1K9fP9w4OHXqlM4XmB9Qz/z58xtz586NsK6gd+/eOkYxJk3QP7j/zTkHc1iqVKmMAQMGhLuHMdei/1ButAPmxN9//91xrjp16ujnKBPmVfRJRPz77796n6ZNm9ZRj2nTpjn95NatWzo3mcfgvJ06ddJxa4J5GmWyfgYwLjB+rGBMZM+eXa+XNGlSXQNwDRPM67169XJczzwG85wJ1ir83vq777//Xj9zXY/sBgUnGxOR4ITBhgGGwWfy4Ycf6sSNBRILyKhRo/TmbdasmeMYTGi4qTF54ZiffvpJJ2ssXuYieuXKFV2gMJFg4Z05c6ZOPLje7Nmz9Vw4Fu9z5sypkw9uVExaYMSIEfqdCW4yLApYNPr27WssWbLEGDp0qE5YmEAABDV8juMgLOJ/fAYwsWDBBpikMamgjlZwbVwTZQUQrHD+mjVr6uQNgadMmTJ6I2NB9wQmw2eeeUYXgS+++EIX+m7duhmFChUKJzihDFhkMJGivJ9//rkuTJUqVXJMlmin3Llz62KOSR9tXq9ePe0D6+SAxQYTJNoKx+CcaAssKN7y22+/aRkhyPTo0UOvB6ENYBFBnVBGnB+LFsbK4MGDHafHgo3fY4KF8AXBLWPGjEbDhg0dx0AAxPnR7xh7ixYt0kkTvxszZozjuBs3buhCgrpPmjRJF2+MQxz3888/O44bNGiQLkJocwhfCxcu1L7Gcb/++qvjOIwpLIIo74YNG4zPPvtM2/C9995zHINFD/0BwQDlmjJligobEKjv37/vsRndjVcIGRh3qBPOhTGJMqFN3IFxgzHwwgsvqOCJ/7dv366H4j7E+XCfQfjDPYdxjEUkX758+kJZcR0IILjOihUrHJeBEInP2rRpo3XEGIEA/Oyzz+oxEAxxPbQPBH78f/bsWf3uxx9/1Dbp0qWLjmU86KBfcC9YhWb0K+6XIUOGGEuXLtX7FAJ+ZIITzoGHMPQh2grX7t69u1GsWDGn++XixYuO+QJ1/OGHH4ydO3dqGb25j8zfY2G1gnZCGXfs2KEfm/cA5hoIExjT6DM8kKDNrA81KGPhwoX1nKtWrdLxiX4yhTBPVKxYUe9XKxhDaL+yZcsaX375pfYlro+yWYWY06dP65wMQQRtgIdRzFG47vr16/WUEPpw/+O306dP1/vVE1gHMMYhxEycOFHnZzzY4Ld4CDIx5xzcLxgHWCtQ3jfeeMNxDMYVfucqOOH+xlpggvbCmMK5UGbUF3XCfWeCNQHjEeMI18O4h7AGIczEneC0f//+cHOEHaHgZGMiEpwAvsciBy5cuKBSvXXxAhAYcDPgaRfgqQLvrU/hWMChATGf3qGlwk1lPQYTJCYCaLqAKTjhpnfF3UKEaw4cONDpUPPpYuPGjY7PsZhj8rZiFZwAtAiYKKxPfBAUMHmbTzSYKDCJWY+5fv26kSFDBkebucNsL1eNUMeOHZ0WAjytYuF+//33nU6zcuVKPc7UzpkTARZLE0zeuI6pycHEgonIVXvYs2dPfVp21ap4wlw0PvjgA6dDNm/erGXC07kVCC2og1kOLHjQfFi1XBg3mzZtcvwMCyzOZf3M7BP81lyY0C6YONFOViCEQdtlCjIoA84HYcgE9cXiDkHBBFrH5557zulca9asMfbs2aOfoZ8h0LVr187pGLQ/BFAIDJ7wNF5NbQlAm0DwgTYgIjA2ISRbgeCEOkIgsTJy5EgVtF01YjVq1NCF3+Svv/4Kp0GAQI1zQttggnsWC5QJ2hjtaF0cAYQM/NYck2h71wUeoA8jE5wgmOEY9IUVaJ/xuamhNQUfCE1WvL2PfBWczAcyk7///lvvMWiGANoc5cN8aIJxhwev0NDQCPsYDzjQdloxxxAERxOMSQjFjRs3dnwG4QL3iVWL++DBA52rSpcu7TjOnBtR74iAUA2h68iRI47DMFZxH6BMAO8hGLoKYOgLlM/EW8EJv8PDgRXMPdu2bXN8hLnG2hYA9zOEPBN3ghPaDGsABDs7Q+fwAAZ2e9O5DhEe8EuBM/nhw4cdr8KFC+txCOkHpsN2165dZfPmzfp7OGTDtmza2tetWyf169fXKD4TnGPWrFnqoG7lxRdf9LoFXX0CXnvtNfUD+P33370+B+jYsaOcO3fOEQoMP4ZffvlFzwe/HPhOrF+/XqpVqyZHjx51tMXFixe1fXbv3u3xer/99ptkyZJF28NKmzZtnN7jOPgxoB2tIJowT548snbtWv04R44c6tA/dOhQtd3D3g+fgBYtWkiaNGn0GBybM2dOPc7adygrfFpQB9S5Tp06Hl8zZ850FMO1T3D+tGnTapSm9fzPPvus1sFM9VC8eHG5evWqvP766xoRiXZEGV544QWnOqJ+rp916NBBf3vo0CE9FmMIEaCIlLGCc58/f17LYYJxV6lSJaf3KAuiikzw/o8//pCPPvpIjh8/rh8j2rJEiRL6/8mTJ+XEiRNaLmsdMW7z5s0bYZ97wtqO8M9ABJW1TL7irl/QB/A5sZa5YsWK2ifoG1CqVClp3769+sScOnVKNmzY4KgP6uyJgwcPyoULF/R81vPDvwX3tnkfYCyjD3D/RDTm3YF7F2McfeHNb13bwNv7yFdc5xrMe2XKlHHMNRkzZpTMmTPLF198IT///LP6LKENXn31VcmQIYPHy+GeuHbtmv7eHdb6wbcKfWcdM7gvmjdv7rj3AfzVOnfurHM0zu0LOB/muQIFCjh+hrH6/fffa7QowHukTsD8gTGFe3T58uU6h5r3ki/gXsT99vbbbzvmDowxBCyZ1KhRQyOhMTdj7K1YsUKDmdAW8OvyBNosXbp0Gp1qZ3z3MCS2AE6AmBQRSQPMCdS6AFnB4gsqV64sEydOlE8//VT/h9MehISPP/7YIVThXK4ToSd8cWx86qmnnE6DCSNr1qx6E/oCygZBY8aMGVK1alW9Ka9cuaLCBYCABKdFhAzj5QoEI0/gxoaA4YrrZ2gj3ORYNFzBQm1OSHACXrx4sbz11lvSrFkzXcgxySBnycsvv+zoO1wXQq470HeYTCJKQWF1lnbtE5QVE3KRIkU8nv/5559XIRCTKcKEEcmYPn16ady4sXzyySfaTyau/QjQHwB9CQED14RQ4K5tANoHEzCAs7Or8641TBkMGTJEhUg4jaI86I9WrVrJBx98oIKAOf6xALkDEUq+4uqAjjL5Ex7url8Qeu2p39GWuL8xNhD5t2zZMh3XGFMYY5FhtgnaydO4Ajg/hAjX8rm7D1xBCgR348H60GXFXRt4cx/5iqcxarYJyrF06VLp0aOHCowoQ7ly5aRXr14qPHkC9yDGgCfHb1cHc+uYgdCChwbXhwlg9ifuCwh43oL6QCiKDAiogwYNUsERD8twxDYf3HwFcxdyn02ZMkWjOHEuzG1YQ0yhEw/yuJ75cA5BE3OYN2A+wPpmZyg4BSjQXuDpp2TJkloDDF6Ap0xoF1yxLgJ46sciCa3CmjVrZNy4cbqYYxLHBIoXtAfRDZ4iIKhZwWcRCTLuwCSHJ3CUG4s8nhiheXjmmWf0cCz4WIgHDx7sdiGNKLoDdf/777/DXdb1CQjHYQK9fv16uAkBx1pTMUALgqdJU1uAMjdq1EgXQjwFou8KFizoMZkeJjhTyIoIT5o7nB/l9ZRXK1WqVHpatBk0Y9DqbN++XcuHNBBI8gcNiCncuHsaND8z+9LTGHI9zlsgLCAtB/p869atMn/+fJ209+/fL0uWLHGMf9wX7h4eIotwigvQL9D2Tp061e3lMY4B0htgIcEDD57kIcRC6wRNQ0SYcwK0C+7yDZnRlugrPHi44s1TP4Qrd8lp8fDiDb7cR8BVcDW1ct6U3XWugWAPLSYeFjZu3CjffvuthsQDT8ITxhHmMHftFRn4Le5ld3NrVO8Lb+ZqCDkQrqARwn2D+Qja7REjRsiAAQPCFTuyNsa4GTNmjP4e7YcHwwkTJmi0944dO7RtME7RvvPmzdMHdFxv1KhR0rdv38iaSX8fkdbPDtBUF4DgKQMZVjFxQuMCzIkR6ndT+MFfCEx4cjUFJ6hLYbqB8AA1Mp5kJ02aJP/995/s3LlTzwUBBFoc10kJC9WCBQui3GJYiK3AnAYtgmluMfHmqR6CE8oMwQ8LQ6dOnRy/xwT19NNP6wJrbQv8j7aLKBUCtCAQcPbt2+dUVteym0LaokWLnI7D+SF4mdoWmACwsOMpD1oPlBvtD20b6m/2HZ40YWq1lhfaBQhmEZXXGzA2Ll26pDlhrOfH0zMmbPMpGUI3JkK0HwRpCFEYZxCo8aRsgvfHjh1zujQmTzxdm9oTtA/CtlEnKxg/MKeif7wFDwgYmxCSsGhhUoZA17p1azVVAAieGOMQ+Kx1hPAAU4G3oebRgTfj1+x31AsaM2uZ0Se49zBGoCmEYIIFB+ZiU/NnmkRdsV4bbYy+xGJmPT8WJbSJOSdgzKOfTPOyieuYd0eFChVUY+X62x9++MHt8a54ex+hzBhfruPOnLNccS07xj7GtjnXwGSOY1BvCH/QOmEugYBujilPwCwW1Z0MUF/cK6754HBfYKy6avkiG0s4H8qL+liBkA3XCoA2guCN+RsuGBBigOsYMseWtY1RTtPNw2Tv3r2ya9cunZcgFEELDE0droM5Gd+hXSFcIW2Op+u5A3Mv5j2r6dGW2NkB60kHzt+1atVSxzs4+CIaCVEKcBiGgy0cX63ACRHRcHA2hpM0HAYRDQenSNOJFuGxeA/HUDjyHjp0yHjttdfUkdcMb0aoNBwOEWWCKBA4oCLCDM6UH330kV7SdA6H46srnpxt4VAIB0VEbcCREM68SJ0A50gTRGsh/B9RQaajsatzuAkiWxCdAgdnhHVbQfQRygsHVzh7IorJTOFgjcRyBakbEEkG51y0GcoAR28z9NaajgBtC0dP1AntjePh4Im+MUNsTWdxONyjHHD2hFOwNQIQ10C9ESWFSCM49cJRHI6bGAOu6Q08YTqH468VtC8cm/FCxA7Oj8ghtB0i7dA2oEGDBuqsDGdhHIOyIEoMY8p0sodzOI4pV66cRmgi7BgOwoh0hIO+CcYexhSia06ePGn8999/GuWDvoJDuAn+t6auMEHbIjrRBOME4wVtjPGDgAJE+iDCyQQRdzg/HMFxPTgAm9FJKKMnPI1Xd9FFiMKKCPQz2gJh6aYDsOkc7pr6AWMckUYIZECABsYQxguCHNC+JhiPcHY+fvy4ngOOt3D6xhiyOuSifapUqaLnNdN8YL5A0AicjTEmcY8jsgm/RaQjwPhAVC0iIDE+UA6MaYw/zBURpSPAb3GvYEyg/dHOcAKGg7w753DXqDjgzX0EkBoC4xWRcoi4xRyGcruLqsN9iWg63M9oW/wW7WA64pvpRRBtiHsS9yaOR309RU6aIMoWDuLW/nQdQyauwQJmoAac5xHQg3sPv8V1ERVnAmdu03ndGgDgCu4tBAUgaAHjAxHROA8CIpAOAKD+OP/bb7+t50L/4l40Ux6YYK7HZ7jvcJ/hfLinEZFojarDnIAgG4w9zBP4i7nLdDRHmXA9ONCb15s6darjema7uXMORyQw6u26ttkNCk42BosmBpH5wmDE5IjoJ0SIuYKb0Az5xrH4HSZhhL2aYCIxI6Nwc+EYhOi75kdBOC0mB/PaEKQw4ZqDPiqCE6KKsEiaZUOIv+sNgokVix+uCwEDeBKcMHHiuObNm7vtRUyAmCxxDK6JF8LKrblE3IGJwBSU8Bu0jykAWQUnTFIQNszz4y9C3xHBY8UM/TfrjagRLLBWDh48qEKX9VzoazOk3Rs8CU4AObSQJ8l6ftTRGtGIRRXCFI4xxwaiYDDZm2DsYBFC9JbZpjge6Q4grFhBqLiZd8k8FoKMNTWAt4ITIsHMPjHLhrJaUzqYuarM71EujCXXaE5XolNwwmJj5spCLiTgSXACaH8IzdZ+QYQVFjET3IsQTHAMxg4WLYxtV8EJoevmYjh8+HD9OYQ/LKDWvsI9gXQgViAEmyk3cCweziCg4Xre5HFCTiD0DwQbhMNDcPRWcPL2PkIZIUCbcxLSCSAKzp3ghPkM/WDecxCkXCNlhw0b5phrcBzGTdeuXZ0icd1hzgVI2WHireAEUDazP3FdM5LYGgGM+RXCsjku0EaeQDnM/G/m8ZjnrII/IpUxh+N6+IsUEmZEqxUIVOacieM6dOigApBVcMJDkdkP5r0GAdZMLwHQHqiXeT30FR66Ub6IBCekIkFqErsTDwW0pSqMqNrSGoEAVWpkGZMBzCpQp+JYOBq6+hWZalqYXOCwB5MdzAWuQL0LvxjYnKHidlUjw7wEPwxXJ0OYl/AbZD0H8EtBlnOcDxFEpjkF53SXARmmCpi8EL0Fp0L4IOA4V18omL9gJnBXBhOYAtEWMHnA78h0Yo4MZDBG3fE7+AbADwgmT5g8XNsK/QQzEkxxZhSjK9i/DGZU1AOmEXd+XVDLQ00OMwqug75zzS4dEWhf9D0cY91lJ8etjnaFmQGmD0TMubY/joGqHW2O8YaxYfUPgk8WrgMzI8YPXjAJeXJWRjua7Q+ToWs0Evw9MF5cf4/2QntYnXxhWsTYgdkFDrbw73Pnr4Zzoq0xPtxd0xXX8QrzKsYgymQ9P3xFUPfIxhCOQRtjnMCE6Hp+V1BOmEkxvjCGXE3XAPcTTE0Y6xgX+A36CPek1ZSL43BfmveOCepjRtHBvOPOURfnhJkZbY+xgb7CuTBWTHOLt8DUjXEOUzbMNZjHcJ/AHOVpmxJv7iOUEfcH2hZjAOPL2g7w84OPG5yh0fcw16N9UGd3cyf6FKYojHuU151TuSuoC64NkzGi14CnPoaJHGZX1yASzEvoD5Qf5XTnl2rO0zBdoWwRBeLgGmbfYZ5zFwyB9sW9iDGJecoc5+Y8bQJzGwJScBz6CmMHZbCOe9ybMAEjQAljEHVwDfDA3Iw6wuyGOcL1eq73Ga6Dc8GXEWZ4W2N3yY4EPp7yg5DAw9Q4EQJg4oKp35r5GkDLB+2CmYgztohI6xqdwMwE7V5keZaI9yAnFawsESWrtQuMqiOEEBIloNmEtgyBKsjbgzQE0HwgEhOBJ96kNAhEunTpotpaaJqgHSb+A60U8tHZMQLWFQpOJMaBmhyTqK+b7BL7gQ2XrTuakycbmMyQs2f69OmaIgHmIOQLQw4fb/ILRTcwtWGu8cbk5g9Y3JHWhUQf/m7mHpvQx4kQQgghxEuYx4kQQgghxEsoOBFCCCGEeAl9nLwAoZcIu0RoZkTbdRBCCCEk8EBKCqTLQYBDZGlgKDh5AYQmd5u5EkIIISR4wLZkkUWDUnDyAjNhGxrU3BCVEEIIIcEBEn9CQeIpQasVCk5eYJrnIDRRcCKEEEKCE2/ccegcTgghhBDiJRScCCGEEEK8hIITIYQQQoiXUHAihBBCCPESWzmHw6v9vffek9OnT8t3332nIYGvvPKK3L59O1x16tatK2+++aZcuXJFWrduHe77zz//XIoXL66fIz/DmDFjZOnSpZI0aVJp0aKFtGnTxssmIoQQYmfuhYVJ6I0bkjF1akkSEmLnopIgwDaC019//SXNmjWT9OnTy59//ukQlmrXri1hYWGOpoYQ9NFHH0np0qX1sxMnTsjq1atl7NixTt7wOI9Jjx49ZPny5fLBBx/oeXv27CmhoaHSp0+fIOhCQkhMwMXY/jx6/FjGzF8oU1aulrthYZI0JEQ61a4pbzd+WRLEp0GFxAy2EZwmTZoknTp1Uu1Rzpw5HbXt0qWLU80hJEH4gTAEzpw5o5k+IQy5w9Rebd++3SFsZcqUSd544w3p3r27aqAIIcSEi3HgAKFp/OKljgJDeDLfv9u0ceBUhAQUthHJR4wYIQMGDIg01Tm0TRCmzMyeEJyQtOrWrVuyb98+uXv3rlMHLFmyRJ5++mmH0AQaN24s9+/fl82bN7vtLHwHs6H1RQh5MjAXYyzCwFyM8TmxD9AIQtPkjqmrVuv3hMQEthGcUqdOHWn9Vq5cKbt27ZL+/fs7joVGCeY6aJ3Kly8vadKkkZEjRzq+R7bv/PnzO507ceLEKmzhO3cMHz5cy2O+uN0KIU8GXIwDB/g0mcKtK3fu/8/niZCYwDaCkzd8/PHH8vrrr8tTTz3lOLxly5aqrTp37pxqnX766SfVXK1atUqPuX79utsU6kirfu3aNbeXxe9v3LjheHkSsAghwQUX48ABjuDwaXJHssQh+j0hMUHACE7Lli2T3bt3O2mbALRM7dq1061QkCodUXiNGjWSefPm6fdZsmSRy5cvh2s7OIdnzZrVbZtCI2Vur8JtVgh5cuBiHDggeg6O4O7oWKsmo+tIjBEwgpOpbXIVdmCqu3r1qlMDZciQQf2TAMx0MO89evTIccy///7r1oRHCHmy4WIcWCB6rmeDeqphAviL9/ickJjCNlF1EQEHbzh+L1q0KNxh77//vhw/fly2bt2qGidol3Dchx9+qMfCERwRd9OmTZPOnTvrZ8jxVKhQIdVWEUKIFXPRhYMxfGWwGEODwcXYfiDlAKLnICwxjxOJLeIZSIxkI86fP68Rc4cOHVLhBpQpU0YqV64sX3zxRbiSHj16VCpVqiQhISFSoUIF2bRpk5QqVUqTXSZM+D+5cObMmZrqoGrVqprKYM+ePbJ48WJ54YUXvKo5tFdwFIffE0x3hJDgh3mcCHly8GWdt53gdOfOHRV04KtkCitz5syRl19+WU1w7oBA8+uvv8qlS5ekZMmSUq9ePadkmGD//v2yZs0azduEpJq5cuXyekT40qCEEEIICSx8WedtJzjZEV8alBBCCCGBhS/rfMA4hxNCCCGExDUUnAghhBBCvISCEyGEEEKIl1BwIoQQQgjxEgpOhBBCCCFeQsGJEEIIIcRLokVwCgsL021NkH/p/v37mlfJZnk1CSGEEEL8xm/BaefOnVK0aFEpV66cvPbaa3Lr1i0ZOXKk1KhRQwUqQgghhJBgwS/BCRvntmrVSurWraub5iZPnlzbpWvXrnLixAmZPXt2sLQTIYSQGARb3JwNDdW/hNgZvzb5hXCEbU7Gjh3rtMVJzpw5pUmTJrJjxw5p3bq1netPCCEkDnn0+LGMmb9QpqxcLXfDwiRpSIh0qv2/TZWxiS8hdsOvUQkN04MHD+Thw4fh6nXt2jXJkiWL3epLCCHERkBoGr94qQpNAH/xHp8TYkf8EpyyZcsmxYsXl+7du6tvE4D5DhvuwkxXv359O9aZEEKIDYBZDpomd0xdtZpmO2JL/DLVAQhIMMelTZtWhaZcuXKpJuq7776TYsWK2bLShBBC4p7QGzccmiZX7twP0+9zZMwY9wUlxILfghP8mTZt2iSHDh2SY8eOSfr06aVEiRIOR3FCCCHEHRlTp1afJtNMZyVZ4hD9nhC74bfgZFK4cGF9EUIIId6Q5P8cweHT5ErHWjX1e0Lshl8+TmDv3r3SqFEjKVCggOTIkcPp9e2339qtvoQQQmwEoud6NqinGiaAv3iPzwmxI/EMP1J8I0t43rx55amnnpJ27dqFi6KD43jBggXtWG+f+O+//yR16tSaET1VqlQ+/ZYQQkjkwFEcPk0wz1HTRGIbX9Z5v0x18GlCHqft27drhB0hhBASFSAs0RGcBAJ+meqgbUqTJo1kzpw5EOpKCCGEEOIXfglOSZMmlXHjxsl7772nW67cvn3b6YXkmIQQQgghwYJfPk6nT5+W3Llze2yLUaNGyTvvvONVW125ckWPPX/+vMycOVMyZMigv0OOKJgDrbRv3143FDaZM2eO/PDDD3L16lXdXHjQoEGSKFEix/ePHz/WsixdulSFvRYtWkiHDh28KhfwxfZJCCGEkMDCl3XeLx8n+DXt27fPY+t4u+XK77//Lq+++qpG5G3dulWdzgH81iEU9e7dW/LkyeO4TpEiRRz/Q2Dq0aOHfPrpp3q9YcOGye7du2Xhwv+frv/111+XtWvXykcffaQasbffflsuX74sffv2DayeJYQQQkic4pfGKbqABqlMmTKqLYJQdO7cOY3Ug6YJwhA0W0i06Qo0Sfny5VMBCMITgMkQQta2bdv0nKdOndJjdu3aJSVLltRjkO28W7duep1kyZJFWg1fJFFCCCGEBBa+rPNR0jjVrVtXmjZtKg0bNpSuXbt6bB0IRC+/HHkujq+++kpSpkwphw8fdjoXBKYECRJoNnLki8qaNatktKTfx2dnzpyRVq1aOX4HrVXlypVV4wTBacmSJVK0aFGH0ASQd6pjx46yefNmqV27drjyQ+Nlar0AGpQQQgghJEqC09NPP61mOvgRRZQtPF26dF61MIQmd0AoguBUqFAhNbFdv35dWrZsKVOmTFFfJWiXIFQhss8KNEzQJgEckz9/fqfvQxD2miOH4xhXhg8fLoMHD/aq7IQQQgh5coiS4ARHaxMIGTFF2bJlNft4zZo11XT3999/S506dWTo0KEyZMgQTUiZIkWKcJeHIHby5En9PKJjIIi5Y8CAAdKnTx/HV9A4QdAihBBCyJNNlASnP/74w6O2xoq/mcNz5cql5j4TnA/+TJMmTVLBCf5PiMZzJTQ0VM16AMcgQWdEx7iSOHFifRFCCCGEWImS4PTjjz/KokWLIm1JpAXwZ8uVixcvSlhYmJNjOExzt27d0mvDBAdt0KFDh5xMhhCUTN8rHDNx4kR59OiRmv0AnM7dmfAIIYQQQiIiSgkwv/76axU8InshDYA/LFiwQLVMZh4nJNScNm2aOn8D5JCqVq2amu5M5s6dq2Y6M89T48aN1dF78uTJjmNGjBihgla5cuX8KR4hhBBCnjCipHFyBQLNwYMHNSoOkXYQVOBnFC9ePL+aE2Y6aLegNYKAhGuA6dOnO86LzOW1atWSZ599Vs1y69atk7Fjxzq2gUmePLn6SeFcEKrgZH7gwAGNtiOEEEII8QW/8zgh2SSSVx49elQTVsK/CIIL/KCQqduawTsyYHabP3++NGvWzJFfCSa2xYsXq9CE/ExIb4CIOivQSK1atUquXbsmVatWlRIlSoS7FMx5SIKJ30LQyp49e2TFceBLfgdCCCGEBBa+rPN+CU5IQIkcSS+++KJGoiF5JRJO3rx5UwUYmNCsW6MEKr40KCGEEEICC1/W+Sj5OJkcP35c/vnnHxk/fryTBge+R0iQ6S6ajRBCCCEkUPFLcMJ2JQ8fPlRzmivIkZQ5c+ZAbRdCCCEByL2wMDkbGqp/CYkJ/HIOR1JKZBF/4403ZMyYMVo+mO/gp4T94LZs2RITZSaEEEKcePT4sYyZv1CmrFwtd8PCJGlIiHSqXVPebvyyJIjvl46AECf8EpwABCT4MWHbE2ifYKZD8sgJEyZoKgFCCCEkpoHQNH7xUsdlIDyZ799t2jimL0+eIPxyDreyb98+OXbsmCaoLFWqlMf95wIRX5zGCCGExC4wy5Xu2VuFJVeSJQ6RXV+NlSQhIbFbKBJQ+LLOR4v+Eua5YsWKSaNGjdR0d+LECbd+T4QQQkh0E3rjhluhCdy5H6bfExJd+C04IalkkyZNtDx79+7VXEvQOCG7NxJhEkIIITFJxtSp1afJHdA44XtCogu/BCfka2rbtq106tRJyzNy5EjN6YRklZcvX5aff/45uspJCCGEuAVmODiCu6NjrZo005FoxS/ncJjk0qVLJw0aNNBCrVmzRrdIwT5wMNvt2bMnWgtLCCGEuAPRc2DqqtVqnoOmCUKT+Tkh0YVfglPWrFnl4sWLcubMGXUMR+6m5557Tst2+vRpKViwYHSVkxBCCPEIUg4geq5ng3rq0wTzHB3CSUzgl+CUKVMmTUUAh3Bs9NulSxdJkSKF5nRCLqcdO3bERJkJIYQQt0BYypExI1uHxBh+CU5gypQp6hyO7OHYgBdg/zqY7ZjHiRBCCCHBRLTlcQpmfMnvQAghhJDAwpd13m+NE0x08+bNU0dw/G+lfv36UqVKlcBqPUIIIYQQD/gtODVv3lzWrl2rjuLXrl2TChUqyM6dOyV16tQUmgghhBASVPiVxwnRdIsWLZL9+/fLpk2bxDAMfX/kyBFJnjy5JsMkhBBCCAkW/BKcbt++rZv7Zs+eXTJnzixhYWFy5coVFZpeeuklNeERQgghhAQLfglOhQoVUr+mVatWaXu88MILsnz5cv0feZ3gbEUIIYQQEiz4JTjFjx9f+vbtK++++662B9IRdOzYUcqVKydz5syRFi1aBEs7EUIIIYRETzoCZAnPlSuXNudPP/2ke9U1bNhQBahgwJcwRUIIIYQEFr6s87bJ4/T48WMZMWKErF+/Xnr16qWpDAAqM2TIENmwYYNmJW/ZsqV07txZ4sWLp99jM2F85sro0aOlRIkS+jHOPWrUKFm6dKkkTZpUNWEdOnTwuld9aVBCCCGEBBa+rPN+pyOAj9PEiRN1c19s+psxY0apVKmSfPDBBw4tVGT8+++/KvycO3dOBaHz58/rTxClB80VspK///77cuDAARk0aJBcvXpV+vXrp8ecPHlSha3x48c7hCmAcpi8/vrrmjLho48+Uof2t99+W68DMyMhhBBCiLf4LTj16NFDli1bplqiIkWKaFQdhChs9vvXX39ptF1krF69WvNALVy4UEqXLu04/NChQ7Jx40bdSBjnadSokR43bNgwh+CElAjZsmWTbt26ub3MqVOnZOrUqbJr1y4pWbKkHpMhQwY9/o033pBkyZJFVjxCCCGEEMUv53AkvMReddiX7r333pMGDRpI+/bt9X3+/PllxowZXjUztEo4FqY4K7lz55bjx487CV83btxQDZQJ/Kty5syparbdu3fLrVu3nM6xZMkS3TvPFJoABDCkTti8ebPb8t2/f1/PZ30RQgghhPglOF24cEHzOEHTZAXRdhUrVnSY3CIjZcqUbg+BP1LevHkd3928eVPNglb/JGickPoAuaQqV64s6dKlk+HDhzt+c/bsWRXirIRg9+wcOdQ06A78HpnPzReOJYQQQgjxS3AqWLCg+hIha7iVO3fuyIoVK1QDFV3AwbtVq1YqqPXv399x2tatW6sjOIQ4aIhmz56t/lW4PoCGylWTBSCsXb9+3W3xBgwYoL8zXxC+CCGEEEKi5OOEbVWQcgDUq1dP6tSpI02bNtWEmHDcnjt3rhQrVkzSpk0bbS0MXypcc8uWLZIoUSLHecuUKaMvE5jhGjduLPPnz9fs5VmyZJHt27eHK0doaKj6S7kjceLE+iKEEEIIsRIlwWnPnj2ybt06x3mQrwlaGVMzgz3qEL2GPetKlSrld4sjkm7x4sXy22+/SaZMmZzOh0i+VKlSqcO3Sfr06dWsB2Cmg3kPflEJEiTQzy5duuTWhEcIIYQQEhFREpxgCsMrNhg3bpx8++23Gl3nbtPgDz/8UAW0bdu2qW8VUhsgOu/jjz/W4kH79Oabb8rkyZM1LQFAvqjChQsHTYJOQgghhMQOURKcYoutW7dqziVEzSF1gJVZs2ZptN3gwYM1bxScw8uXL69aKWjAkCQTYMNhCF6I9oMJEZow5INCtB0hhBBCiC/YJnO4ybx589Q/Co7nMP1ZTYJWXnnlFRWKAMxyEIpggkPaAfg2uYKcUEiCiUi9WrVqqaDlLb5kFCWEEBI33AsLk9AbNyRj6tSSJCQkbgpBAhJf1nnbCU52xJcGJYQQErs8evxYxsxfKFNWrpa7YWGSNCREOtWuKW83flkSxPcreJw8IfiyztvaVEcIIYREBoSm8YuXOg6D8GS+f7dp48h+TohPUBQnhBASsMA8B02TO6auWq3fExKd+CU4TZs2TfeBs4LNd+HUTQghhMQ08GmChskdd+7/z+eJkOjEL8EJ26Fgc1/sVwcQ9v/LL7/opruEEEJITANHcPg0uSNZ4hD9npDoxC/BqUqVKhrW37t3b6levbomqUS+pVy5ckVnGQkhhBC3IHoOjuDu6FirJqPrSLTjl+AEIDzVrl1b0wYgwSSFJkIIIbEJoud6NqinGiaAv3iPzwmJbvyOqkNiyh07dshPP/0k3bp10y1ROnXqFN3lJIQQQtyClAOInoOwxDxOsceTmjfLL8Hpo48+kqVLlzrMc0gqWb9+fcmYMaM0bNgw9nqPEELIEw/MdjkyZnzi2yGmedLzZvmVAPPUqVMSL148J/MctE9FihRxZPUOBnxJjEUIIYQEM6PmznfKm2UCjV+g5s3yZZ33y8cpd+7c4XyaypQpE1RCEyGEEEL+B/NmiX+CEyGEEEKeHJg3i4ITIYQQQryEebMoOBFCCCHES6KaNwsmvrOhoUGxBY5fUXWEEEIIebIw82NhL0Bsa4O8WUg26i5vVjBG4PkVVXf+/Hl56qmndMQ8ePBAdu/eLVevXpVnnnlGUxIEC7542xNCCCFPAt7kcYosAs+bc8QGvqzzUdI4hYWFSfv27WXz5s1y9uxZOXjwoDRp0kSOHDkiiRMnFsMwZNSoUZockxBC7I5dJm9CAonI8mZFFIE3ZdVqefjokfywZp1u0oz9BrF1DrRWSGhqZ6IkOH399dcqNGFDX/Daa69J0aJFZevWrZIqVSr5+eefpWPHjvLcc8+p9okQQuyIt2YEClaE+E5EEXh374fJxGUrHCeF8DR+8VJ9b/dcUFESnFavXq370j3//PMSGhoqf/31l6xcuVLSpEmjlW7durXMnz9fli9fTsGJEGJbxsxf6JTIz3Xy9lawIoSEx4zAw33lLfCbQiJNaLPsSpT0YenSpVNzHEidOrUkS5ZM4ruo1pIkSSL37t2za70JIU843iTyMwUrc+I3BSt8TgiJmIgi8DwBZ3NoquxMlASnFi1ayJgxY3SPupCQEN3ct3fv3nL//n2t6/r161XjVKdOHTvXnRDyBBNZIr9zly/LlJWr3bYQnooheBFCIgY+S9AgIfIO4G+3ui9JkkSJ3P4Q30NTZWeiZKpr0KCBjBgxQh3Cs2bNKgUKFJAVK1bIkiVLJH369PLPP//IkCFD1JTnLVu2bJFPPvlEChYsKF9++aXjZzt37pQvvvhCDh8+rNu5DB48WDJnzux02jlz5sgPP/ygEX01atSQQYMGSSJLpzx+/Fid1bEhcdKkSVXw69Chg7dFI4QEIRGZEcxJ3pOJwXwq5oayhEQMHL3hswThCQ8ruO+giUqYIIHb/e6Q1sDOZjoQJY0T6Nq1q5w7d04++ugjKV++vHTu3Fl69uypQgui6/r06eP1ePrss8+kZs2a6i+1fft2x+/27NkjVatWVWFswIABmv4AwtgNixoPAlO7du2kWrVqGsU3d+5ceeWVV5yuDX+siRMnqsN6w4YN5e2335aRI0d6XT5CSPARWSK/7BkyaKSPOwLhqZgQO2FG4JlCkTtNFIQrd7mg7IZfeZyii8KFC8uMGTPUwRxaK0TngVatWulfROkB5IoqXry4wzQITVK+fPmkb9++0qNHDz0G6RHy5Mkj27ZtUw3VqVOn9Jhdu3ZJyZIl9ZjZs2frOSD4wT8rMnzJ70AICRxM5293ifzwpBxZDhpCiH/YJWLVl3U+yhqniIDmCBqkTZs2edWi0DJByHFl0aJFGqFnAvNb8+bNZeHC/zlm7t27V86cOeMQsECOHDmkcuXKjmMgiCFVgik0gUaNGmkuKqRUcAd8tdCI1hchJPgwzQi7vhorm0cO178wK5h5ZAL5qZiQQMBVExUIRMnHKTJu3bqlSTEhQHlDypQpwx127do1uX37tuTPn9/pO2iPZs6cqZ9BuwQznpkGwQTHQJtkHuN6Dji0Q8Ayj3Fl+PDh6ktFCHky8JTIz5N/BiHkySVGNE4wlUEoadq0aZRb1vRjSpEihdM5IGRdv35dP8Mxrt8DX49xBf5UOLf5goBGCHlyCcSnYkJIzOCXxunhw4fqSwTB4sqVK5IlSxYVmp599lm/S4tzAZw3W7ZsjvNBi4VIPoBj8L0rrsdYHc5NrMe4gm1j8IpN7GLnJYQQQohnoiw4IZoNEXSXLl1SQQaJMCHEXLx4UYoUKSJfffWVRrpFFSTQxAbCf/75pzqEm0AIMk1v+As/pEOHDqmDufUYRP0BHIOIukePHkmCBAn0M5QZwp6rCS8uYGZiQgghJHCIkqlu2bJl0r17dxk4cKAKS6dPn1ZHbaQLwKtevXoa9o+0BP6A9AFjx45Vnylw9OhR+fXXX/VzkDt3bhXOhg4d6rgM0hGcPHlS988DjRs3VmfvyZMnO45BDioIWuXKlfOneNECMxMTQp5koG0/GxrKhKIkYIiSxgnpAd577z0N6XcF5jPkZdq3b58sXrxY3nnnnSg3BnJBYV88aJxKly4tGzZsUEGobt26jnOOGzdOatWqpeZBXHvdunUqbJlJMpMnTy7ffvuttG/fXoUqOJwfOHBAo+3imsi2fLD7fj2EEBJVqG0ngUqUBCdogFwj2VzB9zdv3vSpXZCRHAKSCc4BQQj5naBtQoLNF1980emcEKp2794tq1at0kg8RMSVKFHC6RikMMBna9eu1czhELSyZ8/uU9ligsi2fMD3zExMCAlGIttg2RvoG0rigiglwITPEMxjyKf0zDPPhCs39qlD/iVoi3zZdsWu+JIYyxdw05fu2dvttg7ITIycMnZPPU8IIb7i79xHbRWJbnxZ56OkcYKJbv/+/eojBI0P/IVSpUql/k5//fWXXLhwQbc0CQahKSYxt3wYv3hpuMsEwn49hBASFfzVtkeHtoqQqBIl53CAqLljx45J27ZtNaru3r17kjdvXt3+BDmcYFYjkcPMxIQEBnRijj7MDZbdEdk+gJH5huJ7QmKSKGmcTBDVhg1zSdRhZmJC7A3NQtGPP9p2f7VVhPiLX4ITiT48bflACIlbaBaKGaBtB+42WI4IU1t1141mKTJtFSHRQZScw580fHEaI4QED/46MZPIiUpk3Ki58936hiKFCzZpJsRXfFnno+zjRAghwY43ZiHiH1HZB9Cdb2i3ui/JK5Wep48TiXFoqiOEEA/QLGRPrL6hF69dl1kbN8kPa9bJxGUrJGlIiHSqXVOFKxxHSHTj96jC3nRff/21lguqrpdfflkKFCig6QgIISSQMZ2Y3cGUIXEP+mfOb7+rwGT6POEvUrzAN42QmMAvwQkb51aqVEmzegNstYLtTNq1ayeffvqpbpFCCCGBDFOG2BemJiBxgV+mOmzie/XqVRk1apSWHZnEBw0apILTqVOndBuUqlWrxkW9CPGZqDipkuCHKUPijsjuSW980LhtFYlu/BKcEiVKJA8fPhTDMOTMmTNy8OBBqV69upYxRYoUmhSTELvDPD3EG5gyJPbw9p6kDxqJC/wy1SFTeM6cOaV+/frSrFkzqVixom6eu379evnpp5/UjEeI3THz9NBHghB74O09SR80Ehf4JTjFjx9fFi9eLJkzZ5b8+fPLjz/+qHX4888/pX379tK4MfNpEHtDH4nghNujBC6+3pP0QSOxjV+mOpAnTx6HwGTSr1+/2K4HIVGCPhLBhbcmHmJffL0n6YNGYhu/BKewsDD5+++/PZYZZjtoowixK/SRCC64PUrgE9V7kj5oJLbwy1T3zz//SJkyZdy+4O80c+bM2KoHIVGCPhLBg68mHmJPeE8Su+OXxilXrlwaVWcFkXQrVqzQV/fu3e1ef0KivNkosRe+mniIfYnqBsCExAYxtslv27ZtpV69etKiRYvYqEeM4svmfyRwiSxnDLE33JA3+OA9SWILX9Z5v0x1EXHs2DE5fPhwbNWZEL+JymajxD7QxBN88J4kdsQvU93t27dlwYIFTvUyTXVbt27VLVgIISS2oIknMIhtTVJsX48EN36Z6s6fPx9uS5V48eJJjhw5pGvXrtFipmvQoIHcunUrXC80bNhQ3n77bbl8+bIm33Rl7NixUrJkSf348ePHMmLECFm6dKkkTZpUy9WlSxeve9YXFR4hJO7hQmlPYjtdRGxfjwQuvqzzfmmcnnrqKccGvzFF8+bNNe2BCbZ3QZ4oRO2BEydOyObNm2XixIkqtJlY0yBASNq4caMMHjxYtWT4/ZUrV6R///4xVWxCSBzC0HR7EtvpImL7euTJwC+NU1ywfPlyadKkiQpMWbNmlTlz5kifPn10rzx3nDx5UrOa7969W4oXL66H4DcQpqAxS5YsWaTV8EUSJYQQEp6Ydt531TLG9PVIcOHLOu+3czj2pUMupyRJkqjGx/oaPXp0tLfsxx9/LN26dVOhCUBgwn55165dkx07dsiNGzecrrlkyRIVmEyhCcDMhzQK0FS54/79+9qI1hchhJCo4026iKgAc9youfOldM/e8kLfAfoX7y9eu+bYf9KVO/fDtDyERAW/THV37tyRpk2bSoUKFeSXX36RLFmyOJUBAk10Ah8lZCpfuHCh47QQnI4ePao5pRImTCg3b96UDz/8UD744AM95ty5c5IvXz6nYoSEhGhW87Nnz7ot3vDhw9WsRwghJHqIakbwyPBkjnv46JEkDQlxKzwlSxyi5SEkKvglOMEMBuFp+vTpkj59+qhc3yegbUJSTauAhnxRzz33nGqR4PgNDRM2F3722Welbt26qoFKnjx5uOukSJEinHbKZMCAAWr+M4HmCQ7vhBBCooaZLmL84qXhTtCxVk393lciyhb/49p10q7GizJx2Ypoux4hwC/BqUCBApIyZUq3gkl0s2jRIjlw4IAKRlaeeeYZfZnUr19fBSekSYDgBCHrzz//DFec0NBQh7nPlcSJE+uLEEJI9BHd6SIiM/+9WqWyJEyQQKauWq3mOWiaIKSZ5SCxS7BEu/olOCVIkEBGjRolnTp1kp49e0q6dOmceiFjxozhPosqprbJddNgmOlSp04tmTJlcpw6bdq0Gj0HINxNmDBBfZpgygMXL15UMx2+i2uCZSARQkhkJIgfX95t2lh6NqinQg/MZVHRNJlEZv7LkjZNtF6PRI1gSwvhl3M4/Ifat28vP//8szz//PNSuHBhp9e0adOi1souQHuELOTvvfdeuPN98sknqll69OiRQyjC8VWqVNH30D49ePBAJk2a5PjtsGHDpGjRolK2bNloKV9U8OTQiM8JISSYia6M4N5mi4+u65GoYfqh3f2/1EL4C5MtPg9E/EpHgJxK8P/xBCLt/DV54RowxdWsWVM+//zzcJeCn1WlSpVUo1SuXDnZsmWL5niC8BQ//v/kwrlz50qbNm3UiR2aKGipli1bpu+9wZcwRW9B1Ic7Wz+exPBERgghJHJMbYY7cxw0XCRuCZS0EL6s81ESnCDMWJNNxiSozOLFi6VOnToezX7ILA5B6dKlS5otvEaNGuGKhL3z1q5dqw7k+D5btmxeF9uXBvWGQBlIhBASKNDtwZ6cDQ3VNBGe2DxyuOTImDHOC+/LOh8lHydogFq3bi2tWrXSaDZP9OrVS4/zh1SpUslrr70W4SkQIRfZdZAEEy87EJlDI763w0AihJBAgdni7UlkfmiBmBYiSoITzF5IAQCBBT5OnihWrJg9ezKOCcaBRAghhLgSE2ko4pooCU7WHEdvvPFGXNch4AjGgUQIIYS4I7rTUMQ1UfJxsvL48WNZs2aNHD9+XP+3AiftUqVKxXUd/cYX26e30KGREELIk4Sd/dB8WeejpHEygZM4nLZXr16t4f2IorOC/ErBIDjFBNGdz4QQQgixM8Hih+aX4LR3715ZtWqV7Ny50yl7N/GeYBlIhBAS09hZY0H8x9q/IDqSpMYEfglOyNCNF1IAEEIIITFBsGWeJs649m9C5N+KF8+xUXOn2v/bJscuebn8yg6WM2dOqVatmowYMUI3zEX2busLpjxCCCHEH4It8zRxxrV/Hz5+rEITsGOWcb8EJzBw4EDN6J0mTRrdC876+uKLLzg+CCGERBmYb6asXO3298gWju9J4BJR/1qxU1/7Zaq7c+eOVK5cWTp37iwtW7Z0bKJrkj179sDtTUIIIXEOEwYHNxH1r5U798NskxzaL8HpxIkTkihRIhk9erQkSJAguHs3DqFDJCHkSYUJg4ObiPrXSrLEIbbJMu6XqQ4pCAoVKiQHDx4M7p6NI+Awh82AS/fsrXv94C/e43PiOxBAsW8SVfuEBA5mwmB3dKxVU78ngUtE/WvFTn3tl8bp2rVrUqRIEalVq5a88MILkixZMqeKtmjRQl566aXA7dE4xnSYMzGd5MC7TRvHcekCB0bkEBLYBFvmaeKMa/8mTACdzv+i6qBpgtBkHmMH/MocfvXqVXUO90TTpk2lRo0adqinX/iSUTS6gFYEGiYIS65gIO36aiyftLwEWrrxi5eGO7png3oUQAkJIOi2ENzEZR4nX9Z5vzRO6dKlkwkTJgR3T8YRdIiMHiKLyIHwRFV/YMHF88mFCYODG9f+zZExoy0r7JfgRGIOOkRGDxRAgweaXIMDCr4k0KHgZFNMhzl3JiY7OcnZHQqgwQN9/gIbCr4kWPArqo7ELHCGgykJPk0Af/HeTk5ydocROcFBZCZXRkraH2b/JsECNU42BvvyIHoOwpJdNzsMBBiRE/jQ5BrYRCb40teQBBLRIjitW7dOFixYIIcPH5Zff/1V9u/fL8mTJ5cSJUoEUlvYFjpE+gcF0MCHJtfAxo6CL32tSFTx21Q3atQozeN0/PhxWb9+vTx48EA2btyoW7GcO3cuquUiJNoxBVBG0QUeNLkGNqbg647YzgjNxMLEX/wSnG7duiUffvihbNmyRZYuXSqJEyfW8vTv31/KlSsn06dP97d8ugdetWrVnF4zZsxwOu+sWbOkTp06Ur58ec0rFeaS++jx48cybNgwqVixouaV+u677/wuFyEkdqHPX+BiJ8GXvlbEX/wy1Z06dUpSpUqlQpIrZcuWlUuXLvlVPsMwZP78+fLOO+9I3rx5HecqVqyY4/+pU6fKm2++KcOHD5csWbLIkCFDZM+ePbJkyRLHMV26dFEt2ODBg+X27dvSr18/uXLligp4hJDAIDpMrjTPxB128DWkrxWJDvzKHH7nzh1JkyaNHDp0SAWblClTqjCFxJjYgqVt27bStWvXKJfz4sWLkjVrVjl9+rTkzJkz3HkePXok+fLlUwGoW7du+j3Mg3ny5FEtGIS3kydPSv78+WX37t1SvHhxPWbOnDkqTJ0/fz7cNjHu8CWjKCHEfjAU3j7EpfCKvSqx76cnNo8cbtukiyRm8WWd98tUh73p+vTpI9WrV5eJEyfKw4cPZe7cudKgQQPVNr322mt+1fTMmTOSIEECFcR27dol//zzj9P59u7dK2fPnlVznkn27NnVv2rRokX6ETRPEJhMoQk0bNhQy7p582a35bt//742ovVFCAlcaJ6xD3Hpa2gnXysSuPglOAH4Dn3wwQfy/fffayTd6NGjVXj57bff9L0/mIJToUKFpHbt2pItWzZp3ry5aroAtEsZMmSQ1KlTO10G2i8IVOYx0EpZCQkJ0TKax7gCsx/Oab5y5MjhTzUIIXFIZOaZqOSAwm+gvYjKb0ncYSdfKxK4+C04xY8fXzp27Cjbtm2Ty5cva0oCaJ8yZ87sd6vAd2rKlCmyc+dOCQ0NlQMHDsjWrVvVjwncuHHDrXCWIkUK/c7bY1wZMGCA/t58eRKwCCH2x5tQeG9hRFbgwyAD4i9+OYebmALTvXv3nMpToEABt75J3gK/ptatWzsOL1KkiPTt21cFM2i64AyOa7sCIQu+UQDH/PnnnxEe4wqiA80IQUJIYBOdOaC47UvgEx1BBuTJxm+N0+eff65CStWqVaVRo0ZOLyTF9IcLFy6oc7eVtGnTOkx1EMxu3rwpBw8edDpm+/bt+h3AX2is4NNkAqdzaJHMYwghwUt0mWdiwuRH4o649LUigY1fghO0PTBrDRo0SP+HEGN9vfXWW361zuLFi6VkyZIq6ADkZ0L6gSpVquj7XLlyqWO6aboDs2fP1ig80zG9cePGmpRz0qRJjrJAW1W0aFGNuiOEBD/RYZ6JTpMfISRw8SsdwZEjR6RMmTIqNMHhOrqB6Q9ZyaExgrCEtAe4zoYNG9QEB7C9C47JmDGj+lVt2rRJvvzyS003YIJIvzZt2kiFChU0j9PRo0dl2bJl+t4bfAlTJITYF39C4fHb0j17y12XBLsA2a93fTWWWekJiQH8uW+9xZd13i/BCRQuXFjWrFmjUWoxAbJ+Q8iBOQ75merXry9JkiRxuhR8mlCGa9euqYAFbZIrx44dk7Vr10rSpEk1ezgi9LzFlwYlhAQvo+bOl/GLl4arIDapxYbchJDoIzbzr/myzvslOEF789VXX+nWKq+++qrmW7KCLU5KlSoVfa0YR/jSoISQ4MWcyKeuWi137oeppqljrZpqCoTTMSEk+ojNBxVf1nm/BCckpEQySU/06tXLKSouUPGlQQkhwU9smA4IeZKJbdO4L+u8X+kIEE2HCDZCCHmSMCOyCCExgzfBGDkyZoyZi0dClAQnJLmE0IS96SCleQK+SMyHRAghhBBfiM78a9FNlIzycK4eMWKEbomCTX49vcaPHx/d5SWEEEJIkBNd+ddigihpnObMmSO5c+eW9OnTa4oATyBFACGEEEKIryDoArgLxohL/HIOf1LwxWmMEEIIIdFHbARj+LLOR0njhG1WNm/eHGmrYPNfZO4mhBBCCIkKdgvGiJLgBDMdtlaJjEyZMkV2CCGEEEJIwEBTnRf4osIjhBBCSGDhyzofJY3T1q1bNaIuMpA1vFChQpEdRgghhBASEERJcJo5c6YsXrzYUUFTiEIU3fXr1+X+/fu6d90nn3xCwYkQQgghQUOU8jh9+eWXcvLkSX1hj7oOHTrI5cuX5eLFi3Lnzh2ZNWuWCk4NGjQImoYihBBCCPHLx+n8+fNSsGBB1TIlSpTIqTV79OghTz31lAwcODDgW9kX2ychhBBCAgtf1vkoaZxMoF0yDENfruAzbyLvCCGEEEICBb8Ep/z580u2bNmke/fuqnUCjx8/ll9//VV++OEHqVWrVqC0AyGEEEJIpETJOdwkXrx4MnfuXPVzgmN45syZVYCCtunTTz+VF198MdICEEIIIYQECn4JTqBkyZKyb98+2b59uzqLY/+6Z555RjJkyBAobUAIISSIiI0tOsiTi9+CE0iQIIFUqFBBX4QQQkhc8OjxYxkzf6FMWbla7oaFSdKQEOlUu6ZuFpsgvl+eKYQ48FtwOnLkiOZr2r9/v9y7d8+pad977z1NVUAIIYTENBCaxi9e6rgMhCfz/btNuW8qiR78EpzCwsLUATxlypTSunVryZIli1OpypQpEz2lJIQQQiIA5jlomtwxddVq6dmgnm4WS4i/+CU4HTt2TM6ePau+TTlz5vS3LG6Bszk0Whs2bJAUKVJIy5Yt5fXXX5f4/6d2ReLNZs2ahfvt2LFj1f8KINJvxIgRsnTpUkmaNKm0aNFCunTpEiPlJYTYH/rABB+hN26ohskdd+6H6fc5MmYMvoqTWMcvwSlXrlySNm1azRIeEyA6r2HDhiokQXg6cOCA/oUwNWDAAL3kiRMnZPPmzTJx4kSN8jNBhJ8JhKSNGzfK4MGD5fbt29KvXz+5cuWK9O/fPyaKTQixKfSBCV4ypk6tPk3uhKdkiUP0e0KiA78yh4PvvvtONU59+/ZVISo6OXjwoDz99NNy6dIlyZQpk576+++/l6FDh8rRo0f1/Zw5c6RPnz4eNx1G2ZBvavfu3VK8eHHHbyBMIfN5smTJIi2yLxlFCSH2ZdTc+U4+MCYw49AHJvBh/5Ko4ss671eYATb3hdDy+eefS7p06VTjY32NHj06qnVQ8uTJI6dPn3YITeDq1atqejNBGWAmvHbtmuzYsUNu3LjhdM0lS5aowGQKTQBarIcPH6qmyh3YpBiNaH0RQgKbyHxg8D0JbBA9ByEYGiaAv3iPzwmJLvwy1cEZfM2aNR7L4q/fU5IkSZzOAUEGJrlOnTo5CU7QPsFsmDBhQt3m5cMPP5QPPvhAjzl37pzky5fPqYwhISFqXoR/ljuGDx+uZj1CSPBAH5jgBykHoDmEsIT+hnmODuEkuvFLcIIAUr58+eguk1sePXqkjuFIrAmzoEnbtm3lueeeUy0SHL+hYWrcuLE8++yzUrduXdVAJU+ePNw54Wjuqp0ygf8UNGkmENhy5MgRK/UkhMQM9IF5coCwREdwElNEyVQHH6OffvrJUaZ///1X96az8u6778qMGTOirdyIpINmadGiRZIoUSLHeZGlHFFyEJpA/fr1VXBasGCBvodWDJF3roSGhkrWrFndli9x4sSSKlUqpxchJLDBYopkiO7oWKsmNROEEK+IkuCELVZOnTrluAAcsCEoWTl+/Lg6dUcHiH5buXKlrF69OtxWLhCmILhZgZP6nTt39KMCBQrIzp071afJ5OLFi2qmw3eEkCcH+sAQQvwlSoJTbAIH86lTp6ovFfyYXEF6ApjkYMoDEIqgbapSpYq+h/bpwYMHMmnSJMdPhw0bJkWLFpWyZcvGZlUIIXGM6QOz66uxsnnkcP0Lnxhux0EI8Ra/fJximi1btqg/Exy5O3fu7HS5X3/9Vc1wEJwqVaok2bJlk3LlyulvKlas6HAgR7oBpExo06aNpi5AHidoqZYtWxbTxSeE2BT6wBBCooqtBSdE5U2fPt1t3bDNC0DKgsOHD6uWCabBXr16SY0aNZx+07RpU80ivnbtWvWFwvcQtAghhBBCfCFKCTAR3bZixQpJkyaNXgt5jyC0WNMPwO8I2qB33nnHl/LYEl8SYxFCCCEksPBlnY+SxqlDhw5qHosMmMwIIYQQQoKFKGmcnjR8kUQJIYQQElj4ss7bPqqOEEIIIcQuUHAihBBCCPESCk6EEEIIIV5CwYkQQgghxEsoOBFCCCGEeAkFJ0JIUHMvLEzOhobqX0II8Zco5XEihBC78+jxYxkzf6FMWbla7oaFSdKQEOlUu6Zu9Mu96QghUYWCEyEkKIHQNH7xUkfdIDyZ77GxLyGERAWa6giJBWguil3Q3tA0uWPqqtU02xFCogw1ToTEIDQXxQ2hN26ohskdd+6H6fc5MmaMm8KRcEDQRZ9kTJ1akoSEsIWIW+wyTig4ERKD0FwUN2ABhk+TO+EpWeIQ/Z7EPXywIN5gt3FCUx0hMQTNRXEHtBZwBHdHx1o1qdWwCeaDhSngmn5o+JwQE7uNEwpOhMQQ3piLSMyB6LmeDeqphgngL97jcxL38MGCeIMdxwlNdYTEEDQXxS1IOYDoOQhL9J+xH948WNAPjdhxnFDjREgMQXORPUA/YAGm07G9MB8s3EE/NGJix3FCwYmQGITmIkLcwwcL4g12HCfxDMMwvCn8k8x///0nqVOnlhs3bkiqVKme5KYgUcQuYbSE2AkzWgq5tZAmApomOO8zuzuxEhvjxJd1noKTF/jSoIQQQnyDDxbEG2JynPiyzj8RprpHjx7JJ598IuXLl5eqVavKhAkTvOkjQgghsQD90Ig32GWcPBFRdZ06dZItW7bIp59+Krdv35a+ffvK1atXZeDAgd70FSGEEEKIEvSC0/Hjx2X69Omyd+9eKVq0qFY6TZo00rFjR+ndu7ckT56cQ4EQQgghXhH0prply5ZJ8eLFHUITaNCggZrvNm/e7LaR7t+/r/ZO64sQQgghJOgFp3Pnzkm+fPmcejpRokSSI0cO/c4dw4cP1yg684VjCSGEEEKC3lSHFALuzHH4DN+5Y8CAAdKnTx/HVzguZ86cqoEihBBCSHBhru/eZGgKesEpa9as8scff4Tr4dDQUMmWLZvbnk+cOLG+TMwGpeaJEEIICV5u3ryp1qaICHrBqWDBgvL111/Lw4cPJWHC/1X3woULcvbsWSlQoIBXvQ8BC8enTJlS4sWL59VvXIHwBcEL5wnGJJrBXj/AOgYHwd6PwV4/3ovBg53GKjRNEJo8KVSsBL3g1KhRI+nZs6fmbnrzzTe17kOHDpVixYpJmTJlvBqB8ePHl+zZs3t1bGRAaApGwckk2OsHWMfgINj7Mdjrx3sxeLDLWI1M02QS9IJT0qRJZcqUKfLaa6/Jr7/+qnmcTp8+rdF2hBBCCCG+EPSCE4DWad++fbJu3ToVpKpXry6ZM2f2pZ0IIYQQQp4MwQnkyZNHM4jHFXA2/+ijj5yczoOJYK8fYB2Dg2Dvx2CvH+/F4CFQxyo3+SWEEEII8ZKgT4BJCCGEEBJdUHAihBBCCPESCk6EEEIIIV5CwYkQQgghxEsoOBFCCCGEeAkFJ0IIIYQQL6HgFAMgX1Tp0qXl1q1bHs9+48YN6dy5s+TPn18qVaok06dPj4GSxAz379+XWrVqSe3atSO8wCuvvKLtYH19+eWXMVOoaGDHjh3hyouy//vvv27Pjr2Nhg8fLs8884xu4TNgwADdE9HOzJo1K1wdK1eu7LHII0eODHc8svDbmQcPHki/fv20T4oUKSKDBg3SMeuJw4cPS+PGjSV37txSt25d2bJli52rp1y7dk3nGcwfpUqV0v04PXHnzh2343rDhg22refPP//stsy4H69cueJUbtyH2EYL3xUvXlwGDhwojx49sm3dTD755BO3daxWrVq4so8YMSLcsW3btrV9HTdv3iw1a9bULcuee+453b0jIo4fPy5NmzbVe/Gll16SjRs32rKOT0Tm8Nhkzpw5MmPGDAkLC/N48+JGr1Gjhm4aPHHiRDl58qR069ZNj2/fvn1sFjdKQFjAgM6UKVOEv1+7dq306tVLJzOTQoUKRemasQEW0MuXL8vYsWOdLod+ckf//v21r8eNGydJkiSRd999Vzdxxmd2Zc+ePZIsWTLp06ePo4iJEiXyWNw///xT8ubNK61atXIckzZtWo/H24FXX31Vzp8/r4sNhF4sUKGhofLtt9+GK94///yjDy4QnH788UdZuXKl3pu///67LlR2BHNLlSpVVGjCHpx///23Cu2PHz927MdpBWNy9+7d8tNPPzkl4MUG6HalXLlyKvBaWb16td5buNes4L6bPXu23ocYy++88472//fff2/X6in16tWTokWLOpUR24Ohv1z5448/tL8xtk3SpUtn6/odOnRIH7C7dOmi4xMPJG3atJHkyZPrA4ormHsrVqyoAhPuRez0gf8h4JcvX95edTVItHH9+nUja9asxpAhQwwR0ffuWLp0qZEyZUrjv//+c3w9duxYI3/+/Mbjx4+jrTwxwcGDB43kyZMbH374oZE9e3aPl7hx44a2AY4PFIYOHWpUr17dq+Jeu3bNSJYsmbF69WrH8X///bcRP3584/jx416dIy5o2bKl0atXL68vXbZsWeOrr77y+vi4Zvfu3UaqVKmc7r1Vq1YZ9erVMx48eBCueB988IFRoUIFp89fffVV47XXXovrqnhkxowZRvHixY1Hjx45jpk0aZLRtWtXt79ZuXKlkTp1atvWxxswL5YsWdLo16+f0+GXL182kiZNaqxfv97xOcZAggQJjFOnTnlzatuAOTNt2rTG9OnTw5WpdOnSxoQJE2xTVm8YPXq0UapUKadDcV917NjR7c8x/6Ke1jWwffv2RtOmTb25XKxCU100Ag1E2bJlVdUYEcuXL5c6deo4aTKaNWsmx44dk6NHj0ZjiaIXaMq6du2qWqQSJUpEePIzZ87o99htev369bJ//379vZ1BmXPmzKl9gKcdaCk8AY0bNDcvvvii4xCYhqBRW7FihW2riTrmyJFDNRCbNm2SmzdvRtqP2bJlk23btsnWrVvl3r17tq0b+OWXX+Tll1/W3db37t2rGjNoZ5YsWSIJE4ZXsONebN68uVOdcC/aeRNw1LF169aqoYYmAhonmO3cadSAOa4vXrwoa9as0U3OA4358+erGQfaJSvQRqCv0ccmJUuWVC0ptIeBBDTd0OK3bNkyXLHRh+Z9iHvR7vchgJkcWiRoSE2gCXz66afddgvuRdx78eLFc3yP9/jcbtBUF01AtQ+7PASE27dvR3jWc+fOhTNZ4abAvj0YWHZVoU+ePFkn3w8++ECWLl0a6YIbP358efbZZ1U1ixsfNu558+bZVsWMMsLPCYss+gJmHvj4vPXWW+Hqij7MlSuX1tEK9kREH9oV1PGzzz7TSRp+P/DJglq8fv364YqM79EGPXr00M2x4VeDRWru3Ln6gGBHUL8UKVKo3xkE3//++0+yZs0qCxculMKFC7vtR/SZax+irvANgnBsN1BHzDEoN0xTuCfRHxAu0qdP73bRxW8g2ONB5tSpUw5By7pI2RU8cA0ePFjNkBkyZHAqJu5D+MO41sPu96Er8HkdM2aM+qolSJDA6Wv0Nfy68NCK8Yj/YS7HfYj51a7A5N2wYUP1wYNZEmskzIsdO3Z0W2T0pbt7Efch7kc7uQhQ4xQNwBkVgxoOinCCiww4jbubkLE4YaK3I5cuXVKHWzzVuvoYuAMTORarffv2qQYHflxYyKCVsytDhgxR50UIC/AzmDp1qvoC/fXXX+GK7KkP0a927UOAIAQ8waF+6NPu3bur34Grwy2AUAjtGbRr6D+0CxzJ4RxuV+dbCBGTJk1SPxdMxCgznnA9OdK660fz3oxMGxdXoI6jR4/WsYp+gVCEsvbt29dtkdq1ayczZ86UCxcuqNYGWji8nzZtWlxVwSfgN4p6ok9dCdT70JUvvvhCsmTJ4uTDZALhGNozaIhPnDihYxoPobgP4ddmV9A3uAfv3r2rwjrmRawB0EK5I6J70XZzaqwaBoOUTz/91ChfvrzD5+DQoUMR+ji1adPGeOONN5xa4/79+/qb7du327KVWrRo4WSbnjNnToQ+Tu748ccfjWzZstmyfp4oV66cMWzYsHBf//DDD0ahQoXCfV65cmVj+PDhAVPHhw8fqr/dokWLvCrzhQsXdJweOHDAq+Njm1atWhkvvPCC02X379+vZb506VK44hQoUCCcT8mmTZuMRIkSOfkQ2YkSJUoY3bt3dyrSL7/8YmTMmNHrYvbo0cNo0qSJ18fHFfB3KVq0qDFw4EC3RZgyZYp+78rzzz9vjBw5Mq6K7RNXr15VH7SZM2d6/bszZ87omD5y5IjXv4ltunTpYlSqVEnXNoC+hH8T1kp3wG/vu+++c/rqzz//NOLFi+c4h12gxika+OqrrzSCAKYb+I+Y9nY86SJKwhWYDFxDnuE/Ah8M2ObtBp4S4FexaNEirR9e0LAhIgn///bbb+GKDD8KV00NnuShnbMr6BPX1AOeyow+hE+a9Xg8WaHOdo0cRNtDe2Q1JcMsALOku37BUx6Otz7VmtoYu/Yj7h9Xjahp+nCnJXN3L2IcFChQIJwZ1i54qqMnLSB8oKBpsmL3e9EE0XLQqFmjQK2g/44cOeKkxcD4RvSoXe9Dd9omuGq4+tqZwIyH+9DqI2r3+xDgPmrSpImEhIToe5hTW7RooRpPd/OHp3sR5jvzHLbBLhJcILNt2zbjt99+c7ygWcHTwPLly43z589r1RYsWGCcPn1a/8dneKKdO3euvg8LC9NoLkTz2BGUz1o//P/xxx/rEy7+R/QgIpbwxIRoMzB79mytI572wZ07d4yqVatqVJddgbbIGn21detWrQPqCHbs2GH8/vvvjuK7PvkPHjxYoyrRXnYET3wZMmQw3n33XUfxpk2bptGB//zzj36GCCxoTMG///6r0UnWaJ5Bgwap1tCudYQmLEmSJMbGjRu1DtCoQcNbrFgxfX/v3j3j559/dkS0QtOGKDxTg4Z2yJkzp0a52pV58+ZpP5rRmzdv3lQNizl/INJs1qxZWnfQp08fI3fu3I46Y/7BOB03bpxdq6ggavDpp5/WMWcFWgjcmybQOL355puO94iURMSvuyhKu3HlyhUdf+gvK1g7TG0SxiSida3amP79+xs5cuRw9LEdad68ua5rt27d0uJhzujWrZtRuHBhx3usGWb0OaJfEbG9Z88e/T40NNTImzevMWLECNtVj4JTDLBv3z4nUx0GBtSNUI+bTJ482UicOLGGQuMmxwRx9uzZGChNzIAb3ZqOAAIFhMVvvvnGMem98soresNDNYuJGqpYmHrsyl9//WVkypTJSJ8+vVGxYkXtH0zCJgjNL1iwoOP9zp079XhMBAi7hclrzZo1dq2eMn/+fBWUcuXKZZQpU0ZDuSHoAyxUKVKkcDLhYHFNmDChUaRIEV2g0qVLZ6xdu9bWdfz8889VeIIwAYHhqaeeUqEXQDDEODXrDGESZi+0CQTnNGnSGA0bNjTu3r1r6zrC5IFFBmXGGIQQbz6kIQwcddy1a5fWAYIU7j30NUwnWKgbNWpke8ECAi7KClOWFYSsY740gXsDBEl8hpQF+M26dets3X8m77//vt5X1tQSEHjRVxA8TL744gt9iEEdcS9ijtqwYYOt6whzItwZMI9i3MEcibG6ZcsWLTcebjBOkUrDpHfv3lp3mNsx17z00kvG7du3bVfPeCiRbdRfQQIcNeHIh8zaZgg0wmYRnpk5c2ZHLeHoh0gDfIYkfHaM4PEEHE0PHjwo1atX10NgJli1apU6DyOKzgShs3AQh0kPofsRJVu0AzBPIWQbf5F0DX1mgvpCxWxNxXD16lXta3yOulv7167AxIrxiFsfZbYGNMBkDLMBzM4mMPPAHItoNZihXSOb7AhM56gLyooyIxoQIIoQ4xTZmRGMYQITK1I0wFyA5IuukU12BCYPpFxAX6GOpjkDYxdmEsw/ZrQZ6o1xDWddJF2Ec7Hd2blzpxbRNXLswIEDOt9YE+siuAH3IT7HmI4sOa9dwH2I+801khr9Z7pFmMA1AOsF7sOqVau6jaC0Gxh3SEcD535Et+K+Q/kjWjNgZt21a5e2CeZgd2lE4hoKToQQQgghXkLncEIIIYQQL6HgRAghhBDiJRScCCGEEEK8hIITIYQQQoiXUHAihBBCCPESCk6EEEIIIV5CwYkQQgghxEsoOBFCCCGEeAkFJ0IIIYQQL6HgRAghQQi21Llx44b+jS1GjBgh3333XWxdjpA4gYITIUHM8uXL5aOPPpLBgwfrXmWxRadOnWTx4sWxdblIwX5ZS5YsifQ4cOvWLRk9erT+Bb1795affvrJq9/GFtizbP78+RFebvr06VKyZEmvi4R93jBWnn76aenWrZvH350+fVpGjRol77//vixYsMBJMMuZM6f0799fbt++7fV1CQk0KDgREoRgEWzYsKE0b95cN1nGhqmNGzeWli1bxrgGAhoObHR99+5d27TsokWLvBZ+oKV599139a8pSN25c8c2dQErVqyQqVOnRlimzz//XNq2bevY6Dcy5syZIyNHjpRPPvlE3n77bbeHb968WTe+/vXXX+Xs2bPSoUMHeeWVVxxjqlGjRrrhtbdCKiGBSMJALDQhJGKwcGGRO3HihGTMmFEPPnXqlBQqVEi1Ly+++GKMNWG8ePFk9uzZMXb+2Gby5MkxckkItwkSJIiRc0NY3r9/v2qEvAXj45lnnlFByBMQqtq1aycTJkzQQ6B9ypMnj+zYsUPKli0ryZIlU+EJ2sYWLVp4e2lCAgpqnAgJQs6cOSPp0qWTDBkyOGqXO3duOX/+vDz33HOOz7CwwpyTNGlSKVWqlGoyTB4/fiz9+vWTzJkz6/cQti5cuOD4fuHChVK0aFEJCQmRggULOmkZXnjhBZkxY4bjWAhxFSpU0IUVv5kyZYrju5MnT6pwt3XrVsfiW7NmTS2ryZ49e6Ry5cpajixZsqhmxBPQeA0fPlzrmzJlSunevbsKKVZwbizwadKkkWzZsqnZKSwszO0pcdyXX34pDx8+lMKFC8uPP/7odNyQIUOkVq1a+hmuA7MohIkUKVJI3bp15ejRo47jca4vvvhCateuLUmSJJFLly7pd+PGjdM2RPtUrVpVdu/e7XSNMWPGSN68efWc7du3V61ORGzatEmKFSsm+fPndxwGQap69eraJvgc/kjoY4C2glD0559/al94al+0Lcphkj17dq2HeR6AcXT48OEgvKsI+T8MQkjQcfLkSSNFihRGgwYNjPXr1xsPHz4MV8d169YZqVKlMmbOnGmcPXvWmDhxohESEmL8/fffeuyECROM9OnTG1u3bjXOnDljtG7dWl/gwoULRqJEiYxx48YZly5dMqZMmWJkzZrVuH79un5frFgxY/Lkyfr/gQMH9LzDhw83jh49asyaNctInTq18cMPP+j3x44dg/ey0aJFC2Pv3r3G8uXLjfz58xutWrVylDlfvnxG586djXPnzmm5M2TIYOzYscNtv6FMKPeCBQv0ev3799frvfrqq3p8WFiYUaBAAaNbt27G8ePHjS1btuj13n//ff0e10B58BdUrVrVGDFihP7fq1cv48UXX3Rc9/Hjx0aePHmMr776Sj8bOHCgUbRoUeP33383Tpw4YXTo0MEoWLCgcffuXce5MmbMaHz//ffazigLfpsrVy5jzZo1xqlTp4x3333XyJIli3HlyhX9DdoW5Z89e7a21eDBg/V9/fr1PY7bt99+22jatKnj+9DQUCNNmjRG7969jcOHD2sbZ8+eXc8FTp8+bfTr18+oUKGClsHsR1dQBtT3p59+MrZv3260a9fOKFeunLaDycKFC7V8hAQrFJwICVIgsDRr1sxIkiSJkS5dOuO9995TAcnkpZdeciycJs8//7zRp08ffYuFtHjx4sbt27f1PYQvc4HEognh4q+//nL89sGDB47/rYJTp06djEaNGjm18ujRo40iRYroZ6bgZD0XvscCDR49eqSCFwQ5E+u1XMmbN68KT1ZKlizpEJwguEFQwnlNJk2aZKRNm1bfRiQ47dmzx4gfP76jHTdt2mQkTpxYhRwIQUmTJlWhyeTq1av6/bJly/QjCGEQXqxAgJk7d67jo/v376tgOHXqVP2sRIkSxtChQ51+g36KSHCCwIz+M/n000+NUqVKOZ0D14RwbbYDhEOULyJu3LhhNGnSRPser+TJkxurV692+sm+ffv0OwhrhAQjNNUREqTAiRe+RlevXlU/ne3bt6sPC/xSAPxg8H21atUcLzj8HjhwQL/v2rWrRpbBHPPqq69qlJZpInr22WelWbNmUqZMGTX/DB061GF2cgUmoooVKzp9DFPekSNHnMxjuI4JfLHgYA7ix4+vZqQ33nhDy9+3b189pzvu37+vvjrly5d3+hrlNEG9r127puU2644Qenxmms48UaJECSldurTMnDlTD4E58uWXX1azKPzJ4BAP86Z53iZNmkiiRIkcbQpgajSBA/q5c+fks88+c/wGZjyYvtAPMI2hnVzrA3+iiMA5n3rqKcchnvoA/QtTqbegrvjN8ePHNXJu4sSJ0qBBA9m4caPjFOgDkDhxYm9PS0hAQedwQoIc+AUhog7+NRA8ICxB+IBvEhY9+NRYgd8PgC8LFm04k+OFMPOff/5ZVq1a5XAAx4K8evVqmTVrloaowxcJIelWsIC6+uTAnyhhwoROztFwKvcEhJHWrVvLypUr1fEYghAi5erUqeP0E/Oc5uJtYn2PesOnB/VxBf4/EKAiAqkWvv32W01TgLZE3QHOC3r27OloQxOrr5EVCFWoN6LT4BdlJUeOHPodzhtRfdwBvy2rP5qnPgDwUfKGY8eOyYYNG/Sv6eeEPkEgAgTzKlWq6Gng05U1a1ZtS0KCEWqcCAlCEIbuKlRAWIGGw9QEQPi4fPmyajjMF7Q7pkYEztpwUoaj9rBhw1S7ACEJizg0U3Ash6M3BAg4I+P8v/32W7jWxHWWLl3q9DkWWwhx3kSVQfMEgQ1Oyx07dtT8RdDAIEeVKzhf8eLFnb6D0AbBzwTlgVYKzupmvfPly6ftAufsyEBKBzg/w0kcztpoH5ArVy6HM7553ho1aqhmxoxsdAXXg2YQbWrth3v37kmmTJn0cDhbW+sKbVRkObngaA6h11pnCLxwcDdBH0DAsmqmIsJsG2gwreB98uTJHR8dOnRIr09IsEKNEyFByJtvvqlRaIjqggCFxW3evHkq3GDhBx9//LEKINBIwbSEXE+ItIJgAiEAC+u0adM0SgxaCUR+QTsFAePixYsatv7ee++p2QjaHwhdzz//fLjWxDGI8MJ18RtcZ+zYsU4RfJHx1ltvaRQfTIZ79+6Vbdu2udUYASRxxHWgqYEQhQi+1KlTOy6B9oAwgkg4CH0wPcFUBhOcq/bNHdDIwQQ3YMAANWei3gB/cR5E8UGwhFbrl19+kbVr16pmLlWqVG7Ph9+0atVK/vvvPxVEIWTCFIh2QlTkhx9+qGWGkAahDwJbZGYw1Bv9DVMftFadO3eW8ePHS7169VRjBvPcp59+qpozb4GQhfaH5hJ/oSHDGEEyTmsUHjSC0EQRErQEo+MWIcTQqC04Zj/33HNGtWrVjL59+xoXL150apojR44YXbp0MSpWrKgO3CtWrHB8D2fwkSNH6m8rV65sDBgwQJ2dTeAU/PLLL2skVps2bZycuxEZhug3k/Pnzxs9evQwKlWqpJF5cC43QZkQ8fbff/85Ptu1a5eWywSO7m3bttVroZyLFi2KsIvhjN24cWOjSpUqxrRp04yff/7ZyWEcDtjDhg0zatWqpQ7RcL6+d++enhOO3iiPGdX2ySefaKSYlW3btukxiMpzBRFr+A7tjvZHFJuJu3OB3377TdsFv8Ff1N/K2rVrjVdeeUX74ZtvvjHmzZtnfPbZZx7b4ObNmxrZtmrVKscxiJRDv+AczZs3N1auXOn0e9QR5YsIOOUj+hL9jrbt3r27U/0QFZkyZUp1IickWImHigWtVEgIIU8o77zzjubzQpbv2ALaRWjeoLkkJFih4EQIIUEIhCb4pSFBp2lOJIT4DwUnQgghhBAvYVQdIYQQQoiXUHAihBBCCPESCk6EEEIIIV5CwYkQQgghxEsoOBFCCCGEeAkFJ0IIIYQQL6HgRAghhBDiJRScCCGEEEK8hIITIYQQQoiXUHAihBBCCPGS/wdK/8ZG8TgFQQAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 600x340 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "tg = df[df.treat == 1].groupby('group').agg(sessions=('sessions', 'first'), attendance=('attendance', 'first'),\n",
    "                                             y=('y', 'mean'), x0=('x0', 'mean'))\n",
    "print(tg.sessions.value_counts().sort_index().rename('groups by sessions delivered (of 8)'))\n",
    "print('mean attendance share:', round(tg.attendance.mean(), 3))\n",
    "fig, ax = plt.subplots(figsize=(6, 3.4))\n",
    "ax.scatter(tg.sessions + np.random.default_rng(0).uniform(-0.15, 0.15, len(tg)), tg.y - 0.6 * tg.x0, s=18, color='#1f5f63')\n",
    "ax.set_xlabel('Sessions delivered (of 8)'); ax.set_ylabel('Endline minus 0.6 × baseline')\n",
    "ax.set_title('Descriptive dose–response in treated groups (not causal)')\n",
    "plt.tight_layout(); plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.13.18"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}