{ "cells": [ { "cell_type": "markdown", "id": "05-00", "metadata": {}, "source": [ "# Sequential fits\n", "\n", "*Written by Claude Code, Anthropic's coding agent, for the rietx project.*\n", "\n", "An in-situ experiment measures one specimen many times while something changes: temperature, time, pressure, gas.\n", "rietx fits such a series pattern by pattern, starting each fit from the answer before it.\n", "That warm start makes each fit fast, but it also means a mistake can travel down the chain.\n", "We fit a synthetic heating ramp whose answer we know, then break it twice to see what the series reports.\n", "\n", "**You need** rietx 1.7 or later, which the next cell installs.\n", "Notebook 02 introduces plans and a single refinement.\n", "\n", "**Runtime** is under a minute on a laptop." ] }, { "cell_type": "code", "execution_count": null, "id": "05-01", "metadata": {}, "outputs": [], "source": [ "%pip install rietx" ] }, { "cell_type": "markdown", "id": "05-02", "metadata": {}, "source": [ "## A synthetic heating ramp\n", "\n", "The specimen is LaB₆ on a synchrotron capillary instrument.\n", "Its cell expands by 0.05 % per 100 K from 300 K to 900 K.\n", "We write the structure as a short CIF and make each pattern from the true model with Poisson noise." ] }, { "cell_type": "code", "execution_count": 1, "id": "05-03", "metadata": {}, "outputs": [], "source": [ "import csv\n", "import re\n", "import tempfile\n", "from pathlib import Path\n", "\n", "import numpy as np\n", "\n", "import rietx as rx\n", "\n", "work = Path(tempfile.mkdtemp())\n", "cif = work / \"LaB6.cif\"\n", "cif.write_text(\n", " \"data_LaB6\\n_space_group_name_H-M_alt 'P m -3 m'\\n\"\n", " + \"\".join(f\"_cell_length_{k} 4.1566\\n\" for k in \"abc\")\n", " + \"\".join(f\"_cell_angle_{k} 90\\n\" for k in (\"alpha\", \"beta\", \"gamma\"))\n", " + \"loop_\\n_atom_site_label\\n_atom_site_type_symbol\\n_atom_site_fract_x\\n\"\n", " \"_atom_site_fract_y\\n_atom_site_fract_z\\n_atom_site_B_iso_or_equiv\\n\"\n", " \"La1 La 0 0 0 0.4\\nB1 B 0.1996 0.5 0.5 0.4\\n\",\n", " encoding=\"utf-8\")\n", "lab6 = rx.Structure.from_cif(cif)\n", "\n", "TEMPERATURES = [300, 400, 500, 600, 700, 800, 900]\n", "A0, EXPANSION = 4.1566, 5e-4 # Å at 300 K, and the fraction per step\n", "TRUE_A = [A0 * (1 + EXPANSION * k) for k in range(len(TEMPERATURES))]\n", "\n", "true_instrument = rx.Instrument.debye_scherrer(wavelength=0.4139)\n", "true_instrument.zero_shift.value = 0.008\n", "true_instrument.profile.w.value = 2.5e-4\n", "true_instrument.background = rx.BackgroundChebyshev.with_terms(3)\n", "for coefficient, value in zip(true_instrument.background.coefficients, (40.0, -6.0, 1.5)):\n", " coefficient.value = value\n", "two_theta = np.arange(3.0, 24.0, 0.005)\n", "\n", "\n", "def simulate(a, seed):\n", " \"\"\"One pattern of the true specimen at cell edge a.\"\"\"\n", " structure = lab6.model_copy(deep=True)\n", " for edge in \"abc\":\n", " getattr(structure.phases[0].cell, edge).value = a\n", " structure.phases[0].scale.value = 5e-4\n", " y = rx.Refinement(structure, true_instrument).predict(two_theta)\n", " counts = np.random.default_rng(seed).poisson(np.maximum(y, 1.0)).astype(float)\n", " return rx.PatternData(two_theta=two_theta, intensity=counts)\n", "\n", "\n", "patterns = [simulate(a, seed=11 + k) for k, a in enumerate(TRUE_A)]" ] }, { "cell_type": "markdown", "id": "05-04", "metadata": {}, "source": [ "The starting model is what you would have before the first fit: the cell 0.1 % off, no zero shift, and peaks 50 % too wide." ] }, { "cell_type": "code", "execution_count": 2, "id": "05-05", "metadata": {}, "outputs": [], "source": [ "def starting_model():\n", " structure = lab6.model_copy(deep=True)\n", " for edge in \"abc\":\n", " getattr(structure.phases[0].cell, edge).value = A0 * 1.001\n", " instrument = rx.Instrument.debye_scherrer(wavelength=0.4139)\n", " instrument.profile.w.value = 2.5e-4 * 1.5\n", " instrument.background = rx.BackgroundChebyshev.with_terms(3)\n", " return structure, instrument\n", "\n", "\n", "LABELS = [f\"{t} K\" for t in TEMPERATURES]" ] }, { "cell_type": "markdown", "id": "05-06", "metadata": {}, "source": [ "## A series is a list of patterns\n", "\n", "`refine_sequential` takes the series as an ordinary Python list of `PatternData`, one per measurement, in the order to fit them.\n", "`x` is a list of the same length giving each pattern's place along the series (temperature here), and `labels` names each one.\n", "Nothing else ties them together, so the order of the list is the order of the chain." ] }, { "cell_type": "code", "execution_count": 3, "id": "05-07", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "list 7 PatternData\n", "PatternData of 4200 points, 3–23.995° 2θ, median step 0.005°\n", " intensity: <4200 floats 19…149878>\n", " sigma: None\n", " excluded_regions: none\n", " metadata: none\n" ] } ], "source": [ "print(type(patterns).__name__, len(patterns), type(patterns[0]).__name__)\n", "print(patterns[0])" ] }, { "cell_type": "markdown", "id": "05-08", "metadata": {}, "source": [ "## Loading a real series from files\n", "\n", "A real series arrives as a folder of files, and the work is putting them in the right order with the right temperature.\n", "Here we write our synthetic patterns out as an instrument might, numbered by scan, with the temperatures in a separate log." ] }, { "cell_type": "code", "execution_count": 4, "id": "05-09", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "['scan10.xy', 'scan11.xy', 'scan12.xy', 'scan13.xy', 'scan14.xy', 'scan8.xy', 'scan9.xy']\n" ] } ], "source": [ "folder = work / \"ramp\"\n", "folder.mkdir()\n", "for scan, (pattern, temperature) in enumerate(zip(patterns, TEMPERATURES), start=8):\n", " rows = \"\".join(f\"{t:.4f} {y:.0f}\\n\" for t, y in zip(pattern.two_theta, pattern.intensity))\n", " (folder / f\"scan{scan}.xy\").write_text(rows, encoding=\"utf-8\")\n", "log = \"scan,temperature_K\\n\" + \"\".join(\n", " f\"{scan},{t}\\n\" for scan, t in enumerate(TEMPERATURES, start=8))\n", "(folder / \"temperatures.csv\").write_text(log, encoding=\"utf-8\")\n", "\n", "print([p.name for p in sorted(folder.glob(\"*.xy\"))])" ] }, { "cell_type": "markdown", "id": "05-10", "metadata": {}, "source": [ "Sorting the file names gives the wrong order: `scan10` sorts before `scan8`, because names compare as text, character by character.\n", "Sort by the number in the name instead, and take each temperature from the log rather than assuming the order." ] }, { "cell_type": "code", "execution_count": 5, "id": "05-11", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "['scan8.xy', 'scan9.xy', 'scan10.xy', 'scan11.xy', 'scan12.xy', 'scan13.xy', 'scan14.xy']\n", "[300.0, 400.0, 500.0, 600.0, 700.0, 800.0, 900.0]\n" ] } ], "source": [ "def scan_number(path):\n", " return int(re.search(r\"\\d+\", path.stem).group())\n", "\n", "\n", "files = sorted(folder.glob(\"*.xy\"), key=scan_number)\n", "with open(folder / \"temperatures.csv\", encoding=\"utf-8\", newline=\"\") as handle:\n", " temperature_of = {int(row[\"scan\"]): float(row[\"temperature_K\"]) for row in csv.DictReader(handle)}\n", "\n", "loaded = [rx.read_pattern(path) for path in files]\n", "loaded_x = [temperature_of[scan_number(path)] for path in files]\n", "print([path.name for path in files])\n", "print(loaded_x)" ] }, { "cell_type": "markdown", "id": "05-12", "metadata": {}, "source": [ "Other naming conventions need only a different key.\n", "When the temperature is in the file name itself, read it with a regular expression and sort on it." ] }, { "cell_type": "code", "execution_count": 6, "id": "05-13", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "['LaB6_300K.xye', 'LaB6_650.5K.xye', 'LaB6_1000K.xye']\n" ] } ], "source": [ "names = [\"LaB6_1000K.xye\", \"LaB6_300K.xye\", \"LaB6_650.5K.xye\"]\n", "\n", "\n", "def temperature_in_name(name):\n", " return float(re.search(r\"([\\d.]+)K\", name).group(1))\n", "\n", "\n", "print(sorted(names, key=temperature_in_name))" ] }, { "cell_type": "markdown", "id": "05-14", "metadata": {}, "source": [ "A timestamp in the name sorts correctly as text only when it is written largest unit first, as in `2026-10-07T14-05-00`.\n", "Whatever the convention, print the sorted names and their `x` values once, as above, before fitting.\n", "\n", "## Fit the series\n", "\n", "`refine_sequential` fits the patterns in order, starting each fit from the parameters the previous one finished with.\n", "The `carry` argument says which parameters pass along the chain.\n", "It is a list of parameter-path globs, as in a plan's stages.\n", "A parameter matching one starts each fit from the previous fit's value, and any other starts from the starting model.\n", "`[\"*\"]` carries everything, and it is the default, written out here.\n", "Narrow it only when a parameter must provably not be chained.\n", "For example, `carry=[\"phases.*\", \"instrument.zero_shift\"]` would restart the background and peak widths at every pattern." ] }, { "cell_type": "code", "execution_count": 7, "id": "05-15", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "SeriesResult: 7 pattern(s), rietveld, direction=forward\n", " diagnostics: 2 unresolved\n", " INFO SEQUENTIAL_PERSISTENT_FINDING: CAPILLARY_OFFSET_UNAVAILABLE fired on instrument.geometry.goniometer_radius_mm in 7 of 7 patterns (100% of the series) — read this as a property of the model rather than of any one pattern: CAPILLARY_OFFSET_UNAVAILABLE in most of a series is the same finding repeated, and the thing to change is the model or the plan, not the individual fits. The per-pattern diagnostics carry each occurrence\n", " WARNING SEQUENTIAL_PERSISTENT_FINDING: PATTERN_UNDERSAMPLED fired on the fit in 7 of 7 patterns (100% of the series) — read this as a property of the model rather than of any one pattern: PATTERN_UNDERSAMPLED in most of a series is the same finding repeated, and the thing to change is the model or the plan, not the individual fits. The per-pattern diagnostics carry each occurrence\n", " trajectory (T (K)):\n", " [1/7] 300 K x=300 converged Rwp=0.0408\n", " [2/7] 400 K x=400 converged Rwp=0.0411\n", " [3/7] 500 K x=500 converged Rwp=0.0401\n", " [4/7] 600 K x=600 converged Rwp=0.0408\n", " [5/7] 700 K x=700 converged Rwp=0.0408\n", " [6/7] 800 K x=800 converged Rwp=0.0401\n", " [7/7] 900 K x=900 converged Rwp=0.0405\n" ] } ], "source": [ "series = rx.refine_sequential(loaded, *starting_model(), carry=[\"*\"], plan=\"mccusker_default\",\n", " x=loaded_x, x_label=\"T (K)\", labels=LABELS)\n", "print(series)" ] }, { "cell_type": "markdown", "id": "05-16", "metadata": {}, "source": [ "Every pattern converged.\n", "The two diagnostics are `SEQUENTIAL_PERSISTENT_FINDING`: one finding repeated in every pattern, which is a property of the model and not of any single fit.\n", "Here they are the step size and a capillary correction this model does not use.\n", "\n", "## The trajectory\n", "\n", "A series answers a question about how a parameter changes.\n", "`series.trajectory` returns one parameter across the series, with esds." ] }, { "cell_type": "code", "execution_count": 8, "id": "05-17", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "300 K a = 4.156602 +/- 0.000004 Å true 4.156600 off by +0.6 esd\n", "400 K a = 4.158681 +/- 0.000004 Å true 4.158678 off by +0.6 esd\n", "500 K a = 4.160757 +/- 0.000004 Å true 4.160757 off by +0.1 esd\n", "600 K a = 4.162841 +/- 0.000004 Å true 4.162835 off by +1.5 esd\n", "700 K a = 4.164909 +/- 0.000004 Å true 4.164913 off by -1.1 esd\n", "800 K a = 4.166991 +/- 0.000004 Å true 4.166991 off by -0.1 esd\n", "900 K a = 4.169068 +/- 0.000004 Å true 4.169070 off by -0.5 esd\n", "fitted expansion 4.9979e-06 per K, true 5.0000e-06\n" ] } ], "source": [ "a = series.trajectory(\"phases.0.cell.a\")\n", "for label, value, esd, true in zip(a.labels, a.value, a.stderr, TRUE_A):\n", " print(f\"{label:6} a = {value:.6f} +/- {esd:.6f} Å true {true:.6f} off by {(value - true) / esd:+.1f} esd\")\n", "slope = np.polyfit(TEMPERATURES, a.value, 1)[0]\n", "print(f\"fitted expansion {slope / A0:.4e} per K, true {EXPANSION / 100:.4e}\")" ] }, { "cell_type": "code", "execution_count": 9, "id": "05-18", "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAABjcAAAHLCAMAAACzlzrhAAADAFBMVEX////+/v79/f79/f38/f38/Pz7/P37+/v5+/z6+vr5+fn4+vz4+Pj3+fz2+Pv1+Pv29/j29vb19fX09PTy9vrw9Pnz8/Py8vLx8fHv8vXv7+/u7u7t7e3r8ffp7/bn7vXr7O7k7PTr6+vj6/Tp6eno6Ojn5+fl5eXi6vTf6PLh5enb5fHa5PDj4+PY4+/i4uLg4ODf39/c3+Ld3d3c3Nzb29va2trY2NjX19fV4e7S3+3Q3ezP3ezM2+vJ2OnV1dXT09PS0tLQ0NDOzs7F1ujB0+e90OW7zuTEzNW3zOO3y+PKysrJycnHx8fFxcXDw8PBwcG/v7+9vb28vLy7u7u1yuKxx+Guxd+rw96qwt6owN2kvtu5ubm3t7elutKeutmcudibt9iZtdezs7OwsLCvr6+srKyrq6uqqqqWtNaoqKinp6elpaWioqKfn5+enp6cnJyampqWlpaTstWRsNSOrtOMrNKKq9GHqdCFp8+Cpc5/pM18ocyIma13nsp1nMlymshul8ZrlcWTk5OSkpKQkJCNjY2Li4uJiYmHh4eEhISDg4OBgYF/f39+fn58fHx7e3t5eXloksRkkMNijsJhjsFejMBcir9ZiL5Whr1Vhb1ThLxQgrtNgLlKfrhIfLhGe7dDebZ3d3d1dXVzc3NxcXFubm5ra2toaGhmZmZBd7U9dLM6crI4cbFlZWVjY2NiYmJZYm00brAybK8vaq4saK0rZ6woZasmZKslY6okY6okYqkiYakhYKkgYKggX6gfX6hdXV1bW1taWlpYWFhVVVVUVFRTU1NSUlJPT09MTExKSkpHR0dFRUVDQ0NCQkJBQUE/Pz89PT07Ozs6Ojo4ODg3Nzc1NTUzMzMyMjIwMDAvLy8tLS0rKyspKSknJycmJiYlJSUjIyMhISEfHx8eHh4cHBwbGxsaGhoYGBgXFxcWFhYVFRUTExMSEhIREREQEBAPDw8NDQ0MDAwLCwsKCgoJCQkICAgHBwcGBgYFBQUEBAQDAwMCAgIBAQEAAAAKo9G7AAA16ElEQVR42u3dd3BU9f7/8c8aMIKoFCkDUUYEEkcRM8oFjRIgKCOhg3QGpIuAgAVUBimXJoigREEFEaXYGevlp+ZLFITQRKpICIQmxNBChyTv3ylbztmScL1A9iTPxx/svvdswu5h+bz2fE55KwUAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAKBZWbh9Q+C9i6fZhlhsAQBjbIxOu8G+M6J+yad1HCf/Nj/wur1luAADFKTdqbhVd7ixyAwDIjYKV3C4XZrXvu1ZkErkBAORGgYaJ9NNvU+R0bXIDAFRR2A8eOeSXTeterenNDZe2Q2LDhw+4nxA78cv1m1ZOqOoumyxct3nNkv6lzars0JTfNi7tEHSZaYdsMm5rX/Df4CgzZNmGDV8OqxzkN4XIDb+XAgAolO2L1/40dkAcbeWup2026nM9jHqhmI61M8qZeWZ5xKja/GVWqeUCl5lqioww7/3mDhCPdu6fPV0l8DcFzw2/lwIAKJzcuCir+radlCkn65j1hZyP+7SbdEyyyuv1l7tnDGnXZ8p+ORmtVR1E1g5u22fEqpP6shYXZP+kdh0WnpVVAcvc2oskekb9Y9a/t81FyZrRo92wL7JrBPymELlhfykAgELKDZmp396XJSnu2tgfkZAng/TbaubTyu0Q/YCoJbLBrPW5oog9srOSOxza+C1Tvt0btcx7I0Qq+x6P3C/p5oIKkQG/KURu2F8KAKCQcuNwKePORDlf3qjXmwt22kfnYbJN+/MLWep7qIdII/PeavnCb5nHRJHy5r3BIpbNhH4iLUL+pgL2i5svBQBQSLnxhXnnPpHWRv2qZwQ3Y6D8kI/Xb92+fa/8ZWw0nB1azvOjH8l+bzrs8lvmMVnEvZt8kMh9vscXS4blWX6/KVRu2F4KAKCQcmOK+94lY2LKexzuz/KlftPKvcNa5IRWldH2oZ/bvLB7pJksx7ds2arZlqHvC7ct89Bmp6qY94aKRPkeT5VUy7P8flOI3LC/FABAIeXGRPe9MzIkMDdqHJfDswa2bdlyihi7uyu8mq4P3If0fSDrxOek3zKPASL1zXtT5aIlUdaZO8C9lf03Bc0N/5cCACic3HDvxqgg0icwNybI0ZruzQXPYF2736IjktNKf8ZvLT0S/ZZ5xIv0Ne+lSLrlr02RzZbK/zcFzY0gLwUAUAi54R7Au5v7H/xyY4msMcuPrYN1lX36zo9J3r0SVuYyj5JZssS4E/GXfKys+8vPWg678v9NQXMj+EsBAFzr3MgztxXWS5oKzI1FssOoYk/bBus1slKpuhck2HXOjWVeS+VkTffhVNbT9epelEWWyu83Bc2NEC8FAHCtz984pM0rVVgi5vkafrkxUGSKtlui1T5zx8P04fqxtBH9zshk7XaWXJym7+t2NZnex3/ZoO3bjfMI7zsnW7TH2x93ny7uefwjkcX6GX/xi6oE/Ca/3HD/iN9LAQAUUm5szJX9f54R8/Ap/9yI3CFy/M9DIuvcO6wlO33XMZGd+lkZEV9q10fP+OPgBdF/yL5MO3GjnvF7BufJxbT9WjrV8ZzQYTxearVITnpatkiNgN/klxvuH/F7KQCAwrr+bY807TCmvydFBMsNFfXzJW3pzj79jMG6zdJD+jFPh2dVMp/U93d9ad7O6XX9l3lzQ7XZql236kxKTWXLDeUaqv+1kjYpMuA3Bc8Nv5cCACi866bXSqhf0l03aem+LsgDLe8371Rq0libQqrmOWSqcnxCHcvPV2iQUMd7fK1lWe2WLct4Ho5q0sB7QqDt8cbxFYL+pviWdSw3vh/xeykAgKLQFxYAQG4AAEBuAADIDQDA1eXdDw4AAAAAAAAAAAAAAAAAAAAAAAAAAAAAcIyyiYmxdQAAhSs2MbGsQ3IjUQAA4SCR3AAAFMHciNVeKhuIAFDItC/xsQ7JjToiddjLAwAMxuQGAJAb5AYAkBu8VAAAuQEAIDcAACHFdOnVJYbcAABclrik5BUrViQnxZEbAICCNV++wrS8ObkBAChwa8MTG1pwxJEbAIACJK3wmU1uAADyF5NsyY3kO8kNAEB+qo9eYdWZ3AAAhHBL06fe+HyFXU9yAwAQ6PoHu72y4KcVgdjeAADYue544sW533ty4rs5o/+P/RsAgOAqNn32za88GfHj/DFd7inB8VQAgGBueKjb2A+8h059Nv7Jhjcrzt8AAAThiun40rv/8WTDt2+/0Ow2xfniAIBgKj76XNLXnsj4Yd7ozndfF/ikuNnG9almc30qACjObozrOe4T7xTUp+N6PXJTyOfGdO7Z+U7F9XABoJhy3dXp5fe8uy2+eev5x25jMCY3ACCo2x57/q1vvDu633u5410uBmNyAwCCuemRXuM+9c5MfTK2Z9wNDMbkBgAEc93dnUfP+8ETGV/Pfu7RigzG5AYABHVbsxfe/tYTGf9596WOMS7FYExuAEAQNzfsPf4z70VCPnil20M3KAZjcgMAgihxT9cx83/0ZMZXbz7b9FYGY3IDAIKq/vioOd95IuP7uaOeuMOlGIzJDQAI4lZr04yfFozp+uD1DMbkBgCoAptmfP5G74a3MBiTGwCgQjTNeMfSNGPk49UZjMkNAFDBL01oa5ox2miawWBMbgCAKrhpxiM3MxiTGwCgLqdpxtW+NCG5AQCqqDTN6BSsaQaDMbkBAFrTjIcvv2kGgzG5AUAV96YZ82xNMyo65rWTGwCgin7TjDAfjKNq1/R/qIL9IVdV91qqVtujslZFeoqy5AYAVTSbZowvnKYZ4Z0biTly1FqXbDdz2yU54Hug728X5fTGfvrdjeKxWKsaeYqnyQ0AiqYZxSU3yqRfsOdGrBYEltyI+EKrT5wXidSKVZmGv0X6Gblxyaz7khsAVFFqmjGysJtmhHVuLMx71Z4bdTZNb9XDlxuvSe5rNVXJNmtK+p7SQ/4uZeTGLvZvAChSbgmTphnhnBstcha3tueGroM3N+pelIkBP7RaPlLkBoCipcS9Xce8Hy5NM8I4N8qm76uQb268Jvsi/JfWvCjxntyoWpbcAOB41ZuPmhtWTTPCODc+ymuj8s2NzbJQqUrlbUsny1Zl5sbpLJH9C6PsL9EtkdwAoGiaUdRyo1WudlxUvrlxTPr3Txc5OMsSHekyzJ0bkpet7TU/GGv9WfEhNwCocG+akRyWTTPCNjfK7d1XvoDcyJGPRU7lieys5FnYWk6Z9+svSSyroiaekB3kBgBF04zikBuLtFmqAnNDUmNVueFn9fkqU4p8aXtyi1xpxTwVAEXTjKKfG61yU/SzvfvKsdq1S4XIjROyv4x+O0UOuh+pcsYWE5otMon94gCUQ5pm+C5NGMZNM8I1N6ZappRahsiNP2SZcdtYpIz5yAhJ9/s9q2QGuQFAhX/TjPeWO6VpRrjmxgTzbO9jkpeZmRAiN5bKL8ZtC8mL9CSJ/9bFLvk3uQEgjFV8zGlNM8L8uoae/RvaZQpdAbnRTrKMw2wXS5r5QBO56LnmofvEjkEijckNAIUtpkuvLjGqiDTNcEhuJIjUUMb1cWsPlEPajg/j0NtNktb3gcTFee5jb9USWeP5wa1L+ybEtvsiR9Zy3h+AQhaXpB9Qm5wUp6xNM0Y7tWmGw3Jjl2efxwS9it5rFkvNJ5c9Lj08P7jb/bzNVckNAIWruScfljdXZtOMFxzdNCPMc6NFprmjOz4z05iS2pjp9ozxcPkpWw/uXdXd/eS+mWneCxxWGfrrzoO7V/aP4DojAAp5a8O7WbFiecNHnnR+0wz6/QHAVZW0wuf/ikLTDHIDAK6mmOQVdo5vmkFuAMBVdPNIW2i8UwSaZpAbAHB1XHdP1zHzf7RvbPRktZAbABDE7Y+PnPPdikCdWTXkBgAo/3auvqYZyQvG/58lNpLvZP2QGwCgQrRzfeMpvZ2r9Xiq2awjcgMATNUft7dzbV49yPkbcawncgMAtJkpazvX98d0vbdEfueLg9wAoIr1zNSCn3ztXJ8K1s41brZxfarZbG2QGwCKterNR821tHMdlU8715jOPTuzS5zcAFB83dr0qTe+Cj0zBXIDADyuf7DbK9aZqd7BZqZAbgCA1jTjjidsM1MjH7+dlXKNB2NXgx79TOQGABXuM1Nvememfpw/pss9zEwVwmDcbq+nC5OQGwBUWM9Mea9q+9n43g1vZqUUzmDc+LzsEFm7T2TjanIDgArPmakX3/HOTH07Z2QzZqYKczBOkWVK39IYeCa9CrkBQIXdzNSz1pmp0V3uvo6VUsiD8UFJMHJDTZGF5AaAMHL9Q93GfmCZmXryEWamwmIwviiVVJ6UUup+ySA3AIQH151PvGSZmXr7hWa3sVLCZjA+JpVVtkQrVU7OkRsACl/FR5+d7Z2Z+mH+6M7MTIXZYPyHNFY7RDsEt5WcJDcAFKobtJmpT7zXrf1Um5m6iZUSfoPxQhmhXpX97Vv9KWvJDQCFxRXT8aV3/2OZmXqMmalwHYzjz29TlffoZ2+cbkJuAFCFMjP13OyvvTNT85iZcsJgHLVwx66UeK4zAkAV7szUuF7MTHHRJ3IDgAo5M/WetwvfN289z8wUgzG5AUCFnJlKss5MdbrLxUphMCY3AARzQ1xP28xUz4eZmXLqYDx58mRyA8DV5Lqr48v2mamKrBQn5wbXwwVwFVV87Pm3vvFExvJ5LzMzVQRyY9euXeQGgKvgxod7jvPNTH2izUzdyEphpwG5AUAFnZnq9PI838xUEjNTDMbkBgB1WTNT773ckZkpBmNyA4AKOTP1qW9mamzPuBtYKQzG5AYAFXxmavS8HzyR8XXSc48yM1WEBuMD/sgNAP+L2+wzUy91jGFmqojlhvgjNwD8Qzc90ss2M9XtIWamimJu1PZHbgD47113d2frzNRsZqbYv8FLBaBCzky98Pa3nsj4z7vMTJEbpqjft6+3PRCROGX99lTLE/796++/TKii371/u1sHc1G5EWs2rxoaSW4AqujNTD05npkpciO4lSJHrXX0cX1niG8vet+jxt6R7JLa/UaefSVPm081OkLJ9irkBqCK7szUs8xMcTyVTX9JtedGbO7OD6f7cqPHJdk5rE2fRScjjdzIaGmI0he5tsvB4W0nHpc15AZQRNzWzD4z9cSdzExxPJWfGllb29hzo3Q1pTp4c6P8YVmtb2moysrIDetVrvrK6fu1mxa5kkBuAA4S06VXlxgVbGbqM09kJH/AzBTHU4U4nmrV2fqt7bmh8+XGMMmO8j1uz41fZKlxu04WkhuAY8QlJevBkBSnLDNTXUbP985MffXms02ZmWL/RiiDZJLKNzdSZZmy5sbhL39fv6iVWe2T/sbtRNlEbgBO0dxzHcLlzY369mYj53hnpr5/58Un7mBmitzIR82srRH558YheSZ61oYNi1u4c8OUUlorSuZKvPFoDzlkfYluieQGEI5bG97L165Y3rBhb8vM1IJXuj14PSuI3NDViE8ItWjN2foq/9w4LxP+NqLCmIpqtH/WoHbDNokxQVVVJNp4TkvJVkH3q5AbQNhJWhHEV28+1fRW1g25YXKN2GPsEk/ZWCtw4WBtlir/3HDlSXbmpA6D1okM0coy5qkaUyRX2yFeQ6SmUSbIOXIDcIaYZL/I+H7uKGamyA2riFUiB/XcWCwjAhbWOrolooDcUKflVKwyzvLYaQmjDBmunfQnYixSbSSLeSrAEW4dZQuNd5iZIjcCDJfDLZWeG+1kfcDCKbJPP/c7XS5t3x4fIjf2yirjtrXkRviW/yoztD9PSBv3Zksa+8WBcHf9g91eWfCTfWOjJ6uF3AiwQwYrIzdiJTNg4VTLlFLLELnxq3xs3MaLlPUt/02maH9ulslGtVh+ITeAMOaK6fjSO98H2a/RmXVDbgQ4LTXM3CgjlwIWRpvnfg+Xky1bVgqRG0MlzdjOGCGHfYvrXZDuRu4Yy8pl6uFEbgBhqeKjzyV97cmJH+eP+z9LbCTfyfohNwKcE+3aUXlabkTLyVDP8ezf+Nd274WmLOeLH5LF5ZVqf0wWadWs7vq54332y75S2m2tU7K4tKqwSv4qR24AKizbuX7iDYlPxz/5yM3246lms47IjUB7pJWZG91lR0G5kSD6xol26JW+y+OCtuNjgBEhF+VM2mGR9KrGr7t0IC1bu8iheTrHMyLH/zwll3pwfSpAhfGlCb956/nHblOB52/EsaLIjUCLZF2EyhXlWi/TLzc3dnn2eUwwHm61Tbt7emlV44pUqee14lRKrPtH+2do5a52XEcdUGF6acLl773c6S5XfueLg9yw0y6KvrFfrrRdI0eiQj2nckvzvMDyLVsap2c0bunmOeOjZkJ8Kc+TS9ZNqFfS8sP3JUTTtwkIF7c07P3GZwU1zYibbVyfajZbG+RGcK2yzG2HrBb0+wNUkT7MtuuY93+yXJown6YZMZ17dmaXOLkRUo3p2//O2jmzJn1igaLKdccTL879nhPAQX9xAOoyDrN99s2vPJHx0/tjut5bgpUCcgNAMDfE9RzrO8z2s/G9G97CSsH/PhhHLXYfRjV9cTVyAygqXHd1evk97wG13779QrPbWCm4QoPxVJll3lloXBmE3AAc77bHXnjrG+9htvNGd7rrOlYKruBgvNPT+jsh9Hl/5AbgEDc3tHQAX/HJuJ4P38hKwZUejLOlsnmnshwjNwDnKnFv1zHzf/RExtezn3uUDuC4OoPxGU9uVJGz5AbgSK7qzUfN/c7aAfxODrPF1RuMdxsXrlV6D/B0cgNwnFubWg+zpQM4rsFg/JFsN9pmlNthXM+W3AAc44aHuo39wNvV9fM3OMwW1+j6VCdk1+D4+MFpcjya3ACUUw6z7Wg7zHZks9tZKbhmg3Hb4+b1qY624bw/QDniMNvnfYfZ/jBvdOe7OcwW13gwrj1rZ1bWHzNrcb44EO5ueqTXuE+th9k+chMrBVxnBEAw193TZfR8b5+lr5Oee4zDbEFuAAiu+uOj5ngPs/3Puy91jOEwW5AbAFTww2yfeuNzT2Qka4fZBuuzBBTOYDx58mRyA1BhdZjtKwssh9k+1fRWVgrCajDWDqkiN4Dw4Irp+NK73j5L380Z9Xh1VgrCbzDetWsXuQGoMOiz9FzS197DbOeP7nIPfZbAYExuACr4YbY9x/n6LH06vheH2YLBmNwAgrvu7s6j53kPs/3mrecfo88SGIzJDSC425uNfPtbb5+l917ueBeH2SKcB+PIBj0GDuzRIJLcAArBLQ17v+Hts5T8wVgOs0XYD8YJK0+ZV6c6tTKB3ACupesf7DZmwU+ezPjqzWebcgI4HJAbk/JEcg6lpx/KEcmbRG4A14brzidefMd3mO3cUc2rMzMFZ+RGd5HUdmX0e2XapYqngRO5AairepjtbG+fpR/nj+l6L4fZwkG5sUEW+4rFsp7cAK6mGx/uOdZ3mO1n459seDMrBQ7LjRNiadUULcfJDUBdrT5LnV6et9x3mO0LHGYLZ+bGOSnvK8rLWXIDUFejz9ILb39jOcy2E4fZwrm5kS4DfcUg2U1uAFfWzQ17j//M12dpbM84DrOFs3PjVTk+wP3FxzXghEwlN4ArpsS9Xce8/6PlMNtHOcwWRSA3ymwVObRy1tSps1YeEtlamtwA/lsxXXp1ifF7zHXHE6Pmeg+z/X7ui0/cwcwUisr5G+U+Oidu5z4qx/kbwH8pLknvj5GcFOd9pGLTZ9/0Hmb70/tjuj54PasJRes6I5UHzFq5du3KWQMqc50R4L/V3HOE1PLmyuizNPaDZN9htr0b3sIqAhcLJDcAy9aG98DaFT889/J73urbt19odjurB+QGuQH4SVoR4Id5ozvdfR2rBuQGuQEEivHLjE/G9Xz4RlYLilFuTJ48OcSSiMnT/K95GD9h2gRLGTt08uD73febDJr07+5lzPs1p7nFkxtQRa7P0hxbarzIYbYodrmhHVMVYskUkaPWuuovR7QnH/DWddfqB2TlrdOPNJxxyDg6K2uosaSR52itp8kNqKJ1Avi3/jNUPVkvKHa5sWvXruALGpzba8+NWJH0db7cuO+gnE/9aGmG6L2fdsufqxalZIoMNXMjc6ahMbkBVUQ6gD85/tMVwXRm5YD9G24lt2f29dveGBytOnhzw7VZMurptz0itD/7G395pU1y2GXkxi72b0AVyQ7g2gngnZItsZF8J2sI5IbbVOnX2p4bOl9udJDcJgE/1Nh8FeQGiojbmllmpr5/x30CuPV4qtmsJBS33CgdXzP4gvjzv6h8c2NZsL4dVURqG7mR3m/yhHalyQ0oR89M+TqAL3ilm+8EcMv5G8vjWFMoNrkRtWSG9mfbLMl5PdjikjuOROWfG+kySbWaOLG7LRwGSoby7Rc/bN8tXsctkdyACvOZqS6j51tnppreqvI5XxwoJrkxRRZqFzg8kntQpJ0Kdsncvir/3Dgt/dfp4XCog29p1QMy2MyNtJRFqdnuveTKd+SWB7kB5aiZKX9xs43rU81mawPFKTe2SaJS/WSEmiFrApc2uvCzyj83XNrBVfLb4pRjcj7Bs7DMJkkx7tQwHqq1Tk5UITegnNM0I+TMVKCYzj07s0scxSs3MiVK6y1+obx2dO2+gIWRO7VZqgK2N86LDNAjIk3WemJjnayzzVpVPi6DmKeCE1x3jzYzZWma8ZT/zBRAbqhLol0/fesObbQP0id2uKzXz/b+Qs5MmxYdIjf+ko3GbR85Yz5Qbr2sK2P/PWtlGvvFocJ+ZmrkHN/M1NxRT1SnaQbIjeDbG/GqwvklStWWTBV4DK5PyxC5sUGWGrf1xexWXnmzrPWLDZUq08kNKIe0cy1oZgoo7rmRKitrvK7v/G4h2wIWtjHP9k6RszNn1gmRG1Pdx+G6tzeqbpdU/+NuKxw1d5OTG1DhODM1xjIz9QYzUyA3CtIqR9uW2KNdImSavBbqOZ79G7WmTSsXkBt1z+f10G6q7JRU7abGTllTyvuD3Y2nV1stp6LIDYSf27WZqe9sM1OsE5Abl6HH5r0r62q3aTnxBeVGgkgN/XbYtGnL5IS24yNRGYfq5q1d+HGmnGmkFX9I9mvmJXDv06o92WuXzFp5XGQC5/1BhdvM1BvemamfFoxhZgrkxhXlnxu7PPs8zDh47aJ+f39r/f5J7/6Qtlq1KNu4mzmU88URPkrc09UyM/W5NjNFO1eQG1f6d9adOdW4rT1zprHre8JMtxbm8trPTJ/SPdK4+6pn0cxYvSzXfcT0ie25zgjCxe2Pj7LNTDVnZgr4R4NxjfgE+v2hqLtFm5n63Dsz9f6YrveWYKUA/2gwdo3YY7RsStlYi9xA0VTi3q5j3rfOTDVkZgr454NxxCqRg3puLNauNUJuoMipbp2Z+m7uqMeZmQL+x8F4uBxuqfTcaBfseujkBhQzUwC5YbNDPyVPz43YIOeLkxtQDp6Z+sk3M9WbmSngig3Gp/VDa/XcKCOXyA2oojEzNdc3MzWHmSngCg/G50S7xHmelhvRcpLcgHL6zNRT9pmpe5iZAq74YLxHWpm50V12kBtQjp6ZWsDMFHANBuNFsi5C5YpyrQ96zVpyAw5Qvfmoud/7ZqZGPn476wS4eoNx9HHZ2C9X2q6RI1HkBhznlqaWmakf5zMzBVyDwbhVlnk9qawWnC8ORynxYLdXfDNTn41nZgq4VoNxjenb/87aObMm1xmBc1R/wjIz9e2ckc2YmQKKw2BMbuCfuFWbmfrKMjPV5Z7rWClAIQzGpePZ3kDYuz5gZupmVgpwzQfjqCUztD/bZknO6+QGwpfrDmamgDAZjKfIQu1c8SO52rUN25EbUGE6M/WmZWZqdJe7mZkCCnEw3iZas9d+2rVwZ8gacgMqHGemkn0zU08+wswUUNiDcaZop20svlBeu67hPnIDKqxmpl58xzcz9fYLzW5jpQDhMBhfknJKbd2hX9fwLLkBFS4zU8/6ZqZ+mD+6MzNTQDhtb8SrCueXaE3CuY46VJjMTH1gm5m6iZUChNVgnCora7wufZVqIdvIDShmpgByowCtcrRrjOyJVGqavEZuQIXHzNQ8ZqaAMB6Me2zeu7KudpuWE09uoDBc/1C3sb6ZqU/H92JmCmAwJjegYrr06hLj/6Drzidess1MPcbMFMBgTG5AE5ekb1AkJ8X5Hqr4KDNTgEMHY1eDHv1M5AaukubL3fmwvLle3sDMFODk3Gi3VzzIDVylrQ1PbGjB0azjS+8yMwU4OTcan5cdImv3iWxcTW7g6khaEUibmep0l4t1AzgvN1JkmdK3NAaeSa9CbuCqiEn2y4xPxzEzBTg3Nw5KgpEb5oVxyQ1cca6YMbbQeJOZKcDZuXFRKqk8KaXU/ZJBbuAKq/jYc0lf+21s9GS1AM7OjWNSWWVLtFLl5By5gSvnxod7jvskyH6NFZ1ZN4Czc+MPaax2iHYIbis5SW7ginDd1Wn0PO8hVN8kvfx/lthIvpMVBDg7NxZqLZtelf3tW/0pa8kN/M9ue+yFt77xHnT73ssd9WOmrMdTzWYdAQ7Pjfjz21TlPfrZG6ebkBv4X9z0yJPjP/Xmwydjuz10gwo8fyOOFQUUhYt3RC3csSslnuuM4J+67u4uo+f/4MmGr2c/+2hFlc/54gC46BO5UYzd3mzknG89kfH9Oy8+cWeQs/niZhvXp5rN1gbAYExuFGe3NOz9xufeHd4LXun24PUhnxvTuWdndokDxWUw/lejBv4PRTeqby2r1KvgK6L+ZT3z3F6RG0VFiXu7jnn/J09mfPXGU01vZaUAxSQ3ai/ccTDTFOIZ7fLkqLUuPXBphsgBb+0aliaSt2uYWXXfqRWbW6hgFblRNFRvPmqu99KE388d1bw66wQoRrlRL0sKuB5upf2n7bkRqz0125cbpdaInMvIFInUq0F5cmbfebnQVgVW5IYqAu1cn/I1zfjp/TFd7y3BSgGKWW78KocG31/DFPwZS3Mm2nMj+tcR9Tv4cmORXJhQTjug94uSWhF1UpZVUFFrZX+ZgIrccLbrH+z2ygLvFQo/e6N3w1tYKUBxzI3j0j3/J7TPm9Xanhs6X278K1eG+h6fLOn6ZkflozIooCI3HMuvneuckc1uZ6UAxTY3zkuFfJdXPrC7dL65sVDSLY9vkqnG7RJJCajIDUeq+Ohzs33tXOfTzhUo9rmRLnXzXf7FpQSVb25sk5mq5H11Syr3ZRLNXRlD5M+AyvsS3RLJjTB3w0M9x/ouTUg7V4DcMEyRf+e3uLuWCvnnxnHpN/ywtqP84yilX1RX7jcebSNZ/pWXb0c8uaHC99KEHV9+z3dpwreep2kGQG64lVqfPTD00ioH00rnnxuuXFkpuQeyRfZqwVFDpKbxcIJ+VXZ7RW4oxzTNeN5yacJ5L9POFSA33HYZ0rVthTTzbuBTvtRmqQrY3rgkklJVqX7ZslipqqK38tC0kFP+FfNUygmXJrQ2zfhkXM+Hb2SlAORGsO/9Ic7f6CC/NNIMkRONGpULkRvHzGOm1L/lsFKReRLvnt8KqNgvrsL80oRa0wzfpQmTnvO7NCEAckM18hfwjKmWUGkZIje2yxfmbxMpq9R+6WtUE2Sz8q/IDRXGTTPe9l6a8D/vvtQxhpkpgNz4R4btNuyXnN27G4fIjcWy0rhtJbnaMVUrZYlRpcpHyr8iN8LRzVrTjM+8M1MfWJpmACA3/jHP/o2yjRqVDMiNVnLC2Inxpfyh/TlAsvXjepvkSIuAitwIM9fdozXN+NF7acI3n23KzBRAblyWhEmLFk1KuIzcSBAxrkXyQKNGw+SwNrGlH3qrtKuIDGvRY6XI01oR8afsezphWKasUwEVuRFGbn985JzvrE0z7mBmCiA3LtMDm839F5sfuOzc2OXZ5zFBr6rtMO7nzTCeE3vAqNKiglTkhnJa0wwA5Ebg849ITuqHC1Nz5HB0yA2S3eZu7fq7d1fVb1fvdhtsPBw5LHXn1qWJ7idXmbJh57oJZYJW5EYhK/Fg1zELvE0zPqdpBoB/MBgvkwyjA1ODDFlGv78izHXHE5amGd/NGfU4TTMA/KPBONNzPdzucoTcUEW1acazvqYZP84f0/UemmYA+MeD8QX9rAtdWTlPbqgieGnCbmM/8DXNGN+74c2sFAD/02C8X9ztmmpIBrlRpLhiOr70rq9pxtsvNOPShACuwGA8SyaadybK6+SGKkpNM772Ns2YR9MMAFduMC63+cJUrXNThakXNpcjN4qCG+JsTTPG0TQDwJUdjPXr4eYeOJgrkh7qorjkhnJM04xONM0AcLUH44IviktuKGdcmtDaNOM9mmYAuFqDccEXxSU3VNg3zeg17lNf04yxPeNomgGgSA7G5Ia6Ak0zOlubZsymaQYAcgOh3NaMphkAyA1clpsb9vY1zUimaQYAcgMhlbin6xiaZgAgN3A5qj8+ytI0Y+4ommYAIDegVEyXXl1i/B+8pelTvqYZPy0Y05WmGQDIDejikvSrECYnxXkfuf7Bbq9Ym2b0bngLqwkAuQFTc8/p3subK6NpxovvWJpmjKRpBgByA7atDe9VQlb8v262phmju9A0AwC5QW74S1oR6LPxTz5C0wwADMbkRhAxfpGhNc3g0oQAGIzJDRWqacabttQYQ9MMAAzG5IYK1TRj3Cf+E1Q9WS8AGIzJDZV/0wyrzqwcAAzG5IbKr2lGr2RLbCTfyccZAIMxuaFCNs24we94qtl8mgEwGJMbKv+mGZbzN5bH8WkGwGBMbtz2WP5NM+zniwMAg3Exzo2bHnnyMppmxM02rk81m60NAAzGxTg3rruny+jLbpoR07lnZ3aJA2AwLra5cXuzkXO+9TXNeJGmGQAYjMmNEG5p2NvXNCN5wSvdaJoBgMGY3AiuxL1dx7xvaZrxVNNb+XACYDAmN4Kq3nzUXG871+/mjqJpBgAGY3IjhFu1dq7ephk/vT+m6700zQDAYExuBGW0c/VeJ+Qz2rkCYDAmN0Kxt3P9ds7IZrfzSQRAbpAbKnjTDEs71x+0dq40zQBAbvBSg7vhoW5jP0i2tnO9iQ8gAHKDlxqMK6bjS+/+h3auAMiNYKJ379the6BCj4927PvdXfxrn1sfvUrxVPumBywrKrlR8dHnkr72zkzNG93pLmamAJAbVutEjtr+iouiOeCuGonb03q10VPJiIBlRSA37O1cPx3Xk5kpAORGgGckxZ4bsafXT5loyY30eoZKxraJeb+/nK0SsMzZuWFv5/rNW88/VpFPGwByI9gvPLmhtT03Iksp1cGSG7sCf2iJpIRc5sDc8G/n2vEuLk0IgNwIZX12Xb/c0OWfG2VPSJsikhvB27kCALkRyjAZrgrIjb9Td/+R0sf6DXyIZLhCLXNOblx3d6dQ7VwBgNwIoW72elVQbpjWlfct/V2mhlymvUS3xDDOjYLauQIAuRHchuy6BeXGn5PaN+67SuRn78J6eTnRoZZpxKdOuLdzXRGqnSsAkBvBjJBhqoDciFTu+Sxp4Fm4UNaGXBbWuWG0c/3hMtu5AgC5EeA+fZaqgNzw2G2csWGkxRHpF2pZ+M5T2du5vkM7VwDkxj8wRTL1s73/ktx9+xIKyI1VMsN9r59kRoZaFpb7xWnnCoDcuEKmWqaUWhaQG5tlivKcX/5hyGXhlhu0cwVAblzJwbiqebb3IDler16Z/HOjSY60M+9F59j2ZtiWhVVu+LdzbU47VwDkxpXh2b8Rv29ftYDcWDo4SjvRb2im7I4wH5gmW1SoZeGSG7RzBYBrkRsJIjX027X79h2WS9qOjyFasUfk1OE8kb/jzSe79ssznh/0XxYOuWFv5/o57VwBkBtXPzd2efZ5TNCKNl8e0e4dWlTTuyVyqoLnB/2XFXJu0M4VAK7NYFy2XqxxW7pePWO+qW49tyrm8nJ1q/qeXK1etPVnbcsKMTdo5woA9PtTtHMFAHJD0c4VAMiNwniptHMFAHJDXXY717G0cwUAcsNPTJdeXWIU7VwBgNy4HHFJ+r7u5KQ47yMVaecKAORGKM09GxXLmyvauQIAuVHQ1oZ3LmrFD89b27km0c4VAMiNQEkrAtDOFQDIjVBi/DKDdq4AQG7kp4stNUYxMwUA5Ea+etlyoyf/7gBAbvwX2xud+XcHAHIjXzHJlthIvpN/dwAgN9RlH081m392ACA31GWfv7E8jn92ACA31H91vjgAgNwoaItjtnF9qtlsbQAAuXF5Yjr37MwucQAgNwAADMbkBgAwGJMbAEBuXEOxIol1AACFK1Ek1iG5ob1UAEAYSCQ3AABFMDfKJibG/q/bVolFc5ORt8X74m3xtq7d+4pNTCzLvhzeFm+L98Xb4lMIPtq8Ld4Xb4u3RW6wrnhbvC/eFm+L98W64m3xP5a3xdvifbGueFu8Lf65eFu8L9YVb4u3xfvibfG+yA3eFm+L98Xb4lMIPtq8Ld4Xb4u3RW4AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAIqRSvEJtSxlbOJ9ISuHiW7duobnfrWE+DK+JfbKOWq2dvO0n4xskFDDt9heOe1fK+E+l69IfMClQlTO+Z/V2qu0+UjlJo0r+JbbK0cp06iF5Z8k4l+J0SpE5SgR9VtYx7taifVLhqqKt+gZu/Xu6bufdtd9M7RqT5+gleM+2+ki7vdVJzVHJHthmWCVkwz1tLuvZ5Qlpx4VydvYQAWpHGaw/jnMmmWOo23+1IqDQ9yL7JWDJIqX8c0sKuW8yNmPq5hL7ZWj1NZfumS9GmmWIw5r1Y4WKljlJKVnHNNeesZAd9l4s1ZlTnIFq4q7p0VO/rnngsgU839vnhzbcVzyBgapHOfDvGx3btQ6IDl/7hNJdQVWDsuNzC0G8wvdUm083XlBjtYPUjmK6wuRUzvTLkkjvWp3UbL/+FtkogqsnKS++W+15bDsMDY//pS8NO27zLaygZWjVNuvf9KOi/xslNNEjvxxWs62DFI5SeR6LQv1tzXCDIpsObvzL5GPglTFXptJsdqfNVbKhZr6B+K4fFxalVkmR6sEVI7TMmfxLndupEiG9gW9zwUZHFg5LDe+tFRtJe8ZbZPxD9kQWDnLRMmdoo2f1abqn8bS+yW1koqcJefvC6gcaac5Es2UvxOUanXC/I5mrxxlupzprlSpGSL6VkWDHJkaoaqul90lAypHGS6XtI3aiIlyzvikbZPNNZRrguQlBlZwb6GdEn026t+SoW96ltpvfNDtldOU3bOvvDs3oi5JG/O/6vaAysm5sUaWGV9rc6VBQOUoVU7JdF/VT45X1m+3y8yAyokS5Kz+xSvyqPlFZYQcdPlXzrLe/a17t7ENuNj8ohJ1RjoEVI6yURYbt6tkodInGS/W1avVkhJQwSNDemh/rnP/950pawMqp1mU10a5c2OQ7Dceqid5Uf6V43IjuoVn71zEaTMA1RaZ5F85yzA5aZmt+UK+UOaQujOgcqKl5ljTSs4Ye8crX5LG/pWzpLo/YRuML5QZ7u32lbIooHKU/TLQvfGbZmxUrTeqPpIVUMGtce5ZfZ7qkJEe2v5wY2i1Vw7TKlf78uDOjddllfngcX1otVdOyw19h2T2kmp68YBIeTMi9S0Ne+UsKbJSlYpPqGZW22S4+3v6uQj/yoHKnTQ/ZyNkq/lAmj602ivlsFnFTfo2UvTJ3Hjt7eWKuUNtkr6lYa+cJd09rzJDLpVS6leZYVQ1RWr6V9C+p7Zu3Wd6Vp6xynLc3320/6MBlcP+s+7NKO/NjWWe7z5pMsi/clpuXMpIO6Md86F/eNvIMeX+P7rev3KWbTJlsHYUTu7GeL066P6+ov0fjfKvHGiIZLjM4cj9fWWd/nXdXjlL2S2yfnCHEfuMmcNYkXLK3KpP86+c9u3lD31evpK211877GSr+/uKuiBN/Ctoq0k/TPBAW+PzIHK/8Vh9kUi/ymEW5+rf8dy5sdL9ZUEbn4b5V87SdrA2n1NyRLb8qhU95KByT+Fs9q+c9k1vrZxP258rp/RvKkfd24EVRO7zrxxos0w1bj/ybAf+qs//2iuHKTM9Rz8Ot59+v5HkKfcUzn7/ylla5Mofz3QfsSdb9G2mXZ7twBPSyr+C9s18y5b0M3J0gHa3pLg3MhtJrsuvcpbWucY2hTs3fpZZ5sN/yBD/ypH6Sa52fl8HOWyWE2STf+UsadpB0drGROO98rtWZUo749Gqxvc+e+U89SXHfN0LPbtUU+VV/8pZSi+5lLdfOzT65FCtaCAS4f5QZvhXDjP8ov4d+uxQ4xuKd3Q4rR81Zq9gKj8pN1dfHWfEPOq6jZwIqJykZMbJ7vpZuhkyvbWWfUtkqXLv+errXzmS64S014/xOOdyHxeZ6l85yxY5Udn8pOmHKuyVAcaj2pxHBf/KeT70HFUyVdZ5NkAm+FfOskTS9AnFvqf1KcRoEfMaBc/oZ6nYK6epN3PdhsX1W0qe9kn7zX3CUJk8/euzvYLHl/Kx8f18qDKPb9kRUCmnzb25/WzsxzM3rS5JvH/lTPv1w6aremZufpYP/StnWeaeWiuVIwn6jL/5DbyHsQVlrxynVKb0M+8N8HwDPyLd/StHibwgie5vKFr2ubLdMzcLtYMb/CpnmiLpSj8Kbql7tuVSaf8KHlONb0VL5RdlzrkuCaiUs+beDKdlzxZt9rileQC9Nud6ItK/cqQaxviq9pg76yIPS/+ASjnrOFxzJ2oNY2p0mjtFFhkbTvbKcfpLpvtTVte9v7CJ5ET5V45SR6S2e5tit9JPe5jhPrtxYkDlRJF7jPM3hso+l3mYyY6ACu7PbKU041SednKhkf5ZviStAyoncu/fcB0wdmmU3GqcNWevHCXCPHE/8hf5u6RxhE5GBWNHeHaFgMpZOXg2r435BfaoNjteP8/YpVH3pPGvZ68cZ51v62+ruUtjjXm8m71ykjLnzEyIWG/MtQ2XrJrG0fqX6gZUzlLFPbF4Rt9wjzpn7NKoelAmB1RI2zBpYLuBMw/J2XrmmTx/T2g/Mcv9WbZXzs0NNVhkaff+v8uZ2MDKSSqdWzmhb4+J2lWNjE2Lalmyc3D3D3PMS1XYK2d5TU7N6NE/Jc987cske0r74QdkV0Rg5bTv5paz99vlyao+fVIlp0VgpZx1wGrOsoHth20R40D20nsl45n2086aZ1vbK2dJXza0/ZCNnutTvS4XXu8wZLccqhBYFXub3HsB/mprHiS/y6h2RgWpHJwbapbxRk71CVY5SLnj5r/WRXc0tMoyymWuIJWzNqRSzDe2xHjt5c1PZYaZ6/ZKOe1STlt8xb8v6W/kwrBglaO+l7tHjTzzwMT6B41qXZkglaPsMF75WfdxCiVXGeWRhCAVWs5cuWnd0kHl3GXpYas2rXymVNDKeZZtae++12rphnUzvFvO9spBSvddlLpp1RTvS6/56tqNX3YPXjlLu483rF/s+eodMejn31ZPKB+0cthX8y3WnU2NF61f/1GD4JWTuHosWf976uuNPDkyac1vKd73aa+cpNTgLzeumVrHW/dZtjF1SlTwCgAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAACAKyhq1a81gy4olbI6ntUDAMVOv9V+FtqXL5W1xu3E1as9+VErZfXqvkotkt9YfQBQ7EwVP3/YFsfnitnI+WcRd1/x2AyRRVr78ZpnpA/rDwCKm/oDDFp8rDXvtbMt/lXWKXtuNDos8qpxb7HsZP0BQDHVUuSjIA/fnyv97bnR6rjkDjcfaiLSnjUHAOSGxYdyoowtN7qfkfP9PUvTZA1rDgDIDR/XYVmmrLkx6KKcam/ZN3KhGqsOAMgNr0SRodbcGJErWYm+xa1FnmbVAQC5YT3YqoElN6aJHGxgWVzuonzBqgMAcsNrjWS7fLnxq0h6XdvynbKLVQcA5Iby7fjeq3y5IZJVw758rZxh1QEAueF1TLZZcuO0dr6ffXmKSFXWHQCQGx7nPWf9GbnR44jIEtvyRSLRrDsAIDc8zshGS27E1tdOFV/qsixfIlKLdQcA5IbHYd+VRIzjcOsdEvnCEhwrRcqx7gCA3PDYJgdsuaHuPyCyLMK7fL1kseoAgNzw+lIulbflhrpvv0hKSc/yA7KJVQcA5IbXMJHW9txQdfeJrHQHRy2RWaw6ACA3vOqKTPLLDRW9V2RVpHG3r0hbVh0AkBs+O2S9f26o6HSRNaWUcThVZiSrDgDIDZ8RkhPtnxuq1m6R1NJah/GjQX8GAFAMRA0Y0DjIw+WPyWTzXuKAAZ5d5KqG1hYwVp+myollzQEArKZLRslQy9bJz6wgAIBN5SwZokJNbV18gBUEAFB+h+JmhNj3vV4+ZPUAAPy4+g0I3gs2sv+ASqweAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAFKL/DzZoem+ZHxwSAAAAAElFTkSuQmCC", "text/html": [ "
\"\"\"\"
" ], "text/plain": [ "<_NotebookFigure size 1200x360 with 1 Axes>" ] }, "execution_count": 9, "metadata": { "image/png": { "width": 795 } }, "output_type": "execute_result" } ], "source": [ "series.plot(\"phases.0.cell.a\")" ] }, { "cell_type": "markdown", "id": "05-19", "metadata": {}, "source": [ "## Is the answer path-dependent?\n", "\n", "A warm start means each answer depends on the one before it.\n", "If a fit lands in a slightly wrong place, the next fit starts there.\n", "`direction=\"both\"` runs the chain forward and backward and compares the two.\n", "A parameter the two directions disagree on is flagged `SEQUENTIAL_PATH_DEPENDENT`.\n", "That check is what separates a measured trend from an artefact of the order." ] }, { "cell_type": "code", "execution_count": 10, "id": "05-20", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "['SEQUENTIAL_PERSISTENT_FINDING', 'SEQUENTIAL_PERSISTENT_FINDING']\n" ] } ], "source": [ "both = rx.refine_sequential(patterns, *starting_model(), plan=\"mccusker_default\",\n", " x=TEMPERATURES, x_label=\"T (K)\", labels=LABELS, direction=\"both\")\n", "print([d.code for d in both.diagnostics])" ] }, { "cell_type": "markdown", "id": "05-21", "metadata": {}, "source": [ "No `SEQUENTIAL_PATH_DEPENDENT`: both directions agree, so the trend is not an artefact of fitting from cold to hot.\n", "\n", "## A jump in the series\n", "\n", "Now give the 700 K pattern a cell 2 % larger than its neighbours, as a phase transition might." ] }, { "cell_type": "code", "execution_count": 11, "id": "05-22", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[info] SEQUENTIAL_DISCONTINUITY: phases.0.cell.a steps by 0.08536 between 600 K and 700 K (T (K) 600 → 700), 41× its median step over the series\n", "[4.15660247525093, 4.15868088118438, 4.160757148373907, 4.16284114946175, 4.248201709983883, 4.166991214247255, 4.169067860703045]\n" ] } ], "source": [ "jumped = list(patterns)\n", "jumped[4] = simulate(TRUE_A[4] * 1.02, seed=15)\n", "series_jump = rx.refine_sequential(jumped, *starting_model(), plan=\"mccusker_default\",\n", " x=TEMPERATURES, x_label=\"T (K)\", labels=LABELS)\n", "for d in series_jump.diagnostics:\n", " if d.code == \"SEQUENTIAL_DISCONTINUITY\" and \"cell.a\" in d.message:\n", " print(f\"[{d.level}] {d.code}: {d.message}\")\n", "print(series_jump.trajectory(\"phases.0.cell.a\").value)" ] }, { "cell_type": "markdown", "id": "05-23", "metadata": {}, "source": [ "The 700 K fit found the larger cell.\n", "When a warm start fails, the series escalates: it retries in stages, then from the initial model, and keeps the best attempt.\n", "`SEQUENTIAL_DISCONTINUITY` marks the step.\n", "It cannot tell a real change from a failed fit, so check that pattern's own fit before reading the jump as physics.\n", "\n", "## A frame the chain cannot fit\n", "\n", "Make the jump 5 %, beyond what any starting point here can reach." ] }, { "cell_type": "code", "execution_count": 12, "id": "05-24", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "SeriesResult: 7 pattern(s), rietveld, direction=forward\n", " diagnostics: 3 unresolved\n", " WARNING SEQUENTIAL_RWP_OUTLIER: pattern 4 (700 K) reached Rwp 0.9692 on its best rung, 19.02× the fence at 0.0510 (the reseed factor × the median Rwp of the patterns accepted before it), after every rung the chain tried (warm, warm_staged, cold); against 600 K, the last pattern before it inside the fence: GoF 26.32 to its 1.02, Rexp 0.0368 to its 0.0400 — no starting point brought this fit near its neighbours', so the pattern itself differs from them; open its own fit before reading its values. GoF says which way: well above theirs is a model the pattern has outgrown — a specimen change the model lacks, or a bad frame; like theirs, with a higher Rexp, is fewer counts — a short or blank frame, a weaker beam — where the fit may be sound, so check that each phase is still in the pattern (a scale within a few esds of zero is one it does not contain). The fence follows the median of the accepted patterns, so a run of these that stops is a lasting change the median caught up with, not one that went away. The steps into and out of this pattern are left out of the discontinuity scan\n", " INFO SEQUENTIAL_PERSISTENT_FINDING: CAPILLARY_OFFSET_UNAVAILABLE fired on instrument.geometry.goniometer_radius_mm in 7 of 7 patterns (100% of the series) — read this as a property of the model rather than of any one pattern: CAPILLARY_OFFSET_UNAVAILABLE in most of a series is the same finding repeated, and the thing to change is the model or the plan, not the individual fits. The per-pattern diagnostics carry each occurrence\n", " WARNING SEQUENTIAL_PERSISTENT_FINDING: PATTERN_UNDERSAMPLED fired on the fit in 7 of 7 patterns (100% of the series) — read this as a property of the model rather than of any one pattern: PATTERN_UNDERSAMPLED in most of a series is the same finding repeated, and the thing to change is the model or the plan, not the individual fits. The per-pattern diagnostics carry each occurrence\n", " trajectory (T (K)):\n", " [1/7] 300 K x=300 converged Rwp=0.0408\n", " [2/7] 400 K x=400 converged Rwp=0.0411\n", " [3/7] 500 K x=500 converged Rwp=0.0401\n", " [4/7] 600 K x=600 converged Rwp=0.0408\n", " [5/7] 700 K x=700 converged Rwp=0.9692\n", " [6/7] 800 K x=800 converged Rwp=0.0401\n", " [7/7] 900 K x=900 converged Rwp=0.0405\n" ] } ], "source": [ "broken = list(patterns)\n", "broken[4] = simulate(TRUE_A[4] * 1.05, seed=15)\n", "series_broken = rx.refine_sequential(broken, *starting_model(), plan=\"mccusker_default\",\n", " x=TEMPERATURES, x_label=\"T (K)\", labels=LABELS)\n", "print(series_broken)" ] }, { "cell_type": "markdown", "id": "05-25", "metadata": {}, "source": [ "The 700 K pattern failed on every attempt, and `SEQUENTIAL_RWP_OUTLIER` names it.\n", "The 800 K pattern started from that failure, so its first warm start failed too.\n", "The chain escalated one rung, to a staged warm start, and that answer is good: the entry's `rung` reads `warm_staged`.\n", "\n", "A pattern whose fit diverges on every attempt goes further: it is quarantined (`SEQUENTIAL_UNRECOVERED`), seeds no successor and is left out of the series' statistics.\n", "Either way the rest of the trajectory survives, and the summary says which points to distrust.\n", "\n", "## Checking an agent's work\n", "\n", "When an agent reports a trajectory, check these.\n", "The rules are in the rietx agent skill (`SKILL.md` §4b, the trajectory row, and `references/series.md`).\n", "\n", "- **Was the chain run both ways?** Without `direction=\"both\"`, a trend can be an artefact of the fitting order.\n", "- **Every series diagnostic.** A discontinuity, an outlier or a reseed names the patterns to look at one by one.\n", "- **A finding in most patterns is one finding** about the model, and the thing to change is the model or the plan.\n", "- **Each point's esd, and the precision behind the trend.** A trajectory is only as good as its worst point." ] } ], "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" }, "rietx": { "version": "1.10.0.dev0" } }, "nbformat": 4, "nbformat_minor": 5 }