{
"cells": [
{
"cell_type": "markdown",
"id": "E2PXadNcP-1p",
"metadata": {
"id": "E2PXadNcP-1p"
},
"source": [
"![\"SOS](\"https://indico.in2p3.fr/event/37891/logo-2009395760.png\")
SOS 2026
\n",
"\n",
"# Hands-on session: MCMC with the Metropolis-Hastings algorithm\n",
"\n",
"**IN2P3 School Of Statistics | Florian Ruppin (Université Claude Bernard Lyon 1)**\n",
"\n",
"---\n",
"\n",
"## Context\n",
"\n",
"We analyse a (simulated) **density profile of a galaxy cluster**, in a regime where a secondary structure is merging with the main halo. The cluster is modelled as a *generalised NFW profile* and the merging substructure as a *Gaussian bump*.\n",
"\n",
"The data come with **coloured measurement noise** (we are given its PSD), so the likelihood involves the full noise covariance matrix.\n",
"\n",
"## Goal\n",
"\n",
"Implement the **Metropolis-Hastings algorithm** from scratch and use it to perform Bayesian inference on the 5 free parameters of the model.\n",
"\n",
"## What is provided vs. what you have to code\n",
"\n",
"All the infrastructure is already coded for you (in the section *Provided infrastructure* below):\n",
"\n",
"- the noise covariance matrix from the measured PSD;\n",
"- the proposal function $g(\\vec{\\theta}\\,|\\vec{\\theta}_i)$;\n",
"- the `log_prior`, `log_likelihood`, and `log_probability` functions;\n",
"- the **Gelman-Rubin** convergence test;\n",
"- the integrated **autocorrelation time**.\n",
"\n",
"Your job is to fill in **four blocks**:\n",
"\n",
"1. The **model** $\\rho(r)$.\n",
"2. The **Metropolis-Hastings sampler**.\n",
"3. The **post-processing function** (burn-in, convergence check, thinning).\n",
"4. The **results function** (corner plot and credible intervals).\n",
"\n",
"Estimated duration: **~ 45 minutes**."
]
},
{
"cell_type": "markdown",
"id": "H91Yq9x-P-1p",
"metadata": {
"id": "H91Yq9x-P-1p"
},
"source": [
"---\n",
"\n",
"## 0. Setup\n",
"\n",
"Run the cells below to install `getdist` and import the libraries. Then upload the file `mydata_cluster.npz` provided on the Indico for this hands-on session."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "KopDq-QyP-1q",
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "KopDq-QyP-1q",
"outputId": "adaef4e0-214c-46a0-ef6c-40b45388024f"
},
"outputs": [
{
"output_type": "stream",
"name": "stdout",
"text": [
"\u001b[?25l \u001b[90m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001b[0m \u001b[32m0.0/836.0 kB\u001b[0m \u001b[31m?\u001b[0m eta \u001b[36m-:--:--\u001b[0m\r\u001b[2K \u001b[91m━━━━━━━━━━━━━━━━\u001b[0m\u001b[91m╸\u001b[0m\u001b[90m━━━━━━━━━━━━━━━━━━━━━━━\u001b[0m \u001b[32m348.2/836.0 kB\u001b[0m \u001b[31m9.4 MB/s\u001b[0m eta \u001b[36m0:00:01\u001b[0m\r\u001b[2K \u001b[91m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001b[0m\u001b[91m╸\u001b[0m \u001b[32m829.4/836.0 kB\u001b[0m \u001b[31m16.7 MB/s\u001b[0m eta \u001b[36m0:00:01\u001b[0m\r\u001b[2K \u001b[90m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001b[0m \u001b[32m836.0/836.0 kB\u001b[0m \u001b[31m10.6 MB/s\u001b[0m eta \u001b[36m0:00:00\u001b[0m\n",
"\u001b[?25h"
]
}
],
"source": [
"!pip install -q getdist"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "efReClkdP-1q",
"metadata": {
"id": "efReClkdP-1q"
},
"outputs": [],
"source": [
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"from scipy.linalg import toeplitz\n",
"from getdist import MCSamples, plots\n",
"\n",
"rng = np.random.default_rng(42) # global RNG, fixed seed for reproducibility"
]
},
{
"cell_type": "markdown",
"id": "kYolyx_NP-1q",
"metadata": {
"id": "kYolyx_NP-1q"
},
"source": [
"### Load the data\n",
"\n",
"On Google Colab, upload `mydata_cluster.npz` by dragging it into the file panel on the left."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "C81Wc3fcP-1r",
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "C81Wc3fcP-1r",
"outputId": "d64f2ec4-2dbe-4bde-d067-e01ecbfd1ba3"
},
"outputs": [
{
"output_type": "stream",
"name": "stdout",
"text": [
"Number of data points : 300\n",
"Radius range : 1.0 - 3000.0 kpc\n"
]
}
],
"source": [
"data = np.load('mydata_cluster.npz')\n",
"r_data = data['r'] # radii (kpc)\n",
"y_data = data['y'] # measured density profile\n",
"f_data = data['f'] # frequencies of the noise PSD\n",
"psd = data['psd'] # noise power spectral density\n",
"sigma_n = 1e-3 # noise standard deviation\n",
"\n",
"print(f'Number of data points : {len(r_data)}')\n",
"print(f'Radius range : {r_data.min():.1f} - {r_data.max():.1f} kpc')"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "lCEh1cvqP-1r",
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/",
"height": 506
},
"id": "lCEh1cvqP-1r",
"outputId": "de84072d-5ce4-403e-d71b-d3ee2fd4921d"
},
"outputs": [
{
"output_type": "display_data",
"data": {
"text/plain": [
""
],
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAxYAAAHpCAYAAAAf5apCAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjAsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvlHJYcgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAdONJREFUeJzt3Xl8VOXZ//HvmT3bJIRshB0MIIJEQRYVcaFStSLWBdEqLnWrgEprq9a6t/q0xaq4/bTuGz5oRR5qqYhQqeDGoqKyqCyyJARCMllnPb8/hhmykzBJJoHP+/XKK86Zc85c58wJnuvc93XfhmmapgAAAAAgBpZ4BwAAAACg8yOxAAAAABAzEgsAAAAAMSOxAAAAABAzEgsAAAAAMSOxAAAAABAzEgsAAAAAMSOxAAAAABAzEgsAAAAAMSOxAIBO7OSTT9bJJ5/c4fbVkW3evFmGYeiFF16IdygyDEN33313vMNodS+//LIGDRoku92utLQ0SfWvr470PQBoHSQWADqMF154QYZhyDAM/fe//633vmma6tmzpwzD0M9+9rM4RIhY/OlPf9K8efPiHUaHtnz5ct19990qKSmJdygHbd26dbr88svVv39/PfPMM3r66afjHRKAdmKLdwAAUJfL5dJrr72mE088sdby//znP9q2bZucTmecIkMs/vSnP+n888/XpEmT4h1Kh1FVVSWbbf//ipcvX6577rlHl19+efRJf2ezdOlShUIhPfLIIzriiCOiy9977704RgWgPdBiAaDDOfPMMzV37lwFAoFay1977TUNHz5cOTk5cYqs9VRUVMQ7hENCdXW1QqFQvMM4aC6Xq1Zi0RG19FrdtWuXJNVLjBwOhxwOR2uFBaADIrEA0OFMmTJFe/bs0aJFi6LLfD6f3nzzTV188cUNbhMKhfTwww/rqKOOksvlUnZ2tq699lrt3bu31nrvvPOOzjrrLOXm5srpdKp///667777FAwGa623ceNGnXfeecrJyZHL5VKPHj100UUXqbS0VFLT/cPr9pu/++67ZRiGvvnmG1188cXq0qVLrdaYV155RcOHD1dCQoLS09N10UUX6ccff6y336efflr9+/dXQkKCRo4cqWXLlh3wXNb0yiuvaOTIkUpMTFSXLl100kknNfkUOdI1bfPmzbWWL126VIZhaOnSpdFlBzpfhmGooqJCL774YrS72+WXXx7dfvv27bryyiuVnZ0tp9Opo446Ss8991yDnztnzhzdcccd6t69uxITE+XxeBo9hpKSEl1++eVKTU1VWlqapk6d2mg3o3Xr1un8889Xenq6XC6XRowYofnz5zd4Tj766CPNnDlTmZmZSkpK0rnnnquioqJa637++eeaMGGCMjIylJCQoL59++rKK6+stU7Na+Xuu+/WLbfcIknq27dv9Dxt3rxZ48aN07BhwxqMe+DAgZowYUKj50CS+vTpo5/97Gd67733lJ+fL5fLpcGDB+sf//hHg8f3n//8R7/61a+UlZWlHj16RN9/4okndNRRR8npdCo3N1c33HBDrfPZp08f3XXXXZKkzMzMWsfX3Bqe5nwPADqmjv2YBMBhqU+fPhozZoxef/11nXHGGZKkf/3rXyotLdVFF12kRx99tN421157rV544QVdccUVmjFjhjZt2qTHHntMq1ev1kcffSS73S4pfOOUnJysmTNnKjk5WR988IHuvPNOeTwe/eUvf5EUTmImTJggr9er6dOnKycnR9u3b9eCBQtUUlKi1NTUgzquCy64QHl5efrTn/4k0zQlSX/84x/1hz/8QRdeeKF++ctfqqioSLNnz9ZJJ52k1atXR5/6Pvvss7r22mt1/PHH66abbtIPP/ygiRMnKj09XT179jzgZ99zzz26++67dfzxx+vee++Vw+HQJ598og8++ECnn376QR1PRHPO18svv6xf/vKXGjlypK655hpJUv/+/SVJhYWFGj16tAzD0LRp05SZmal//etfuuqqq+TxeHTTTTfV+rz77rtPDodDv/nNb+T1eht9Cm6aps455xz997//1XXXXacjjzxSb7/9tqZOnVpv3a+//lonnHCCunfvrltvvVVJSUn63//9X02aNElvvfWWzj333FrrT58+XV26dNFdd92lzZs36+GHH9a0adP0xhtvSAo/tT/99NOVmZmpW2+9VWlpadq8eXO9G/mafv7zn2vDhg16/fXX9be//U0ZGRmSwjfol156qa6++mqtXbtWQ4YMiW7z2WefacOGDbrjjjsO8C2Fk7/Jkyfruuuu09SpU/X888/rggsu0MKFC/WTn/yk1rq/+tWvlJmZqTvvvDPaYnH33Xfrnnvu0fjx43X99ddr/fr1evLJJ/XZZ59F/8YefvhhvfTSS3r77bf15JNPKjk5WUcfffQBY4to6fcAoIMxAaCDeP75501J5meffWY+9thjZkpKillZWWmapmlecMEF5imnnGKapmn27t3bPOuss6LbLVu2zJRkvvrqq7X2t3DhwnrLI/ur6dprrzUTExPN6upq0zRNc/Xq1aYkc+7cuY3GumnTJlOS+fzzz9d7T5J51113RV/fddddpiRzypQptdbbvHmzabVazT/+8Y+1ln/11VemzWaLLvf5fGZWVpaZn59ver3e6HpPP/20KckcN25co3Gapmlu3LjRtFgs5rnnnmsGg8Fa74VCoeh/jxs3rta+It/Hpk2bam2zZMkSU5K5ZMkS0zSbd75M0zSTkpLMqVOn1lt+1VVXmd26dTN3795da/lFF11kpqamRr+zyOf269evwe+xrnnz5pmSzD//+c/RZYFAwBw7dmy97+60004zhw4dGr0GTDN8bo4//ngzLy8vuixyTsaPH1/r3N18882m1Wo1S0pKTNM0zbfffjt6LTel7rXyl7/8pcFzXlJSYrpcLvN3v/tdreUzZswwk5KSzPLy8iY/p3fv3qYk86233oouKy0tNbt162Yec8wx9Y7vxBNPNAOBQHT5rl27TIfDYZ5++um1rqHHHnvMlGQ+99xz0WWR672oqKhWDHWvr4b+hpr7PQDomOgKBaBDuvDCC1VVVaUFCxaorKxMCxYsaLQb1Ny5c5Wamqqf/OQn2r17d/Rn+PDhSk5O1pIlS6LrJiQkRP+7rKxMu3fv1tixY1VZWal169ZJUrRF4t///rcqKytb7Ziuu+66Wq//8Y9/KBQK6cILL6wVd05OjvLy8qJxf/7559q1a5euu+66Wk/nI118DmTevHkKhUK68847ZbHU/mffMIyYjyuW82Wapt566y2dffbZMk2z1nmYMGGCSktLtWrVqlrbTJ06tdb32Jh3331XNptN119/fXSZ1WrV9OnTa61XXFysDz74QBdeeGH0mti9e7f27NmjCRMmaOPGjdq+fXutba655ppa527s2LEKBoPasmWLpP31BQsWLJDf72/ROWlIamqqzjnnHL3++uvR1q5gMKg33nhDkyZNUlJS0gH3kZubW+uJv9vt1mWXXabVq1eroKCg1rpXX321rFZr9PX7778vn8+nm266qdY1dPXVV8vtduuf//xnrId4UN8DgI6FxAJAh5SZmanx48frtdde0z/+8Q8Fg0Gdf/75Da67ceNGlZaWKisrS5mZmbV+ysvLo8WkUrirxbnnnqvU1FS53W5lZmbqF7/4hSRF6wH69u2rmTNn6u9//7syMjI0YcIEPf7449H3D1bfvn3rxW2apvLy8urF/e2330bjjtys5uXl1drebrerX79+B/zc77//XhaLRYMHD44p/sbEcr6KiopUUlKip59+ut45uOKKKySp1vcX+bzm2LJli7p166bk5ORaywcOHFjr9XfffSfTNPWHP/yhXgyReoG6MfTq1avW6y5dukhStKZn3LhxOu+883TPPfcoIyND55xzjp5//nl5vd5mxd6Qyy67TFu3bo3W1rz//vsqLCzUpZde2qztjzjiiHqJ5IABAySpXh1N3XMcuQbrnjuHw6F+/fpF34/FwXwPADoWaiwAdFgXX3yxrr76ahUUFOiMM85odPjNUCikrKwsvfrqqw2+n5mZKSlcyDtu3Di53W7de++96t+/v1wul1atWqXf/e53tUYXmjVrli6//HK98847eu+99zRjxgw98MAD+vjjj9WjR49Gn/TXLQKvqe5T9lAoJMMw9K9//avW0+GIujfE7a0lx3ig89WYyDn/xS9+0WDtg6R6ffSb01rREpEYfvOb3zRaBF1z2FRJDX5fkqKtCYZh6M0339THH3+s//u//9O///1vXXnllZo1a5Y+/vjjg/puJ0yYoOzsbL3yyis66aST9MorrygnJ0fjx49v8b4OpLXPcXMczPcAoGMhsQDQYZ177rm69tpr9fHHH0eLYhvSv39/vf/++zrhhBOavCFaunSp9uzZo3/84x866aSToss3bdrU4PpDhw7V0KFDdccdd2j58uU64YQT9NRTT+n++++PPqGuO8JQS57c9u/fX6Zpqm/fvtEnxw3p3bu3pHALx6mnnhpd7vf7tWnTpkZHC6r5OaFQSN98843y8/ObHV9Lj7Gp8yU1nKhkZmYqJSVFwWCw1W+Qe/furcWLF6u8vLzWjfz69etrrRdp9bHb7a0ew+jRozV69Gj98Y9/1GuvvaZLLrlEc+bM0S9/+csG12+qa5rVatXFF1+sF154Qf/zP/+jefPm1euy1JRIi0DNz9iwYYOk8IAJTYlcg+vXr6/VSubz+bRp06ZWOW9t+T0AaB90hQLQYSUnJ+vJJ5/U3XffrbPPPrvR9S688EIFg0Hdd9999d4LBALRG+PIDVjkqbIUvjF64oknam3j8XjqzaExdOhQWSyWaFcWt9utjIwMffjhh7XWq7uvpvz85z+X1WrVPffcUyumSIx79uyRJI0YMUKZmZl66qmn5PP5ouu88MILzZqhedKkSbJYLLr33nvrzflQ93NriozaVPMYg8FgvZmUm3O+JCkpKalevFarVeedd57eeustrV27tl4MdYdwbYkzzzxTgUBATz75ZK34Z8+eXWu9rKwsnXzyyfp//+//aefOna0Sw969e+ud20hS11R3qEitRGPf66WXXqq9e/fq2muvVXl5ebQbX3Ps2LFDb7/9dvS1x+PRSy+9pPz8/APODTN+/Hg5HA49+uijtY7r2WefVWlpqc4666xmx9GYtvgeALQvWiwAdGiNdY+pady4cbr22mv1wAMPaM2aNTr99NNlt9u1ceNGzZ07V4888ojOP/98HX/88erSpYumTp2qGTNmyDAMvfzyy/VuAD/44ANNmzZNF1xwgQYMGKBAIKCXX345ehMc8ctf/lIPPvigfvnLX2rEiBH68MMPo0+Am6N///66//77ddttt2nz5s2aNGmSUlJStGnTJr399tu65ppr9Jvf/EZ2u13333+/rr32Wp166qmaPHmyNm3apOeff75ZNRZHHHGEfv/73+u+++7T2LFj9fOf/1xOp1OfffaZcnNz9cADDzS43VFHHaXRo0frtttuU3FxsdLT0zVnzpx6SURzz9fw4cP1/vvv66GHHlJubq769u2rUaNG6cEHH9SSJUs0atQoXX311Ro8eLCKi4u1atUqvf/++youLm72Oa3p7LPP1gknnKBbb71Vmzdvjs7b0FDtx+OPP64TTzxRQ4cO1dVXX61+/fqpsLBQK1as0LZt2/TFF1+06LNffPFFPfHEEzr33HPVv39/lZWV6ZlnnpHb7daZZ57Z6HbDhw+XJP3+97/XRRddJLvdrrPPPjuacBxzzDEaMmSI5s6dqyOPPFLHHntss2MaMGCArrrqKn322WfKzs7Wc889p8LCQj3//PMH3DYzM1O33Xab7rnnHv30pz/VxIkTtX79ej3xxBM67rjjWpTgNKW1vwcA7az9B6ICgIbVHG62KXWHm414+umnzeHDh5sJCQlmSkqKOXToUPO3v/2tuWPHjug6H330kTl69GgzISHBzM3NNX/729+a//73v2sNn/rDDz+YV155pdm/f3/T5XKZ6enp5imnnGK+//77tT6vsrLSvOqqq8zU1FQzJSXFvPDCC81du3Y1Otxs3eE3I9566y3zxBNPNJOSksykpCRz0KBB5g033GCuX7++1npPPPGE2bdvX9PpdJojRowwP/zww3pDeDblueeeM4855hjT6XSaXbp0MceNG2cuWrQo+n5D+/r+++/N8ePHm06n08zOzjZvv/12c9GiRQd1vtatW2eedNJJZkJCgimp1tCzhYWF5g033GD27NnTtNvtZk5OjnnaaaeZTz/9dHSdyHCzBxrWtqY9e/aYl156qel2u83U1FTz0ksvjQ6PW3eo4O+//9687LLLzJycHNNut5vdu3c3f/azn5lvvvlmdJ3GrtG6Q/CuWrXKnDJlitmrVy/T6XSaWVlZ5s9+9jPz888/r7Vd3WvFNE3zvvvuM7t3725aLJYGh57985//bEoy//SnPzX7PET+Zv7973+bRx99tOl0Os1BgwbVO5cH+ht87LHHzEGDBpl2u93Mzs42r7/+enPv3r211olluFnTbN73AKBjMkyziXZwAADQoTzyyCO6+eabtXnz5nqjUzWmT58+GjJkiBYsWNDG0QE4nFFjAQBAJ2Gapp599lmNGzeu2UkFALQXaiwAAOjgKioqNH/+fC1ZskRfffWV3nnnnXiHBAD1kFgAANDBFRUV6eKLL1ZaWppuv/12TZw4Md4hAUA91FgAAAAAiBk1FgAAAABiRleoBoRCIe3YsUMpKSlNzoIKAAAAHMpM01RZWZlyc3NlsTTdJkFi0YAdO3aoZ8+e8Q4DAAAA6BB+/PFH9ejRo8l1SCwakJKSIil8At1ud5yjAQAAAOLD4/GoZ8+e0fvjppBYNCDS/cntdpNYAAAA4LDXnPIAircBAAAAxIzEAgAAAEDMSCwAAAAAxIwaCwAAAHR6wWBQfr8/3mF0Ona7XVartVX2RWIBAACATss0TRUUFKikpCTeoXRaaWlpysnJiXn+NhILAAAAdFqRpCIrK0uJiYlMbtwCpmmqsrJSu3btkiR169Ytpv2RWAAAAKBTCgaD0aSia9eu8Q6nU0pISJAk7dq1S1lZWTF1i6J4GwAAAJ1SpKYiMTExzpF0bpHzF2uNCokFAAAAOjW6P8Wmtc4fiQUAAAAgyR8M6Z012+UPhuIdSqdEYgEAAABICgRNzV+zQ4GgGe9QOiUSCwAAAKCDOPnkk3XTTTfFO4yDEvfE4vHHH1efPn3kcrk0atQoffrpp02uP3fuXA0aNEgul0tDhw7Vu+++W+v98vJyTZs2TT169FBCQoIGDx6sp556qi0PAQAAAGh3S5culWEYHWYOj7gmFm+88YZmzpypu+66S6tWrdKwYcM0YcKE6Fi6dS1fvlxTpkzRVVddpdWrV2vSpEmaNGmS1q5dG11n5syZWrhwoV555RV9++23uummmzRt2jTNnz+/vQ4LAAAAndDeCp/KqgPaW+GLdyidUlwTi4ceekhXX321rrjiimjLQmJiop577rkG13/kkUf005/+VLfccouOPPJI3XfffTr22GP12GOPRddZvny5pk6dqpNPPll9+vTRNddco2HDhjXZEuL1euXxeGr9AAAA4PCxbGORZsxZrW92lmrGnNVatrGozT+zoqJCl112mZKTk9WtWzfNmjWr1vsvv/yyRowYoZSUFOXk5Ojiiy+OPoDfvHmzTjnlFElSly5dZBiGLr/8cknSwoULdeKJJyotLU1du3bVz372M33//fdtfjxxSyx8Pp9Wrlyp8ePH7w/GYtH48eO1YsWKBrdZsWJFrfUlacKECbXWP/744zV//nxt375dpmlqyZIl2rBhg04//fRGY3nggQeUmpoa/enZs2eMRwcAAIDOorjCp9mLv1NJpV92i0UllX7N/uA7Fbdxy8Utt9yi//znP3rnnXf03nvvaenSpVq1alX0fb/fr/vuu09ffPGF5s2bp82bN0eTh549e+qtt96SJK1fv147d+7UI488IimcsMycOVOff/65Fi9eLIvFonPPPVehUNuOdhW3mbd3796tYDCo7OzsWsuzs7O1bt26BrcpKChocP2CgoLo69mzZ+uaa65Rjx49ZLPZZLFY9Mwzz+ikk05qNJbbbrtNM2fOjL72eDxxTS78wZDe/WqnzhzaTXZr3MtgAAAADmmFnmqVVfuVmexQSaUv/LvKr11l1UpPcrTJZ5aXl+vZZ5/VK6+8otNOO02S9OKLL6pHjx7Rda688srof/fr10+PPvqojjvuOJWXlys5OVnp6emSpKysLKWlpUXXPe+882p91nPPPafMzEx98803GjJkSJscj9QBirdb2+zZs/Xxxx9r/vz5WrlypWbNmqUbbrhB77//fqPbOJ1Oud3uWj/xVOTx6tWPt6rI441rHAAAAIeDbLdLKS67isp9CoZMFZX7lJJgV1aKq80+8/vvv5fP59OoUaOiy9LT0zVw4MDo65UrV+rss89Wr169lJKSonHjxkmStm7d2uS+N27cqClTpqhfv35yu93q06dPs7aLVdwSi4yMDFmtVhUWFtZaXlhYqJycnAa3ycnJaXL9qqoq3X777XrooYd09tln6+ijj9a0adM0efJk/fWvf22bA2ll8ejfBwAAcDhLT3Jo+mlHKC3RLn8opLREu6afekSbtVY0R0VFhSZMmCC3261XX31Vn332md5++21J4ZKCppx99tkqLi7WM888o08++USffPJJs7aLVdwSC4fDoeHDh2vx4sXRZaFQSIsXL9aYMWMa3GbMmDG11pekRYsWRdf3+/3y+/2yWGofltVqbfM+Za0hXv37AAAADndj8zL16EXH6KhuqXr0omM0Ni+zTT+vf//+stvt0Zt+Sdq7d682bNggSVq3bp327NmjBx98UGPHjtWgQYPqjZzqcIQTn2AwGF22Z88erV+/XnfccYdOO+00HXnkkdq7d2+bHktE3GospPDQsFOnTtWIESM0cuRIPfzww6qoqNAVV1whSbrsssvUvXt3PfDAA5KkG2+8UePGjdOsWbN01llnac6cOfr888/19NNPS5LcbrfGjRunW265RQkJCerdu7f+85//6KWXXtJDDz0Ut+Nsrnj07wMAAEBYlySHkl02dWmH+67k5GRdddVVuuWWW9S1a1dlZWXp97//ffQBea9eveRwODR79mxdd911Wrt2re67775a++jdu7cMw9CCBQt05plnKiEhQV26dFHXrl319NNPq1u3btq6datuvfXWNj8eKc41FpEuSnfeeafy8/O1Zs0aLVy4MFqgvXXrVu3cuTO6/vHHH6/XXntNTz/9tIYNG6Y333xT8+bNq1WEMmfOHB133HG65JJLNHjwYD344IP64x//qOuuu67dj6+l4tG/DwAAAPHxl7/8RWPHjtXZZ5+t8ePH68QTT9Tw4cMlSZmZmXrhhRc0d+7c6D1t3a793bt31z333KNbb71V2dnZmjZtmiwWi+bMmaOVK1dqyJAhuvnmm/WXv/ylXY7HME3TbJdP6kQ8Ho9SU1NVWlra7oXcyzYW6ZH3N2pdgUeDcty6cXxemzfFAQAAdEbV1dXatGmT+vbtK5cr9gexVb6gpr22So9dfKwSHNZWiLBzaOo8tuS+OK5doVDf2LxM9c9I1s1vrNHfJucrt0tCvEMCAAA4LNishibm58pmNeIdSqdEYtEBZbqdunh0L2W6nfEOBQAA4LBht1p0Tn73eIfRaZFYdEBc1AAAAOhsDrkJ8gAAAAC0PxILAAAAdGqdYb6yjqy1zh9doQAAANApORwOWSwW7dixQ5mZmXI4HDIMCq+byzRN+Xw+FRUVyWKxRCfcO1gkFgAAAOiULBaL+vbtq507d2rHjh3xDqfTSkxMVK9evaKT8x0sEgsAAAB0Wg6HQ7169VIgEFAwGIx3OJ2O1WqVzWZrlZYeEgsAAAB0aoZhyG63y263xzuUwxrF2wAAAABiRmIBAAAAIGYkFgAAAABiRmIBAAAAIGYkFgAAAABiRmIBAAAAIGYkFgAAAABiRmIBAAAAIGYkFgAAAABiRmIBAAAAIGYkFgAAAABiRmIBAAAAIGYkFgAAAABiRmIBAAAAIGYkFgAAAABiRmIBAAAAIGYkFgAAAABiRmIBAAAAIGYkFgAAAABiRmIBAAAAIGYkFgAAAABiRmIBAAAAIGYkFgAAAABiRmIBAAAAIGYkFgAAAABiRmIBAAAAIGYkFgAAAABiRmIBAAAAIGYkFgAAAABiRmIBAAAAIGYkFgAAAABiRmIBAAAAIGYkFgAAAABiRmIBAAAAIGYkFgAAAABiRmIBAAAAIGYkFgAAAABiRmIBAAAAIGYkFgAAAABiRmIBAAAAIGYkFgAAAABiRmIBAAAAIGYkFgAAAABiRmIBAAAAIGYkFgAAAABiRmIBAAAAIGYkFgAAAABiRmIBAAAAIGYkFgAAAABiRmIBAAAAIGYkFgAAAABiRmIBAAAAIGYkFgAAAABiRmIBAAAAIGYkFgAAAABiRmIBAAAAIGYkFgAAAABiRmIBAAAAIGYkFgAAAABiRmIBAAAAIGYkFgAAAABiRmIBAAAAIGYkFgAAAABiRmIBAAAAIGYkFgAAAABiRmIBAAAAIGYkFgAAAABiFvfE4vHHH1efPn3kcrk0atQoffrpp02uP3fuXA0aNEgul0tDhw7Vu+++W2+db7/9VhMnTlRqaqqSkpJ03HHHaevWrW11CAAAAMBhL66JxRtvvKGZM2fqrrvu0qpVqzRs2DBNmDBBu3btanD95cuXa8qUKbrqqqu0evVqTZo0SZMmTdLatWuj63z//fc68cQTNWjQIC1dulRffvml/vCHP8jlcrXXYQEAAACHHcM0TTNeHz5q1Cgdd9xxeuyxxyRJoVBIPXv21PTp03XrrbfWW3/y5MmqqKjQggULostGjx6t/Px8PfXUU5Kkiy66SHa7XS+//HKz4/B6vfJ6vdHXHo9HPXv2VGlpqdxu98EeHgAAANCpeTwepaamNuu+OG4tFj6fTytXrtT48eP3B2OxaPz48VqxYkWD26xYsaLW+pI0YcKE6PqhUEj//Oc/NWDAAE2YMEFZWVkaNWqU5s2b12QsDzzwgFJTU6M/PXv2jO3gAAAAgMNM3BKL3bt3KxgMKjs7u9by7OxsFRQUNLhNQUFBk+vv2rVL5eXlevDBB/XTn/5U7733ns4991z9/Oc/13/+859GY7nttttUWloa/fnxxx9jPDoAAADg8GKLdwCtKRQKSZLOOecc3XzzzZKk/Px8LV++XE899ZTGjRvX4HZOp1NOp7Pd4gQAAAAONXFrscjIyJDValVhYWGt5YWFhcrJyWlwm5ycnCbXz8jIkM1m0+DBg2utc+SRRzIqFAAAANCG4pZYOBwODR8+XIsXL44uC4VCWrx4scaMGdPgNmPGjKm1viQtWrQour7D4dBxxx2n9evX11pnw4YN6t27dysfAQAAAICIuHaFmjlzpqZOnaoRI0Zo5MiRevjhh1VRUaErrrhCknTZZZepe/fueuCBByRJN954o8aNG6dZs2bprLPO0pw5c/T555/r6aefju7zlltu0eTJk3XSSSfplFNO0cKFC/V///d/Wrp0aTwOEQAAADgsxDWxmDx5soqKinTnnXeqoKBA+fn5WrhwYbRAe+vWrbJY9jeqHH/88Xrttdd0xx136Pbbb1deXp7mzZunIUOGRNc599xz9dRTT+mBBx7QjBkzNHDgQL311ls68cQT2/34AAAAgMNFXOex6KhaMl4vAAAAcKjqFPNYAAAAADh0kFgAAAAAiBmJBQAAAICYkVgAAAAAiBmJBQAAAICYkVgAAAAAiBmJBQAAAICYkVgAAAAAiBmJBQAAAICYkVgAAAAAiBmJBQAAAICYkVgAAAAAiBmJBQAAAICYkVgAAAAAiBmJBQAAAICYkVgAAAAAiBmJBQAAAICYkVgAAAAAiBmJBQAAAICYkVgAAAAAiBmJBQAAAICYkVgAAAAAiBmJBQAAAICYkVgAAAAAiBmJBQAAAICYkVgAAAAAiBmJBQAAAICYkVgAAAAAiBmJBQAAAICYkVgAAAAAiBmJBQAAAICYkVgAAAAAiBmJBQAAAICYkVgAAAAAiBmJBQAAAICYkVgAAAAAiBmJBQAAAICYkVgAAAAAiBmJBQAAAICYkVgAAAAAiBmJBQAAAICYkVgAAAAAiBmJBQAAAICYkVgAAAAAiBmJBQAAAICYkVgAAAAAiBmJBQAAAICYkVgAAAAAiBmJBQAAAICYkVgAAAAAiBmJBQAAAICYkVgAAAAAiBmJBQAAAICYkVgAAAAAiBmJBQAAAICYHVRiUVJSor///e+67bbbVFxcLElatWqVtm/f3qrBAQAAAOgcbC3d4Msvv9T48eOVmpqqzZs36+qrr1Z6err+8Y9/aOvWrXrppZfaIk4AAAAAHViLWyxmzpypyy+/XBs3bpTL5YouP/PMM/Xhhx+2anAAAAAAOocWJxafffaZrr322nrLu3fvroKCglYJCgAAAEDn0uLEwul0yuPx1Fu+YcMGZWZmtkpQAAAAADqXFicWEydO1L333iu/3y9JMgxDW7du1e9+9zudd955rR4gAAAAgI6vxYnFrFmzVF5erqysLFVVVWncuHE64ogjlJKSoj/+8Y9tESMAAACADq7Fo0KlpqZq0aJF+u9//6svv/xS5eXlOvbYYzV+/Pi2iA8AAABAJ2CYpmnGO4iOxuPxKDU1VaWlpXK73fEOBwAAAIiLltwXN6vF4tFHH232h8+YMaPZ6wIAAAA4NDSrxaJv3761XhcVFamyslJpaWmSwjNxJyYmKisrSz/88EObBNqeaLEAAAAAWnZf3Kzi7U2bNkV//vjHPyo/P1/ffvutiouLVVxcrG+//VbHHnus7rvvvlY5AAAAAACdS4trLPr3768333xTxxxzTK3lK1eu1Pnnn69Nmza1aoDxQIsFAAAA0AYtFjXt3LlTgUCg3vJgMKjCwsKW7g4AAADAIaDFicVpp52ma6+9VqtWrYouW7lypa6//nqGnAUAAAAOUy1OLJ577jnl5ORoxIgRcjqdcjqdGjlypLKzs/X3v/+9LWIEAAAA0MG1eIK8zMxMvfvuu9qwYYPWrVsnSRo0aJAGDBjQ6sEBAAAA6BxanFhEDBgwgGQCAAAAgKSDSCyuvPLKJt9/7rnnDjoYAAAAAJ1TixOLvXv31nrt9/u1du1alZSU6NRTT221wAAAAAB0Hi1OLN5+++16y0KhkK6//nr179+/VYICAAAA0Lm0eFSoBndisWjmzJn629/+dlDbP/744+rTp49cLpdGjRqlTz/9tMn1586dq0GDBsnlcmno0KF69913G133uuuuk2EYevjhhw8qNgAAAAAH1iqJhSR9//33DU6cdyBvvPGGZs6cqbvuukurVq3SsGHDNGHCBO3atavB9ZcvX64pU6boqquu0urVqzVp0iRNmjRJa9eurbfu22+/rY8//li5ubktjgsAAABA8xmmaZot2WDmzJm1XpumqZ07d+qf//ynpk6dqscee6xFAYwaNUrHHXdcdLtQKKSePXtq+vTpuvXWW+utP3nyZFVUVGjBggXRZaNHj1Z+fr6eeuqp6LLt27dr1KhR+ve//62zzjpLN910k2666aYGY/B6vfJ6vdHXHo9HPXv2bNbU5QAAAMChyuPxKDU1tVn3xS2usVi9enWt1xaLRZmZmZo1a9YBR4yqy+fzaeXKlbrttttq7W/8+PFasWJFg9usWLGiXnIzYcIEzZs3L/o6FArp0ksv1S233KKjjjrqgHE88MADuueee1oUOwAAAID9WpxYLFmypNU+fPfu3QoGg8rOzq61PDs7Ozr5Xl0FBQUNrl9QUBB9/T//8z+y2WyaMWNGs+K47bbbaiUrkRYLAAAAAM3T4hqLU089VSUlJfWWezyeDjHc7MqVK/XII4/ohRdekGEYzdrG6XTK7XbX+gEAAADQfC1OLJYuXSqfz1dveXV1tZYtW9aifWVkZMhqtaqwsLDW8sLCQuXk5DS4TU5OTpPrL1u2TLt27VKvXr1ks9lks9m0ZcsW/frXv1afPn1aFB8AAACA5ml2V6gvv/wy+t/ffPNNra5HwWBQCxcuVPfu3Vv04Q6HQ8OHD9fixYs1adIkSeH6iMWLF2vatGkNbjNmzBgtXry4ViH2okWLNGbMGEnSpZdeqvHjx9faZsKECbr00kt1xRVXtCg+AAAAAM3T7MQiPz9fhmHIMIwGuzwlJCRo9uzZLQ5g5syZmjp1qkaMGKGRI0fq4YcfVkVFRTQJuOyyy9S9e3c98MADkqQbb7xR48aN06xZs3TWWWdpzpw5+vzzz/X0009Lkrp27aquXbvW+gy73a6cnBwNHDiwxfEBAAAAOLBmJxabNm2SaZrq16+fPv30U2VmZkbfczgcysrKktVqbXEAkydPVlFRke68804VFBQoPz9fCxcujBZob926VRbL/h5bxx9/vF577TXdcccduv3225WXl6d58+ZpyJAhLf5sAAAAAK2jxfNYHA5aMl4vAAAAcKhq9Xks5s+frzPOOEN2u13z589vct2JEyc2P1IAAAAAh4RmtVhYLBYVFBQoKyurVrekejszDAWDwVYNMB5osQAAAADaoMUiFAo1+N8AAAAAIB3EPBYAAAAAUFezWiweffTRZu9wxowZBx0MAAAAgM6pWTUWffv2bd7ODEM//PBDzEHFGzUWAAAAQBvUWGzatKlVAgOA5vIHQ3r3q506c2g32a302gQAoKOL6f/WpmmKaTAAtIUij1evfrxVRR5vvEMBAADNcFCJxbPPPqshQ4bI5XLJ5XJpyJAh+vvf/97asQE4TC3bWKQZc1brm52lmjFntZZtLIp3SAAA4ACa1RWqpjvvvFMPPfSQpk+frjFjxkiSVqxYoZtvvllbt27Vvffe2+pBAjh8FFf4NHvxdyqp9Mtusaik0q/ZH3yno3JTlZ7kiHd4AACgES1OLJ588kk988wzmjJlSnTZxIkTdfTRR2v69OkkFgBiUuipVlm1X5nJDpVU+sK/q/zaVVZNYgEAQAfW4q5Qfr9fI0aMqLd8+PDhCgQCrRIUgMNXttulFJddReU+BUOmisp9SkmwKyvFFe/QAABAE1qcWFx66aV68skn6y1/+umndckll7RKUAAOX+lJDk0/7QilJdrlD4WUlmjX9FOPoLUCAIAOrsVdoaRw8fZ7772n0aNHS5I++eQTbd26VZdddplmzpwZXe+hhx5qnSgBHFbG5mWqf0aybn5jjf42OV+5XRLiHRIAADiAFicWa9eu1bHHHitJ+v777yVJGRkZysjI0Nq1a6PrGYbRSiECOBxlup26eHQvZbqd8Q4FAAA0Q7Nm3j7cMPM2AAAA0LL7YqazBQAAABCzFneFqq6u1uzZs7VkyRLt2rVLoVCo1vurVq1qteAAAAAAdA4tTiyuuuoqvffeezr//PM1cuRIaikAAAAAtDyxWLBggd59912dcMIJbREPAAAAgE6oxTUW3bt3V0pKSlvEAgAAAKCTanFiMWvWLP3ud7/Tli1b2iIeAAAAAJ1Qi7tCjRgxQtXV1erXr58SExNlt9trvV9cXNxqwQEAAADoHFqcWEyZMkXbt2/Xn/70J2VnZ1O8DQAAAKDlicXy5cu1YsUKDRs2rC3iAQAAANAJtbjGYtCgQaqqqmqLWAAAAAB0Ui1OLB588EH9+te/1tKlS7Vnzx55PJ5aPwAAAAAOP4ZpmmZLNrBYwrlI3doK0zRlGIaCwWDrRRcnHo9HqampKi0tldvtjnc4AAAAQFy05L64xTUWS5YsafS9r776qqW7AwAAAHAIaHGLRV1lZWV6/fXX9fe//10rV66kxQIAAAA4RLTkvrjFNRYRH374oaZOnapu3brpr3/9q0499VR9/PHHB7s7AAAAAJ1Yi7pCFRQU6IUXXtCzzz4rj8ejCy+8UF6vV/PmzdPgwYPbKkYAAAAAHVyzWyzOPvtsDRw4UF9++aUefvhh7dixQ7Nnz27L2AAAAAB0Es1usfjXv/6lGTNm6Prrr1deXl5bxgQAAACgk2l2i8V///tflZWVafjw4Ro1apQee+wx7d69uy1jAwAAANBJNDuxGD16tJ555hnt3LlT1157rebMmaPc3FyFQiEtWrRIZWVlbRknAAAAgA4spuFm169fr2effVYvv/yySkpK9JOf/ETz589vzfjiguFmAQAAgHYablaSBg4cqD//+c/atm2bXn/99Vh2BQAAAKATi3mCvEMRLRYAAABAO7ZYAAAAAIBEYgEAAACgFZBYAAAAAIgZiQUAAACAmJFYAAAAAIgZiQUAAACAmJFYAAAAAIgZiQUAAACAmJFYAAAAAIgZiQUAAACAmJFYAAAAAIgZiQUAAACAmJFYAAAAAIgZiQUAAACAmJFYAEA78AdDemfNdvmDoXiHAgBAmyCxAIB2EAiamr9mhwJBM96hAADQJkgsAKAd7K3wqaw6oL0VvniHAgCHFFqEOw4SCwBoY8s2FmnGnNX6ZmepZsxZrWUbi+IdEgAcMoo8Xr368VYVebzxDuWwR2IBAG2ouMKn2Yu/U0mlX3aLRSWVfs3+4DsV03IBADHjwU3HQmIBAG2o0FOtsmq/MpMdsloMZSY7VFbl166y6niHBgCdGg9uOh4SCwBoQ9lul1JcdhWV+xQMmSoq9yklwa6sFFe8QwOATo0HNx0PiQUAtKH0JIemn3aE0hLt8odCSku0a/qpRyg9yRHv0ACgU+PBTcdDYgEAbWxsXqYevegYHdUtVY9edIzG5mXGOyQA6PR4cNPx2OIdAAAcDjLdTl08upcy3c54hwIAh4yxeZnqn5Gsm99Yo79Nzldul4R4h3RYI7EAgHZgt1p0Tn73eIcBAIccHtx0HCQWAAAA6LR4cNNxUGMBAIc4ZqUFALQHEgsAOMQFgqbmr9mhQNCMdygA0G54qNL+SCwA4BC3t8KnsuqA9jJpFIBDUN0EIvK6yhfkoUo7I7EAgEPYso1FmjFntb7ZWaoZc1Zr2caieIcEAK2qbqts5HVRmZeHKu2MxAIADlHFFT7NXvydSir9slssKqn0a/YH36mY/8kCOAT4gyH9Y9U2vfrJFnn2zbj9zprt2lVWrW17q3TTG2t4qNLOGBUKAA5RhZ5qlVX7lZnsUEmlL/x73/98mUAKQGcXCJp64aPN2rynQsGQqV/8/ROFTFNul10bd5XJYbXIZjWiD1WOyk3l3742RosFAByist0upbjsKir3KRgyVVTuU0qCXVkprlb7DIojAcSDPxjSc//9Qd8XlcsbCMk0TW3bW6XtJdX6rqg8vJIhVftDSku0qWzfQxW0LRILADhEpSc5NP20I5SWaJc/FFJaol3TTz2iVZ/YFXm8evXjrSryeFttnwDQmMjDjHfWbNcjH2xUhS8obyCkSn9IkRJtQ6ZCpiRTCpmmNu2uVLLL1qoPVdAwukIBwCFsbF6m+mck6+Y31uhvk/OV2yWh1fa9bGORHnl/o9YVeDRjzmrdOD5PY/MyW23/AFBXkcerR9/fqK17K+VvZLSnQEgyJIUkmaZU7Q/q4pG96AbVDmixAIBDXKbbqYtH91Km29lq+6QwHEB7W7axSNe/ulI/7K5oNKmQpJApGYaUk+KU026RKenVT7ZSwN0OSCwA4BBnt1p0Tn532a2t909+zcJwq8VQZrKDPswA2kzkYcaOkio1Z1aKJIdFqYl2pSbY5bDy8KO9dIjE4vHHH1efPn3kcrk0atQoffrpp02uP3fuXA0aNEgul0tDhw7Vu+++G33P7/frd7/7nYYOHaqkpCTl5ubqsssu044dO9r6MADgsNEeheEAEFHoqdbeSq8q/cFmrV/mDWl7SZX6dU2SjYcf7SbuicUbb7yhmTNn6q677tKqVas0bNgwTZgwQbt27Wpw/eXLl2vKlCm66qqrtHr1ak2aNEmTJk3S2rVrJUmVlZVatWqV/vCHP2jVqlX6xz/+ofXr12vixInteVgAcEhrj8JwAIjIdrvktFnlD5gymrlNcYVf3+8u5+FHOzJM04zrPOejRo3Scccdp8cee0ySFAqF1LNnT02fPl233nprvfUnT56siooKLViwILps9OjRys/P11NPPdXgZ3z22WcaOXKktmzZol69etV73+v1yuvdP6KJx+NRz549VVpaKrfbHeshAsAha8feqjYpDAeAut79aodm/u8XqvbvH97akJrsGuWyWWS1GDqym5sBJg6Sx+NRampqs+6L49pi4fP5tHLlSo0fPz66zGKxaPz48VqxYkWD26xYsaLW+pI0YcKERteXpNLSUhmGobS0tAbff+CBB5Samhr96dmzZ8sPBgAOQ12SHEp22dQlTi0VzKMBHD5G98uo1+JgtxpNtmCkJznULyNJj150DElFO4hrYrF7924Fg0FlZ2fXWp6dna2CgoIGtykoKGjR+tXV1frd736nKVOmNJpl3XbbbSotLY3+/PjjjwdxNABw+LFZDU3Mz5XN2tzOCa0rEDQ1f80OBZoYIQbAoeFfa3fWq5EwDCkrxam+XRMb3MZTHdCIvumtOioeGhf3Gou25Pf7deGFF8o0TT355JONrud0OuV2u2v9AAAOrC1GnGqJvRU+lVUHtJeRXoBDWnGFT29+vk1mqO5DBENV/qBmXZivYT1SlZ3i3LdUshpSisumDzcUMRpUO4lrYpGRkSGr1arCwsJaywsLC5WTk9PgNjk5Oc1aP5JUbNmyRYsWLSJZAIBDzLKNRZoxZ7W+2VmqGXNWN2uMerpOAZ3T9pJKFXqqlZJgr7U8N9WlZJdNFov0i9G91b2LS26XTUlOq5x2q1x2iwo81SooZTSo9hDXxMLhcGj48OFavHhxdFkoFNLixYs1ZsyYBrcZM2ZMrfUladGiRbXWjyQVGzdu1Pvvv6+uXbu2zQEAAOLiYCfoo+sU0Dl1TXSqzBtQWZVfFiPcImExwl2henRJVK/0JE06prvOG95TednJSnTY5PMHtbOkSpL0341FPFBoB3HvCjVz5kw988wzevHFF/Xtt9/q+uuvV0VFha644gpJ0mWXXabbbrstuv6NN96ohQsXatasWVq3bp3uvvtuff7555o2bZqkcFJx/vnn6/PPP9err76qYDCogoICFRQUyOejGQwADgUHO0EfXaeAzunTTcXyB0PyB02ZZri+y2GzKD3JGR3q2m616JJRvXXT+AHKcbsUMCVfwJQ/ENKTS7/Xgi+Y06yt2eIdwOTJk1VUVKQ777xTBQUFys/P18KFC6MF2lu3bpXFsj//Of744/Xaa6/pjjvu0O233668vDzNmzdPQ4YMkSRt375d8+fPlyTl5+fX+qwlS5bo5JNPbpfjAgC0nYYm6EtPdjQ5Rv2yjUV65P2NWlfg0Yw5qxl6Eugk3v1qh+54Z628/pAMhZOKvMwkuew2PXrRMQ0OdR15eBCS5Aua8geD+sM7X8tpt+rMod3a9wAOI3Gfx6Ijasl4vQCA+KiZKAzKaXqM+uIKn657eaWKK3zaU+5V12Sn0pMdeuoXw5nUD+jAiit8uuzZT7ShsFzBUEimJIfVoh5dEnTe8B765dh+tQaPKK7w6aoXPtPXO0rlq9Pl0W41NCgnRS9eOYq/+xboNPNYAABwsMbmZerRi47RUd1SDzhG/cF2nQIQX4WeankDQTntFoUUrq0ImaYSHFZNPq5XvRHpCj3VKqn0ydLACNh2q0XVgRB/922IxAIA0Gllup26eHSvA45R31DXqZQEe5NdpwDEX7bbpS6JTiXarTJNKWRKDptV1598RIOtDtlul9ISHVID0+a5XXZ1SWy6yyRiQ2IBAOi0mjuPRnqSQ9NPO0JpiXb5QyGlJdqjBZ8AOq70JIdOHpSpCl9AUngUqF+M6tVonUR6kkMzTx+gnumJ0VYLQ5LTFu4+xd992yKxAADEVXvNLdGSrlMAOobiCp+WritSRrJTViOcIHy+ZW+TQ0uPzcvUi1eM1NDuqTq6h1s3nnaEju7B3317ILEAAMRVe84t0dyuUwA6hkh9VLbbqQSHVblpCSr3Bg5YJ5Hpdmrq8X009fi+uvLEfuFuULRUtLm4DzcLADi81ZxbIsFRf9jI1hTpOgWgc4jUR+0p98lusciQIXcz6qPsVot+fmwPSeFW0Yn5ubJZG6joRquixQIAEDfLNhZpxpzV+mZnqWbMWa1lG4viHRKADsIfDGnZxiJdf0r/mOqjmluLhdhxhgEAcVFc4dPsxd+ppNIvu8Wikkq/Zn/wXZN9pwEcPoo8Xr368VYNzEqJqT6qveq4QGIBAIgT5pYA0JiarZnTXl+l55dvUqLTelB1Eu1Zx3W4I7EAAMQFc0sAaEjd1szSKr/+9/MfNW5A5kHVSdSs40LbIrEAAMQFc0sAaEjd1kx3gk2+QEj9M5NbXCdBHVf7IrEAAMTN4TC3BP27gZap2ZpZ7Q/q2x1l8gZCmvXehhYlBtRxtT8SCwBAXHXUuSVaKyGIFKAWebytFBlwaIu0ZiY5rar0ByVJiQ6ryr2BFiUG1HG1PxILAEBcddShIFsjIaAbBnBwxuZl6tc/GSinzaoju7klGeqabG9RYkAdV/vrWP+KAwDQAbRGQkA3DCA2A7ulKDfNpeIKn6p8ARV6vC1KDKjjan8kFgAA1NBaCQHdMIDYZLtdOm94D+2p8CpoSrvLvTp5YGaLEoPDoY6rIyGxAACghpYkBE3VYdANA4hNcYVPS9cVKSPZKashZSQ7tXR9UYuT/I5ax3UoIrEAAKCGliQETdVh0A0DiE0kyc92O5XgsCnb7TyoVr+OWsd1KOIMAwBQQ3MTgubUYdANAzh42W6Xklw2bd5TKZvF0J5y/0G3+vmDIf1j1Ta9tWobQz+3IRILAADqOFBC0JI6jC5JDiW7bOpCSwXQIulJDl13Un9V+4PyB2Nr9SvyePXaJ1v11sptCgTNNogWkmSLdwAAAHRETfXLrlmHUVLpC//e10Wj7k2PzWpoYn6ubFajvUIHDhmDu7nVp2uSJOnRi45RbpeEFu9j2cYiPfL+Rn2zo1QWw9DCtTt17rE9WjtUiMQCAIAGRfplN6ShOoz0ZEeDXTSa2g+AxkUSgu92lctmMfTNTk+LE4tI6+KOkip5gyGFQtIf3vlaTrtVZw7t1kaRH77oCgUAQAtRmA20rVrdDa0WuRxW/b8Pvz+oYZ/3VnpV5g0oZEoWQ/IFgnpyKXPKtAUSCwAADgKF2UDbqdnd0GYx1K9rksqrAy0eESrb7ZLTZlW1PyjTlExJTrtV1YEQc8q0ARILAAAOEuPjA22jteaBSU9y6ORBmQqGwgXbIVNyWC1KS7Br5Za9jBDVykgsAAA4SIyPD7SN1upuWFzh06c/7FVmcjj5NyR5AyFdMrK3Pvh2FyNEtTL+JQQA4BDQ1CzgQGfUGt0NI12qemckymZIyU6rst1OdU12qKw6oL3UWbQqEgsAAA4BTc0CDnRWsXY3jHSp2lPuk8NukWEYqqgO6PZ5X+nrJia3xMEhsQAAoJNrzizgQGcSmSl7/hc7dObQbgfd3TDSpapLokMypUAwpJIqv4o8XoVCpvZW+hqd3BItR2IBAEAn1pJZwIGOKtKVr9IX0DtrtmtHSZVeXL5ZLy3fHHMr3Ni8TP1iVG9ZLBb5AiFV7/up8odkmtK2vZXaUVLVSkdyeCOxAACgE6s5LKfVYigz2aGyfbOAA51FIGhq/podqvaF9PxHm3XlC5/pq+2l+mp7qaY+/2lMrXDFFT698skW+QJBBerUav+wu0LF5V5uiFsJ5xEAgE6stYblBOJpb4VPZdUBfbmtRBsLy7Rtb7gFwTSlH/dW6m+LNhx0K1yhp1olVb7wkFAN8AZM3f72WroQtgISCwAA4qQ1RnJiFnB0dpEaoS+3l+j6V1aqwheM/k1YLYYMSXsrD74VLtvtkstmldnIn5kp6fvd5brrna9V6KGlLxYkFgAAxElrjeTELODorCI1QrvLvfIFQtGEImSGf4IhU6akLokH3wqXnuTQr07pL6fd2ug6ZdUBbS+p1IaCsoP6DISRWAAAEAetPZJTlySHkl02daGlAp1ApLVue0mlPNV++YOmQqaU4LTK0P5eS4Yh9UpP1M0/GRBTK9yZQ3N1/zlDlNBEcuENmLpr/td0iYoBiQUAAO2sLUZyslkNTczPlc3aSEdyoAOJFGt3TXTKZbOqyBOuqajwBmUq3D3JkNSzS4JeuHxkq7TCnTWsm351cj85mvgb2VFaFVM9x+GOxAIAgHbWFiM52a0WnZPf/aDH+wfaU6RY+9NNxdpVXi3fvvqH0L5RmyxGuL5ib5W/1T7TbrVoYn53pTgbbrWwtUI9x+GOf30AAGhn8R7JqTWKxoGDFekG+PWOEv1+3lcqKvPWG7DJkGS1SIl2m4orW6f1YNnGIv3yxc9VXBlo8P3WqOc43JFYAADQzuI9klNrFY0DLVWzG6DFsMgbCEVv6GsKmpIvYCojyaHctIRW+dyH3tugH4srG12nteo5DmckFgAAxEG8RnJq7aJxoCVqdgO0Ww057ZZo96e6TEll3oC+3lHaKp9bUumTxTBktYS7WtVtJbFbDV0/rj+jqsWAxAIAgDhp75Gc2qJoHGiJSDfAXeVeVQdCclottW7wrcb+m/5wnYVa5RrNdruUluhQyDSjQ9lG8hljXwD+oKkXl2/h7yEGJBYAAMRJe4/k1NKicWox0Noi3QCTHFZV+YLyBkLKSnHKIslukRxWi0wznFQkOW3KTnHFPLBB5HNnnj5APdMTJYUTF5vFkGFIVoU/z2YxVFLto3A7BiQWAADESXuP5NTSonFqMdBWzH3tFMGQqS5JDiU6rerVNUnZqS4lOKzqkuSQ1TBadWCDsXmZevGKkRraPVWDu6VoQHaybBZDISmczFgMpSc6KNyOAYkFAACHiZYUjbd2LQatH5D2d8er9AXkshsKmqY276mQbV9/JMMw1LdronqlJyrQBgMbZLqd+sXo3spyuzTzJwPVp2vSvs+lcLs1kFgAANABNXUjHstNemNF4zX32Ra1GLR+QNrfHS8r2SmX3SabxVAwaMobDCk1wS63y6bLT+irv03Ob5OBDexWiyYd012TjumucQP3t2AM7Z7aahPxHc5ILAAA6ICauhGPzFocCDYynE4NDSUhmW6nLh7dS5luZ4P7bO0J/BiJChGR7nhbiivlqfaryh9SIGQqK8WpO84aLKfNqhP6Zyg3LaHeNdpaanZBzHQ7NfX4Prrs+D5t8lmHGxILAAA6mAPdiEdmLd7bjBaEhpKQhmo7au6zNSfwYyQq1JSe5NDlJ/SWpzqgmg1uO0qq9eu5X0Sv+Y9/2NMu9Ud2q0U/P7aHzju2B7PWtwLOIAAAHciBbsRb+vS/OUlI3X1+vaO01Sbwa+3WD3R+fTOSw9dDjaFlAyFT20uqZDMMks9OjMQCAIAOpKkb8ZY+/W9OEtLYPo/KTW2VCfxas/UDh4b0JId8QVOGYchQeD4JY1+SIYPkszMjsQAAoANp6ka8JU//m5uENLXPuhP4HUzReEtGomoPdY+B0araXzBoymW3yGkLz7odmVNCkmSSfHZmJBYAAHQgTd2It+Tpf3OTkKb2WXcCv7r1Gs29SW9sJKp4qHsMLSmER+wirWiFnmplJjvlslvkcljVu2tSeIhZ04x78tlcJKX1kVgAANDBNHYj3pKn/81NQpraZ90i77r1Gi25SW9oJKp4qHsMLSmER2zqtqJZrYZsVkN5Wcl68YqReuHykXFLPit9Af1h3lpV+gLN3oaktD4SCwAAOqDGbsSb+/S/JUlIc/bZUL1GS27S23uW8YbUPYYnln7HMLjtqG4rWkayQ8lOuxw2i7okOeKafO4sqdb/fbFdO0uaX9dBUlqfYZomaVYdHo9HqampKi0tldvtjnc4AADU4g+G9O5XO3Xm0G4HvFHfsbdKN7+xRn+bnK/cLgkHtc/iCp+ue3mliit82lPuVddkp6xWQwk2q74rKtOgHLdOPTJLH3y7S+sKPBqU49aN4/NifurckuM8kLrHkJboUFGFV1nJDu2t8KtrslPpyQ499YvhHb4LTmfV0HXUJcmuScd014UjesYt6Vy2sUh//tc6fb3Do6Ny3frtGYMOeO0uWb9L98z/WjtLq9Ut1aW7Jh6lUwZmtVPE7asl98W0WAAA0Mm05Ol/c58CN7XPuk+a0xJs2rqnUsWVXtktFu0p9+mJpd+ruMJbq1C80FPd4tnDay5vza4mdY8hxWWVPxBSisserUHxVPn1xmdb6TPfQs2tNWioFW3GaXm6ZFTvuCUVxRU+3b/gW23cVa6QpI27ynX/P7+tNchB3eMrrvDp8X3Xty8Q0ra9lXp40QaGxxWJBQAAh7TW6IJUt16jsCw8G3hWirPRm/SyKr82FJS1ePbwmjOOt2ZXk7rHUFYdlN1mUVm1P1qD4rRZtOjrwgbjRVhDSURLEsCOVMgvSRsKy7S1uFLGvkGpDEPauqdSz3/0Q/QYa16TUjhJ3VFSpWp/SKakQEj6ZqdHC9fubNXYOmNxOIkFAABoUt0nzelJDvXqmqCSykCjN+khSbPe29Ci2cNr1kBMff5TTX3+01arf6h7DF2THfrVyf2VnuSUPxSSYUiVvqA27Cqj3qIJdW+ypebXGkRulJNdtlrDGHcEphSdU8M0Q1r09S4VebwN1hZ9ta1EO0urVTON8gdNPbPsh1ZttWjoXHd0JBYAAOCAaj5pfvKS4brjrMGN3qQnOa2STJV7A82ePbzmiEFWw9CPxZX6cW9lq87EXPdp+a9OPkKPXnSMBmSlKMlhk2mqWRMPHq4a+t6WrN+lS579RF83IwGMtGxIqjWMcTwNyE5Rr64JkhlOKoJBUzIM/bi3Ute/ulL3//ObWnPBPLRog577aLMainxnSbU2Fpa1SlzNmdyyIyKxAAAAzVKzXqOxm/SjuqXq1z8ZKIuMFs0evqGwLFoDEZ6FOTwrc2vNxBx5Wp6WZK9Vc5Lpdmr8UdmqDgSbNfFgR9Ne3WUa+t4eWrRBjyzaqEJPdYO1NZW+QK3YIi0b5dWBuI8QFpGe5NAdZw3WgOwUGZIMi6EuiQ45rBYVV/i0dU+V0hJt0euipMKnSl9QTnvt2PdNGq4GM44Wqnuu91b6dPf8r1Xo6fjXY/y/UQAA0CnUrdeoWxgeeT2wW0qLZw83DEW3MU0pZJrhriatNBNz5Gm5IUPn5HeXJL2zZrsk6aLjesndjDk/OmKf9/bqLtPQ91ZS4dOeCq9C+7qS1a2tWbdzf41NR34CPzYvU7MvPkZZKQ6Zpqkct0tWi6HslPB1vavMG70u0pIc6prskNtlr7UPm9VQr65JystKiTmeuufanWDTjpIqbShondaQtkRiAQAADkrdRCPyOtvtavHs4XlZKdFtgqapnumJrToTc1OT+zU250eKy1Yrkehofd7b82a9oe/NYjW0p8Kn6oApT5VfW/ZWRmtrvtxeokuf/URfbi9psEtRR+tq9sPuClktFvmDpr7e6VF1IKiSqoB6dU1UeqIzel3M/MkAzfzJAPVMT5TTakQbKI7MSdEdZx3ZKkMV1zzX1f6gvt1RJm8gpFnvbehQCVlDSCwAAECrO5jZw2tu0xozMUdaGJas33XAyf0airdmIlH3Jn7J+l2t2nrR0taQxrqUNfdm3R8M6R+rtumtVdua9Zl1v7ckp1U2i6EUp02SFDKl0kq/AsGQtpdUyh8MyRcMyR8IaXe5t16Xoo7U1ay4wqcH/7V/iFlfIKQKb1DJTpvuOOtIPXHJsTqqW6pmXThMxRU+je7XVY9edIyO7pmmaaf0V7+MRM2++NhWGeHKHwxp6fpdyu+VpkSHRRW+oCQp0WFVuTfQ4RKyumzxDgAAAByauiQ5Ghz9Z2xepvpnJDc4cV/d7lWxzMQcCJp68/NtqvIHa92A3//Pb6OT+82Yszo6mV/NeJdtLNIj72/UugKPrn91paoDQQWD+4u7H//gOyXYrTp9cI7s1oM/RzVjnb9mR7P3V7O7TEmlL/x73816c56aB4Km3ly5Tb5ASGP6dm1y8sSImt/bDaccod/P+0oV/vCNr2FICXZr+KbcF5RpSgkOqyp8gXA3szKvCj1eBUKmfthToV7piTF1bWtNGwrLtHVPlWTsHxlKkq45qa/G5mXKU+VXuS+gan9Qr368Vcf1Tlem26kpI3spEDJ1w6l5yk078PlrjiKPVy8u3ywpXOC+dodH/TOTVeipbvF3HA+0WAAAgDZhsxqNjv5TM4Go+bS+Zveqmv99MPUNeyt82lvpi954NzS5X80n/ZF4PdX+Wq0BDRXxRvYbafFoKr7mxN7SOTuy3S4luWz6YU+FAk3UhTTm32sL9OW2Un2z06Npr6/S/yxc16xzG/neslOd2lvplzeSWMhQmdevQk+1Kn1Bhczw8L0ypXJvQBaLoZBpyhcIqsIb0MUje3W4m+NIUmFIslqk1MRwfHsrfdpT5tVVL3weHf3qv9/tlqfKr7mf/agxfbvWKkRv6PtuzjWwZP0u/fzJ5fpqe6m+2Faqvyxcr5ApfV9UrupAsFVqjdoaiQUAAGgTTU3OV/O95tQutLS+IdJ1adPuchV4vNqyt7LByf1qdsuJxFRc4atVPFu3iHdLcaUKPV79sLsi2rUqEt+OkqroDWTkZnJHSVWTsbe0VsIfDGnZxiJdfWI/VfuD8gdDLapDeferHbrjnbWq8Abk9Qe1fW+VXlmxRRv3FQc3dRMcOUchU+qSZJfTFv5uXXZDpgyZkiJ5ZMg0ZTEM7SrzygyZ2l3mVdAMd5t69ZOtHaZeYEB2SrieZ98Ef6Ykm9Wi0ip/+Lt5fbV2lXu1o7RKUjjR+O2bX+ov/16vr3d66n1nDV2rdZdFznFk5KxCT7VmLVyvQk+1QvuaTLxBU5Z98VT6wl2zYq01amskFgAAIG6ac1Pd0hvvmvUHDqtVbpdNnqqAfMH6k/s19BS4bqFyzSJeXzAoT3VA7gSbHFZL9CbzuldX6pudpZr26io9/sF3evXjLdpRUqUnlnyvaa81HLs/GNKrn2zRI4s31quViAzZ2tDNfeQmNTPFqT5dk3REVnKz6lAin/fYku/kC4RkMcI3+UXlXpV7A/rNm19Gk6RXVmzRqx9vabQFZtXWvcp2u5SR7JTVkFITHLIYksNmkcNmUYrTqtzUBFmthkIhU/6gKW/QlC8QkmR2qALu9CSHbhx/hBz7kiSLEV72zLIf9NCiDdEYTVPhRC4QUlGZV5W+oLz+oHaUVkWPZdnGIk17fZW+2Faiaa+v0rKNRQ1ev5HvcO22Uj36/ka98dlWbSwqU925y0OS+nRNlNNm0a9PHxD3mcoPhBoLAAAQF40VIB+Vmxp9KtucdaTwze67X+3UmUO71as/6J2eqF3lXiXarXrykuH6fne5Hnl/owo8VQ0+6U9x2TS8Txd9uqlYBZ4qpSbaNbJvui46rqemvbpKlf6gDEPy+kNKdFi0qagiPIeBKW3cVS5fMKT/WbhOLy7foi3FFXJYLbJZjXqxF3m8eu3jLdpd4VOfronaW+mTjPAIS5EhW4/rnR6tf/AHQ3po0QZ98sMefbPTo0uf/UTBkCmnzaJvdnqU6XZGz0FDrURFHq9e/2SrqrzheRjKq0P7ZpoO30xX+wPR+pMNhR59u9Oj1AS7zj22R63zW+Tx6p3VO/SLkb318idbtLO0Wq59hSGV3qBMSU6bJXqOat4sm5K8AVNpibZoS1FHeALfNyNZOakubd9bKZfDpkSHVXsrfDIMQy67JdqKEAxJW/eGWy7MfT9l1QHtrfRp464yzV4cTjACwZCKK3x6aNEGSVJZVUA2i0U/7q3UfQu+UaLdpnWFHl3xwmeq9AX12AffKRiqm1aEFZV75bRZ1T8zuR3ORGxosQAAAHHR2JwWNUcLas46Uu2uJg0NjZqe5NBFoxqe3K/uU+Aij1crN+/V7888Ukd1S9XvzzxSKzfvlSRdNKqXkpw2bdtbpSpfQOsLyhWSVFLpU9W+hCNkStX+kLaVVIZ3aIRf17yZjjzF/qGoXEVlXm3eHR5JaVtxlXyBoGa9t6FeK0ckMSiu8MkfDKnaH+5u5bBb9P8+/F4bC8oa7XJVs2tYUblPdqsRToYUvhl0WC1KS7Rr655KbS+tlDcYLsL+wztf692vdkaLy/+zfv/T91c+2aJfjOqlo7q5leyyyWW3RodfDYVCshiG6uY3hiQzZGpXmbdD1Qtku13qkuiQxWJRaN/3kOSwKtlp1Y/FVY3Oexfa14rhsllkmtLO0irtrvApaEq7K3zasbdKJZX7r98qX1A/FlepqLxagWBIXn84EfMFTVkbmYm8pMKvkf3SD3oQg/ZEYgEAAOKisTktmuqW1NA6dbuafLGtRMP7dFFqjSFtZ5yWp0tG9W50cr+G9nX/P79Rptup+//5jb7ZWapf/+8XKq3yq8oXlNcfUqU/tK+OIDxEadCUAvu6DpmS/Pv67BsK1xtEbqbtVotmL/5OO0qq5AuZCpnh+o2y6oC8gZB2lFZrV3l1veLy73aVyxsIhVsHzMinSLmpCSoordZv3vyywS5XDXUN8wVMOaxGtNuPPxjSzhKvAsGgqvaN6hQ+rqCeXPqdNhWVq7jCp8c+qN169NqnW3VCXka0XiOSrBiGIW8gpL4ZSbXOr9USjrpLoqND1QukJzl0w6lHyO2yq9wXlC8Yki9oanifLjKM2hNq1739t1stuv7kI5SZsq+gfd8QsV5fUOW+gJKdNhWVR5LBgAKhkFIS7AqFpJqNFL5Awy0WVouhnSVVKqsOtO5BtwESCwAAEBdNzWnR3HUa6ir18KIN+uT7PdEWh4ZaJRoqLK+7r70Vfi3dUKSSSp/sFov2lPv0xNLvFQyFlOgIP5132q1Kcljl3NcVKFCjJMGQJHP/zWN6olPTTz1Clb6AthRXqMwbCN/Aa/8QpwkOq0KhcPeamrNZL1y7U7MWrZc3ENSmonKZCu/XVLi1ZG+VX9X+QINzWmwvqdS2kkp1TbbLajHUOz1R2W6n+mel6JSBmaryBxU0pZJKr8x9CUFkhCSLEZ4E7zdvfql1BR6tL/TIbjOirUfl1QG5XTYVlFaryh9OriLHHh45yy9njSfxNsOQ3SrddfZRHa5eYFiPNGW5nbIa4fk5QqappeuK1DcjSXbb/mOoeftvSLr9jCP1k8HZ+vfXBeqSaI/WajhsFqUnOTRlZE8ZRnh0LF9QCgZNFZVWK6RwDcWBOO0WVQdCHWbej6aQWAAAgLg5ULekA61Tt6uU3WpoQ2G5vi0o0/3//EbH9unS7C4kdfeV4rLKHwgpxWWv99phs8hpt8gfCMkwDGUmO+W0GXLYLLIo/LTfZbdIhtQrPUF5WSl64pJjNbpfV324vkieqvC8CPtv4MPD81othpx2i7z+YHQYWZfDqrkrt6nCG1Si3aqQuT+pMPf9dlotykxxymIxonUakRtRqwx5qgIq9HhrdQ07Jz9XpdWB8HaSQjIUCJoKBMOtKDKlQCik0spw0uKwhrv7fFdYLn8wpB92V6jSF9Q/vypQoqP25Bu5aQkKhUIqq/LLHzRls4RvzE3DkD8o3f1/X3eYUaEiCj3VqvAGZLVIFoshd4JNhZ5qHdktRWaNbMJSo8kiwW7RqYOyFAiaWvHdHrnsVvn3ZZL+kKlEh02j+nVVkiM8VLHVCCcTeyr9jcbhdtU+ly67VV0SHR2m21hTSCwAAEBcNTaRXnPWqdlVyh8MadPuCpky5bSGb6ZXbtnb7C4kdbtdlVUHZbdZVFbtr/dappSe6JDDZlUwZCo9yaG8rBQNyE6Rw25RitMWnmvCYdNfLximX53SX5lup/6zvkhPf/iDQqGQgvu6QRmSbJZwUmGVlOy0yWIYCoRMpSXadf7wHqr2BZWZ7AiPumS3yJCUaLfIabPIZQsnIoUeb7ROw2mzaOWWvVqyfpfunP+1AqHwSEZVgaCCpqlfnXKExuZlqrw63GoSUniSO4skW427w2Ao3KXLabfKZrXInWBXMGTKGwipwheQKVPl1X71Tk+UzbK/m9CPe6vkD0nVQVOmpO77JpAzjPBPSaWvw4wKFbFpd7kKPV55A6ZKK/36clupqv0h/fPLguhQtNL+FiiLIbkcVu0qr9Ybn21VWbVfoZAZvbm2SDJNU/NWb1eZN5xIhHTgm+8eXRJrvfZU+XVSXkaH6TbWFBILAAAQV01NpHegdSIjOKUm2uULhmQYUt+MJNmslkYLvRtTt9tV12SHfnVyf6UnORt8nZuWoPvOOUpDclP1t8n5uvqkfrrxtCPULTVBATOcbPTPStaAbLfOye+usuqAHvvgO3mD4TqJzGSnDCPc1aV316TwXAqmqcxkh3LTEpS3bxjZM4Z0q5E8mfLvGyrWabcqPdmhrXur5AuaKiitVll1eFjdCm9Qr3+yVY/vq4lIsFmVmRIeGjbRYdHArBRlu11KdtpU6KmKdtsyjHA3MWPfOU92WmUxpE27K1TpC6iozLuvG1Z4CNmiMq8KPF7t8FTLabNGi7WDITNab2FK2lG6/zMshqGsFGeLvpu2Vlzh0wsfbZHbZQsX4CucVDlsUnUg1GiXJV/A1CffF+uv723QuoIy7SitUvcuCbJZpEHdUlRa6dd7X++UYYRbgyJJXFO2l1TVWWJq0TeFHSoJawyJBQAAiKumJtI70Do1R3A6Msetgdlu+QNmo4XeB1K329WvTj6iydc/G5ari0f3Um5ags7J767xg3P06lWjNKRbONm44oQ+slkN+YMhPfivb7Wu0COfPzziUlqiXS6bRX0zkvXiFSP1wuUj1SMtUaYMFZZVy1Md0PpdZbUSnpAZksNmjbYcFJeHbzZdtvDkdDKkJKdNgVBIPxRVaE+FV5nJDlkshqr9IVX5w8tnzFmtr3eU6hejestqCQ8NGwyFWxcq9w0R2yUp3CLTJyNRvkB428gQq75guNXCGwjJYQu3DknhpM5uDXfnshrhlphIS0YkJXTZLSqpDHSoUaEi3eB6pycqyWmLxms1LI2OCOWyW+W0GXpk8UZVeAPhAv6Qqe0lVbJbLfq+qFy7yrz6YXelSit89eaoaEzkXEYEQtKGwnItXLszhiNsHyQWAACgU6o5gtN9C75RSNK0U5suBm+OuiNG1e2GVfN1QwlPZPtIsmG3WrSxoEz//HJnuK/+vpv4LcWVykxxqUuiXV2SHLJaDflDIVX6AnJYLLIahp5c+r2KK3zRhGdIbpruP2eI+mcmy1+jhcZiCY/wZJjhfe+uCA9/W1Dq1brCMnkDARV6qmUxDDms4QLv+//5rV5Yvlmmacq+rytWpJjcbjVUWuVXcN+wttrXZauhG0d/MNzK0i8jSX+9YFi429O+OpBgyJTNatGArBT1y0zSgOwUSTro76at1OoGFwxGE6gqf6jRxMIwTFX6Q/KHwi1IMsIjOMmU/IGQyqqD0dHBilrQ2lBvkjxTMmXqzZXbOnyrBYkFAADoVBqasdpT7VdZtV9De6QesBj8QOomC3W7YR2o61bd7ZdtLNJv3vwiOs+FoX21BjL0y7F9deFxPWWzGiqu8MlhtSgr2SmrxVDXZLu27a3Ujn1dYyIJy1nDuunRi46p1UJjmvsnuvP6wzNCWwzJ7bKptMqvCm+4tcFmDbdspCXYtHVPpYorvXLYrEpy2fa1dliVmmhX/8xkGTLkDZqy7ysqd9gs9W56DYU/L9FpU1qiQwOy3bpv0hD1TA/XCRiG1Cs9Ub/+6UBNPy1Psy+O7btpK5FWoSSnVdX7hn2NFGmbRv0hZiWp0huS2xmevyOk8Dr+YLguppGRYw9aN7dLVb5gh+k61pgOkVg8/vjj6tOnj1wul0aNGqVPP/20yfXnzp2rQYMGyeVyaejQoXr33XdrvW+apu68805169ZNCQkJGj9+vDZu3NiWhwAAANpJZKK4kgpfdASnrGSnHFaL9lb6mlUM3hJ1E4XmdN2KiAxhW+ULP9UOhsJJgN1qaEB2ss46Oje6r7rF44Uer8qqAtGbtZqfm+l2auoJfXTjT/KUlmhX0DTVMz1R3bskRGe+TnLaZN83kpPNIln3fX6FN6CdnvANalaKUzaLoX5dk2QoXDBumuEb5AHZyeqe5pLNYlHQNBWqMzO0Zd9kgA6bVdeO66fzR/SQzWpobF6mXrxipIZ2T9XQ7ql64fKROmVgls7J767ctIQG5w/pCMbmZerXPxkoh82iJKdVyU6bEu3h85eRXP9ashpSSoJdiXZruHbCDNee+PzBZn+mwyKlJ9pkbSR5idheWqWQ1GG6jjUm7onFG2+8oZkzZ+quu+7SqlWrNGzYME2YMEG7du1qcP3ly5drypQpuuqqq7R69WpNmjRJkyZN0tq1a6Pr/PnPf9ajjz6qp556Sp988omSkpI0YcIEVVd37CwPAAA0reYM0gUer7bsraxXT9GcYvD2Eum7n+N2KtERfrptSrJbrZp+al6jc3ZUBYLaXe6VPxTSnfPrD80aSTJOGZgVbaF58YqRevKS4Upy2JSR7FQgGNK2kqrwcLYWIzyC1L7PT3ba1KtrgkoqAwqGTO2t8is1wa6sFGe0G9lVY/sqZIaTowS7NdzNRwrXhXRNkmvfDfh95xylicO610q2Mt1O/WJ0b6UnOWoleC1JyuKhpNonmVKlN6jy6oAq/OFS67rRhkeEsmncgAxV+PbXRCTYLSrzNn8iu0SnTb0zkuSyW6ND8jYs/tdyc8T9W33ooYd09dVX64orrtDgwYP11FNPKTExUc8991yD6z/yyCP66U9/qltuuUVHHnmk7rvvPh177LF67LHHJIVbKx5++GHdcccdOuecc3T00UfrpZde0o4dOzRv3rwG9+n1euXxeGr9AACAjqWhGaQ9VeFRkGr22e9IN681WyGshqGuSQ4lOa26flw/jRvY8Jwd9048SnaLRRnJTiXYrPUmvKurZk1I38xk9c9KlqFwEXZEMCR5AyGlJdiV5LDpbxfm646zBkfrUZIcVqUm2PXw5GM0ZF9XpbysFFkMqdIflMNm0ZG5KXLZLBqYk6K/Tx2ho7qn6jenD9TPhuXWi8lutWjSMd016ZjuHSLBa47iCp9eWbFV6UmO6MhQEaXVASXZrdHuUYkOqy4/obe+2uZRRnJ4HhCrRarwheRvzqx3+3iqAioo9SoQDEXnv2hIn64Jskh0hWqKz+fTypUrNX78+Ogyi8Wi8ePHa8WKFQ1us2LFilrrS9KECROi62/atEkFBQW11klNTdWoUaMa3ecDDzyg1NTU6E/Pnj1jPTQAANDK6k5gF5lBul9GUofrsx9Rfwhbpy4d00fXjOvfaOITlKmUBJuy3c7oDNeeKr/e+Gyr/MH6d601Eymb1dAFI3rItm/I2EgXm5AZnjMj2WlT/6xk9c1MjhaER0ai2l5SWWtSwUirUHl1QJ4qvwpKvXLYLTr3mO7qlubSL0b31iWjezd6HB0pwWuOmiNDJeyb8C+SEkUSjSMykpTstOqlK0fprKG5Kqv2Kz3JLlPh5K2lQpJ2l3tlsRg1RqLa/34kkanwBjvUKFqNies3vXv3bgWDQWVnZ9danp2drYKCgga3KSgoaHL9yO+W7PO2225TaWlp9OfHH388qOMBAABtp24NQlG5T10SHUpLdLRaPUVbqDmE7eyLjtHvfjqoyZvt7mmJ6pGWqD3l/uhxOm0WLfq6UEUeb5OfZbdadGyvLuFJAu2W6KzeVkNKsFv14HlDo0PgSqo1EpXdYolOKri1uDI6r4PFEq4fKKsO6OfH9tAlo3sr0WHrVElDc9S8vmpWqUdmR5dMnTwoS73Sk3REVnJ0/cJSX3QUqYPhD5rh0aeM8BwfDptl/0zl+3a6t8KvkwdmdphRtBpz6FwNMXA6nXK73bV+AABAx1L36X9aol03nHpEtGi4I6s7hG1T6h6nYUiFHq/W7yrTjDmr69Vb1JXtdik90an0xH1deszwZHrh4V7dtRKCuiNRRSYVXF9QFn1673bZlewMt6BMPq7nIZVM1FTzvJsKj4Zls+ybH8SUenVN0lVj++nacf2U4LAqPcmh60/pr0qfv9F9tuSqjAwrGwiaOiIzPJO5LdxwomAopHfW7GC42aZkZGTIarWqsLCw1vLCwkLl5OQ0uE1OTk6T60d+t2SfAACgc6g7gV1ktKGOfrPb0m5BkeMckJUiu9VQSZVPNsM4YL2FtP8GOTctQU5ruMD6nrOH6OqT+tVLwBpqBUpJsGtgTkp0uWlKic7wjXRH74oTq8h5z8tKUc8uCerVNTxsbp+MJN1x1pHKdrtqfY8Ds1Lka6IPlL3O+bYd4OtPdFgVMk1NHtlLFsOQdd/05YYhbd1TqY2FZTEcXduL61+hw+HQ8OHDtXjx4uiyUCikxYsXa8yYMQ1uM2bMmFrrS9KiRYui6/ft21c5OTm11vF4PPrkk08a3ScAAOg8Wns42Y6qS5JDFb6gfthdEZ2t224zVFblP2ARbzQB21dgPfGY3AYTm4ZagaafeoT6ZyZHlwdCIfXskqgZp+V1+K44rWH9rjJ5qgPaUVqtnfvmEGlsxKZ5q7er0td4YtE12alEu0WJdkNH5bqVkdx0i1WFNyhDUpU3KKvFEs4opBq/W3w47Sru6f3MmTP1zDPP6MUXX9S3336r66+/XhUVFbriiiskSZdddpluu+226Po33nijFi5cqFmzZmndunW6++679fnnn2vatGmSJMMwdNNNN+n+++/X/Pnz9dVXX+myyy5Tbm6uJk2aFI9DBAAAragjDSfbljzVfnkDwXDn/X1zUGzaXaEEh7XJlgN/MKR31mxXWpL9gAXWUv1WoEgRfGPLD2XFFT49ueT78OzZwZC8gfDM2+W+QL2Wone/2qHHl37XaG2Fsa+2JTwhoUOj+6arrDogRxPXrSkpYErzv9ihrsn2/YUbptSra6LyslJa6Ujbhi3eAUyePFlFRUW68847VVBQoPz8fC1cuDBafL1161ZZLPv/GI4//ni99tpruuOOO3T77bcrLy9P8+bN05AhQ6Lr/Pa3v1VFRYWuueYalZSU6MQTT9TChQvlch3azXcAABwOIt2KDnWR+SP6ZiRp465yhWTKkKHzh/dosuUgEDQ1f80OnT44p9nnqbFWoMOldSgiMjJUissavam3GFKqyx5tKUpPcqi4wqcnlnwvf9CsUdxdm2GG92dYLNpT6ddrn2xVdSAkp90iBfdvEZloMLqdwkP8JjntcjvtWldQpgHZybrlp4M6fItR3BMLSZo2bVq0xaGupUuX1lt2wQUX6IILLmh0f4Zh6N5779W9997bWiECAAC0q0j9Q3GFT4kOq3yBkAZkJ+unQ7o1ud3eCp/KqgPaW+FTgiOhWZ/VWCvQ4dI6FBGtOSkLj74VMsMtD2XVfmWkuKItRYWeankDQTntFgVCIcmsnVzYLYb8IVOhkGRTSMHQ/taL4L6kwmkLDzFrtVhU4ds/W7fFEp4RvcIb1G8mDNQf//mtZl98rHp3TWqns3Dw4t4VCgAAAPXVrH8IhsLDx9adrbuuyMzk3+wsbdYIUhGNFZd3trkoYhU55zarsW+UpvCQrzarJToBoxROQLokOpXitNUqf3DUGMbXZgm3RkRaNJx2iyyGZLeF5xixWy26dEwf3XPOUXLa9iduoZC0o7RaKQl2HdU9VTPG5yk3rXkJYrwdHlcJAABAJxSpcxic69bN4wc0OFt3RM2Zye0WS7NGkEJ9R+WmKslhk9NuUbLTKqfNogS7VUflpkbXqTvyVoLDoox9Ezda99XDREZ0iszEXeELt2zsn329v37700Eam5epbHeC7DXuyj1VAV1+fJ96o1B1dB2iKxQAAAAa1iXJIbfLrsnH9WryBrPmzOQllb7w7xp1AWieQk+1qvxBuWxWWYzwELHV/lC98zg2L1P9M5I1Y85q+YMhVfmDcifYVVBaJYvFompfQDL3t1hEfqcnOfSzYbnR2deLK3wyTVNWiyU6b0mK06Z+mR2/61NdJBYAAAAdWHPrHBqakyI9+dCfe6K1ZbtdSkmwqaTKIq8/KKfdKneCrcHzmOl26uJRvVTmDej/vtihkgq/HFarEhxWBYKmAqGgEh1WVfuDOjInRd/vrtBvTh+osQP2tzw5bBbtrfQruG+abUOGyryBTtNKUVPnixgAAOAw0tw6h8bmpKC1omXSkxy68bQ89eySqEDIbHIOD7vVop8f20NTx/TRjaflRc99blqCbj9zkJKdNmWkOOW0W1UVCC8fkFN7yFhfIKQuSXY5982el+CwqEuSXf4mJt7rqGixAAAAOEREuufc/MYa/W1yvnK7dI6i347mYM5j3W26JDn01qrtshjSnnKvuiQ6dOP4PGW7a7d8ZLtd6uZOkNUwVOWrVEayU1luV6dsaaLFAgAA4BCS6Xbq4tG9lOluepZnNO1gzmPNbWxWQ1ec0Ed/m5zf5CSDkZamrslO2awWdU12dtqWJlosAAAADiGHywSCbe1gzmPdbc7J764qX/CAkwweKi1NtFgAAAAAbaS5xfeHQksTLRYAAABAG2luy8eh0NJEiwUAAACAmJFYAAAAAIgZiQUAAACAmJFYAAAAAIgZiQUAAACAmJFYAAAAAIgZiQUAAACAmJFYAAAAAIgZiQUAAACAmJFYAAAAAIgZiQUAAACAmJFYAAAAAIgZiQUAAACAmJFYAAAAAIgZiQUAAACAmNniHUBHZJqmJMnj8cQ5EgAAACB+IvfDkfvjppBYNKCsrEyS1LNnzzhHAgAAAMRfWVmZUlNTm1zHMJuTfhxmQqGQduzYoZSUFBmGIUk67rjj9Nlnnx1w29ZYz+PxqGfPnvrxxx/ldrtbFnwn0dzz1FljaK19x7Kflm7bkvWbs+6B1uE67/yfH+/r/GC24zpvXfG+xts6hnhf4wezbWtf4wda71C/xiWuc9M0VVZWptzcXFksTVdR0GLRAIvFoh49etRaZrVam/UH05rrud3uQ/aPtLnnqbPG0Fr7jmU/Ld22Jes3Z93m7o/rvPN+fryv84PZjuu8dcX7Gm/rGOJ9jR/Mtq19jTd3vUP1Gpe4ziUdsKUiguLtZrrhhhvist6hqiMcf1vG0Fr7jmU/Ld22Jes3Z92O8B3HW7zPQVt/fryv84PZjuu8dXWE4+ff8oNfn3uW5ukIx98ZrnOJrlAdksfjUWpqqkpLS+OeIQNtheschwOucxzquMZREy0WHZDT6dRdd90lp9MZ71CANsN1jsMB1zkOdVzjqIkWCwAAAAAxo8UCAAAAQMxILAAAAADEjMQCAAAAQMxILAAAAADEjMQCAAAAQMxILDqhBQsWaODAgcrLy9Pf//73eIcDtLpzzz1XXbp00fnnnx/vUIA28eOPP+rkk0/W4MGDdfTRR2vu3LnxDglodSUlJRoxYoTy8/M1ZMgQPfPMM/EOCW2M4WY7mUAgoMGDB2vJkiVKTU3V8OHDtXz5cnXt2jXeoQGtZunSpSorK9OLL76oN998M97hAK1u586dKiwsVH5+vgoKCjR8+HBt2LBBSUlJ8Q4NaDXBYFBer1eJiYmqqKjQkCFD9Pnnn3PPcgijxaKT+fTTT3XUUUepe/fuSk5O1hlnnKH33nsv3mEBrerkk09WSkpKvMMA2ky3bt2Un58vScrJyVFGRoaKi4vjGxTQyqxWqxITEyVJXq9XpmmK59mHNhKLdvbhhx/q7LPPVm5urgzD0Lx58+qt8/jjj6tPnz5yuVwaNWqUPv300+h7O3bsUPfu3aOvu3fvru3bt7dH6ECzxHqNA51Ba17nK1euVDAYVM+ePds4aqBlWuM6Lykp0bBhw9SjRw/dcsstysjIaKfoEQ8kFu2soqJCw4YN0+OPP97g+2+88YZmzpypu+66S6tWrdKwYcM0YcIE7dq1q50jBQ4O1zgOB611nRcXF+uyyy7T008/3R5hAy3SGtd5WlqavvjiC23atEmvvfaaCgsL2yt8xIOJuJFkvv3227WWjRw50rzhhhuir4PBoJmbm2s+8MADpmma5kcffWROmjQp+v6NN95ovvrqq+0SL9BSB3ONRyxZssQ877zz2iNMICYHe51XV1ebY8eONV966aX2ChU4aLH8ex5x/fXXm3Pnzm3LMBFntFh0ID6fTytXrtT48eOjyywWi8aPH68VK1ZIkkaOHKm1a9dq+/btKi8v17/+9S9NmDAhXiEDLdKcaxzo7JpznZumqcsvv1ynnnqqLr300niFChy05lznhYWFKisrkySVlpbqww8/1MCBA+MSL9qHLd4BYL/du3crGAwqOzu71vLs7GytW7dOkmSz2TRr1iydcsopCoVC+u1vf8voCug0mnONS9L48eP1xRdfqKKiQj169NDcuXM1ZsyY9g4XOCjNuc4/+ugjvfHGGzr66KOj/dZffvllDR06tL3DBQ5Kc67zLVu26JprrokWbU+fPp1r/BBHYtEJTZw4URMnTox3GECbef/99+MdAtCmTjzxRIVCoXiHAbSpkSNHas2aNfEOA+2IrlAdSEZGhqxWa73CpsLCQuXk5MQpKqD1cI3jcMB1jsMB1zkaQmLRgTgcDg0fPlyLFy+OLguFQlq8eDHdQHBI4BrH4YDrHIcDrnM0hK5Q7ay8vFzfffdd9PWmTZu0Zs0apaenq1evXpo5c6amTp2qESNGaOTIkXr44YdVUVGhK664Io5RA83HNY7DAdc5Dgdc52ixOI9KddhZsmSJKanez9SpU6PrzJ492+zVq5fpcDjMkSNHmh9//HH8AgZaiGschwOucxwOuM7RUoZpMrc6AAAAgNhQYwEAAAAgZiQWAAAAAGJGYgEAAAAgZiQWAAAAAGJGYgEAAAAgZiQWAAAAAGJGYgEAAAAgZiQWAAAAAGJGYgEAAAAgZiQWAIB2c/nll8swDBmGoXnz5kmSNm/eLMMwtGbNmjb73BdeeCH6uTfddFObfQ4AHM5ILAAA7eqnP/2pdu7cqTPOOKPdPnPy5MnauXOnxowZ026fCQCHG1u8AwAAHHp8Pp8cDkeD7zmdTuXk5LRrPAkJCUpISGg0JgBA7GixAADE7OSTT9a0adN00003KSMjQxMmTDjofQWDQV155ZUaNGiQtm7dKkkyDENPPvmkzjjjDCUkJKhfv3568803a223bds2TZkyRenp6UpKStKIESP0ySefxHRcAIDmI7EAALSKF198UQ6HQx999JGeeuqpg9qH1+vVBRdcoDVr1mjZsmXq1atX9L0//OEPOu+88/TFF1/okksu0UUXXaRvv/1WklReXq5x48Zp+/btmj9/vr744gv99re/VSgUapVjAwAcGF2hAACtIi8vT3/+858Pevvy8nKdddZZ8nq9WrJkiVJTU2u9f8EFF+iXv/ylJOm+++7TokWLNHv2bD3xxBN67bXXVFRUpM8++0zp6emSpCOOOOLgDwYA0GK0WAAAWsXw4cNj2n7KlCmqqKjQe++9Vy+pkFSv8HrMmDHRFos1a9bomGOOiSYVAID2R2IBAGgVSUlJMW1/5pln6ssvv9SKFStavG1CQkJMnw0AiB2JBQCgQ7j++uv14IMPauLEifrPf/5T7/2PP/643usjjzxSknT00UdrzZo1Ki4ubpdYAQD1kVgAADqM6dOn6/7779fPfvYz/fe//6313ty5c/Xcc89pw4YNuuuuu/Tpp59q2rRpksLdqHJycjRp0iR99NFH+uGHH/TWW28dVOsHAODgkFgAADqUm266Sffcc4/OPPNMLV++PLr8nnvu0Zw5c3T00UfrpZde0uuvv67BgwdLkhwOh9577z1lZWXpzDPP1NChQ/Xggw/KarXG6zAA4LDDqFAAgJgtXbr0oLft06ePTNOstWzmzJmaOXNmrWW5ubl67733Gt1P7969681tAQBoP7RYAADa1YIFC5ScnKwFCxa022e++uqrSk5O1rJly9rtMwHgcGOYdR8TAQDQRnbt2iWPxyNJ6tatW7NHkjIMQ2+//bYmTZp0UJ9bVlamwsJCSVJaWpoyMjIOaj8AgMaRWAAAAACIGV2hAAAAAMSMxAIAAABAzEgsAAAAAMSMxAIAAABAzEgsAAAAAMSMxAIAAABAzEgsAAAAAMSMxAIAAABAzP4/M+nk6KGUoxoAAAAASUVORK5CYII=\n"
},
"metadata": {}
}
],
"source": [
"# Quick visual look at the data\n",
"fig, ax = plt.subplots(figsize=(8, 5))\n",
"ax.errorbar(r_data, y_data, yerr=sigma_n, fmt='o', ms=3, lw=0.8, capsize=0,\n",
" color='tab:blue', alpha=0.7, label='data')\n",
"ax.set_xscale('log')\n",
"ax.set_xlabel('r [kpc]')\n",
"ax.set_ylabel('Amplitude')\n",
"ax.set_title('Measured cluster density profile')\n",
"ax.legend()\n",
"plt.tight_layout()\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"id": "pxwhIefuP-1r",
"metadata": {
"id": "pxwhIefuP-1r"
},
"source": [
"---\n",
"\n",
"## 1. Provided infrastructure\n",
"\n",
"All the functions in this section are already written for you. **Read them carefully** - you will use them in your own code below."
]
},
{
"cell_type": "markdown",
"id": "jyGEW4YFP-1r",
"metadata": {
"id": "jyGEW4YFP-1r"
},
"source": [
"### 1.1. Noise covariance matrix\n",
"\n",
"The noise is **coloured**, so the log-likelihood is\n",
"$$\\log\\mathcal{L}(\\vec{\\theta}) = -\\tfrac{1}{2}\\,\\vec{r}^{\\top}\\mathbf{C}^{-1}\\vec{r}\\qquad\\text{with}\\qquad \\vec{r} = \\vec{y} - \\vec{m}(\\vec{\\theta}).$$\n",
"\n",
"By the Wiener-Khinchin theorem, the autocovariance of a stationary process is the inverse Fourier transform of its PSD. We then build $\\mathbf{C}$ exploiting its Toeplitz structure $C_{ij} = c(|i-j|)$."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "3bRXV9tFP-1s",
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "3bRXV9tFP-1s",
"outputId": "77904d7d-3f2c-4fb7-d1a4-3df3d816a133"
},
"outputs": [
{
"output_type": "stream",
"name": "stdout",
"text": [
"Covariance matrix shape : (300, 300)\n",
"C[0, 0] = 1.000e-06 (expected sigma^2 = 1.000e-06)\n"
]
}
],
"source": [
"def covariance_matrix(psd, sigma):\n",
" \"\"\"Noise covariance matrix from a PSD (Wiener-Khinchin + Toeplitz).\n",
"\n",
" Parameters\n",
" ----------\n",
" psd : ndarray, shape (N,)\n",
" Power spectral density of the noise.\n",
" sigma : float\n",
" Standard deviation of the noise (sets the overall amplitude).\n",
"\n",
" Returns\n",
" -------\n",
" C : ndarray, shape (N, N)\n",
" \"\"\"\n",
" psd_norm = psd / psd.mean() # normalise so mean(PSD) = 1\n",
" autocov = np.real(np.fft.ifft(psd_norm)) * sigma**2 # autocovariance via Wiener-Khinchin\n",
" return toeplitz(autocov)\n",
"\n",
"C = covariance_matrix(psd, sigma_n)\n",
"Cinv = np.linalg.inv(C)\n",
"\n",
"print(f'Covariance matrix shape : {C.shape}')\n",
"print(f'C[0, 0] = {C[0, 0]:.3e} (expected sigma^2 = {sigma_n**2:.3e})')"
]
},
{
"cell_type": "markdown",
"id": "TyVPScoAP-1s",
"metadata": {
"id": "TyVPScoAP-1s"
},
"source": [
"### 1.2. Proposal function $g(\\vec{\\theta}\\,|\\vec{\\theta}_i)$\n",
"\n",
"Symmetric Gaussian proposal centred on the current point, with per-parameter step sizes $\\delta\\vec{\\theta}$."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "rdeP0ILHP-1s",
"metadata": {
"id": "rdeP0ILHP-1s"
},
"outputs": [],
"source": [
"def proposal(theta, step, rng):\n",
" \"\"\"Symmetric Gaussian proposal centred on `theta`.\n",
"\n",
" Parameters\n",
" ----------\n",
" theta : ndarray, shape (n_params,)\n",
" step : ndarray, shape (n_params,)\n",
" rng : numpy.random.Generator\n",
"\n",
" Returns\n",
" -------\n",
" theta_prop : ndarray, shape (n_params,)\n",
" \"\"\"\n",
" return theta + rng.normal(0.0, step)"
]
},
{
"cell_type": "markdown",
"id": "0N5p566AP-1s",
"metadata": {
"id": "0N5p566AP-1s"
},
"source": [
"### 1.3. Priors, likelihood, posterior\n",
"\n",
"Uniform priors on physically reasonable boxes; Gaussian likelihood using the full inverse covariance (equivalent to $-\\tfrac12\\chi^2$)."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "ViBGZXTIP-1s",
"metadata": {
"id": "ViBGZXTIP-1s"
},
"outputs": [],
"source": [
"# Uniform-prior bounds (lower, upper) for each parameter\n",
"# theta = [rho0, rp, A, mu, sigma_g]\n",
"PRIOR_BOUNDS = np.array([\n",
" [1e-4, 1.0 ], # rho0 amplitude of the NFW profile\n",
" [10.0, 2000.0], # rp characteristic radius (kpc)\n",
" [1e-4, 0.2 ], # A amplitude of the Gaussian bump\n",
" [500., 2500. ], # mu position of the Gaussian bump (kpc)\n",
" [10., 500. ], # sigma_g width of the Gaussian bump (kpc)\n",
"])\n",
"\n",
"def log_prior(theta):\n",
" \"\"\"Uniform log-prior: 0 inside the box, -inf outside.\"\"\"\n",
" for v, (lo, hi) in zip(theta, PRIOR_BOUNDS):\n",
" if not (lo < v < hi):\n",
" return -np.inf\n",
" return 0.0\n",
"\n",
"def log_likelihood(theta, model_func, r, y, Cinv):\n",
" \"\"\"Gaussian log-likelihood with full inverse covariance (i.e. -chi^2 / 2).\"\"\"\n",
" residual = y - model_func(r, *theta)\n",
" return -0.5 * residual @ Cinv @ residual\n",
"\n",
"def log_probability(theta, model_func, r, y, Cinv):\n",
" \"\"\"Log of the (unnormalised) posterior: log prior + log likelihood.\"\"\"\n",
" lp = log_prior(theta)\n",
" if not np.isfinite(lp):\n",
" return -np.inf\n",
" return lp + log_likelihood(theta, model_func, r, y, Cinv)"
]
},
{
"cell_type": "markdown",
"id": "_SQG4CD2P-1t",
"metadata": {
"id": "_SQG4CD2P-1t"
},
"source": [
"### 1.4. Gelman-Rubin convergence test\n",
"\n",
"Given $M$ chains of length $n$, compares the within-chain variance $W$ to the between-chain variance $B$. Convergence is declared when $\\hat R$ is close to 1 (we will use $\\hat R < 1.03$ for every parameter)."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "uWhToxJDP-1t",
"metadata": {
"id": "uWhToxJDP-1t"
},
"outputs": [],
"source": [
"def gelman_rubin(chains):\n",
" \"\"\"Gelman-Rubin statistic R for each parameter.\n",
"\n",
" Parameters\n",
" ----------\n",
" chains : ndarray, shape (n_chains, n_steps, n_params)\n",
" Chains AFTER burn-in removal.\n",
"\n",
" Returns\n",
" -------\n",
" R : ndarray, shape (n_params,)\n",
" \"\"\"\n",
" M, n, p = chains.shape\n",
" means = chains.mean(axis=1) # per-chain means, shape (M, p)\n",
" grand = means.mean(axis=0) # global mean, shape (p,)\n",
" B = n / (M - 1) * ((means - grand) ** 2).sum(axis=0)\n",
" W = chains.var(axis=1, ddof=1).mean(axis=0)\n",
" var_hat = (1.0 - 1.0 / n) * W + B / n\n",
" return np.sqrt(var_hat / W)"
]
},
{
"cell_type": "markdown",
"id": "e-NxeFxXP-1t",
"metadata": {
"id": "e-NxeFxXP-1t"
},
"source": [
"### 1.5. Autocorrelation time\n",
"\n",
"The integrated autocorrelation time $\\tau$ tells you how many MCMC steps separate two *effectively independent* samples."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "gnXUeyzYP-1t",
"metadata": {
"id": "gnXUeyzYP-1t"
},
"outputs": [],
"source": [
"def autocorr_function_1d(x):\n",
" \"\"\"Normalised autocorrelation function of a 1D chain (FFT-based).\"\"\"\n",
" x = np.asarray(x, dtype=float) - np.mean(x)\n",
" n = 1\n",
" while n < len(x):\n",
" n *= 2\n",
" f = np.fft.fft(x, n=2 * n)\n",
" acf = np.real(np.fft.ifft(f * np.conjugate(f))[: len(x)])\n",
" acf /= acf[0]\n",
" return acf\n",
"\n",
"def autocorr_time(chain_1d, threshold=0.05):\n",
" \"\"\"Integrated autocorrelation time (summed until ACF drops below `threshold`).\"\"\"\n",
" acf = autocorr_function_1d(chain_1d)\n",
" tau = 1.0\n",
" for t in range(1, len(acf)):\n",
" if acf[t] < threshold:\n",
" break\n",
" tau += 2.0 * acf[t]\n",
" return tau"
]
},
{
"cell_type": "markdown",
"id": "VlDCiYDfP-1t",
"metadata": {
"id": "VlDCiYDfP-1t"
},
"source": [
"---\n",
"\n",
"## Block 1 - Define the model\n",
"\n",
"The model is a **generalised NFW profile** plus a **Gaussian bump** that accounts for the merging substructure:\n",
"\n",
"$$\\rho(r) \\;=\\; \\frac{\\rho_0}{\\left(\\dfrac{r}{r_p}\\right)^{c}\\left[1+\\left(\\dfrac{r}{r_p}\\right)^{a}\\right]^{(b-c)/a}} \\;+\\; A \\, \\exp\\!\\left[-\\frac{(r-\\mu)^2}{2\\sigma_g^{\\,2}}\\right]$$\n",
"\n",
"with the **shape parameters fixed** to $a = 1.1$, $b = 5.5$, $c = 0.31$ and the **free parameters** $\\vec{\\theta} = (\\rho_0, r_p, A, \\mu, \\sigma_g)$.\n",
"\n",
"**Your task**: complete the function `model(r, rho0, rp, A, mu, sigma_g)` below."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "4UXCK8zzP-1t",
"metadata": {
"id": "4UXCK8zzP-1t"
},
"outputs": [],
"source": [
"# Fixed shape parameters of the generalised NFW profile\n",
"a_NFW, b_NFW, c_NFW = 1.1, 5.5, 0.31\n",
"\n",
"def model(r, rho0, rp, A, mu, sigma_g):\n",
" \"\"\"Generalised NFW profile + Gaussian bump.\n",
"\n",
" Parameters\n",
" ----------\n",
" r : ndarray\n",
" Radii (kpc).\n",
" rho0, rp, A, mu, sigma_g : float\n",
" Free parameters of the model.\n",
"\n",
" Returns\n",
" -------\n",
" rho : ndarray, same shape as r\n",
" \"\"\"\n",
" # =================== TODO ===================\n",
" # 1) Compute the NFW component using a_NFW, b_NFW, c_NFW.\n",
" # 2) Compute the Gaussian component.\n",
" # 3) Return their sum.\n",
" raise NotImplementedError\n",
" # ============================================"
]
},
{
"cell_type": "markdown",
"id": "JlgmAimPP-1t",
"metadata": {
"id": "JlgmAimPP-1t"
},
"source": [
"**Sanity check.** Run the cell below to plot the data together with the model evaluated at *guess* parameters. Tweak the values until the model roughly matches the data - this gives a good starting point for the MCMC."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "MfuLq1efP-1t",
"metadata": {
"id": "MfuLq1efP-1t"
},
"outputs": [],
"source": [
"theta_guess = np.array([0.01, 600.0, 0.02, 1500.0, 200.0])\n",
" # rho0, rp, A, mu, sigma_g\n",
"\n",
"r_smooth = np.logspace(np.log10(r_data.min()), np.log10(r_data.max()), 500)\n",
"\n",
"fig, ax = plt.subplots(figsize=(8, 5))\n",
"ax.errorbar(r_data, y_data, yerr=sigma_n, fmt='o', ms=3, lw=0.8,\n",
" color='tab:blue', alpha=0.6, label='data')\n",
"ax.plot(r_smooth, model(r_smooth, *theta_guess), '-',\n",
" color='crimson', lw=2, label='guess model')\n",
"ax.set_xscale('log')\n",
"ax.set_xlabel('r [kpc]'); ax.set_ylabel('Amplitude')\n",
"ax.legend(); plt.tight_layout(); plt.show()"
]
},
{
"cell_type": "markdown",
"id": "7Cnem6E7P-1t",
"metadata": {
"id": "7Cnem6E7P-1t"
},
"source": [
"---\n",
"\n",
"## Block 2 - Implement Metropolis-Hastings\n",
"\n",
"Using **only** the helper functions already provided (`proposal`, `log_probability`), write a function that runs Metropolis-Hastings for `n_steps` iterations starting from `theta0`.\n",
"\n",
"Recall the algorithm:\n",
"\n",
"1. Propose $\\vec{\\theta}\\, \\sim g(\\vec{\\theta}\\,|\\vec{\\theta}_i)$.\n",
"2. Compute the log acceptance ratio \n",
"$\\quad \\log\\alpha = \\log p(\\vec{\\theta}\\,\\,|\\,D) - \\log p(\\vec{\\theta}_i\\,|\\,D)$ \n",
"(symmetric proposal: the proposal terms cancel).\n",
"3. Draw $u \\sim \\mathcal{U}(0,1)$. If $\\log u < \\log\\alpha$ --> accept, else stay at $\\vec{\\theta}_i$.\n",
"4. Store the current state in the chain; iterate.\n",
"\n",
"> Trick: Compare $\\log u$ to $\\log\\alpha$ (not $u$ to $\\alpha$) - this avoids overflow when $\\log\\alpha$ is very negative."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "4taZl2jTP-1t",
"metadata": {
"id": "4taZl2jTP-1t"
},
"outputs": [],
"source": [
"def metropolis_hastings(theta0, step, n_steps, model_func, r, y, Cinv, rng):\n",
" \"\"\"Run Metropolis-Hastings for n_steps iterations.\n",
"\n",
" Parameters\n",
" ----------\n",
" theta0 : ndarray, shape (n_params,)\n",
" Starting point of the chain.\n",
" step : ndarray, shape (n_params,)\n",
" Proposal standard deviations (one per parameter).\n",
" n_steps : int\n",
" model_func : callable\n",
" Model function with signature ``model_func(r, *theta)``.\n",
" r, y : ndarray\n",
" Data arrays.\n",
" Cinv : ndarray\n",
" Inverse noise covariance matrix.\n",
" rng : numpy.random.Generator\n",
"\n",
" Returns\n",
" -------\n",
" chain : ndarray, shape (n_steps, n_params)\n",
" acceptance_rate : float\n",
" \"\"\"\n",
" # =================== TODO ===================\n",
" # 1) Allocate `chain` as a 2D numpy array of shape (n_steps, len(theta0)).\n",
" # 2) Initialise theta = theta0 (as a float array) and compute its log_probability.\n",
" # 3) Initialise an `n_accept` counter.\n",
" # 4) Loop n_steps times:\n",
" # a. Propose a new theta via the `proposal` function.\n",
" # b. Compute its log_probability.\n",
" # c. Compute log_alpha = (new log_probability) - (current log_probability).\n",
" # d. Accept if np.log(rng.random()) < log_alpha — then update theta, lp, n_accept.\n",
" # e. Store the *current* theta into chain[i].\n",
" # 5) Return (chain, n_accept / n_steps).\n",
" raise NotImplementedError\n",
" # ============================================"
]
},
{
"cell_type": "markdown",
"id": "CfQTV2Q3P-1u",
"metadata": {
"id": "CfQTV2Q3P-1u"
},
"source": [
"**Run several chains** with random starting points drawn near the guess.\n",
"\n",
"We use `step = 1 % × init`, as a first guess for the step. Aim for an acceptance rate between roughly **20 % and 50 %** - if it is much lower, your `theta_guess` is probably far from the true mode."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "X0HE7ymqP-1u",
"metadata": {
"id": "X0HE7ymqP-1u"
},
"outputs": [],
"source": [
"n_chains = 4\n",
"n_steps = 30000 # ~ 6 s on Colab for the four chains\n",
"n_params = 5\n",
"\n",
"chains_raw = np.zeros((n_chains, n_steps, n_params))\n",
"acceptances = np.zeros(n_chains)\n",
"\n",
"for c in range(n_chains):\n",
" # Each chain starts by jittering `theta_guess` by ±15 %\n",
" init = theta_guess * rng.uniform(0.85, 1.15, size=n_params)\n",
" step = init * 0.01 # 1 % of starting values\n",
" chain, acc = metropolis_hastings(\n",
" init, step, n_steps, model, r_data, y_data, Cinv, rng\n",
" )\n",
" chains_raw[c] = chain\n",
" acceptances[c] = acc\n",
" print(f'Chain {c+1}: acceptance = {acc:.2%}')"
]
},
{
"cell_type": "markdown",
"id": "A2qPWkmJP-1u",
"metadata": {
"id": "A2qPWkmJP-1u"
},
"source": [
"**Trace plot** (provided). Inspect how each chain wanders in parameter space - you should see an initial transient (the *burn-in*) followed by stationary fluctuations around the mode."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "XeAIVTvvP-1u",
"metadata": {
"id": "XeAIVTvvP-1u"
},
"outputs": [],
"source": [
"param_labels = [r'$\\rho_0$', r'$r_p$', r'$A$', r'$\\mu$', r'$\\sigma_g$']\n",
"\n",
"fig, axes = plt.subplots(n_params, 1, figsize=(9, 1.7 * n_params), sharex=True)\n",
"for j, ax in enumerate(axes):\n",
" for c in range(n_chains):\n",
" ax.plot(chains_raw[c, :, j], lw=0.5, alpha=0.8)\n",
" ax.set_ylabel(param_labels[j])\n",
"axes[-1].set_xlabel('MCMC step')\n",
"fig.suptitle('Trace plots of the Markov chains', y=1.0)\n",
"plt.tight_layout(); plt.show()\n",
"\n",
"# 2D walk in (rho0, rp) space\n",
"fig, ax = plt.subplots(figsize=(6.5, 5))\n",
"for c in range(n_chains):\n",
" ax.plot(chains_raw[c, :, 0], chains_raw[c, :, 1],\n",
" lw=0.4, alpha=0.7, label=f'chain {c+1}')\n",
"ax.set_xlabel(r'$\\rho_0$'); ax.set_ylabel(r'$r_p$ [kpc]')\n",
"ax.set_title(r'Markov chains in the $(\\rho_0, r_p)$ plane')\n",
"ax.legend(fontsize=8); plt.tight_layout(); plt.show()"
]
},
{
"cell_type": "markdown",
"id": "XucleS4yP-1u",
"metadata": {
"id": "XucleS4yP-1u"
},
"source": [
"---\n",
"\n",
"## Block 3 - Burn-in, convergence, and thinning\n",
"\n",
"Now use the **Gelman-Rubin** test and the **autocorrelation time** (both provided in Section 1) to turn the raw chains into a set of approximately independent samples.\n",
"\n",
"Your function should:\n",
"\n",
"1. Drop the first $\\lfloor \\text{burnin}_\\text{frac} \\times n_\\text{steps} \\rfloor$ samples of each chain.\n",
"2. Compute $\\hat R$ on the remaining chains; record whether $\\hat R < \\text{R}_\\text{threshold}$ for **every** parameter.\n",
"3. Compute the per-parameter integrated autocorrelation time $\\tau_j$ on the first chain, and define $\\tau_{\\max} = \\max_j \\tau_j$.\n",
"4. **Thin** each chain by keeping one sample every $\\lceil \\tau_{\\max} \\rceil$ steps.\n",
"5. Stack the thinned chains into a single flat array of final samples."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "fCLXZw7lP-1u",
"metadata": {
"id": "fCLXZw7lP-1u"
},
"outputs": [],
"source": [
"def postprocess(chains_raw, burnin_frac=0.2, R_threshold=1.03, verbose=True):\n",
" \"\"\"Return final independent samples + diagnostics.\n",
"\n",
" Parameters\n",
" ----------\n",
" chains_raw : ndarray, shape (n_chains, n_steps, n_params)\n",
" burnin_frac : float\n",
" Fraction of each chain to discard as burn-in.\n",
" R_threshold : float\n",
" Gelman-Rubin convergence criterion.\n",
" verbose : bool\n",
" Print diagnostics if True.\n",
"\n",
" Returns\n",
" -------\n",
" samples_flat : ndarray, shape (n_samples, n_params)\n",
" Concatenated, burn-in-removed, thinned samples from all chains.\n",
" info : dict\n",
" Contains 'R' (per parameter), 'tau' (per parameter),\n",
" 'converged' (bool), 'tau_max' (int).\n",
" \"\"\"\n",
" # =================== TODO ===================\n",
" # 1) Compute n_burn = int(burnin_frac * n_steps) and drop the first n_burn\n",
" # samples of each chain to obtain `chains_clean` of shape (M, n', p).\n",
" # 2) Call `gelman_rubin(chains_clean)` and store the result as R.\n",
" # Set `converged = bool(np.all(R < R_threshold))`.\n",
" # 3) For each parameter j, compute the autocorrelation time on the FIRST chain\n",
" # via `autocorr_time(chains_clean[0, :, j])`. Take tau_max = max over j\n",
" # (cast to int via np.ceil); guard against degenerate values with `max(tau_max, 1)`.\n",
" # 4) Thin each chain by `chains_clean[:, ::tau_max, :]`.\n",
" # 5) Reshape to (-1, n_params) to get `samples_flat`.\n",
" # 6) If `verbose`, print R, tau, tau_max and the number of final samples.\n",
" # 7) Return samples_flat and the info dictionary.\n",
" raise NotImplementedError\n",
" # ============================================"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "Ay74EI6FP-1v",
"metadata": {
"id": "Ay74EI6FP-1v"
},
"outputs": [],
"source": [
"samples_flat, info = postprocess(chains_raw, burnin_frac=0.2,\n",
" R_threshold=1.03, verbose=True)"
]
},
{
"cell_type": "markdown",
"id": "fCLl6YEOP-1v",
"metadata": {
"id": "fCLl6YEOP-1v"
},
"source": [
"> **If the chains have not converged** ($\\hat R > 1.03$ for some parameter, typically $\\rho_0$ and $r_p$), it is because these two parameters are **strongly anti-correlated** - look at the 2-D trace in the $(\\rho_0, r_p)$ plane just above! A diagonal Gaussian proposal explores correlated directions slowly.\n",
">\n",
"> Two ways forward:\n",
">\n",
"> - re-run the previous cell with `n_steps = 60000` - usually enough to get $\\hat R < 1.03$;\n",
"> - or accept the partial convergence and proceed, keeping in mind the limitation. This is exactly the kind of failure mode that motivates the more elaborate algorithms in `emcee` (covered later in the course)."
]
},
{
"cell_type": "markdown",
"id": "hv3GPYK-P-1v",
"metadata": {
"id": "hv3GPYK-P-1v"
},
"source": [
"---\n",
"\n",
"## Block 4 - Corner plot and credible intervals\n",
"\n",
"Write a function that, given the flat samples, produces:\n",
"\n",
"1. A **corner plot** with `getdist.MCSamples` and `getdist.plots`.\n",
"2. The **68 % credible intervals** (16th, 50th, 84th percentiles) for each parameter, printed in a readable form (`median + (upper-median) / - (median-lower)`).\n",
"\n",
"Then call it on `samples_flat`."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "YvTE2YdKP-1v",
"metadata": {
"id": "YvTE2YdKP-1v"
},
"outputs": [],
"source": [
"PARAM_NAMES = ['rho0', 'rp', 'A', 'mu', 'sigma_g']\n",
"PARAM_LATEX = [r'\\rho_0', r'r_p', r'A', r'\\mu', r'\\sigma_g']\n",
"\n",
"def report_results(samples_flat,\n",
" param_names=PARAM_NAMES,\n",
" param_latex=PARAM_LATEX):\n",
" \"\"\"Corner plot + 68 % credible intervals.\n",
"\n",
" Parameters\n",
" ----------\n",
" samples_flat : ndarray, shape (n_samples, n_params)\n",
" param_names, param_latex : lists of str\n",
"\n",
" Returns\n",
" -------\n",
" intervals : dict\n",
" Maps each parameter name to a tuple (median, lower16, upper84).\n",
" \"\"\"\n",
" # =================== TODO ===================\n",
" # 1) Build an MCSamples object:\n",
" # mc = MCSamples(samples=samples_flat,\n",
" # names=param_names,\n",
" # labels=param_latex)\n",
" # 2) Draw the corner plot:\n",
" # g = plots.get_subplot_plotter()\n",
" # g.triangle_plot([mc], filled=True)\n",
" # 3) For each parameter j, compute the 16th, 50th, 84th percentiles via\n",
" # np.percentile, store them in the `intervals` dictionary.\n",
" # 4) Print the result as e.g.:\n",
" # rho0 = 0.01 (+0.001 / -0.001)\n",
" # 5) Return `intervals`.\n",
" raise NotImplementedError\n",
" # ============================================"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "hE2WWGX5P-1v",
"metadata": {
"id": "hE2WWGX5P-1v"
},
"outputs": [],
"source": [
"intervals = report_results(samples_flat)"
]
},
{
"cell_type": "markdown",
"id": "CG9tiFqnP-1v",
"metadata": {
"id": "CG9tiFqnP-1v"
},
"source": [
"---\n",
"\n",
"## Going further (if time allows)\n",
"\n",
"- Re-run the MCMC with a **larger step** (e.g. `5 % × init`) and a **smaller step** (e.g. `0.1 % × init`). Look at the trace plots and the acceptance rate. Which step gives the best mixing?\n",
"- Trade off **number of chains** vs **number of steps**, keeping `n_chains × n_steps` constant. Does the Gelman-Rubin diagnostic verdict change?\n",
"- Plot the **best-fit model** (using the medians of the marginals) on top of the data. Are the residuals consistent with the assumed noise level?\n",
"- Look at the 2-D marginal in the $(\\rho_0, r_p)$ plane in your corner plot. The strong correlation you see there is what made the Gelman-Rubin test slow to pass. Algorithms like the *affine-invariant ensemble sampler* implemented in `emcee` are explicitly designed to handle this kind of degeneracy efficiently."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "z4Voj7YeTca0",
"metadata": {
"id": "z4Voj7YeTca0"
},
"outputs": [],
"source": []
}
],
"metadata": {
"colab": {
"provenance": []
},
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"name": "python",
"version": "3.11"
}
},
"nbformat": 4,
"nbformat_minor": 5
}