{ "cells": [ { "cell_type": "markdown", "id": "b689393c", "metadata": {}, "source": [ "# The Half-Life of History: Reproducible Analysis (v2)\n", "\n", "This notebook is the generator of record for `results.json`. Pipeline: `data.json` (input) -> this notebook -> `results.json` -> site components. Re-running this notebook end to end regenerates every published statistic.\n", "\n", "**Spec (canonical, v2):** per-metric z-scores standardized on the 32 fit periods (the partial 21st century is excluded from both standardization and fitting), with `ddof=1` (the ddof choice rescales every metric by the same constant and leaves the half-life unchanged); composite = arithmetic mean of the four z-columns; x = years before present (2026) at each period midpoint; primary model = exponential with intercept, y = a*exp(-lambda*t) + c, fit on the 32 points.\n", "\n", "**Lineage:** the originally published results.json was produced by a standalone script that standardized each metric on century-length rows only, including the partial 21st century. As of June 2026 this notebook is the sole generator. The standardization-variants table in the sensitivity section reproduces the originally published headline figure (84.0 years) under that old specification.\n", "\n", "**Changelog:** v2, June 2026: intercept-class primary model, canonical standardization, moving-block bootstrap, leave-one-out panel; headline half-life corrected from 84 to roughly 95 years.\n" ] }, { "cell_type": "code", "execution_count": 1, "id": "9b6d71ee", "metadata": { "execution": { "iopub.execute_input": "2026-06-09T21:29:18.589737Z", "iopub.status.busy": "2026-06-09T21:29:18.589603Z", "iopub.status.idle": "2026-06-09T21:29:20.174079Z", "shell.execute_reply": "2026-06-09T21:29:20.173624Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Loaded 33 periods\n" ] } ], "source": [ "import json\n", "import numpy as np\n", "import pandas as pd\n", "from scipy.optimize import curve_fit\n", "from scipy.stats import spearmanr, shapiro\n", "import statsmodels.stats.stattools as smtools\n", "import statsmodels.stats.diagnostic as smdiag\n", "import matplotlib.pyplot as plt\n", "import warnings\n", "warnings.filterwarnings('ignore', category=RuntimeWarning)\n", "\n", "CURRENT_YEAR = 2026\n", "np.random.seed(42)\n", "\n", "# data.json is the only input; results.json is OUTPUT (written by this notebook)\n", "with open('data.json') as f:\n", " raw = json.load(f)\n", "\n", "periods = raw['periods']\n", "print(f\"Loaded {len(periods)} periods\")\n" ] }, { "cell_type": "code", "execution_count": 2, "id": "7da88541", "metadata": { "execution": { "iopub.execute_input": "2026-06-09T21:29:20.175201Z", "iopub.status.busy": "2026-06-09T21:29:20.175104Z", "iopub.status.idle": "2026-06-09T21:29:20.186828Z", "shell.execute_reply": "2026-06-09T21:29:20.186434Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "=== Metric completeness ===\n", " openalex_work_count: 33/33\n", " named_individuals: 33/33\n", " source_proxy: 33/33\n", " ngram_discourse: 33/33\n", " films_set_in_period: 31/33\n" ] }, { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
labelstart_yearend_yearmid_yearspan_yearsopenalex_work_countopenalex_work_count_tiernamed_individualsnamed_individuals_tiersource_proxysource_proxy_tierngram_discoursengram_discourse_tierfilms_set_in_periodfilms_set_in_period_tieryears_before_present
03rd millennium BC-3000-2000-2500.0100010156T1385T150000T36.540000e-08T1NaNNaN4526.0
12nd millennium BC-2000-1000-1500.0100012471T1201T1130000T27.510000e-08T1NaNNaN3526.0
210th century BC-1000-900-950.0100526T111T1500T31.340000e-08T19.0T12976.0
39th century BC-900-800-850.0100619T117T1600T32.130000e-08T11.0T12876.0
48th century BC-800-700-750.01001713T137T12000T36.200000e-08T15.0T12776.0
\n", "
" ], "text/plain": [ " label start_year end_year mid_year span_years \\\n", "0 3rd millennium BC -3000 -2000 -2500.0 1000 \n", "1 2nd millennium BC -2000 -1000 -1500.0 1000 \n", "2 10th century BC -1000 -900 -950.0 100 \n", "3 9th century BC -900 -800 -850.0 100 \n", "4 8th century BC -800 -700 -750.0 100 \n", "\n", " openalex_work_count openalex_work_count_tier named_individuals \\\n", "0 10156 T1 385 \n", "1 12471 T1 201 \n", "2 526 T1 11 \n", "3 619 T1 17 \n", "4 1713 T1 37 \n", "\n", " named_individuals_tier source_proxy source_proxy_tier ngram_discourse \\\n", "0 T1 50000 T3 6.540000e-08 \n", "1 T1 130000 T2 7.510000e-08 \n", "2 T1 500 T3 1.340000e-08 \n", "3 T1 600 T3 2.130000e-08 \n", "4 T1 2000 T3 6.200000e-08 \n", "\n", " ngram_discourse_tier films_set_in_period films_set_in_period_tier \\\n", "0 T1 NaN NaN \n", "1 T1 NaN NaN \n", "2 T1 9.0 T1 \n", "3 T1 1.0 T1 \n", "4 T1 5.0 T1 \n", "\n", " years_before_present \n", "0 4526.0 \n", "1 3526.0 \n", "2 2976.0 \n", "3 2876.0 \n", "4 2776.0 " ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Build DataFrame with metric values and tier info\n", "rows = []\n", "for p in periods:\n", " m = p['metrics']\n", " row = {\n", " 'label': p['label'],\n", " 'start_year': p['start_year'],\n", " 'end_year': p['end_year'],\n", " 'mid_year': (p['start_year'] + p['end_year']) / 2,\n", " 'span_years': p['span_years'],\n", " }\n", " for key in ['openalex_work_count', 'named_individuals', 'source_proxy', 'ngram_discourse', 'films_set_in_period']:\n", " if m.get(key) and isinstance(m[key], dict):\n", " row[key] = m[key]['value']\n", " row[f'{key}_tier'] = m[key].get('tier')\n", " else:\n", " row[key] = np.nan\n", " row[f'{key}_tier'] = None\n", " rows.append(row)\n", "\n", "df = pd.DataFrame(rows)\n", "df['years_before_present'] = CURRENT_YEAR - df['mid_year']\n", "\n", "print(\"=== Metric completeness ===\")\n", "for col in ['openalex_work_count', 'named_individuals', 'source_proxy', 'ngram_discourse', 'films_set_in_period']:\n", " print(f\" {col}: {df[col].notna().sum()}/{len(df)}\")\n", "df.head()\n" ] }, { "cell_type": "code", "execution_count": 3, "id": "6d76cf4f", "metadata": { "execution": { "iopub.execute_input": "2026-06-09T21:29:20.187698Z", "iopub.status.busy": "2026-06-09T21:29:20.187642Z", "iopub.status.idle": "2026-06-09T21:29:20.191713Z", "shell.execute_reply": "2026-06-09T21:29:20.191413Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "=== Canonical z-scores (standardized on the 32 fit periods, ddof=1) ===\n" ] }, { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
labelz_openalex_work_countz_named_individualsz_source_proxyz_ngram_discourse
03rd millennium BC-0.286597-0.230595-0.224250-0.549383
12nd millennium BC-0.222242-0.231343-0.218923-0.547469
210th century BC-0.554300-0.232116-0.227546-0.559648
39th century BC-0.551715-0.232091-0.227539-0.558088
48th century BC-0.521303-0.232010-0.227446-0.550054
\n", "
" ], "text/plain": [ " label z_openalex_work_count z_named_individuals \\\n", "0 3rd millennium BC -0.286597 -0.230595 \n", "1 2nd millennium BC -0.222242 -0.231343 \n", "2 10th century BC -0.554300 -0.232116 \n", "3 9th century BC -0.551715 -0.232091 \n", "4 8th century BC -0.521303 -0.232010 \n", "\n", " z_source_proxy z_ngram_discourse \n", "0 -0.224250 -0.549383 \n", "1 -0.218923 -0.547469 \n", "2 -0.227546 -0.559648 \n", "3 -0.227539 -0.558088 \n", "4 -0.227446 -0.550054 " ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "active_metrics = ['openalex_work_count', 'named_individuals', 'source_proxy', 'ngram_discourse']\n", "z_cols = [f'z_{m}' for m in active_metrics]\n", "\n", "# The partial 21st century (index 32, only 26 of 100 years elapsed) is excluded\n", "# from BOTH standardization and fitting.\n", "fit_mask = df.index.values < 32\n", "\n", "def zscores_on_basis(frame, basis_mask, ddof=1, transform=None):\n", " # Standardize each active metric using mean/std computed on basis_mask rows.\n", " # ddof rescales all metrics by the same factor, so the half-life is invariant to it.\n", " out = pd.DataFrame({'label': frame['label']})\n", " for col in active_metrics:\n", " v = frame[col].values.astype(float)\n", " if transform is not None:\n", " v = transform(v)\n", " mu = v[basis_mask].mean()\n", " sd = v[basis_mask].std(ddof=ddof)\n", " out[f'z_{col}'] = (v - mu) / sd\n", " return out\n", "\n", "z_scores = zscores_on_basis(df, fit_mask, ddof=1)\n", "print(\"=== Canonical z-scores (standardized on the 32 fit periods, ddof=1) ===\")\n", "z_scores.head()\n" ] }, { "cell_type": "code", "execution_count": 4, "id": "d62670a2", "metadata": { "execution": { "iopub.execute_input": "2026-06-09T21:29:20.192617Z", "iopub.status.busy": "2026-06-09T21:29:20.192559Z", "iopub.status.idle": "2026-06-09T21:29:20.252230Z", "shell.execute_reply": "2026-06-09T21:29:20.251876Z" } }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAA90AAAHqCAYAAAAZLi26AAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjksIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvJkbTWQAAAAlwSFlzAAAPYQAAD2EBqD+naQAAYTVJREFUeJzt3QeY1FTb//F76b0XRZoKUhXBijyiCIJgw/LYULGXBxWwd8SG/RUVe8HHhgWwN0SqFBGkilhAQenSe9n8r99538x/dnZ2mVkm7M7k+7muYZkkkzlJTjK5c1qW53meAQAAAACAlCuW+lUCAAAAAACCbgAAAAAAAkRJNwAAAAAAASHoBgAAAAAgIATdAAAAAAAEhKAbAAAAAICAEHQDAAAAABAQgm4AAAAAAAJC0A0AAAAAQEAIuoE01rBhQ7vooousKMnKyrJ77rknkHX/8ccfbv2DBw+2TBPEtmldWqfWjcKj80HHIV3pGqNrTWHYsGGD1apVy9566y3LNEFfz0aPHu3Wr78omr+r8c6tIH9D87MnvvfII4+0m2++OdDvAIoqgm6giPADpB9++CHu/GOPPdZatmy529/z+eefF8oP+p7k32x+8MEHOaZv27bNTjrpJCtWrJi9+uqrhZa+MBk+fLh17drVatSoYaVKlbI6derYWWedZd9++62F2YMPPmgffvhhyterG3jl8UwxcOBAq1ixop1zzjmFnZTQGzlypF1yySV2wAEHWLly5Wy//fazyy67zJYsWVLkzwsUjd//W265xQYNGmRLly7lkCB0CLqBNDZv3jx76aWXkv7R7d+/f2Bp2rx5s915551W1Gzfvt3OPPNMt/3aZ7p5RHA8z7OLL77YTj/9dFu2bJldf/319vzzz1uvXr1s/vz51rFjR5swYUIoDoHOB50X0QguEjtnFXQrsCtevLhlmgYNGrh8ccEFF1g6UMCkB5qnnXaaPfXUU+5ByHvvvWetW7dOWRDFeRHsb2h+v/974rf71FNPtUqVKtmzzz4b6PcARVGJwk4AgIIrXbp0kdh92dnZrhS5TJky7lUUb95Vuvrpp5/aCy+8YJdeemlhJynjPf744672Rp8+feyJJ57IUb36jjvusDfeeMNKlAjHT5C2Myzbmko6X1esWOHO3Uykc6IoXi/zovP4X//6l6sp5DvhhBPsmGOOsWeeecbuv//+Qk1fUbRlyxZXwyd6n+1KYeWJPfG92g96+P3f//7XBf/p3OwGSBYl3UAGtT1TcKkfssaNG7sf0OrVq7ubpBEjRrj5WlZVu0Q/dv7Lt3HjRrvhhhusXr16LqBv0qSJPfbYY67UMpo+c80117h2li1atHDLfvnll5F5sdXX/v77bxfoqmqxlt13333t6quvdoG6rFq1ym688UY78MADrUKFCu5JuKokz5gxY7f30Y4dO1yJzEcffWTPPfecXX755XHb2/72229u/1SpUsUqV67sSmk3bdqUa1333Xef7b///m47tP9vv/1227p1a2QZlehqv0fvs2uvvdZ9h0qHfCr91TSlKT8///yzu0mpVq2aO6aHHnqoffzxx7mWmzNnjh133HFWtmxZq1u3rrsB1sOQWJqmbdaxUBXRDh062E8//RS3HeOaNWtc0Oznh0aNGtnDDz8cd72xJSYDBgywpk2buvwT78ZKpXuHH3545L1Kv//973+77VS61Pbvs88+i9tsQKVryuf77LOPq3qs/bN27Vp3HJRetQFWPtIxjD42sXlX+Vv79JBDDrGxY8fmSuOPP/7o8qHyo9an0vlJkyblWGZX51y8Nt36v861119/PXIORu97nS+qiVG7dm2333WOFbQ5hN9uWMfhxRdfjOTdww47zKZMmZJreVXtVTMWbYv+qnlAPMoDTz75pEubllVar7zySlu9enVkmX79+rmbbFVLjnbFFVe4QGRX57fSonypNMc7LxSM16xZ0+V5HUs9zPH9+eef9p///MdN13wdF+Wv2P4N/GY93333nTt3tb7y5cu70lwF/LFUQudf83QOqeaGzpN4TYF0Xun8Un5WXn3kkUcSatOdqm1LhJrgKA1jxozJNU8PKDVv9uzZ7n379u1zBY+apnN27ty5u/yuX3/91c444wzba6+9XJ7RdUrXZp27uzovgjieukbrOql0+NdCXUdjJfr75F+fhgwZ4kqMdcy13nXr1iV1bkX/hvp5JK+Xb9y4cW5/1K9f3+VNXbP79u2bo4bNrn7/4/12J3INTPYcOv74493xnD59etztBzIVj96BIkY3ICtXrsw1XTf3u6IfTAU7qo6pgEY/9mojPm3aNPdDp5vixYsXu4BAJY2xNyCnnHKKjRo1ygXIBx98sH311Vd20003uSDgf/7nf3Isrza5Cn4UwKi9bl4dLen7lBbdmOpmW4GY1qebPQW1uvlWwKUbEt00KCBXQKobPpWg6MZVN7cFoSD53HPPdTc3utnQ9udFN7n6bu0/7a+XX37ZBW8KMn3ar7ohVJCnhxOTJ092y+uG07+BOvroo92+0s2b3wZfN0S6WdXf6667LjLNv2nNi9bRrl07d/N26623uhsZ7fPu3bvb0KFD3U2NqGqnbhi1vf5yCrB0cxrrtttuczf/J598snXp0sXdOOqvSmSi6dho/+tYab/pZk7VwfV5teFUwJWX8ePHuxtVBcCJVAvW8T7qqKPcd2r/6IZa+1n5UfnE306f9rm2TduqhyVPP/20lSxZ0u1jBX06D3RjqJtBHdO77747x+cVYLz77rvuu3SDqkBKJXbff/995Jhp3+tY6mZTHf9o/cqTCqj0+SOOOCKhcy4enXv+8jonxA8stS/0wMF/OKAb2C+++MKdk1q39mlBvP3227Z+/Xp3LLVu5QFV/de5p22Tr7/+2gVFzZs3d9v0zz//uAcXCkpiaT3av5qv/bhgwQJX2qmbdN18a50KPD755BOX9lmzZrkHJLqmqHmHHl61atUq3zQrv7Vp0ybX9JkzZ7pjo+/Q/tO15/fff3ff9cADD7hl9EBBn1dQp/QreNEDLh0/XVMUDEXTg7GqVau6BwVaVvlb+1/5xKdjrQcsnTp1cg8N1bxH69R3+dvsUz5UntI+1rVF+VjVsxW4KYjJSxDblp8TTzzRBVO6ruh8j6Zt1wOG/PoSUUd3euk3ID96wKrrjB6CaV8r8Na1RbUZ9NugB535nRdBHE9dFxR0d+vWzb10znbu3DnyMNiX7O+T8rZ+1xSoa3v1/2TOrWg6/2N/q3UvoIBa6/W9//777vqpfKnrp65lui7+9ddfbp7k9/sfT6LXwGT2ueghp+icUdMEIDQ8AEXCa6+9pqLRfF8tWrTI8ZkGDRp4PXv2jLxv1aqVd+KJJ+b7Pb169XLrivXhhx+66ffff3+O6WeeeaaXlZXl/fbbb5FpWq5YsWLenDlzcq1H8/r16xd5f+GFF7plp0yZkmvZ7Oxs93fLli3ezp07c8xbsGCBV7p0ae/ee+/NMU3r177Kz6hRo9xy2j/6O2jQoDyXVVq1zCWXXJJj+mmnneZVr1498n769OluucsuuyzHcjfeeKOb/u2337r3y5cvd++fffZZ937NmjVu+//97397tWvXjnzuuuuu86pVqxbZB/G2rWPHjt6BBx7o9k/0PjvqqKO8xo0bR6b16dPHfXby5MmRaUpH5cqV3XStW5YuXeqVKFHC6969e45tuOeee9xy0Xnpvvvu88qXL+/98ssvOZa99dZbveLFi3sLFy7Mc58OHDjQrW/48OFeIvz0jxs3LjJt/fr13r777us1bNgwkjf849qyZUtv27ZtkWXPPfdcl0e7du2aY71t27Z1eSCafy798MMPkWl//vmnV6ZMGXfMfdpHpUqV8n7//ffItMWLF3sVK1b02rdvn9Q55+exaNq30fvbd+mll3p77723t3LlyhzTzznnHHc8N23alO93aXuj0+PnK+XlVatWRaZ/9NFHbvonn3wSmXbwwQe771ae9X399deRc8mn46Rpb731Vo7v/vLLL3NNnzVrltuPOm9Wr17t7bPPPt6hhx7qbd++Pd/t0Hwd0xtuuCHXPO1/HQcdt2j+uSTx9tPEiRNd+v773//muu526tQpx+f79u3r8rm/L3Q+aTs6d+6c41r1zDPPuM+/+uqrkWnHHHNMru/ZunWrt9dee3lnnHFGZFq8cz6V2+afL/qbH50/tWrV8nbs2BGZtmTJEnfdir7+xqPrhL5j5MiR+S73448/uuXef//9fJfL67wI6njqXIle7vbbb891LUz098nf3/vtt1+u9CZ6bsX7DY31n//8x22L/5uT1/4ZMGCAO4ei81Jev//xvjfRa2Ci+zya1nv11VfnuY1AJqJ6OVDEqERWT6JjXwcddNAuP6uq0Xo6rWp8BelgRaWSfkmsTyW6+j1WaVs0PeXXU/v8qAqqSghUqqpq0bH8qm0qbfSrLe7cudOVAqj0RVUJVfpQUCqRUFtalU7sylVXXZXjvZ7wKx1+1UDtH1H1udj9I35VaJVMqDTfr66sp/nar6oxoPT4x0Yl3aqGnFebNpUUqzaBSslUQqnaD3opTSox0npUUuSnTaWj0dW1lY4ePXrkWKeq+ao0XNU0Y0soYql0RPtAJRf+d+ulUj4do3jVsX3+PlPJZiKUfqVd+8On46/SLpWYqDQp2oUXXpijVFElLsqjsZ3jafqiRYvcNkdr27ZtpLRFVIqvDn5UCqtt00slU6pRoB6afXvvvbedd955riTf38bdOediaRtUg0Hni/4fvd91zFULpqDnw9lnn+2OpU/H1i/FE9VeUHXPnj17ulJHn0rrY89z5Q0to3nRadQ+1XFTbRmfSklVOqyaI9oGLadaDLtq4678r30QnWZRdVXlPR1rHbdo0edSdC0PlQzqvFHzCB2vePtQeS3689o/ygeqBivffPONKwFVTYPoKtZqrqKSwNimENoP559/fuS9SiWVx/39HU9Q25ZI3li+fHmOocVUMq/rt+blRWnVsdU1Sk1b8uPnKZ1jsc12EhHU8fSb/vji1SRJ9vdJ51B0epM5t3ZFbaFVM0c1VVS7Kd7+URV9nWeqPaRzSLVPkpXMNTDRfR7N/10BwoSgGyhidGOmwCb2FXvzGc+9997rquppSBdVY1Sgp+qKidAPo6rJxQZKzZo1i8yPlkggq5tI/TDvaqgz3dypSrbaxeoGR1UVFTQq7X57v4LQjYluXlUdXMFvfmJvcv397bdR1fbrxks3etFUTVI3ftH7RzcbfvVx/dUDB73U9lHvtU9UrdsPfOJRtWndMN11111uX0S/VH1PdKPsp037LpZuCqP5aYzdBqUrNn8piFQ7/djvVl6M/u54FISIHhYkQumKTWt+eS/2WPk3smrHGDtdeSs2D8XbVzpnFAwoz+ql/+eVJq1TwfzunnOx9L1al5oGxO53VUXd1X7f3fwtieQj5Q3tUzW/iE2nqhrHplH7RFXJVeVVeTeZQCO2Pwk/aN3VNUVtWVV92O+PwL+maP/Gu6Ykun9i94WCaQUlsXlU1YZjH6hpndFt3mMFtW27omrwOleiqwHr/2pipHwdj9qdq9mH0qoHKrui3ws9sNSySq8ewOgBc6LpDep4xuZ3rTP2Wpjs71Psb2My51Z+FLjr4bCaTMU+/F24cKFrs61ruR4IKH1+c4GC5IlkroGJ7vPY85pO1BA2tOkGMojaB6v9nzoN01Nq3eDoZkFDNamtXCrFay+8O8PEKLhUCY/aw+nGQQGuSh121WlXfvRUXrUEVIKqtotqh5ZXO9K82h7H60RuV/R9areqm2gF2Qqu9TlN13s93NB25Rd0+9utdoG6QY0nNnhOJX2/SmLUli+evG7GRSX9ona8KilJtbyOVaLHsKiec/4xVwmpSsXiSaTGS9D7RulUwK3O6OLRDX80nQd+TQDliUToGqBzJr8gNT8qwXzttdfcNUQ1GxRUan1qExzvmpLqvBNkXkx223ZFgaTOU/VLoVJU1cjRQ0pdl+NRsKW2z/pe1VJJtEaLRjRQYOifK6pVpfbN6n9hV22bC/N4Jvv7lMrfRp/OA7UJ13U39iGHSpN1rVbtEPUboOuv+vVQTSjt7935DU1GMvtcD0t21Q8AkGkIuoEMoxsClYrppVInBQXqAMgPAPIKGjVmrKrcqXQy+iZKJRr+/GTp5lulnn7vt3lRVUZVlXvllVdS/sOsUihVadRTfwWvCnrjlTjsirZfNy8KHvwSWNENqtIZvX/8YFoBvzoAUodfomOhzn8UdOumKLqKc7x0i6pR+6XL+aUtXvVmdfQUu5xfih5dGqPqkrHBjTowUv7Z1XfHo4cLKuV45513XO/uu+pMTemKTevu5r38xNtXv/zyi+uMyQ8Y9f+80qQb7uhS9V2dc/HEOw/13Tr3dBNdkP2+O/x9nEg+Ut7QtUKd/O0qwNA5oxt/XQcUpCiAUc0TdTCWH1U/1/eog7Z450Ui1xQ9uFCg51NngbE9jSe7f7Qvoqvbqoqy0piK41VY2yaqRq5q/2qCoo4hFSjFq1qua4UCbnUQpmX1YDMZqg2ilzrZU8doykN6QOUPN5bX71NQx1P5Pfp4qoQ39lq4u79PyZxbeZ1Daiqk79N5F9tpnB5k6fql46emN77oERR8iZYu61qUzDUwGXoYoPMm+ncUCAOqlwMZRDdE0VTNTKWh0cMmKdiT2JsV9d6qm331QBxNpXb6oc6vx9286IdZJSjqeVc9Ouf1BFxBWezTcLUb9dss7y7d5KnNpQIilQgUZL3aPxLba7fGrhWVpPsU0KrHce07tT/UjaUfjKtUVDdxaoOdX7tWlSSql1j1Fqs2gbGih2JR2lRapOq70fNjSyI13Iu+M3aYsthjLmqnOXHiRPfAIpbyTmw76Wi6WVOJi27e9TdeScebb74ZSa/Sr//r+6LbJaqatXpvTrbd467oe6LbYqrkTqVvCiaUF/XS/zUtekgiPWBRL+B6qOBXoU/knItH52HsOajvVWmW2nXHC7ziDb+TKgqeVJ1YN+7R1VF14x7bpl55Q9cKlfrFUr6I3i6dHwqudCy1vNqZqoflRNpzqkQz9rqhYEAPNTSEmqrURovOZ/GuKerNWekuCAXVqkquYf+i16tATPsr+vwvqMLaNn/79PBI1cr1UjOn2GrSOid1rur6qRLuZB5eqklN7DVD12X9RsT+PsULpIM4nnqgqXVErzfeqAy7+/uUzLkVj9rN6zqsh5jxmnX5DzWj06j/Dxw4MNeyef3+x1tnotfAZE2dOtX91bUACBNKuoEMouBEgZpKUHUDpRtWBXgatsPnl66qap9KfvXjqip66rxJT/M1Hqx+ZFUNW1UA9aOrEqp4Y+UmQiVbWo9KmtXRip5uK4jUTYs6Y1F76JNOOsm1jVVJoX6I9eReAWN0CcTu0g38sGHD3HYq8FaJt4ZWSZT2h0paFDzohkXbo0BRN1J6sBDdqY0fYGu8Vt1Y+m3bNPyRbnpUKqHOaHZFbR51c6N1qMMm7Q/d9Cho1FAw/jixqgKuIWDUNrN3796RIcNUwhLdvlhjKWu+Sos0HJeW1zrUSZ5KbKJLQdQOV+OB69iopFL5RjfdOjbKU8oj+ZXy6PPqYEzfpY61VLqp9u8a3kyd62nfKRgT1QTQDaUe7ChfKu9qv6oEUQFo7NjAu0vtUJX3o4cM829ufSp585smqOM5PazQAxAFCNHjLSdyzsWj5VVqpaBUNR90M62O3x566CG3v/R/HXOtX9VG9ZBAy+v/QVFVXwWP2mZVpdV3KSjRsFF6YOVT3tfwQ1pe7Ux1c64ARiV5Oq91s6/jrYcuqpar/KPzTjTMmAIQ7VMNU5UfdW6nfK3zJbo5gwJfpVHnk64p2nfKj3qw5o/9q3yrz6oasvahzhntv2TO+diAWMPlKY/ovNH5o1JA5R2NeR7dadruKIxtEx0/1T7QNUvnucZ1j6XSVp23yhs6ttFjc+thU35NSdQppM4JDbulY6kAXNvgP2ja1XkRxPFU0x3lYa1bDxPU4Zh/LYyWit+nRM+tWPouPazSwxj1laCHldGU71SdXL/P2h49CFAwrOtmvKYZef3+x5PoNTBZWqfafzNcGEKnsLtPB5Bz2I14Q2v5w9DsasgwDfd1+OGHe1WqVPHKli3rNW3a1HvggQdyDK+kYWGuvfZar2bNmm44kejLgIZp0jAfderU8UqWLOmGpXr00UdzDAMi+oyGHokn3nAnGrJEQ4fpOzXMioZU0ec1jI4/JIuGBtKQKkp3u3bt3HAw2ma9CjpkWLzhad599103FM5hhx3mrVu3LjKc04oVK+IeD3+4LX8Yo/79+7uhrLR/6tWr59122205hvTyaZgyfT52WBQNqxJviJ28tk3DtWjfabghfaeGXDrppJO8Dz74IMdyM2fOdPtKQ19pGQ3l88orr+TaBh3/u+66y61P+/q4447z5s6d64aUuuqqq3KsU/lB29eoUSM3xEuNGjXccGWPPfZYjjyVH6VTwyxpeDQNV6ZjfPbZZ3ujR4/OtZ0ank55V9ugfPzpp58mdFzzOnfiHVs/77755psufys/tm7dOu6wStOmTfO6dOniVahQwStXrpzXoUMHb8KECTmWSeScizdk2M8//+yG3dFnYocoWrZsmUuj8peOuY6Vho978cUXd7m/8xoyTOdxIufq0KFDvWbNmrn90rx5c2/YsGEubbHDGonSc8ghh7ht0DBCGt7u5ptvdsMKKZ/pHKtbt26uIYP8IeV0LuZH1wflOeXlWLNnz3ZDvPn5pUmTJi5f+zQ82cUXX+w+r+On46h9HnvNzCvv5DXcloYI0zHWcdEQgDq/9V27ulZL7H7M65xP1bYlOmSYb8SIEW55/S4sWrQo13x/CMZ4r3j5I9r8+fPdsIz777+/2yZdD3Q+ffPNNwmdF0EcTw0Dpuu5/7tz7LHHun0fu85Ef5/y+91J5tyKPi/9deb18v3000/ut0X7Rvvo8ssv92bMmJErf+X3+x/vepDINTDZfa79eOedd8bdR0Amy9I/hR34AwAKj0ruVRqvkg3VdMhkKs3v1atX3Cr1KHpUyqcOtFSKvqu+AQAUbarlpFpeamaVbH8AQLqjTTcAhIiG3onlt2NUNWmgKOnbt6+rfqtqzwDS28MPP+yaGRBwI4xo0w0AIaJOktSuVm0Y1Q5T7erVnlrtcv0O34CiQnm0oGOTAyhaojvLBMKGoBsAQkTjPKtDHHWEox6F/c7V/CF7AAAAkFq06QYAAAAAICC06QYAAAAAICAE3QAAAAAABCSt23RnZ2fb4sWLrWLFim4YGAAAAAAA9gSNvr1+/XqrU6eOFStWLDODbgXc9erVK+xkAAAAAABCatGiRVa3bt3MDLpVwu1vZKVKlQo7OUiDmhErVqywmjVr5vskCkh35HWEAfkcYUA+Rxhkp/E9ukaCUSGwH5dmZNDtVylXwE3QjURO6C1btri8km4nNJAM8jrCgHyOMCCfIwyyM+AefVdNndNzqwAAAAAASAME3QAAAAAABISgGwAAAACAgBB0AwAAAAAQEIJuAAAAAAACQtANAAAAAEBACLoBAAAAAAgIQTcAAAAAAAEh6AYAAAAAICAE3QAAAAAABISgGwAAAACAgBB0AwAAAAAQkBJBrRgAAAAAkNk2b9thwyYtsM+nLbRVG7ZYtQplrFub+nb6kfta2VKEm8JeAAAAAAAUKOC+8fWJ9vuydeZ5/ztt5fot9sbYX2zCvKX2WM+2BN5ULwcAAAAAFIRKuKMDbp/ea7rmgzbdAAAAAIACUJXy2IDbp+maD4JuAAAAAEABqA337swPC3ovBwAAAAAkTZ2m7c78sCDoBgAAAAAkTb2UZ2XFn6fpmg+CbgAAAABAAWhYsP1rV8oVeOu9pms+GDIMAAAAAFAAGodbw4IxTnf+GKcbAAAAAFDgwLtH+8buhfho0w0AAAAAQEAIugEAAAAACAhBNwAAAAAAASHoBgAAAAAgIATdAAAAAAAEhKAbAAAAAICAEHQDAAAAABAQgm4AAAAAAAJC0A0AAAAAQEAIugEAAAAACAhBNwAAAAAAASHoBgAAAAAgIATdAAAAAAAEhKAbAAAAAICAEHQDAAAAABAQgm4AAAAAAAJC0A0AAAAAQEAIugEAAAAACAhBNwAAAAAAASHoBgAAAAAgIATdAAAAAAAEhKAbAAAAAICAEHQDAAAAABAQgm4AAAAAAAJC0A0AAAAAQKYH3Q899JBlZWVZnz59CjspAAAAAABkTtA9ZcoUe+GFF+yggw4q7KQAAAAAAJA5QfeGDRusR48e9tJLL1nVqlULOzkAAAAAAGRO0N2rVy878cQTrVOnToWdFAAAAAAAUqqEFaIhQ4bYtGnTXPXyRGzdutW9fOvWrXN/s7Oz3QvIj/KI53nkFWQ88jrCgHyOMCCfIwyy0/gePdE0F1rQvWjRIuvdu7eNGDHCypQpk9BnBgwYYP379881fcWKFbZly5YAUolMopNi7dq17qQuVqzQK3kAgSGvIwzI5wgD8jnCIDuN79HXr1+f0HJZnrauEHz44Yd22mmnWfHixSPTdu7c6Xow185WiXb0vLxKuuvVq2erV6+2SpUq7dH0Iz1PaD2gqVmzZtqd0EAyyOsIA/I5woB8jjDITuN7dMWj6pdMDw3yi0cLraS7Y8eONmvWrBzTLr74YmvatKndcsstuQJuKV26tHvF0sFJtwOEwuE/1CG/INOR1xEG5HOEAfkcYZCVpvfoiaa30ILuihUrWsuWLXNMK1++vFWvXj3XdAAAAABAZtu8bYcNm7TAPp+20FZt2GLVKpSxbm3q2+lH7mtlSxVqd2S7JX1TDgAAAADImID75jcm2+/L1pnfAHrl+i32xthfbMK8pfZYz7ZpG3gXqVSPHj26sJMAAAAAANjDhk/+I0fA7dN7TVcJeI/2jdPyuKRXpXkAAAAAQMb54seFuQJun6arynm6IugGAAAAABSqVRu27mJ++g4RTdANAAAAAChU1SqU3sX8MpauCLoBAAAAAIWqa+v6lpUVf56mqxfzdEXQDQAAAAAoVKcd0dD2r10pV+Ct95quYcPSVZHqvRwAAAAAED5lS5Vww4IxTjcAAAAAAAEF3j3aN07bocHyQvVyAAAAAAACQtANAAAAAEBACLoBAAAAAAgIQTcAAAAAAAEh6AYAAAAAICAE3QAAAAAABISgGwAAAACAgBB0AwAAAAAQEIJuAAAAAAACQtANAAAAAEBACLoBAAAAAAgIQTcAAAAAAAEh6AYAAAAAICAE3QAAAAAABISgGwAAAACAgBB0AwAAAAAQEIJuAAAAAAACQtANAAAAAEBACLoBAAAAAAgIQTcAAAAAAAEh6AYAAAAAICAE3QAAAAAABISgGwAAAACAgBB0AwAAAAAQkBLJLDx37lwbMmSIjRs3zv7880/btGmT1axZ01q3bm1dunSxM844w0qXLh1UWgEAAAAAyLyS7mnTplmnTp1ccD1+/Hg74ogjrE+fPnbffffZ+eefb57n2R133GF16tSxhx9+2LZu3Rp8ygEAAAAAyISSbpVg33TTTfbBBx9YlSpV8lxu4sSJNnDgQHv88cft9ttvT2U6AQAAAADIzKD7l19+sZIlS+5yubZt27rX9u3bU5E2AAAAAAAyv3p5XgH3li1bkloeAAAAAIAwSbr38uzsbNeWe5999rEKFSrY/Pnz3fS77rrLXnnllSDSCAAAAABAOILu+++/3wYPHmyPPPKIlSpVKjK9ZcuW9vLLL6c6fQAAAAAAhCfo/u9//2svvvii9ejRw4oXLx6Z3qpVK/v5559TnT4AAAAAAMITdP/999/WqFGjuNXO6UANAAAAAIDdCLqbN29u48aNyzVdw4lpHG8AAAAAAJDEkGHR7r77buvZs6cr8Vbp9rBhw2zevHmu2vmnn36a7OoAAAAAAMhYSZd0n3rqqfbJJ5/YN998Y+XLl3dB+Ny5c920448/PphUAgAAAACQ6SXdO3bssAcffNAuueQSGzFiRHCpAgAAAAAgbCXdJUqUcEOFKfgGAAAAAAAprl7esWNHGzNmTLIfAwAAAAAgdJLuSK1r165266232qxZs+yQQw5x7bqjnXLKKalMHwAAAAAA4Qm6//Of/7i/TzzxRK55WVlZtnPnztSkDAAAAACAsAXdGiYMAAAAAAAE0KYbAAAAAAAEGHSrI7WTTz7ZGjVq5F5qxz1u3LiCrAoAAAAAgIyVdND95ptvWqdOnaxcuXJ23XXXuVfZsmVdr+Zvv/12MKkEAAAAACAMbbofeOABN1Z33759I9MUeKtjtfvuu8/OO++8VKcRAAAAAIBwlHTPnz/fVS2PpSrmCxYsSFW6AAAAAAAIX9Bdr149GzlyZK7p33zzjZsHAAAAAAAKWL38hhtucNXJp0+fbkcddZSb9t1339ngwYNt4MCBya4OAAAAAICMlXTQffXVV9tee+1ljz/+uL333ntuWrNmzezdd9+1U089NYg0AgAAAAAQjqBbTjvtNPcCAAAAAAApbNM9ZcoUmzx5cq7pmvbDDz8kuzoAAAAAADJW0kF3r169bNGiRbmm//33324eAAAAAAAoYND9008/WZs2bXJNb926tZsHAAAAAAAKGHSXLl3ali1blmv6kiVLrESJAjURBwAAAAAgIyUddHfu3Nluu+02W7t2bWTamjVr7Pbbb7fjjz8+1ekDAAAAACBtJV00/dhjj1n79u2tQYMGrkq5aMzu2rVr2xtvvBFEGgEAAAAACEfQvc8++9jMmTPtrbfeshkzZljZsmXt4osvtnPPPddKliwZTCoBAAAAAEhDBWqEXb58ebviiit2+8ufe+459/rjjz/c+xYtWtjdd99tXbt23e11AwAAAACQdm26X3/9dfvss88i72+++WarUqWKHXXUUfbnn38mta66devaQw89ZFOnTnVjfB933HF26qmn2pw5c5JNFgAAAAAA6R90P/jgg65KuUycONGeeeYZe+SRR6xGjRrWt2/fpNZ18sknW7du3axx48Z2wAEH2AMPPGAVKlSwSZMmJZssAAAAAADSv3r5okWLrFGjRu7/H374oZ155pmuqnm7du3s2GOPLXBCdu7cae+//75t3LjR2rZtG3eZrVu3updv3bp17m92drZ7AflRHvE8j7yCjEdeRxiQzxEG5HOEQXYa36Mnmuakg26VRP/zzz9Wv359+/rrr+36669308uUKWObN29OOqGzZs1yQfaWLVvcuocPH27NmzePu+yAAQOsf//+uaavWLHCfR7Y1Umhoe50UhcrlnQlDyBtkNcRBuRzhAH5HGGQncb36OvXrw8m6NZY3JdddpkbLuyXX35x1cNF7bAbNmyYdEKbNGnihhzTjv7ggw+sZ8+eNmbMmLiBt8YH94N8v6S7Xr16VrNmTatUqVLS343wndBZWVkuv6TbCQ0kg7yOMCCfIwzI5wiD7DS+R1fBcyBB96BBg+zOO+901cyHDh1q1atXd9PVGZqGDUtWqVKlItXVDznkEJsyZYoNHDjQXnjhhVzLli5d2r1i6eCk2wFC4dAJTX5BGJDXEQbkc4QB+RxhkJWm9+iJpjfpoFs9lavztFjxqn0X9ElHdLttAAAAAADS1W49SjjwwANdiXdBqbr42LFj3Tjdatut96NHj7YePXrsTrIAAAAAACgSki7pjqZgefv27QX+/PLly+3CCy+0JUuWWOXKle2ggw6yr776yrUbBwAAAAAg1EH37nrllVcK8+sBAAAAACi61cuPPvpoK1u2bOpSAwAAAABABtmtku7PP/88dSkBAAAAACDsJd3Fixe3Dh062KpVq3JMX7ZsmZsHAAAAAAAKGHR7nueG9Dr00ENtzpw5ueYBAAAAAIACBt0auHzo0KF28sknW9u2be2jjz7KMQ8AAAAAAOxGSbeqkQ8cONAee+wxO/vss+3++++nlBsAAAAAgFR2pHbFFVdY48aN7d///reNHTt2d1YFAAAAAEDGSbqku0GDBjk6TFOnapMmTbJFixalOm0AAAAAAISrpHvBggW5pjVq1Mh+/PFH14M5AAAAAAAoYEl3XsqUKeNKwQEAAAAAQIqDbgAAAAAAkBNBNwAAAAAAASHoBgAAAAAgIATdAAAAAAAUld7Ld+7caYMHD7aRI0fa8uXLLTs7O8f8b7/9NpXpAwAAAAAgPEF37969XdB94oknWsuWLS0rKyuYlAEAAAAAELage8iQIfbee+9Zt27dgkkRAAAAAABhbdNdqlQpa9SoUTCpAQAAAAAgzEH3DTfcYAMHDjTP84JJEQAAAAAAYa1ePn78eBs1apR98cUX1qJFCytZsmSO+cOGDUtl+gAAAAAACE/QXaVKFTvttNOCSQ0AAAAAAGEOul977bVgUgIAAAAAQNjbdAMAAAAAgIBKuuWDDz5ww4YtXLjQtm3blmPetGnTCrJKAAAAAAAyTtIl3U899ZRdfPHFVrt2bfvxxx/t8MMPt+rVq9v8+fOta9euwaQSAAAAAIAwBN3PPvusvfjii/b000+7MbtvvvlmGzFihF133XW2du3aYFIJAAAAAEAYgm5VKT/qqKPc/8uWLWvr1693/7/gggvsnXfeSX0KAQAAAAAIS9C911572apVq9z/69evb5MmTXL/X7BggXmel/oUAgAAAAAQlqD7uOOOs48//tj9X227+/bta8cff7ydffbZjN8NAAAAAMDu9F6u9tzZ2dnu/7169XKdqE2YMMFOOeUUu/LKK5NdHQAAAAAAGSvpoLtYsWLu5TvnnHPcCwAAAAAA7Gb1chk3bpydf/751rZtW/v777/dtDfeeMPGjx9fkNUBAAAAAJCRkg66hw4dal26dHE9l2uc7q1bt7rpGi7swQcfDCKNAAAAAACEI+i+//777fnnn7eXXnrJSpYsGZnerl07mzZtWqrTBwAAAABAeILuefPmWfv27XNNr1y5sq1ZsyZV6QIAAAAAIJzjdP/222+5pqs993777ZeqdAEAAAAAEL6g+/LLL7fevXvb5MmTLSsryxYvXmxvvfWW3XjjjXb11VcHk0oAAAAAAMIwZNitt97qxunu2LGjbdq0yVU1L126tAu6r7322mBSCQAAAABAGko66Fbp9h133GE33XSTq2a+YcMGa968uVWoUCGYFAIAAAAAEJag21eqVCkXbAMAAAAAgBQF3Vu2bLGnn37aRo0aZcuXL3dVzaMxbBgAAAAAAAUMui+99FL7+uuv7cwzz7TDDz/cVTcHAAAAAAApCLo//fRT+/zzz61du3bJfhQAAAAAgFBJesiwffbZxypWrBhMagAAAAAACHPQ/fjjj9stt9xif/75ZzApAgAAAAAgrNXLDz30UNeZ2n777WflypWzkiVL5pi/atWqVKYPAAAAAIDwBN3nnnuu/f333/bggw9a7dq16UgNAAAAAIBUBd0TJkywiRMnWqtWrZL9KAAAAAAAoZJ0m+6mTZva5s2bg0kNAAAAAABhDrofeughu+GGG2z06NH2zz//2Lp163K8AAAAAABAAauXn3DCCe5vx44dc0z3PM+17965c2eyqwQAAAAAICMlHXSPGjUqmJQAAAAAABD2oPuYY44JJiUAAAAAAISxTffChQuTWqmGFAMAAAAAIOwSCroPO+wwu/LKK23KlCl5LrN27Vp76aWXrGXLljZ06NBUphEAAAAAgMytXv7TTz/ZAw88YMcff7yVKVPGDjnkEKtTp477/+rVq938OXPmWJs2beyRRx6xbt26BZ9yAAAAAAAyoaS7evXq9sQTT9iSJUvsmWeescaNG9vKlSvt119/dfN79OhhU6dOtYkTJxJwAwAAAABQkI7UypYta2eeeaZ7AQAAAACAFPdeDgAAAABIX5u37bBhkxbY59MW2qoNW6xahTLWrU19O/3Ifa1sKULEVGOPAgAAAECIAu4bX59ovy9bZ573v9NWrt9ib4z9xSbMW2qP9WxL4F0YbboBAAAAAOlPJdzRAbdP7zVd85FaBN0AAAAAEBKqUh4bcPs0XfNRhINuL6+jBwAAAAAodGrDvTvzsQeC7osuusg2btyYa/off/xh7du3L0ASAAAAAAB7gjpN25352ANB94wZM+yggw5yY3L7Xn/9dWvVqpXVqFEjqXUNGDDADjvsMKtYsaLVqlXLunfvbvPmzUs2SQAAAACABKiX8qys+PM0XfNRyEH3999/b6effrode+yxdvvtt9tZZ51l11xzjT322GM2fPjwpNY1ZswY69Wrl02aNMlGjBhh27dvt86dO8ctSQcAAAAA7B4NC7Z/7Uq5Am+913TNRyEPGVayZEl79NFHrVy5cnbfffdZiRIlXPDctm3bpL/8yy+/zPF+8ODBrsR76tSpVFUHAAAAgBTTONwaFoxxuotw0K3S6FtvvdUGDRpkt912m40fP96VfL/yyivWrVu33UrM2rVr3d9q1art1noAAAAAAHkH3j3aN3YvFMGg+9BDD7VNmzbZ6NGj7cgjj3Q9lj/yyCMu8L7kkkvs2WefLVBCsrOzrU+fPtauXTtr2bJl3GW2bt3qXr5169ZFPqsXsKs8pvxKXkGmI68jDMjnCAPyOcIgO43v0RNNc4GC7qeeesrKly/v3mdlZdktt9zi2mJfcMEFVlBq2z179mxXcp5fx2v9+/fPNX3FihW2ZQtd22PXJ4VqU+ikLlaMIeqRucjrCAPyOcKAfI4wyE7je/T169cntFyWl8LBtVUKXbp06aQ/p47YPvroIxs7dqztu2/eDffjlXTXq1fPVq9ebZUqVSpwuhGeE1oPaGrWrJl2JzSQDPI6woB8jjAgnyMMstP4Hl3xaNWqVd1Dg/zi0aRLuuWNN96w559/3hYsWOCGDmvQoIE9+eSTLmA+9dRTE16P4v1rr73W9Xqu6ur5BdyigD5eUK+Dk24HCIVDNTPILwgD8jrCgHyOMCCfIwyy0vQePdH0Jr1Vzz33nF1//fWu07Q1a9bYzp073fQqVaq4wDvZKuVvvvmmvf32226s7qVLl7rX5s2bk00WAAAAAABFTtJB99NPP20vvfSS3XHHHVa8ePEcbb1nzZqVdACvoniN+b333ntHXu+++26yyQIAAAAAoMhJunq5qpS3bt0613RV+964cWNS60phc3IAAAAAANK/pFvtrqdPn55r+pdffmnNmjVLVboAAAAAAAhfSbfac6sttoboUkn1999/b++8844bzuvll18OJpUAAAAAAIQh6L7sssusbNmyduedd9qmTZvsvPPOszp16tjAgQPtnHPOCSaVAAAAAACkoQINGdajRw/3UtC9YcMGq1WrVupTBgAAAABAGINuX7ly5dwLAAAAAAAUMOhWb+UasDwR06ZNS2g5AAAAAAAyXUJBd/fu3SP/Vwdqzz77rDVv3tzatm3rpk2aNMnmzJlj//nPf4JLKQAAAAAAmRh09+vXL0dHatddd53dd999uZZZtGhR6lMIAAAAAEBYxul+//337cILL8w1/fzzz7ehQ4emKl0AAAAAAIQv6NZwYd99912u6ZpWpkyZVKULAAAAAIDw9V7ep08fu/rqq12HaYcffribNnnyZHv11VftrrvuCiKNAAAAAACEI+i+9dZbbb/99rOBAwfam2++6aY1a9bMXnvtNTvrrLOCSCMAAAAAAOEZp1vBNQE2AAAAAAABBN2ybds2W758uWVnZ+eYXr9+/YKuEgAAAACAcAfdv/76q11yySU2YcKEHNM9z7OsrCzbuXNnKtMHAAAAAEB4gu6LLrrISpQoYZ9++qntvffeLtAGAAAAAAApCLqnT59uU6dOtaZNmyb7UQAAAAAAQiXpcbqbN29uK1euDCY1AAAAAACEOeh++OGH7eabb7bRo0fbP//8Y+vWrcvxAgAAAAAABaxe3qlTJ/e3Y8eOOabTkRoAAAAAALsZdI8aNSrZjwAAAAAAEEpJB93HHHNMMCkBAAAAACCsQffMmTMTWu6ggw7anfQAAAAAABC+oPvggw92Y3Kr7XZeNH/nzp2pShsAAAAAAOEIuhcsWBBsSgAAAAAACGvQ3aBBg2BTAgAAAABA2MfpBgAAAAAAiSHoBgAAAAAgIATdAAAAAAAEhKAbAAAAAICiFHTv2LHDvvnmG3vhhRds/fr1btrixYttw4YNqU4fAAAAAACZ33u5788//7QTTjjBFi5caFu3brXjjz/eKlasaA8//LB7//zzzweTUgAAAAAAMr2ku3fv3nbooYfa6tWrrWzZspHpp512mo0cOTLV6QMAAAAAIDwl3ePGjbMJEyZYqVKlckxv2LCh/f3336lMGwAAAAAA4Srpzs7Otp07d+aa/tdff7lq5gAAAAAAoIBBd+fOne3JJ5+MvM/KynIdqPXr18+6deuW7OoAAAAAAMhYSVcvf/zxx61Lly7WvHlz27Jli5133nn266+/Wo0aNeydd94JJpUAAAAAAIQh6K5bt67NmDHD3n33XfdXpdyXXnqp9ejRI0fHagAAAAAAhF3SQffYsWPtqKOOckG2XtFjd2te+/btU51GAAAAAADC0aa7Q4cOtmrVqlzT165d6+YBAAAAAIACBt2e57nO02L9888/Vr58+WRXBwAAAABAxkq4evnpp5/u/irgvuiii6x06dKReRpCbObMma7aOQAAAAAASDLorly5cqSkW+NxR3eaVqpUKTvyyCPt8ssvT3R1AAAAAABkvISD7tdee839bdiwod14441UJQcAAAAAINW9l/fr1y/ZjwAAAAAAEEoJBd1t2rSxkSNHWtWqVa1169ZxO1LzTZs2LZXpAwAAAAAgs4PuU089NdJxWvfu3YNOEwAAAAAA4Qm6o6uUU70cAAAAAICAxuletGiR/fXXX5H333//vfXp08defPHFZFcFAAAAAEBGSzroPu+882zUqFHu/0uXLrVOnTq5wPuOO+6we++9N4g0AgAAAAAQjqB79uzZdvjhh7v/v/fee3bggQfahAkT7K233rLBgwcHkUYAAAAAAMIRdG/fvj3Sqdo333xjp5xyivt/06ZNbcmSJalPIQAAAAAAYQm6W7RoYc8//7yNGzfORowYYSeccIKbvnjxYqtevXoQaQQAAAAAIBxB98MPP2wvvPCCHXvssXbuuedaq1at3PSPP/44Uu0cAAAAAAAkOGRYNAXbK1eutHXr1lnVqlUj06+44gorV64c+xQAAAAAgIIG3VK8eHHbsWOHjR8/3r1v0qSJNWzYsCCrAgAAAAAgYyVdvXzjxo12ySWX2N57723t27d3rzp16till15qmzZtCiaVAAAAAACEIei+/vrrbcyYMfbJJ5/YmjVr3Oujjz5y02644YZgUgkAAAAAQBiqlw8dOtQ++OAD17bb161bNytbtqydddZZ9txzz6U6jQAAAAAAhKOkW1XIa9eunWt6rVq1qF4OAAAAAMDuBN1t27a1fv362ZYtWyLTNm/ebP3793fzAAAAAABAAauXDxw40Lp06WJ169aNjNE9Y8YMK1OmjH311VfJrg4AAAAAgIyVdNDdsmVL+/XXX+2tt96yn3/+2U0799xzrUePHq5dNwAAAAAA2I1xusuVK2eXX355QT4KAAAAAEBoFCjonjdvnj399NM2d+5c975Zs2Z2zTXXWNOmTVOdPgAAAAAAwtORmoYMUxXzqVOnujbdek2bNs0OPPBANy8ZY8eOtZNPPtnq1KljWVlZ9uGHHyabHAAAAAAAMqek++abb7bbbrvN7r333hzT1aO55p1xxhkJr2vjxo0uaL/kkkvs9NNPTzYpAAAAAABkVtC9ZMkSu/DCC3NNP//88+3RRx9Nal1du3Z1LwAAAAAAMlHSQfexxx5r48aNs0aNGuWYPn78eDv66KMtSFu3bnUv37p169zf7Oxs9wLyozzieR55BRmPvI4wIJ8jDMjnCIPsNL5HTzTNSQfdp5xyit1yyy2uTfeRRx7ppk2aNMnef/9969+/v3388cc5lk2lAQMGuO+ItWLFCtuyZUtKvwuZRyfF2rVr3UldrFjS3RkAaYO8jjAgnyMMyOcIg+w0vkdfv359Qstledq6JCS6I9Qx2s6dOxNer5YfPny4de/ePamS7nr16tnq1autUqVKCX8XwntC6wFNzZo10+6EBpJBXkcYkM8RBuRzhEF2Gt+jKx6tWrWqe2iQXzyadEl3YRb7ly5d2r1i6eCk2wFC4dDDHfILwoC8jjAgnyMMyOcIg6w0vUdPNL3ptVUAAAAAAKSRpEu6ZcqUKTZq1Chbvnx5rpLvJ554IuH1bNiwwX777bfI+wULFtj06dOtWrVqVr9+/YIkDQAAAACA9A26H3zwQbvzzjutSZMmVrt2bVcVwBf9/0T88MMP1qFDh8j766+/3v3t2bOnDR48ONmkAQAAAACQ3kH3wIED7dVXX7WLLrpot79cw48l2Y8bAAAAAABpo1hBGou3a9cumNQAAAAAABDmoLtv3742aNCgYFIDAAAAAECYq5ffeOONduKJJ9r+++9vzZs3t5IlS+aYP2zYsFSmDwAAAACA8ATd1113neu5XB2gVa9ePenO0wAAAAAACIukg+7XX3/dhg4d6kq7AQAAAABACtt0awxtVS0HAAAAAAApDrrvuece69evn23atCnZjwIAAAAAECpJVy9/6qmn7Pfff7fatWtbw4YNc3WkNm3atFSmDwAAAACA8ATd3bt3DyYlAAAAAACEPehW1XIAAAAAABBA0O2bOnWqzZ071/2/RYsW1rp164KuCgAAAACAjJR00L18+XI755xzbPTo0ValShU3bc2aNW7c7iFDhljNmjWDSCcAAAAAAJnfe/m1115r69evtzlz5tiqVavca/bs2bZu3Tq77rrrgkklAAAAAABhKOn+8ssv7ZtvvrFmzZpFpjVv3twGDRpknTt3TnX6AAAAAAAIT0l3dnZ2rmHCRNM0DwAAAAAAFDDoPu6446x37962ePHiyLS///7b+vbtax07dkx2dQAAAAAAZKykg+5nnnnGtd9u2LCh7b///u617777umlPP/10MKkEAAAAACAMbbrr1atn06ZNc+26f/75ZzdN7bs7deoURPoAAAAAAAjXON1ZWVl2/PHHuxcAAAAAANjN6uXffvut66Vc1chjrV271lq0aGHjxo1LdHUAAAAAgCRs3rbD3hr7q/V4cqR1vf8z91fvNR0ZEHQ/+eSTdvnll1ulSpVyzatcubJdeeWV9sQTT6Q6fQAAAAAQegqsb3x9or0x9hdbuX6LZXvm/uq9phN4Z0DQPWPGDDvhhBPynK8xuqdOnZqqdAEAAAAA/s+wSQvs92XrzPNy7hK913TNR5oH3cuWLYs7PrevRIkStmLFilSlCwAAAADwfz6ftjBXwO3TdM1Hmgfd++yzj82ePTvP+TNnzrS99947VekCAAAAAPyfVRu27NZ8pEHQ3a1bN7vrrrtsy5bcB3Pz5s3Wr18/O+mkk1KdPgAAAAAIvWoVyuzWfKTBkGF33nmnDRs2zA444AC75pprrEmTJm66xuoeNGiQ7dy50+64444g0woAAAAAodStTX3XaVq8KuZZWf87H2kedNeuXdsmTJhgV199td12223m/d/R1pjdXbp0cYG3lgEAAAAApNbpR+5rE+YtzdWZmgLu/WtXcvOR5kG3NGjQwD7//HNbvXq1/fbbby7wbty4sVWtWjW4FAIAAABAyJUtVcIe69nW9VKuTtPUhltVylXCrYBb81E0FejIKMg+7LDDUp8aAAAAAEBcCqx7tG/sXsjAjtQAAAAAAEByCLoBAAAAAAgIQTcAAAAAAAEh6AYAAAAAICAE3QAAAAAABISgGwAAAACAgBB0AwAAAAAQEIJuAAAAAAACQtANAAAAAEBACLoBAAAAAAgIQTcAAAAAAAEh6AYAAAAAICAE3QAAAAAABISgGwAAAACAgBB0AwAAAAAQEIJuAAAAAAACQtANAAAAAEBASgS1YgAAAAAIu83bdtiwSQvs82kLbdWGLVatQhnr1qa+nX7kvla2FOFYGHCUAQAAACCggPvG1yfa78vWmef977SV67fYG2N/sQnzltpjPdsSeIcA1csBAAAAIAAq4Y4OuH16r+maj8xH0A0AAAAAAVCV8tiA26fpmo/MR9ANAAAAAAFQG+7dmY/MQNANAAAAAAFQp2m7Mx+ZgaAbAAAAAAKgXsqzsuLP03TNR+Yj6AYAAACAAGhYsP1rV8oVeOu9pms+Mh9DhgEAAABAADQOt4YFY5zucCPoBgAAAIAAA+8e7Ru7F8KJoBsAAAAAYmzetoMSaqQEQTcAAAAAxATcN74+0X5fti4yzvbK9VvsjbG/2IR5S12VcZVgA4mgIzUAAAAAiKI22NEBt0/vNV3zgUQRdAMAAABAlM+nLcwVcPs0XfOBRBF0AwAAAECUVRu27NZ8IBpBNwAAAABEqVahzG7NB6IRdAMAAABAlG5t6ltWVvxdoumaDySKoBsAAAAAopx+5L62f+1KuQJvvdd0zQfSKugeNGiQNWzY0MqUKWNHHHGEff/994WdJAAAAAAhpeHANCzYBe0PsBoVy1ixLHN/9Z7hwpCsQh9c7t1337Xrr7/enn/+eRdwP/nkk9alSxebN2+e1apVq7CTBwAAACCkgXeP9o3dC0jroPuJJ56wyy+/3C6++GL3XsH3Z599Zq+++qrdeuuthZ08AAAAAGli87YdbgxtDemlHsbV4ZnaX6s6uIJoIHTVy7dt22ZTp061Tp06/f8EFSvm3k+cOLEwkwYAAAAgzQLuG1+faG+M/cVWrt9i2Z65v3qv6ZoPFIZCfdyzcuVK27lzp9WuXTvHdL3/+eefcy2/detW9/KtW7fO/c3OznYvID/KI57nkVeQ8cjrCAPyOcKAfJ6coRPn2+/L1pnn5Zyu95qu+ecd3SiVhwghz+fZCaY5repYDBgwwPr3759r+jvvvGNly5a14447znXCtmHDBqtataq1bNnSxo0b55Zp2rSp2ym//PKLe3/MMcfY9OnTbe3atVapUiVr06aNjR492s1r3LixlShRwubOneve/+tf/3L//+eff6x8+fJ25JFH2siRI928/fbbz8qVK2ezZ89279u2bWu///67LV++3KXp6KOPtq+//trNU2dxlStXthkzZrj3asP+559/2tKlS61kyZLWsWNH++qrr1ymq1u3rmvTPm3aNLfsIYccYsuWLbO//vrLihcvbscff7xLw/bt223vvfd2y0+ZMsUte/DBB9uqVats4cKF7v0JJ5xgo0aNcg8s9EBDafZrEhx00EFuf82fP9+992sZbNy40apXr+7223fffefmNW/e3NVO+O2339z7Dh062A8//GDr16+3KlWquHWNHTvWzWvSpIn7q7b50r59e5s5c6atWbPGKlasaIceeqhLkzRq1MhKlSplP/30k3vfrl0799DF39/ap998801kf1eoUMGty9/fSrv2TenSpV2avvzySzevfv36Vq1aNXecRd+ptK9evdp9n/b3iBEj3IMf7T/tG9W8EOUHHUPt76ysLNfPgL+/99prL2vQoIFNnjzZLduqVSuXj/744w/3vnPnzi7fbd682R3D/fffP7K/lSc3bdoU2d9Kw6RJkyL7u1mzZjZ+/Hg3T//fsWOH/frrr+79scce6/KDHjYpH+k4jxkzxs074IADXC0R/2GV8p3ypLZV++vwww+3b7/9NrK/ta/mzJkT2d86L1asWOHyst5rv/h5VueHv7+j86y/v/08q/2tbfjxxx8j+3vx4sXu5edZf3/XqVPHvZR/pHXr1u54K8/6+9vPs7H7W/lM+8Df31qv8qj2a82aNd2+8PNsixYt3Dr8PBuWa4TywmGHHebSoDzENSKxa4T2mfbfkiVLIvuba0TRvUbomqZrZ8296tgf26ra6HkrbfPOYlapTHFrVbu41d6x2EoW4xrBfUR6XyNq1Kjh1qs0aT3cR+R/jfhg0gbzvOIWj65DH3z3szWvuIH7iCIWa4wfP97dr2hZfV86xRpabyKyPO31QqIATj/gH3zwgXXv3j0yvWfPnm6HffTRR7ss6a5Xr54LLPSjD+RHAZVuGhWYKSgBMhV5HWHJ54sWL7VHv1xg85fnLNnyh/R55IIjaMOJtMb1PDknPviFq1KeF/VA/tntXXf3sCDFstP4Hl3xqApy/EKaIlnSrScDeqqipyh+0K2drvfXXHNNruX1dEGvWDo46XaAUDj0lJj8gjAgryMMvp69IlfAHV2V9MPv/6TXYaQ9rueJU6dpasOd33xihqIpK03v0RNNb6FvlYYLe+mll+z111931TOvvvpqV93A780cAAAgnlE/r8gVcPs0Xb0XAwgP9VKumi7xaLrmA4Wh0Nt0n3322a46wd133+3agKmNgOrJx3auBgAAEG3Npu357hANFwQgPMN8afkJ85bm6kzNb3Ki+UAog25RVfJ41ckBAADyUqVcSVu9Me/AWzfuAIp2sOwP8xUdKPvDfCmAfqxn24TXpeW0PON0o6gpEkE3AABAsjo0rWnDpy2OW8WcqqRAsFIVLCtAzm+YL83v0b5xwunSd2r5ZD4DBK3Q23QDAAAURJcDa7kqo7FtOKlKCgQvkWA5ESolp28GZDqCbgAAkJbKlCzuhgW7oP0BVqNiGTcckP7qfTJVUgEkL1XB8q76XqBvBmQCfo0AAEDaoiopUDhSFSwnMswXkO4IugEAAICQSFVP4akKlvXdagdO3wzIZFQvBwAAAELU+ZmCXAXM2d7/7/xM0zV/T4+JrWCfvhmQ6Qi6AQAAgBBIVednqQyW/WG+6JsBmYzq5QAAAEAIJNL5WaJDbaVyTGz6ZkCmI+gGAADI8Pa3QBA9hRMsA4nhag0AAFCE299GVwf2299OmLc06WHRCODTl47dh9//WWQ6PwOQHIJuAACAIhiUJtL+NtGqwKkM4FO5n4raPi+Kadqyfafd98bklBw7egoHCgdBNwAACL2iWKqcyva3qQrgUx28F7V9XhTT9NWs5Sl7+KLv1XbEri/Zzs8AJIfeywEAQOilslfnVA3LlMr2t4kE8Ht6PxXFfV4U0zTq5xUpOXZCT+FA4SDoBrDbdOPw1thfrceTI63r/Z+5v3qfzHifAFCYUhWUpjJw21X72mTa36YqgE/lfiqK+7wopmnNpu2BdH72Vp+O9sWdJ7q/ek/HfEBwCLoDRjCCTJeqJ/kAUJiKYqmyqiHHjoHs03TNT1SqAvhU7qeiuM+LYpqqlCuZ73w6PwOKPoLuABGMIAxSWRUPAApLUSxVVvtatbONDbwL0v42VQF8KvdTUdznRTFNHZrWTNnDFwCFg6A7QAQjCINUVsUDEB5FrSZYUSxVTmX721QF8KncT0VxnxfFNHU5sFbKHr4AKBz0Xh6gVPY6ChRVqayKB6S7ojiUUlFNUyp7iE6FVPbqnMphmfz2t7t7v+AH8Lt7/FK5n4riPi+KaSpTsrg9csERKRmnG0DhyPK8vMLCom/dunVWuXJlW7t2rVWqVMmKGj25V/vWvOiJtTqwwJ6RnZ1ty5cvt1q1almxYlTySBWVTulmOS8qlVEnLdhzwpLXi1owGS+QjL5Z392hlAqyrqKYJlGJdn7BiEpydxVkBpHPg3yoUNB9VRQVtXMviLxeVNIUlus5wi07jfN5ovEoQXeACEaKlnQ+oYuyVNw8I5i8XrFKtZSVjBS1UteiGEym8lxI1bqKYppS9ftY1K/pqQxMkb77fHfTVNTzOZAK2SEIurnqByiV1cuAoiqVVfGQOlu277T73pickuq7qaoKnMoqxYn0mZFoAJiqdaWySVGq1lUU0xSWZimpqhaO9N7nRTFNAPa89HqUkGZS2esoUFSlsqMfpM5Xs5anrFf5VHUKmcrOJYviWMFFcSilopimVPcQDQBAUUfQHSCCEYSF/yRf1UHVT4H+6j0Bd+EZ9fOKIheUpjJQLorBZFEcSqkopinVPUQDAFDUEXQHjGAEQGFYs2l7Rpe6FsVgsigOpVQU0yTUBAMAhAlBNwBkoCrlSmZ0qWtRDCZTGUimal1FMU1CTTAAQJjQ2BIAMlCHpjVt+LTFKenIMVWdQqayc8miOFZwqsZBTuW6imKaotdHB1MAgDBgyDCERjoPRwAkm9cX/r3EHv1yfpEaUivVYxcXxbGCsedwTUcYkM8RBtlpfI/OON1ABp3QQDLCME43wDUdYUA+Rxhkp/E9OkE3kEEnNJAM8jrCgHyOMCCfIwyyQxB0p9dWAQAAAACQRgi6AQAAAAAICEE3AAAAAAABIegGAAAAACAgBN0AAAAAAASEoBsAAAAAgIAQdAMAAAAAEBCCbgAAAAAAAkLQDQAAAABAQAi6AQAAAAAICEE3AAAAAAABIegGAAAAACAgJSyNeZ7n/q5bt66wk4I0kJ2dbevXr7cyZcpYsWI8b0LmIq8jDMjnCAPyOcIgO43v0f041I9LMzLo1sGRevXqFXZSAAAAAAAhtH79eqtcuXKe87O8XYXlRfypyOLFi61ixYqWlZVV2MlBGjyJ0gOaRYsWWaVKlQo7OUBgyOsIA/I5woB8jjBYl8b36AqlFXDXqVMn31L6tC7p1obVrVu3sJOBNKOTOd1OaKAgyOsIA/I5woB8jjColKb36PmVcPvSq9I8AAAAAABphKAbAAAAAICAEHQjNEqXLm39+vVzf4FMRl5HGJDPEQbkc4RB6RDco6d1R2oAAAAAABRllHQDAAAAABAQgm4AAAAAAAJC0A0AAAAAQEAIupHWHnroIcvKyrI+ffpEpm3ZssV69epl1atXtwoVKtgZZ5xhy5Yty/G5hQsX2oknnmjlypWzWrVq2U033WQ7duzIsczo0aOtTZs2rlOHRo0a2eDBg/fYdiHcnnvuOTvooIMi41W2bdvWvvjii8h88jgywYABA+ywww6zihUruutw9+7dbd68eTmWIa8jE4wdO9ZOPvlkq1Onjrtn+fDDD3PMV/dKd999t+29995WtmxZ69Spk/366685llm1apX16NHD/SZUqVLFLr30UtuwYUOOZWbOnGlHH320lSlTxurVq2ePPPLIHtk+YHcNGjTIGjZs6PLuEUccYd9//33G7VSCbqStKVOm2AsvvOCCk2h9+/a1Tz75xN5//30bM2aMLV682E4//fTI/J07d7qAe9u2bTZhwgR7/fXXXUCtHzzfggUL3DIdOnSw6dOnu6D+sssus6+++mqPbiPCqW7duu6B0tSpU+2HH36w4447zk499VSbM2eOm08eRybQ9VkPSCdNmmQjRoyw7du3W+fOnW3jxo2RZcjryATK061atXKBRTwKjp966il7/vnnbfLkyVa+fHnr0qWLe+jkU8Ct3wCdK59++qkL5K+44orI/HXr1rnzp0GDBu6349FHH7V77rnHXnzxxT2yjUBBvfvuu3b99de73sunTZvmzhXl/+XLl2fWTlXv5UC6Wb9+vde4cWNvxIgR3jHHHOP17t3bTV+zZo1XsmRJ7/33348sO3fuXPXQ702cONG9//zzz71ixYp5S5cujSzz3HPPeZUqVfK2bt3q3t98881eixYtcnzn2Wef7XXp0mUPbSGQU9WqVb2XX36ZPI6MtXz5cnetHjNmjHvP9RyZSHl8+PDhkffZ2dneXnvt5T366KORacr7pUuX9t555x33/qeffnKfmzJlSmSZL774wsvKyvL+/vtv9/7ZZ591vxP+fYzccsstXpMmTfbQlgEFc/jhh3u9evWKvN+5c6dXp04db8CAARm1SynpRlpS6YhKolUFK5qe7qq0JHp606ZNrX79+jZx4kT3Xn8PPPBAq127dmQZPVHTU2K/JFHLxK5by/jrAPYU1cwYMmSIKylRNXPyODLV2rVr3d9q1aq5v+R1hIFq1i1dujTHPUflypVdFdvo+xZVKT/00EMjy2j5YsWKuZJxf5n27dtbqVKlcty3qMnG6tWr9+g2AYnatm2bu9ZH53/la73PtHvuEoWdACBZCkBU/UTVy2Pph0s/OPpxiqYAW/P8ZaIDbn++Py+/ZRSYb9682bW5AoI0a9YsF2SreqH6Jhg+fLg1b97cNXcgjyPTZGdnu2Y87dq1s5YtW7ppXM8RBv59R7x7juh7EvV7EK1EiRLuAVX0Mvvuu2+udfjzqlatGuh2AAWxcuVKV7gQL////PPPGbVTCbqRVhYtWmS9e/d2bZrU2QKQqZo0aeICbJX+ffDBB9azZ0/XBhbI1NpLs2fPtvHjxxd2UgAASDmqlyOtqAqKOlZQr+J6yquXAhF1QKL/68mYqqqsWbMmx+fUe/lee+3l/q+/sb2Z++93tYx6DaWUG3uCSrPVa/4hhxzienlWxyIDBw50eZM8jkxyzTXXuI6hRo0a5ToR9JHXEQb+fUe8e47oe5LYTqU04op6NE/m3gYoamrUqGHFixfPN/9nCoJupJWOHTu6arcqAfRfauOkXj39/5csWdJGjhwZ+YzaM2mIMFXVFf3VOqJ/wFRyroBa1Xf9ZaLX4S/jrwMojOq3W7dudUE4eRyZQH1KKeBW04lvv/02V9VY8jrCQPlewUX0PYeasqmtdvR9iwoTVPDg0zmj3wW1/faXUY/m6tcm+r5FtaaoWo6iXMBwyCGH5Mj/ytd6n3H33IXdkxuwu6J7L5errrrKq1+/vvftt996P/zwg9e2bVv38u3YscNr2bKl17lzZ2/69Onel19+6dWsWdO77bbbIsvMnz/fK1eunHfTTTe53s8HDRrkFS9e3C0LBO3WW291PTgvWLDAmzlzpnuvXmq//vpr8jgyxtVXX+1VrlzZGz16tLdkyZLIa9OmTZFluJ4jU0Zc+fHHH91Lt95PPPGE+/+ff/7p5j/00ENelSpVvI8++shd80899VRv33339TZv3hxZxwknnOC1bt3amzx5sjd+/Hg3gsu5556bo8fz2rVrexdccIE3e/Zsb8iQIe4+5oUXXiiUbQYSNWTIENdb/+DBg11P/VdccYU7H6JHGcoEBN3IuKBbP1L/+c9/3NAZ+sE57bTT3I1ctD/++MPr2rWrV7ZsWa9GjRreDTfc4G3fvj3HMqNGjfIOPvhgr1SpUt5+++3nvfbaa3tsmxBul1xyidegQQOX9/RAqGPHjpGAW8jjyAQKPuK9oq+15HVkAt1PxMvrPXv2jAwbdtddd7mgWcGHrvnz5s3LsY5//vnHBdkVKlRwQ5xefPHFLpiPNmPGDO9f//qXW8c+++zjgnkgHTz99NOuwEz3PRpCbNKkSV6mydI/hV3aDgAAAABAJqJNNwAAAAAAASHoBgAAAAAgIATdAAAAAAAEhKAbAAAAAICAEHQDAAAAABAQgm4AAAAAAAJC0A0AAAAAQEAIugEAAAAACAhBNwAAAAAAASHoBgBgD/E8z6644gqrVq2aZWVl2fTp09n3AWjfvr29/fbbRXLfHnnkkTZ06NDCTgYAYA8i6AYAFGoQ2qlTJ+vSpUuuec8++6xVqVLF/vrrLytsxx57rAuS/Vft2rXt3//+t/35559JrefLL7+0wYMH26effmpLliyxli1bWmFtR5kyZax58+ZuPxd12mfKC4n4+OOPbdmyZXbOOee496tWrbJrr73WmjRpYmXLlrX69evbddddZ2vXrs3xuYULF9qJJ55o5cqVs1q1atlNN91kO3bsiMwfNmyYHX/88VazZk2rVKmStW3b1r766qtc3//333/b+eefb9WrV3ffd+CBB9oPP/wQmX/nnXfarbfeatnZ2buxRwAA6YSgGwBQaBT8vfbaazZ58mR74YUXItMXLFhgN998sz399NNWt27dlH7n9u3bC/S5yy+/3AXKixcvto8++sgWLVrkgqtk/P7777b33nvbUUcdZXvttZeVKFGiQA8qooPBgm7HTz/9ZGeddZb16tXL3nnnnbjLbtu2zdLNU089ZRdffLEVK/a/tzg6Xno99thjNnv2bBfA6+HHpZdeGvnMzp07XcCt7Z0wYYK9/vrrbrm77747sszYsWNd0P3555/b1KlTrUOHDnbyySfbjz/+GFlm9erV1q5dOytZsqR98cUXbh8//vjjVrVq1cgyXbt2tfXr17v5AICQ8AAAKGSDBw/2KlSo4M2fP9/Lzs72OnTo4J122mnerFmzvBNOOMErX768V6tWLe/888/3VqxYEfncF1984bVr186rXLmyV61aNe/EE0/0fvvtt8j8BQsWePqpGzJkiNe+fXuvdOnS3muvveb98ccf3kknneRVqVLFK1eunNe8eXPvs88+yzN9xxxzjNe7d+8c09544w332Wj5pbdnz54uLf6rQYMGbvqWLVu8a6+91qtZs6ZLn7bn+++/j6xz1KhRbvnPP//ca9OmjVeyZEk3befOnd6DDz7oNWzY0CtTpox30EEHee+//36++znedjRu3Ng755xzIvN79erllqlevbp37LHH7nK7RN/bsmVLlw4dh44dO3obNmyIzH/ppZe8pk2buu1r0qSJN2jQoFzHaOjQoe77ypYt67ZlwoQJObY/+tWvX7+427d8+XIvKyvLmz17dr774b333vNKlSrlbd++3b3Xvi1WrJi3dOnSyDLPPfecV6lSJW/r1q15rkf5pn///pH3t9xyi/evf/3L25WLL77Y7UMAQDhQ0g0AKHQ9e/a0jh072iWXXGLPPPOMK5FUyfdxxx1nrVu3dtVzVTqpasMqnfVt3LjRrr/+ejd/5MiRrnTztNNOy1V1V9V5e/fubXPnznVV2VW6u3XrVld6OWvWLHv44YetQoUKCadXVZbfe+89O+KIIyLT1qxZk296Bw4caPfee68ruVdJ85QpU9x0leirja9KV6dNm2aNGjVyadR3xG7DQw895LbhoIMOsgEDBth///tfe/75523OnDnWt29fV/I+ZsyYpPa9qkBHl2grHaVKlbLvvvvOrXtX26VtOffcc92xU9pGjx5tp59+uiuRl7feesuVGD/wwANu/oMPPmh33XWX+55od9xxh914442unfsBBxzg1qkSfdUKePLJJ12Vbn2XXlounvHjx7vq4c2aNct3m1W1XOvzaxpMnDjRVQNXswGfjsG6devcvo1HeUwl1mqfH121/dBDD3VND1RFXfvspZdeyvXZww8/3MaNG5dvGgEAGaSwo34AAGTZsmVejRo1XInj8OHDvfvuu8/r3Llzjp2zaNEiV9I5b968uDtNpa+ar5LZ6FLUJ598MsdyBx54oHfPPfckvONVAqwSZpX0qnRb6zzggAPc+n2JpPd//ud/IiXcotJgrfett96KTNu2bZtXp04d75FHHslR0vvhhx9GllHpuNLhlwb7Lr30Uu/cc89NqKR7x44drrRe637mmWci81u3bp3jM7varqlTp7r/q/ZAPPvvv7/39ttv51pn27Ztcxyjl19+OTJ/zpw5btrcuXPde9VOUG2GXdH+3W+//fJdRnmkfv363u233x6Zdvnll+faxo0bN0ZqGMTz8MMPe1WrVnX51qeSfL1uu+02b9q0ad4LL7zgSv9VkyPaRx995PK5aisAADIfJd0AgCJBJYNXXnmlK6Xs3r27zZgxw0aNGuVKoP1X06ZNI22j5ddff3Ulovvtt58ruWzYsGGkU6xoKn2Mpo607r//ftf+tl+/fjZz5sxdpq9Hjx6uFFbpUomqSqQ7d+7sSjslkfTG0nS1MVc6fGoPrJJQlQrntQ2//fabbdq0ybUxjv4+lXzn9V0+dZymZVXCrfbdKiG/+uqrI/MPOeSQHMvvartatWrlaimopFglvCrZVdtmvyaCllH76ejPa9/HplOl9z61e5fly5dbMjZv3uw6iMuLSq7VdlsdyN1zzz1WUOoZvX///q62g/JtdOl3mzZtXGm+SrnVU732sWoMRNO+17KqbQEAyHzJ9+ACAEBAVN3Xr/K7YcMG11GVqn7H8oMyzW/QoIEL9OrUqeMCGfUIHtsBWPny5XO8v+yyy1z14c8++8y+/vprV1VbHV6pl+u8VK5c2QXaor+vvPKKS8e7777r1pdIendH9Dbou0Tp32effXIsV7p06V0+PFBVbgV+Spff4Vi87/G/K7/tKl68uI0YMcJ1QKZ9qc7vtH51jqeq3qLjE10VX/S5aHrYEN3BniTbw3eNGjUiAX8sPRw54YQTrGLFijZ8+PAc36dO7b7//vscy6sKvT8v2pAhQ9zxfv/9913P+7H7QwF9ND1Eih0iTE0HtJ91DAAAmY+gGwBQJKnEUMGKSq/j9fL9zz//2Lx581xAd/TRR7tpKoFOVL169eyqq65yr9tuu82tJ7+gO5YfNKp0NZH0xrP//vtH2k/r4YGo5Fvtvfv06ZPn5xTYKbhWif4xxxxjyYh+eJCIRLZLQbJK6/VS+21tiwJbtbfXw5D58+e7YL+gtI/Uw/iuqHR56dKlLvCO7jFcJdx6yKJ9pnbXsaXhGv5Lbc5Vsu6XXOtBgmpPRAfR6uVdbdcVeKvEPJa2X3ky2i+//BI5tj71WaC0AgDCgerlAIAiSZ2dqURQ1ccVhKo6ssZF1nBQCsAUVGks5BdffNFVt/72229dkJcIBbRal4YmU+dlqj69q863VJ1bAZ1eqnKtKtkK3lTFPJH0xqPSTq1HY0KrgzINMaXqyPqu6CGtYqm0Vp2JqWq4OiTTd2k7VMoc20HZ7trVdqlEW9Wp1cmaHgJoPOsVK1ZE9qeqYasmgYbyUgCqjus0TNwTTzyRcBoU8KvEXZ3lrVy50u2feBTIqrRbDzGiA24dI1V1V+0EvfePo39cNF/B9QUXXOCOrbZP42lr2/2aA6pSfuGFF7oaESq199cRPd63jsekSZPc/lCe1GeUP7WeaOpEzc83AIAQKOxG5QAA+DQUVKtWrSLvf/nlFzd0mIb20lBSGnaqT58+blgxGTFihNesWTPXeZWGmRo9erTr/EodsUV30vXjjz/m2MnXXHON6+BLn9NQXRdccIG3cuXKPA+EOhiLHrJKHWhp2rfffptjuV2lN7YjNdm8ebMbMkydyOU3ZNjq1atzfE7rVAdxGoJLnbFpO7p06eKNGTMmqSHDEpmf33b99NNP7nv9Ic/UwdzTTz+d4/PqKO7ggw92w3Rp32n4tmHDhuV5jLStmqZt91111VVuGLP8hgyTm2++OTIEWvT+i/eK7ghPHcF17drVbZ+OxQ033BAZUszfN/HWoaHgon3yySdu+DTtC+2nF198Mcf8v/76yx0vdUYHAAiHLP1T2IE/AABAKqj0uUWLFq7kP7Zad1Fwyy23uOrvKgEHAIQD1csBAEDGUMdnqkYe24N9UaE24/fdd19hJwMAsAdR0g0AAAAAQEAo6QYAAAAAICAE3QAAAAAABISgGwAAAACAgBB0AwAAAAAQEIJuAAAAAAACQtANAAAAAEBACLoBAAAAAAgIQTcAAAAAAAEh6AYAAAAAwILx/wA2eZJ+OmLjXwAAAABJRU5ErkJggg==", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "Ancient-tail mean (first 12 periods): -0.362\n", "The composite never approaches zero in the ancient tail; a no-intercept\n", "exponential is therefore structurally unable to follow it.\n" ] } ], "source": [ "composite = z_scores[z_cols].mean(axis=1).values\n", "ybp = df['years_before_present'].values\n", "\n", "fig, ax = plt.subplots(figsize=(10, 5))\n", "ax.scatter(ybp, composite, c='steelblue', s=30, zorder=3)\n", "ax.set_xlabel(f'Years Before Present ({CURRENT_YEAR})')\n", "ax.set_ylabel('Composite Index (mean z-score)')\n", "ax.set_title('Historical Knowledge Composite Index (canonical v2 standardization)')\n", "ax.invert_xaxis()\n", "ax.axhline(0, color='gray', linewidth=0.5, linestyle='--')\n", "ax.grid(True, alpha=0.3)\n", "plt.tight_layout()\n", "plt.show()\n", "\n", "print(f\"Ancient-tail mean (first 12 periods): {composite[:12].mean():.3f}\")\n", "print(\"The composite never approaches zero in the ancient tail; a no-intercept\")\n", "print(\"exponential is therefore structurally unable to follow it.\")\n" ] }, { "cell_type": "code", "execution_count": 5, "id": "ed3501ba", "metadata": { "execution": { "iopub.execute_input": "2026-06-09T21:29:20.253207Z", "iopub.status.busy": "2026-06-09T21:29:20.253146Z", "iopub.status.idle": "2026-06-09T21:29:20.270581Z", "shell.execute_reply": "2026-06-09T21:29:20.270213Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "=== Model Comparison (AICc; k = number of curve parameters) ===\n", "Model R2 AICc RSS\n", "----------------------------------------------------------\n", "biexponential_intercept 0.9939 -155.93 0.1667\n", "exponential_intercept 0.9902 -146.26 0.2674\n", "power_law_intercept 0.9813 -125.61 0.5098\n", "exponential 0.9092 -77.53 2.4719\n", "power_law 0.8609 -63.89 3.7858\n", "logarithmic 0.6435 -33.77 9.7039\n", "linear 0.2457 -9.79 20.5304\n", "\n", "Best by AICc: biexponential_intercept\n", "Primary model (parsimony, 2 fewer params): exponential_intercept\n", "exponential_intercept: a=8.2827 lam=0.007277 c=-0.2880\n", "Half-life (primary) = 95.25 years\n", "Biexponential components: fast 83.9 yr, slow 817.2 yr\n", "Legacy no-intercept exponential on this composite: 78.46 years\n" ] } ], "source": [ "def exp_model(x, a, lam):\n", " return a * np.exp(-lam * x)\n", "\n", "def exp_c_model(x, a, lam, c):\n", " return a * np.exp(-lam * x) + c\n", "\n", "def biexp_c_model(x, a1, lam1, a2, lam2, c):\n", " return a1 * np.exp(-lam1 * x) + a2 * np.exp(-lam2 * x) + c\n", "\n", "def power_model(x, a, b):\n", " return a * x ** (-b)\n", "\n", "def power_c_model(x, a, b, c):\n", " return a * x ** (-b) + c\n", "\n", "def log_model(x, a, b):\n", " return a - b * np.log(x)\n", "\n", "def linear_model(x, a, b):\n", " return a - b * x\n", "\n", "def aicc(rss, n, k):\n", " # Corrected AIC. Convention: k = number of curve parameters.\n", " aic = n * np.log(rss / n) + 2 * k\n", " if n - k - 1 > 0:\n", " return aic + (2 * k * (k + 1)) / (n - k - 1)\n", " return np.inf\n", "\n", "x_fit = ybp[fit_mask]\n", "y_fit = composite[fit_mask]\n", "n_fit = len(x_fit)\n", "\n", "def fit_curve(func, x, y, p0, bounds=None, maxfev=300000):\n", " kwargs = {'p0': p0, 'maxfev': maxfev}\n", " if bounds is not None:\n", " kwargs['bounds'] = bounds\n", " popt, _ = curve_fit(func, x, y, **kwargs)\n", " pred = func(x, *popt)\n", " res = y - pred\n", " rss = float(np.sum(res ** 2))\n", " r2 = float(1 - rss / np.sum((y - np.mean(y)) ** 2))\n", " return popt, res, rss, r2\n", "\n", "# Multistart for the biexponential. Without it, many starts collapse to a\n", "# degenerate single-exponential local minimum (RSS 0.267 vs the global 0.167).\n", "BIEXP_BOUNDS = ([0, 1e-5, 0, 1e-7, -5], [100, 1, 100, 0.05, 5])\n", "BIEXP_P0S = [\n", " [3, 0.01, 0.5, 0.0005, -0.4],\n", " [5, 0.008, 1, 0.001, 0],\n", " [8, 0.012, 0.3, 0.0003, -0.5],\n", " [2, 0.02, 2, 0.002, -0.3],\n", "]\n", "\n", "def fit_biexp(x, y):\n", " best = None\n", " for p0 in BIEXP_P0S:\n", " try:\n", " cand = fit_curve(biexp_c_model, x, y, p0, bounds=BIEXP_BOUNDS)\n", " except Exception:\n", " continue\n", " if best is None or cand[2] < best[2]:\n", " best = cand\n", " popt, res, rss, r2 = best\n", " a1, l1, a2, l2, c = popt\n", " if l1 < l2: # order components: fast component first\n", " a1, l1, a2, l2 = a2, l2, a1, l1\n", " return np.array([a1, l1, a2, l2, c]), res, rss, r2\n", "\n", "MODEL_SPECS = {\n", " 'exponential': (exp_model, [5.0, 0.005], ['a', 'lam']),\n", " 'exponential_intercept': (exp_c_model, [5.0, 0.005, -0.4], ['a', 'lam', 'c']),\n", " 'power_law': (power_model, [100.0, 1.0], ['a', 'b']),\n", " 'power_law_intercept': (power_c_model, [100.0, 1.0, -0.4], ['a', 'b', 'c']),\n", " 'logarithmic': (log_model, [5.0, 0.7], ['a', 'b']),\n", " 'linear': (linear_model, [0.5, 0.0004], ['a', 'b']),\n", "}\n", "\n", "fit_results = {}\n", "for name, (func, p0, pnames) in MODEL_SPECS.items():\n", " popt, res, rss, r2 = fit_curve(func, x_fit, y_fit, p0)\n", " fit_results[name] = {'func': func, 'popt': popt, 'param_names': pnames,\n", " 'residuals': res, 'rss': rss, 'r_squared': r2,\n", " 'aicc': aicc(rss, n_fit, len(popt))}\n", "\n", "popt, res, rss, r2 = fit_biexp(x_fit, y_fit)\n", "fit_results['biexponential_intercept'] = {\n", " 'func': biexp_c_model, 'popt': popt,\n", " 'param_names': ['a1', 'lam1', 'a2', 'lam2', 'c'],\n", " 'residuals': res, 'rss': rss, 'r_squared': r2, 'aicc': aicc(rss, n_fit, len(popt))}\n", "\n", "half_lives = {\n", " 'exponential': float(np.log(2) / fit_results['exponential']['popt'][1]),\n", " 'exponential_intercept': float(np.log(2) / fit_results['exponential_intercept']['popt'][1]),\n", " 'biexponential_intercept': {\n", " 'fast': float(np.log(2) / fit_results['biexponential_intercept']['popt'][1]),\n", " 'slow': float(np.log(2) / fit_results['biexponential_intercept']['popt'][3]),\n", " },\n", "}\n", "\n", "PRIMARY_MODEL = 'exponential_intercept'\n", "aicc_order = sorted(fit_results, key=lambda m: fit_results[m]['aicc'])\n", "BEST_BY_AICC = aicc_order[0]\n", "\n", "print(\"=== Model Comparison (AICc; k = number of curve parameters) ===\")\n", "print(f\"{'Model':<26} {'R2':>8} {'AICc':>10} {'RSS':>10}\")\n", "print('-' * 58)\n", "for name in aicc_order:\n", " r = fit_results[name]\n", " print(f\"{name:<26} {r['r_squared']:>8.4f} {r['aicc']:>10.2f} {r['rss']:>10.4f}\")\n", "\n", "print()\n", "print(f\"Best by AICc: {BEST_BY_AICC}\")\n", "print(f\"Primary model (parsimony, 2 fewer params): {PRIMARY_MODEL}\")\n", "p = fit_results['exponential_intercept']['popt']\n", "print(f\"exponential_intercept: a={p[0]:.4f} lam={p[1]:.6f} c={p[2]:.4f}\")\n", "print(f\"Half-life (primary) = {half_lives['exponential_intercept']:.2f} years\")\n", "print(f\"Biexponential components: fast {half_lives['biexponential_intercept']['fast']:.1f} yr, \"\n", " f\"slow {half_lives['biexponential_intercept']['slow']:.1f} yr\")\n", "print(f\"Legacy no-intercept exponential on this composite: {half_lives['exponential']:.2f} years\")\n" ] }, { "cell_type": "code", "execution_count": 6, "id": "ef047c74", "metadata": { "execution": { "iopub.execute_input": "2026-06-09T21:29:20.271638Z", "iopub.status.busy": "2026-06-09T21:29:20.271572Z", "iopub.status.idle": "2026-06-09T21:29:22.242096Z", "shell.execute_reply": "2026-06-09T21:29:22.241742Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "block b=8 (primary): CI [86.99220459111012, 107.21449004810914], median 95.5, filtered 0\n", "block b=5: CI [88.33252889809108, 103.89239742816571], median 95.4, filtered 0\n", "iid residual: CI [88.32114433757764, 102.7121580098244], median 95.1, filtered 0\n", "Point estimate: 95.3 years\n", "Median sits on the point estimate: no leverage bias, unlike the v1 iid-pairs\n", "bootstrap whose median (70.6) sat far below its point estimate (84).\n", "Curve band computed on 185 grid points.\n" ] } ], "source": [ "# Moving-block residual bootstrap on the primary (exponential_intercept) fit.\n", "# Residuals are kept in chronological order; blocks preserve local autocorrelation.\n", "primary = fit_results['exponential_intercept']\n", "yhat_fit = exp_c_model(x_fit, *primary['popt'])\n", "resid_fit = y_fit - yhat_fit\n", "N_BOOT = 10000\n", "\n", "def block_bootstrap(b, seed, collect_popts=False):\n", " np.random.seed(seed)\n", " n = len(resid_fit)\n", " nblocks = int(np.ceil(n / b))\n", " hls, popts = [], []\n", " n_filtered = 0\n", " for _ in range(N_BOOT):\n", " starts = np.random.randint(0, n - b + 1, nblocks)\n", " idx = (starts[:, None] + np.arange(b)[None, :]).ravel()[:n]\n", " ystar = yhat_fit + resid_fit[idx]\n", " try:\n", " popt, _ = curve_fit(exp_c_model, x_fit, ystar, p0=[5.0, 0.005, -0.4], maxfev=20000)\n", " except Exception:\n", " n_filtered += 1\n", " continue\n", " if popt[1] <= 0:\n", " n_filtered += 1\n", " continue\n", " hls.append(np.log(2) / popt[1])\n", " if collect_popts:\n", " popts.append(popt)\n", " return np.array(hls), popts, n_filtered\n", "\n", "def iid_bootstrap(seed):\n", " np.random.seed(seed)\n", " n = len(resid_fit)\n", " hls = []\n", " n_filtered = 0\n", " for _ in range(N_BOOT):\n", " idx = np.random.randint(0, n, n)\n", " ystar = yhat_fit + resid_fit[idx]\n", " try:\n", " popt, _ = curve_fit(exp_c_model, x_fit, ystar, p0=[5.0, 0.005, -0.4], maxfev=20000)\n", " except Exception:\n", " n_filtered += 1\n", " continue\n", " if popt[1] <= 0:\n", " n_filtered += 1\n", " continue\n", " hls.append(np.log(2) / popt[1])\n", " return np.array(hls), n_filtered\n", "\n", "def ci95(arr):\n", " return [float(np.percentile(arr, 2.5)), float(np.percentile(arr, 97.5))]\n", "\n", "hl_b5, _, filt5 = block_bootstrap(5, seed=41)\n", "hl_b8, popts_b8, filt8 = block_bootstrap(8, seed=42, collect_popts=True)\n", "hl_iid, filt_iid = iid_bootstrap(seed=43)\n", "\n", "print(f\"block b=8 (primary): CI {ci95(hl_b8)}, median {np.median(hl_b8):.1f}, filtered {filt8}\")\n", "print(f\"block b=5: CI {ci95(hl_b5)}, median {np.median(hl_b5):.1f}, filtered {filt5}\")\n", "print(f\"iid residual: CI {ci95(hl_iid)}, median {np.median(hl_iid):.1f}, filtered {filt_iid}\")\n", "print(f\"Point estimate: {half_lives['exponential_intercept']:.1f} years\")\n", "print(\"Median sits on the point estimate: no leverage bias, unlike the v1 iid-pairs\")\n", "print(\"bootstrap whose median (70.6) sat far below its point estimate (84).\")\n", "\n", "# Pointwise 95% envelope of the fitted curve across b=8 bootstrap fits\n", "t_grid = np.arange(0, 4601, 25)\n", "curves = np.array([exp_c_model(t_grid, *p) for p in popts_b8])\n", "band_lower = np.percentile(curves, 2.5, axis=0)\n", "band_upper = np.percentile(curves, 97.5, axis=0)\n", "print(f\"Curve band computed on {len(t_grid)} grid points.\")\n" ] }, { "cell_type": "code", "execution_count": 7, "id": "3cf9fe64", "metadata": { "execution": { "iopub.execute_input": "2026-06-09T21:29:22.243319Z", "iopub.status.busy": "2026-06-09T21:29:22.243244Z", "iopub.status.idle": "2026-06-09T21:29:22.249749Z", "shell.execute_reply": "2026-06-09T21:29:22.249143Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Dropped period Half-life\n", "----------------------------------\n", "3rd millennium BC 95.2\n", "2nd millennium BC 95.2\n", "10th century BC 95.0\n", "9th century BC 95.0\n", "8th century BC 95.0\n", "7th century BC 95.0\n", "6th century BC 95.0\n", "5th century BC 95.1\n", "4th century BC 95.1\n", "3rd century BC 95.0\n", "2nd century BC 95.0\n", "1st century BC 95.1\n", "1st century 95.9\n", "2nd century 95.3\n", "3rd century 95.3\n", "4th century 95.4\n", "5th century 95.3\n", "6th century 95.3\n", "7th century 95.3\n", "8th century 95.3\n", "9th century 95.3\n", "10th century 95.3\n", "11th century 95.3\n", "12th century 95.4\n", "13th century 95.3\n", "14th century 95.0\n", "15th century 94.5\n", "16th century 93.6\n", "17th century 95.3\n", "18th century 103.6\n", "19th century 92.5\n", "20th century 96.0\n", "\n", "LOO range: 92.5 (drop 19th century) to 103.6 (drop 18th century)\n" ] } ], "source": [ "# Leave-one-period-out half-life panel (primary model, drop each of the 32 fit points)\n", "loo_panel = []\n", "fit_labels = df['label'].values[fit_mask]\n", "for i in range(n_fit):\n", " keep = np.ones(n_fit, dtype=bool)\n", " keep[i] = False\n", " popt_i, _, _, _ = fit_curve(exp_c_model, x_fit[keep], y_fit[keep], [5.0, 0.005, -0.4])\n", " loo_panel.append({'label': str(fit_labels[i]), 'half_life': float(np.log(2) / popt_i[1])})\n", "\n", "loo_hls = [e['half_life'] for e in loo_panel]\n", "i_min, i_max = int(np.argmin(loo_hls)), int(np.argmax(loo_hls))\n", "print(f\"{'Dropped period':<22} {'Half-life':>10}\")\n", "print('-' * 34)\n", "for e in loo_panel:\n", " print(f\"{e['label']:<22} {e['half_life']:>10.1f}\")\n", "print()\n", "print(f\"LOO range: {loo_hls[i_min]:.1f} (drop {loo_panel[i_min]['label']}) \"\n", " f\"to {loo_hls[i_max]:.1f} (drop {loo_panel[i_max]['label']})\")\n" ] }, { "cell_type": "code", "execution_count": 8, "id": "90d126bf", "metadata": { "execution": { "iopub.execute_input": "2026-06-09T21:29:22.251739Z", "iopub.status.busy": "2026-06-09T21:29:22.251597Z", "iopub.status.idle": "2026-06-09T21:29:22.262229Z", "shell.execute_reply": "2026-06-09T21:29:22.261870Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "=== Spearman rho (BH-adjusted p in parens) ===\n", " openalex_work_count vs named_individuals: rho=0.933 (p_adj=1.65e-14)\n", " openalex_work_count vs source_proxy: rho=0.915 (p_adj=2.59e-13)\n", " openalex_work_count vs ngram_discourse: rho=0.794 (p_adj=4.17e-08)\n", " named_individuals vs source_proxy: rho=0.905 (p_adj=9.62e-13)\n", " named_individuals vs ngram_discourse: rho=0.872 (p_adj=5.63e-11)\n", " source_proxy vs ngram_discourse: rho=0.781 (p_adj=8.01e-08)\n", "rho range: 0.78 to 0.93\n" ] } ], "source": [ "# Spearman rank correlations between the 4 active metrics (raw values, all 33\n", "# periods; rank-based, so independent of any standardization choice), with\n", "# Benjamini-Hochberg correction across the 6 pairwise tests.\n", "n_m = len(active_metrics)\n", "rho = np.eye(n_m)\n", "pairs, pvals = [], []\n", "for i in range(n_m):\n", " for j in range(i + 1, n_m):\n", " r_ij, p_ij = spearmanr(df[active_metrics[i]], df[active_metrics[j]], nan_policy='omit')\n", " rho[i, j] = rho[j, i] = r_ij\n", " pairs.append((i, j))\n", " pvals.append(p_ij)\n", "\n", "pvals = np.array(pvals)\n", "m_tests = len(pvals)\n", "order = np.argsort(pvals)\n", "p_adj_flat = np.empty(m_tests)\n", "running = 1.0\n", "for rank_idx in range(m_tests - 1, -1, -1):\n", " k = order[rank_idx]\n", " running = min(running, pvals[k] * m_tests / (rank_idx + 1))\n", " p_adj_flat[k] = running\n", "\n", "p_adjusted = np.zeros((n_m, n_m))\n", "significant = np.zeros((n_m, n_m), dtype=bool)\n", "for (i, j), p_adj in zip(pairs, p_adj_flat):\n", " p_adjusted[i, j] = p_adjusted[j, i] = p_adj\n", " significant[i, j] = significant[j, i] = bool(p_adj < 0.05)\n", "np.fill_diagonal(significant, True)\n", "\n", "print(\"=== Spearman rho (BH-adjusted p in parens) ===\")\n", "for i in range(n_m):\n", " for j in range(i + 1, n_m):\n", " print(f\" {active_metrics[i]} vs {active_metrics[j]}: rho={rho[i,j]:.3f} (p_adj={p_adjusted[i,j]:.2e})\")\n", "off_diag = rho[~np.eye(n_m, dtype=bool)]\n", "print(f\"rho range: {off_diag.min():.2f} to {off_diag.max():.2f}\")\n" ] }, { "cell_type": "code", "execution_count": 9, "id": "1fc19f2a", "metadata": { "execution": { "iopub.execute_input": "2026-06-09T21:29:22.263171Z", "iopub.status.busy": "2026-06-09T21:29:22.263102Z", "iopub.status.idle": "2026-06-09T21:29:22.267529Z", "shell.execute_reply": "2026-06-09T21:29:22.266793Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Most overstudied (fit periods): 19th century (+2.20)\n", "Most understudied (fit periods): 20th century (-1.27)\n", "21st century (excluded from extremes): -12.25 (standardization artifact)\n" ] } ], "source": [ "# Attention vs sources: z_attention = mean of the 3 attention metrics,\n", "# bias = z_attention - z_source (positive = more studied than its surviving\n", "# sources would predict).\n", "att_cols = ['z_openalex_work_count', 'z_named_individuals', 'z_ngram_discourse']\n", "z_attention = z_scores[att_cols].mean(axis=1).values\n", "z_source = z_scores['z_source_proxy'].values\n", "bias_vals = z_attention - z_source\n", "\n", "bias_rows = []\n", "for i in range(len(df)):\n", " bias_rows.append({\n", " 'label': str(df['label'].values[i]),\n", " 'z_source': float(z_source[i]),\n", " 'z_attention': float(z_attention[i]),\n", " 'bias_residual': float(bias_vals[i]),\n", " })\n", "\n", "# Extremes are reported over the 32 fit periods only. The 21st century is\n", "# scored against a basis that excludes it, so its z_source (about +13) is a\n", "# standardization artifact, not evidence about study patterns.\n", "bias_fit = bias_vals[fit_mask]\n", "labels_fit = df['label'].values[fit_mask]\n", "bias_i_max = int(np.argmax(bias_fit))\n", "bias_i_min = int(np.argmin(bias_fit))\n", "fit_extremes = {\n", " 'max_label': str(labels_fit[bias_i_max]), 'max_value': float(bias_fit[bias_i_max]),\n", " 'min_label': str(labels_fit[bias_i_min]), 'min_value': float(bias_fit[bias_i_min]),\n", "}\n", "print(f\"Most overstudied (fit periods): {fit_extremes['max_label']} ({fit_extremes['max_value']:+.2f})\")\n", "print(f\"Most understudied (fit periods): {fit_extremes['min_label']} ({fit_extremes['min_value']:+.2f})\")\n", "print(f\"21st century (excluded from extremes): {bias_vals[32]:+.2f} (standardization artifact)\")\n" ] }, { "cell_type": "code", "execution_count": 10, "id": "2432bd0a", "metadata": { "execution": { "iopub.execute_input": "2026-06-09T21:29:22.268769Z", "iopub.status.busy": "2026-06-09T21:29:22.268689Z", "iopub.status.idle": "2026-06-09T21:29:22.282587Z", "shell.execute_reply": "2026-06-09T21:29:22.282083Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "=== Sensitivity Analysis (all refits use the primary exponential_intercept model) ===\n", " drop openalex_work_count HL= 74.1 R2=0.990\n", " drop named_individuals HL= 118.5 R2=0.984\n", " drop source_proxy HL= 123.0 R2=0.982\n", " drop ngram_discourse HL= 70.9 R2=0.994\n", " exclude millennia (n=30) HL= 95.3 R2=0.990\n", " exclude 1st century (n=31) HL= 96.0 R2=0.992\n", " median composite HL= 95.7 R2=0.991\n", " log-scale composite HL= 817.7 R2=0.859\n", "\n", " Standardization variants (half-life in years):\n", " basis no intercept intercept\n", " all33_ddof0 85.2 109.1\n", " century31_ddof1 84.0 109.3\n", " fit32_ddof1 78.5 95.3\n", "\n", " The originally published 84.0 = century-basis standardization with no\n", " intercept: 84.01. That cell of the grid is the v1 result.\n" ] } ], "source": [ "print(\"=== Sensitivity Analysis (all refits use the primary exponential_intercept model) ===\")\n", "sens = {}\n", "\n", "# Leave one METRIC out (per-metric z-scores are unchanged by dropping another metric)\n", "lomo = []\n", "for m_drop in active_metrics:\n", " kept = [f'z_{c}' for c in active_metrics if c != m_drop]\n", " comp_k = z_scores[kept].mean(axis=1).values\n", " popt_k, _, _, r2_k = fit_curve(exp_c_model, x_fit, comp_k[fit_mask], [5.0, 0.005, -0.4])\n", " lomo.append({'metric_dropped': m_drop,\n", " 'half_life': float(np.log(2) / popt_k[1]), 'r_squared': r2_k})\n", " print(f\" drop {m_drop:<24} HL={lomo[-1]['half_life']:>6.1f} R2={r2_k:.3f}\")\n", "sens['leave_one_metric_out'] = lomo\n", "\n", "# Exclude the two millennium bins (re-standardized on the surviving 30 century fit rows)\n", "cent_mask = (df['span_years'].values == 100) & fit_mask\n", "zz = zscores_on_basis(df, cent_mask, ddof=1)\n", "comp_c = zz[z_cols].mean(axis=1).values\n", "popt_c, _, _, r2_c = fit_curve(exp_c_model, ybp[cent_mask], comp_c[cent_mask], [5.0, 0.005, -0.4])\n", "sens['exclude_millennia'] = {'half_life': float(np.log(2) / popt_c[1]), 'r_squared': r2_c,\n", " 'n': int(cent_mask.sum())}\n", "print(f\" exclude millennia (n={cent_mask.sum()}) HL={sens['exclude_millennia']['half_life']:>6.1f} R2={r2_c:.3f}\")\n", "\n", "# Exclude the 1st century CE (re-standardized on the remaining 31 fit rows)\n", "mask1 = fit_mask & (df['label'].values != '1st century')\n", "zz1 = zscores_on_basis(df, mask1, ddof=1)\n", "comp_1 = zz1[z_cols].mean(axis=1).values\n", "popt_1, _, _, r2_1 = fit_curve(exp_c_model, ybp[mask1], comp_1[mask1], [5.0, 0.005, -0.4])\n", "sens['exclude_first_century'] = {'half_life': float(np.log(2) / popt_1[1]), 'r_squared': r2_1,\n", " 'n': int(mask1.sum())}\n", "print(f\" exclude 1st century (n={mask1.sum()}) HL={sens['exclude_first_century']['half_life']:>6.1f} R2={r2_1:.3f}\")\n", "\n", "# Median composite\n", "comp_med = z_scores[z_cols].median(axis=1).values\n", "popt_m, _, _, r2_m = fit_curve(exp_c_model, x_fit, comp_med[fit_mask], [5.0, 0.005, -0.4])\n", "sens['median_composite'] = {'half_life': float(np.log(2) / popt_m[1]), 'r_squared': r2_m}\n", "print(f\" median composite HL={sens['median_composite']['half_life']:>6.1f} R2={r2_m:.3f}\")\n", "\n", "# Log-scale composite (z of log10 values; measures the slow cultural layer)\n", "zz_log = zscores_on_basis(df, fit_mask, ddof=1, transform=np.log10)\n", "comp_log = zz_log[z_cols].mean(axis=1).values\n", "popt_l, _, _, r2_l = fit_curve(exp_c_model, x_fit, comp_log[fit_mask], [10.0, 0.001, -5.0])\n", "sens['log_scale'] = {'half_life': float(np.log(2) / popt_l[1]), 'r_squared': r2_l}\n", "print(f\" log-scale composite HL={sens['log_scale']['half_life']:>6.1f} R2={r2_l:.3f}\")\n", "\n", "# Standardization variants grid: how the half-life depends on the z basis and\n", "# the intercept, fit always on the first 32 periods. The century-basis\n", "# no-intercept cell reproduces the originally published 84.0.\n", "variants = {\n", " 'all33_ddof0': (np.ones(len(df), dtype=bool), 0),\n", " 'century31_ddof1': (df['span_years'].values == 100, 1),\n", " 'fit32_ddof1': (fit_mask, 1),\n", "}\n", "grid = {'no_intercept': {}, 'intercept': {}}\n", "for vname, (bmask, dd) in variants.items():\n", " zz_v = zscores_on_basis(df, bmask, ddof=dd)\n", " comp_v = zz_v[z_cols].mean(axis=1).values\n", " p_no, _, _, _ = fit_curve(exp_model, x_fit, comp_v[fit_mask], [5.0, 0.005])\n", " p_yes, _, _, _ = fit_curve(exp_c_model, x_fit, comp_v[fit_mask], [5.0, 0.005, -0.4])\n", " grid['no_intercept'][vname] = float(np.log(2) / p_no[1])\n", " grid['intercept'][vname] = float(np.log(2) / p_yes[1])\n", "sens['standardization_variants'] = grid\n", "\n", "print()\n", "print(\" Standardization variants (half-life in years):\")\n", "print(f\" {'basis':<18} {'no intercept':>13} {'intercept':>10}\")\n", "for vname in variants:\n", " print(f\" {vname:<18} {grid['no_intercept'][vname]:>13.1f} {grid['intercept'][vname]:>10.1f}\")\n", "print()\n", "print(f\" The originally published 84.0 = century-basis standardization with no\")\n", "print(f\" intercept: {grid['no_intercept']['century31_ddof1']:.2f}. That cell of the grid is the v1 result.\")\n" ] }, { "cell_type": "code", "execution_count": 11, "id": "29b13979", "metadata": { "execution": { "iopub.execute_input": "2026-06-09T21:29:22.283887Z", "iopub.status.busy": "2026-06-09T21:29:22.283780Z", "iopub.status.idle": "2026-06-09T21:29:22.289307Z", "shell.execute_reply": "2026-06-09T21:29:22.288840Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Durbin-Watson (exponential_intercept): 1.166\n", "Durbin-Watson (no-intercept exponential, same composite): 0.158\n", "The no-intercept misfit is what drove the v1 autocorrelation alarm; the\n", "intercept resolves most of it.\n", "Shapiro-Wilk p: 0.366 (normality not rejected)\n", "Goldfeld-Quandt ratio: 1.24 (p=0.344)\n", "Anomalies (|standardized residual| > 2): [('1st century', 2.45), ('18th century', -2.38)]\n" ] } ], "source": [ "# Residual diagnostics for the primary model\n", "res_c = fit_results['exponential_intercept']['residuals']\n", "dw = float(smtools.durbin_watson(res_c))\n", "sw_stat, sw_p = shapiro(res_c)\n", "X_design = np.column_stack([np.ones(n_fit), x_fit])\n", "gq_stat, gq_p, _ = smdiag.het_goldfeldquandt(res_c, X_design)\n", "legacy_dw = float(smtools.durbin_watson(fit_results['exponential']['residuals']))\n", "\n", "std_res = res_c / res_c.std(ddof=1)\n", "anomalies = []\n", "for i in range(n_fit):\n", " if abs(std_res[i]) > 2:\n", " anomalies.append({'label': str(fit_labels[i]), 'std_residual': float(std_res[i])})\n", "\n", "print(f\"Durbin-Watson (exponential_intercept): {dw:.3f}\")\n", "print(f\"Durbin-Watson (no-intercept exponential, same composite): {legacy_dw:.3f}\")\n", "print(\"The no-intercept misfit is what drove the v1 autocorrelation alarm; the\")\n", "print(\"intercept resolves most of it.\")\n", "print(f\"Shapiro-Wilk p: {sw_p:.3f} (normality not rejected)\")\n", "print(f\"Goldfeld-Quandt ratio: {float(gq_stat):.2f} (p={float(gq_p):.3f})\")\n", "print(f\"Anomalies (|standardized residual| > 2): {[(a['label'], round(a['std_residual'], 2)) for a in anomalies]}\")\n" ] }, { "cell_type": "code", "execution_count": 12, "id": "af580bf2", "metadata": { "execution": { "iopub.execute_input": "2026-06-09T21:29:22.290468Z", "iopub.status.busy": "2026-06-09T21:29:22.290396Z", "iopub.status.idle": "2026-06-09T21:29:22.293346Z", "shell.execute_reply": "2026-06-09T21:29:22.292968Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Cronbach's alpha (32 fit rows): 0.954\n", "PCA PC1 variance explained (32 fit rows): 0.878\n", "Originally published v1 values (old standardization basis, all 33 rows):\n", "alpha 0.879, PC1 0.748. The v2 values use the canonical fit-row basis.\n" ] } ], "source": [ "# Scale reliability on the 32 fit rows (numpy only; no sklearn dependency)\n", "Z_fit = z_scores[z_cols].values[fit_mask]\n", "k_items = Z_fit.shape[1]\n", "item_vars = Z_fit.var(axis=0, ddof=1)\n", "total_var = Z_fit.sum(axis=1).var(ddof=1)\n", "cronbach_alpha = float(k_items / (k_items - 1) * (1 - item_vars.sum() / total_var))\n", "\n", "cov = np.cov(Z_fit, rowvar=False, ddof=1)\n", "evals = np.linalg.eigh(cov)[0]\n", "pca_pc1 = float(evals[-1] / evals.sum())\n", "\n", "print(f\"Cronbach's alpha (32 fit rows): {cronbach_alpha:.3f}\")\n", "print(f\"PCA PC1 variance explained (32 fit rows): {pca_pc1:.3f}\")\n", "print(\"Originally published v1 values (old standardization basis, all 33 rows):\")\n", "print(\"alpha 0.879, PC1 0.748. The v2 values use the canonical fit-row basis.\")\n" ] }, { "cell_type": "code", "execution_count": 13, "id": "4b48ebc3", "metadata": { "execution": { "iopub.execute_input": "2026-06-09T21:29:22.294907Z", "iopub.status.busy": "2026-06-09T21:29:22.294819Z", "iopub.status.idle": "2026-06-09T21:29:22.301726Z", "shell.execute_reply": "2026-06-09T21:29:22.301395Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Wrote results.json\n", "Headline half-life: 95.25 years, CI 87.0 to 107.2 (block b=8)\n" ] } ], "source": [ "# Assemble and WRITE results.json (this notebook is the generator of record)\n", "def py6(o):\n", " # Convert numpy types and round floats to 6 significant digits\n", " if isinstance(o, dict):\n", " return {k: py6(v) for k, v in o.items()}\n", " if isinstance(o, (list, tuple, np.ndarray)):\n", " return [py6(v) for v in o]\n", " if isinstance(o, (bool, np.bool_)):\n", " return bool(o)\n", " if isinstance(o, (np.floating, float)):\n", " return float(f'{float(o):.6g}')\n", " if isinstance(o, (np.integer, int)):\n", " return int(o)\n", " return o\n", "\n", "model_fits_out = {}\n", "for name, r in fit_results.items():\n", " entry = {\n", " 'converged': True,\n", " 'params': dict(zip(r['param_names'], [float(v) for v in r['popt']])),\n", " 'aicc': r['aicc'],\n", " 'rss': r['rss'],\n", " 'n': n_fit,\n", " 'r_squared': r['r_squared'],\n", " }\n", " if name == 'exponential_intercept':\n", " entry['residuals'] = list(r['residuals'])\n", " model_fits_out[name] = entry\n", "\n", "results_out = {\n", " 'metadata': {\n", " 'version': '2.0',\n", " 'generated': '2026-06-09',\n", " 'n_periods': int(len(df)),\n", " 'n_fit': int(n_fit),\n", " 'current_year': CURRENT_YEAR,\n", " 'spec': ('Per-metric z-scores standardized on the 32 fit periods (partial 21st '\n", " 'century excluded from standardization and fit), ddof=1; composite = '\n", " 'arithmetic mean of 4 z-columns; x = years before 2026 at period '\n", " 'midpoint; primary model = exponential with intercept.'),\n", " 'lineage': ('v1 results were produced by a standalone script standardizing on '\n", " 'century rows only (including the partial 21st century). This '\n", " 'notebook is the sole generator as of v2; the standardization_variants '\n", " 'grid reproduces the published v1 figure (84.0) under the old spec.'),\n", " 'license': 'CC-BY-4.0',\n", " },\n", " 'composite': {\n", " 'labels': [str(l) for l in df['label'].values],\n", " 'years_before_present': list(ybp),\n", " 'values': list(composite),\n", " },\n", " 'model_fits': model_fits_out,\n", " 'primary_model': PRIMARY_MODEL,\n", " 'best_model_by_aicc': BEST_BY_AICC,\n", " 'half_lives': half_lives,\n", " 'bootstrap_ci': {\n", " 'model': 'exponential_intercept',\n", " 'method': 'moving-block residual bootstrap',\n", " 'n_boot': N_BOOT,\n", " 'block8': {'half_life_ci': ci95(hl_b8), 'median': float(np.median(hl_b8)),\n", " 'n_filtered': int(filt8), 'seed': 42},\n", " 'block5': {'half_life_ci': ci95(hl_b5), 'median': float(np.median(hl_b5)),\n", " 'n_filtered': int(filt5), 'seed': 41},\n", " 'iid_residual': {'half_life_ci': ci95(hl_iid), 'median': float(np.median(hl_iid)),\n", " 'n_filtered': int(filt_iid), 'seed': 43},\n", " 'curve_band': {'t': list(t_grid), 'lower': list(band_lower), 'upper': list(band_upper)},\n", " },\n", " 'loo_half_life': {\n", " 'panel': loo_panel,\n", " 'min': float(loo_hls[i_min]), 'min_label': loo_panel[i_min]['label'],\n", " 'max': float(loo_hls[i_max]), 'max_label': loo_panel[i_max]['label'],\n", " },\n", " 'sensitivity_analysis': sens,\n", " 'residual_diagnostics': {\n", " 'model': 'exponential_intercept',\n", " 'durbin_watson': dw,\n", " 'shapiro_wilk_p': float(sw_p),\n", " 'goldfeld_quandt_ratio': float(gq_stat),\n", " 'goldfeld_quandt_p': float(gq_p),\n", " 'legacy_no_intercept_durbin_watson': legacy_dw,\n", " 'anomalies': anomalies,\n", " },\n", " 'scale_reliability': {\n", " 'cronbach_alpha': cronbach_alpha,\n", " 'pca_variance_explained_pc1': pca_pc1,\n", " 'basis': '32 fit periods',\n", " 'published_v1': {'cronbach_alpha': 0.879, 'pca_variance_explained_pc1': 0.748},\n", " },\n", " 'bias_residual': {'rows': bias_rows, 'fit_extremes': fit_extremes},\n", " 'correlation_matrix': {\n", " 'metrics': active_metrics,\n", " 'r': rho.tolist(),\n", " 'p_adjusted': p_adjusted.tolist(),\n", " 'significant': significant.tolist(),\n", " },\n", "}\n", "\n", "with open('results.json', 'w') as f:\n", " json.dump(py6(results_out), f, indent=2)\n", "print(\"Wrote results.json\")\n", "print(f\"Headline half-life: {half_lives['exponential_intercept']:.2f} years, \"\n", " f\"CI {ci95(hl_b8)[0]:.1f} to {ci95(hl_b8)[1]:.1f} (block b=8)\")\n" ] }, { "cell_type": "code", "execution_count": 14, "id": "d166d994", "metadata": { "execution": { "iopub.execute_input": "2026-06-09T21:29:22.303014Z", "iopub.status.busy": "2026-06-09T21:29:22.302945Z", "iopub.status.idle": "2026-06-09T21:29:22.956291Z", "shell.execute_reply": "2026-06-09T21:29:22.955890Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " OK: Composite index (max diff 2.44e-06)\n", " OK: Half-life (primary) (95.25 vs 95.2538)\n", " OK: R-squared (primary) (0.9902 vs 0.990177)\n", " OK: half_life = ln2/lam consistency \n", " OK: AICc ordering (biexponential_intercept > exponential_intercept > power_law_intercept ...)\n", " OK: Primary/best model labels \n", " OK: Durbin-Watson (1.166)\n", " OK: Cronbach's alpha (0.954)\n", " OK: Spearman correlations \n", " OK: Published v1 figure (84.0) reproduces under v1 spec (84.01)\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ " OK: Bootstrap b=8 CI ([87.0, 107.2] vs [86.9922, 107.214])\n", "\n", "ALL CHECKS PASSED. Results are reproducible.\n" ] } ], "source": [ "# Self-contained verification: re-read data.json and the just-written\n", "# results.json from disk and recompute every headline quantity from scratch.\n", "checks = []\n", "\n", "def check(name, ok, detail=''):\n", " checks.append((name, bool(ok)))\n", " status = 'OK' if ok else ('FA' + 'ILED')\n", " print(f\" {status}: {name} {detail}\")\n", "\n", "with open('data.json') as f:\n", " raw_v = json.load(f)\n", "with open('results.json') as f:\n", " res_v = json.load(f)\n", "\n", "rows_v = []\n", "for p in raw_v['periods']:\n", " m = p['metrics']\n", " row = {'label': p['label'], 'span_years': p['span_years'],\n", " 'mid_year': (p['start_year'] + p['end_year']) / 2}\n", " for key in active_metrics:\n", " row[key] = m[key]['value'] if (m.get(key) and isinstance(m[key], dict)) else np.nan\n", " rows_v.append(row)\n", "df_v = pd.DataFrame(rows_v)\n", "ybp_v = CURRENT_YEAR - df_v['mid_year'].values\n", "mask_v = df_v.index.values < 32\n", "\n", "zz_v = pd.DataFrame()\n", "for col in active_metrics:\n", " v = df_v[col].values.astype(float)\n", " zz_v[f'z_{col}'] = (v - v[mask_v].mean()) / v[mask_v].std(ddof=1)\n", "comp_v = zz_v.mean(axis=1).values\n", "\n", "check('Composite index',\n", " np.allclose(comp_v, res_v['composite']['values'], atol=1e-6),\n", " f\"(max diff {np.max(np.abs(comp_v - np.array(res_v['composite']['values']))):.2e})\")\n", "\n", "popt_v, res_resid_v, rss_v, r2_v = fit_curve(exp_c_model, ybp_v[mask_v], comp_v[mask_v], [5.0, 0.005, -0.4])\n", "hl_v = float(np.log(2) / popt_v[1])\n", "check('Half-life (primary)', abs(hl_v - res_v['half_lives']['exponential_intercept']) < 0.1,\n", " f\"({hl_v:.2f} vs {res_v['half_lives']['exponential_intercept']})\")\n", "check('R-squared (primary)', abs(r2_v - res_v['model_fits']['exponential_intercept']['r_squared']) < 0.005,\n", " f\"({r2_v:.4f} vs {res_v['model_fits']['exponential_intercept']['r_squared']})\")\n", "check('half_life = ln2/lam consistency',\n", " abs(res_v['half_lives']['exponential_intercept']\n", " - np.log(2) / res_v['model_fits']['exponential_intercept']['params']['lam']) < 0.05)\n", "\n", "order_file = sorted(res_v['model_fits'], key=lambda m_: res_v['model_fits'][m_]['aicc'])\n", "order_recomputed = sorted(fit_results, key=lambda m_: fit_results[m_]['aicc'])\n", "check('AICc ordering', order_file == order_recomputed, f\"({' > '.join(order_file[:3])} ...)\")\n", "check('Primary/best model labels',\n", " res_v['primary_model'] == 'exponential_intercept'\n", " and res_v['best_model_by_aicc'] == 'biexponential_intercept')\n", "\n", "dw_v = float(smtools.durbin_watson(res_resid_v))\n", "check('Durbin-Watson', abs(dw_v - res_v['residual_diagnostics']['durbin_watson']) < 0.01,\n", " f\"({dw_v:.3f})\")\n", "\n", "Zf_v = zz_v.values[mask_v]\n", "alpha_v = Zf_v.shape[1] / (Zf_v.shape[1] - 1) * (1 - Zf_v.var(axis=0, ddof=1).sum() / Zf_v.sum(axis=1).var(ddof=1))\n", "check(\"Cronbach's alpha\", abs(alpha_v - res_v['scale_reliability']['cronbach_alpha']) < 0.005,\n", " f\"({alpha_v:.3f})\")\n", "\n", "rho_v = np.eye(len(active_metrics))\n", "for i in range(len(active_metrics)):\n", " for j in range(i + 1, len(active_metrics)):\n", " r_ij, _ = spearmanr(df_v[active_metrics[i]], df_v[active_metrics[j]], nan_policy='omit')\n", " rho_v[i, j] = rho_v[j, i] = r_ij\n", "check('Spearman correlations', np.allclose(rho_v, res_v['correlation_matrix']['r'], atol=1e-5))\n", "\n", "# The published v1 figure must reproduce under the v1 spec (century-basis z, no intercept)\n", "cmask_v = df_v['span_years'].values == 100\n", "zz_old = pd.DataFrame()\n", "for col in active_metrics:\n", " v = df_v[col].values.astype(float)\n", " zz_old[f'z_{col}'] = (v - v[cmask_v].mean()) / v[cmask_v].std(ddof=1)\n", "comp_old = zz_old.mean(axis=1).values\n", "popt_old, _, _, _ = fit_curve(exp_model, ybp_v[mask_v], comp_old[mask_v], [5.0, 0.005])\n", "hl_old = float(np.log(2) / popt_old[1])\n", "check('Published v1 figure (84.0) reproduces under v1 spec', abs(hl_old - 84.0) < 0.1,\n", " f\"({hl_old:.2f})\")\n", "\n", "# Bootstrap reproduction: re-seed 42 and re-run the b=8 block bootstrap.\n", "# Tolerance 1.0 year exists for cross-environment RNG/library differences.\n", "hl_b8_v, _, _ = block_bootstrap(8, seed=42)\n", "ci_v = ci95(hl_b8_v)\n", "ci_file = res_v['bootstrap_ci']['block8']['half_life_ci']\n", "check('Bootstrap b=8 CI', max(abs(ci_v[0] - ci_file[0]), abs(ci_v[1] - ci_file[1])) < 1.0,\n", " f\"([{ci_v[0]:.1f}, {ci_v[1]:.1f}] vs {ci_file})\")\n", "\n", "n_fail = sum(1 for _, ok in checks if not ok)\n", "print()\n", "if n_fail:\n", " raise RuntimeError(f\"{n_fail} verification check(s) did not pass\")\n", "print(\"ALL CHECKS PASSED. Results are reproducible.\")\n" ] }, { "cell_type": "markdown", "id": "e4f2ff97", "metadata": {}, "source": [ "## Changelog and license\n", "\n", "**v2 (June 2026).** Corrected analysis. Two defects in the v1 pipeline were found and fixed: (1) the primary model had no intercept, forcing the curve to decay to zero while the z-score composite sits on a floor near -0.33 in the ancient tail (the misfit was visible as Durbin-Watson 0.17); (2) metrics were standardized on century-length rows only, including the partial 21st century. The corrected specification (intercept model, standardization on the 32 fit periods) moves the headline half-life from 84 to roughly 95 years. The standardization-variants grid above reproduces the published v1 figure exactly under the v1 specification.\n", "\n", "**v1 (April 2026).** Original release.\n", "\n", "Data, code, and results are released under [CC BY 4.0](https://creativecommons.org/licenses/by/4.0/). Mark Pavlyukovskyy.\n" ] } ], "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.12" } }, "nbformat": 4, "nbformat_minor": 5 }