Model & Analysis Replication Package/OIC-OO v7.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Important:** The model estimation code is intended to work with an experimental parallelised Vensim engine. With appropriate modifications to the main function calls (but not the analytical procedure), the same analysis can be run on regular commercially available Vensim DSS, though it will take *much* longer. Please contact [Tom Fiddaman](mailto:tom@ventanasystems.com) for information on the experimental Vensim engine.\n",
    "\n",
    "For more information on the model estimation procedure, see S4 of the Supplementary Materials of the paper.\n",
    "\n",
    "**Note:** if running in Jupyter, the `keyboard` module may need to be installed directly in the Notebook; see [here](https://stackoverflow.com/questions/38368318/installing-a-pip-package-from-within-a-jupyter-notebook-not-working) for example. The keyboard module is *only* used to bypass Vengine error messages; if not running Vengine (e.g. using normal Vensim DSS), it is not necessary, and you can safely remove the import statement and all press commands in the code.\n",
    "\n",
    "---\n",
    "\n",
    "This notebook contains code cells for the following:\n",
    "1. Class & function definitions\n",
    "2. Controlfile input & parsing\n",
    "3. Main iterative estimation process & MCMC\n",
    "4. Baseline projections & sensitivity analysis of input projection assumptions\n",
    "5. Policy scenarios w/ uncertainty quantification\n",
    "6. Loop knockout (structural sensitivity) analyses\n",
    "7. Parametric assumption elasticity analyses\n",
    "8. Alternative data condition analyses\n",
    "9. Holdout data (out-of-sample validation) analysis\n",
    "10. Synthetic data generation & analysis\n",
    "11. PSRF statistic parsing\n",
    "\n",
    "The first three of these (definitions, controlfile input, main estimation) should always be run first and in order (using controlfile bypass settings to skip main estimation if not needed). Remaining analyses can be run as needed, and do not need to be in order."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import subprocess\n",
    "import regex\n",
    "import json\n",
    "import time\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "from keyboard import press\n",
    "from shutil import copy\n",
    "from distutils.dir_util import copy_tree\n",
    "from statsmodels.tsa.tsatools import detrend"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class Script(object):\n",
    "    \"\"\"Master object for holding and modifying .cmd script settings, \n",
    "    creating .cmd files, and running them through Vensim/Vengine\"\"\"\n",
    "    def __init__(self, controlfile):\n",
    "        print(\"Initialising\", self)\n",
    "        for k, v in controlfile['simsettings'].items():\n",
    "            self.__setattr__(k, v if isinstance(v, str) else v.copy())\n",
    "        self.setvals = []\n",
    "        self.runcmd = \"SIMULATE>REPORT|1\\nMENU>RUN_OPTIMIZE|o\\n\"\n",
    "        self.savecmd = f\"MENU>VDF2TAB|!|!|{self.savelist}|\\n\"\n",
    "        self.basename = controlfile['baserunname']\n",
    "        self.cmdtext = []\n",
    "        \n",
    "    def copy_model_files(self, dirname):\n",
    "        \"\"\"Create subdirectory and copy relevant model files to it,\n",
    "        then change working directory to subdirectory\"\"\"\n",
    "        os.makedirs(dirname, exist_ok=True)\n",
    "        os.chdir(f\"./{dirname}\")\n",
    "\n",
    "        # Copy needed files from the working directory into the sub-directory\n",
    "        for s in ['model', 'payoff', 'optparm', 'sensitivity', 'savelist', 'senssavelist']:\n",
    "            if getattr(self, s):\n",
    "                copy(f\"../{getattr(self, s)}\", \"./\")\n",
    "        for slist in ['data', 'changes']:\n",
    "            for file in getattr(self, slist):\n",
    "                copy(f\"../{file}\", \"./\")\n",
    "        \n",
    "    def add_suffixes(self, settingsfxs):\n",
    "        \"\"\"Modify mdl, voc, etc. with suffixes specified as dict\"\"\"\n",
    "        for s, sfx in settingsfxs.items():\n",
    "            if hasattr(self, s):\n",
    "                self.__setattr__(s, getattr(self, s)[:-4] + sfx + getattr(self, s)[-4:])\n",
    "   \n",
    "    def update_changes(self, chglist, setvals=[]):\n",
    "        \"\"\"Combines and flattens list of paired names & suffixes for \n",
    "        changes, and appends to changes setting; also updates setvals\"\"\"\n",
    "        flat = [i for s in \n",
    "                [[f\"{self.basename}_{n}_{sfx}.out\" for n in name] if isinstance(name, list) \n",
    "                 else [f\"{self.basename}_{name}_{sfx}.out\"] for name, sfx in chglist] for i in s]\n",
    "        self.changes.extend(flat)\n",
    "        self.setvals = setvals\n",
    "        \n",
    "    def write_script(self, scriptname):\n",
    "        \"\"\"Write actual .cmd file based on controlfile attributes\"\"\"\n",
    "        self.cmdtext.extend([\"SPECIAL>NOINTERACTION\\n\", \n",
    "                             f\"SPECIAL>LOADMODEL|{self.model}\\n\"])\n",
    "        \n",
    "        for s in ['payoff', 'sensitivity', 'optparm', 'savelist', 'senssavelist']:\n",
    "            if hasattr(self, s):\n",
    "                self.cmdtext.append(f\"SIMULATE>{s}|{getattr(self, s)}\\n\")\n",
    "        \n",
    "        if hasattr(self, 'data'):\n",
    "            datatext = ','.join(self.data)\n",
    "            self.cmdtext.append(f\"SIMULATE>DATA|\\\"{','.join(self.data)}\\\"\\n\")\n",
    "\n",
    "        if hasattr(self, 'changes'):\n",
    "            self.cmdtext.append(f\"SIMULATE>READCIN|{self.changes[0]}\\n\")\n",
    "            for file in self.changes[1:]:\n",
    "                self.cmdtext.append(f\"SIMULATE>ADDCIN|{file}\\n\")\n",
    "        \n",
    "        self.cmdtext.extend([\"\\n\", f\"SIMULATE>RUNNAME|{scriptname}\\n\"])\n",
    "        \n",
    "        if hasattr(self, 'setvals'):\n",
    "            for var, val in self.setvals:\n",
    "                self.cmdtext.append(f\"SIMULATE>SETVAL|{var}={val}\\n\")\n",
    "        \n",
    "        self.cmdtext.extend([self.runcmd, self.savecmd, \n",
    "                             \"SPECIAL>CLEARRUNS\\n\", \"MENU>EXIT\\n\"])\n",
    "        \n",
    "        with open(f\"{scriptname}.cmd\", 'w') as scriptfile:\n",
    "            scriptfile.writelines(self.cmdtext)\n",
    "\n",
    "    def run_script(self, scriptname, controlfile, subdir, logfile):\n",
    "        \"\"\"Run the compiled .cmd file, calling Vengine by default\"\"\"\n",
    "        return run_vengine_script(scriptname, controlfile['vensimpath'], \n",
    "                                  controlfile['timelimit'], '.log', check_opt, logfile)\n",
    "        \n",
    "\n",
    "class SubdirScript(Script):\n",
    "    \"\"\"Script subclass for optimization runs in subdirectory\"\"\"\n",
    "    def prep_subdir(self, scriptname, controlfile, subdir):\n",
    "        self.copy_model_files(subdir)\n",
    "        copy(f\"../{scriptname}.cmd\", \"./\")\n",
    "\n",
    "    def run_script(self, scriptname, controlfile, subdir, logfile):\n",
    "        \"\"\"Copy necessary files to subdirectory and run there, copying \n",
    "        .out file back to parent directory when complete\"\"\"\n",
    "        self.prep_subdir(scriptname, controlfile, subdir)\n",
    "        payoff = run_vengine_script(scriptname, controlfile['vensimpath'], \n",
    "                                    controlfile['timelimit'], '.log', check_opt, logfile)\n",
    "        copy(f\"./{scriptname}.out\", \"..\") # Copy the .out file to parent directory\n",
    "        os.chdir(\"..\")\n",
    "        return payoff\n",
    "\n",
    "\n",
    "class MCScript(Script):\n",
    "    \"\"\"Script subclass for MCMC optimizations\"\"\"\n",
    "    def run_script(self, scriptname, controlfile, subdir, logfile):\n",
    "        run_vengine_script(scriptname, controlfile['vensimpath'], \n",
    "                           controlfile['timelimit'], '_MCMC_points.tab', check_MC, logfile)\n",
    "        \n",
    "        downsample(scriptname, controlfile['samplefrac'], remove=True) # Create MCMC subsample\n",
    "\n",
    "        \n",
    "class LongScript(Script):\n",
    "    \"\"\"Script subclass for long calibration runs e.g. all-params\"\"\"\n",
    "    def run_script(self, scriptname, controlfile, subdir, logfile):\n",
    "        return run_vengine_script(scriptname, controlfile['vensimpath'], \n",
    "                                  controlfile['timelimit']*5, '.log', check_opt, logfile)\n",
    "\n",
    "    \n",
    "class RunScript(Script):\n",
    "    \"\"\"Script subclass for simple single runs (not optimzations)\"\"\"\n",
    "    def __init__(self, controlfile):\n",
    "        super().__init__(controlfile)\n",
    "        self.runcmd = \"MENU>RUN|o\\n\" # Change .cmd run command\n",
    "\n",
    "    def run_script(self, scriptname, controlfile, subdir, logfile):\n",
    "        \"\"\"Run the compiled .cmd file with Vensim instead of Vengine\"\"\"\n",
    "        return run_vensim_script(scriptname, controlfile['vensim7path'], 10, logfile)\n",
    "            \n",
    "\n",
    "class ScenScript(Script):\n",
    "    \"\"\"Script subclass for scenario analysis with .cin files\"\"\"\n",
    "    def update_changes(self, chglist, setvals=None):\n",
    "        scen = []\n",
    "        try: # Pull out all str items from end of chglist\n",
    "            while type(chglist[-1]) == str:\n",
    "                scen.append(chglist.pop())\n",
    "        except IndexError: pass\n",
    "        super().update_changes(chglist, setvals)\n",
    "        scen.reverse() # Remember to reverse to maintain scen order!\n",
    "        self.changes.extend(scen) # Then add str items to self.changes\n",
    "        chglist.extend(scen) # And add them back to chglist as well\n",
    "        \n",
    "    def run_script(self, scriptname, controlfile, subdir, logfile):\n",
    "        return run_vengine_script(scriptname, controlfile['vensimpath'], \n",
    "                                  controlfile['timelimit']/5, '.vdf', check_run, logfile)\n",
    "    \n",
    "\n",
    "class ScenRunScript(ScenScript):\n",
    "    \"\"\"Script subclass for scenario analysis runs (not optimizations)\"\"\"\n",
    "    def __init__(self, controlfile):\n",
    "        super().__init__(controlfile)\n",
    "        self.runcmd = \"MENU>RUN|o\\n\"\n",
    "\n",
    "\n",
    "class ScenSensScript(ScenScript):\n",
    "    \"\"\"Script subclass for scenario sensitivity analysis\"\"\"\n",
    "    def __init__(self, controlfile):\n",
    "        super().__init__(controlfile)\n",
    "        self.sensitivity = self.basename + '_full.vsc'\n",
    "        self.runcmd = \"MENU>RUN_SENSITIVITY|o\\n\"\n",
    "        self.savecmd = f\"MENU>SENS2FILE|!|!|%#T\\n\" # Exports only percentiles\n",
    "\n",
    "        \n",
    "class FullSensScript(ScenScript):\n",
    "    \"\"\"Script subclass for scenario sensitivity w/ full simulation export\"\"\"\n",
    "    def __init__(self, controlfile):\n",
    "        super().__init__(controlfile)\n",
    "        self.sensitivity = self.basename + '_full.vsc'\n",
    "        self.runcmd = \"MENU>RUN_SENSITIVITY|o\\n\"\n",
    "        self.savecmd = f\"MENU>SENS2FILE|!|!|#T[\\n\" # Exports full sample (VERY LARGE!)\n",
    "        \n",
    "    def run_script(self, scriptname, controlfile, subdir, logfile):\n",
    "        \"\"\"Run with extended timelimit to allow time for vdf export\"\"\"\n",
    "        return run_vengine_script(scriptname, controlfile['vensimpath'], \n",
    "                                  controlfile['timelimit']*5, '.log', check_run, logfile)\n",
    "\n",
    "\n",
    "def compile_script(controlfile, scriptclass, name, namesfx, settingsfxs, \n",
    "                   logfile, chglist=[], setvals=[], subdir=None):\n",
    "    \"\"\"Master function for assembling & running .cmd script\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    controlfile : JSON object\n",
    "        Master control file specifying sim settings, runname, etc.\n",
    "    scriptclass : Script object\n",
    "        Type of script object to instantiate, depending on run type\n",
    "    name : str\n",
    "    namesfx : str\n",
    "        Along with `name`, specifies name added to baserunname for run\n",
    "    settingsfxs : dict of str\n",
    "        Dict of suffixes to append to filenames in simsettings; use to \n",
    "        distinguish versions of e.g. .mdl, .voc, .vpd etc. files\n",
    "    logfile : str of filename/path\n",
    "    chglist : list of tuples of (str or list, str)\n",
    "        Specifies changes files to be used in script; specify as tuples \n",
    "        corresponding to `name`, `namesfx` of previous run .out to use; \n",
    "        tuples can also take a list of `names` as first element, taking \n",
    "        each with the same second element; if used with ScenScript run, \n",
    "        `chglist` can also take non-tuple strs as final elements, which \n",
    "        will be added directly (e.g. .cin files for scenarios)\n",
    "    setvals : list of tuples of (str, int or float, <str>)\n",
    "        Specifies variables and values to change for a given run using \n",
    "        Vensim's SETVAL script command; by default all SETVAL commands \n",
    "        will be implemented together for main run, but if `scriptclass` \n",
    "        is MultiScript, each SETVAL command will be implemented and run \n",
    "        separately in sequence; if used with MultiScript, each tuple in \n",
    "        `setvals` will require a third str element specifying the suffix \n",
    "        with which to save the run\n",
    "    subdir : str, optional\n",
    "        Name of subdirectory to create/use for run, if applicable\n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "    float\n",
    "        Payoff value of the script run, if applicable, else 0\n",
    "    \"\"\"\n",
    "    mainscript = scriptclass(controlfile)\n",
    "    mainscript.add_suffixes(settingsfxs)\n",
    "    mainscript.update_changes(chglist, setvals)\n",
    "    scriptname = f\"{mainscript.basename}_{name}_{namesfx}\"    \n",
    "    mainscript.write_script(scriptname)\n",
    "    return mainscript.run_script(scriptname, controlfile, subdir, logfile)\n",
    "\n",
    "\n",
    "def write_log(string, logfile):\n",
    "    \"\"\"Writes printed script output to a logfile\"\"\"\n",
    "    with open(logfile,'a') as f:\n",
    "        f.write(string + \"\\n\")\n",
    "    print(string)\n",
    "    \n",
    "\n",
    "def check_opt(scriptname, logfile):\n",
    "    \"\"\"Check function for use with run_vengine_script for optimizations\"\"\"\n",
    "    if check_zeroes(scriptname):\n",
    "        write_log(f\"Help! {scriptname} is being repressed!\", logfile)\n",
    "    return not check_zeroes(scriptname)\n",
    "\n",
    "def check_MC(scriptname, logfile, threshold=0.1):\n",
    "    \"\"\"Check function for use with run_vengine_script for MCMC\"\"\"\n",
    "     # Check for discrepancy between .out and .rep payoff\n",
    "    if abs(compare_payoff(scriptname)) >= threshold:\n",
    "        write_log(f\"{scriptname} is a self-perpetuating autocracy! re-running MC...\", logfile)\n",
    "        return False\n",
    "    return True\n",
    "\n",
    "def check_run(scriptname, logfile):\n",
    "    \"\"\"Check function for use with run_vengine_script for normal & sens runs\"\"\"\n",
    "    if not os.path.exists(f\"./{scriptname}.vdf\"): # Check that output VDF exists\n",
    "        write_log(f\"Help! {scriptname} is being repressed!\", logfile)\n",
    "    return os.path.exists(f\"./{scriptname}.vdf\")\n",
    "\n",
    "\n",
    "def run_vengine_script(scriptname, vensimpath, timelimit, checkfile, check_func, logfile):\n",
    "    \"\"\"Call Vengine with command script using subprocess; monitor output \n",
    "    file for changes to see if Vengine has stalled out, and restart if \n",
    "    it does, or otherwise bugs out; return payoff if applicable\"\"\"\n",
    "\n",
    "    write_log(f\"Initialising {scriptname}!\", logfile)\n",
    "\n",
    "    while True:\n",
    "        proc = subprocess.Popen(f\"{vensimpath} \\\"./{scriptname}.cmd\\\"\")\n",
    "        time.sleep(2)\n",
    "        press('enter') # Necessary to bypass the popup message in Vengine\n",
    "        while True:\n",
    "            try: # See if run completes within timelimit\n",
    "                proc.wait(timeout=timelimit)\n",
    "                break\n",
    "            except subprocess.TimeoutExpired: # If timelimit reached, check run status\n",
    "                try: # If Vengine gives error popup on exit, attempt to bypass it\n",
    "                    print(\"Attempting bypass...\")\n",
    "                    press('enter')\n",
    "                    proc.wait(3) # Process should complete if bypass successful\n",
    "                    break\n",
    "                except subprocess.TimeoutExpired:\n",
    "                    try: # If bypass unsuccessful, check if run still going\n",
    "                        write_log(f\"Checking for {scriptname}{checkfile}...\", logfile)\n",
    "                        timelag = time.time() - os.path.getmtime(f\"./{scriptname}{checkfile}\")                        \n",
    "                        if timelag < (timelimit): # Compare time since last output with timelimit\n",
    "                            write_log(f\"At {time.ctime()}, {round(timelag,3)}s since last output, \"\n",
    "                                      \"continuing...\", logfile)\n",
    "                            continue\n",
    "                        else: # If run seems to have stalled out, kill and restart\n",
    "                            proc.kill()\n",
    "                            write_log(f\"At {time.ctime()}, {round(timelag,3)}s since last output. \"\n",
    "                                      \"Calibration timed out!\", logfile)\n",
    "                            break\n",
    "                    except FileNotFoundError: # If check fails, kill and restart\n",
    "                        proc.kill()\n",
    "                        write_log(\"Calibration timed out!\", logfile)\n",
    "                        break\n",
    "        # Check if process successfully completed or bugged out / was killed\n",
    "        if proc.returncode != 1: # Note that Vengine returns 1 on MENU>EXIT, not 0!\n",
    "            write_log(f\"Return code is {proc.returncode}\", logfile)\n",
    "            write_log(\"Vensim! Trying again...\", logfile)\n",
    "            continue\n",
    "        try: # If process completed successfully, run final check for errors in output\n",
    "            if check_func(scriptname, logfile):\n",
    "                break # NOTE: this is the only successful completion outcome!\n",
    "        except FileNotFoundError: # Catch output error and restart run\n",
    "            write_log(\"Outfile not found! That's it, I'm dead.\", logfile)\n",
    "            pass\n",
    "    \n",
    "    time.sleep(2)\n",
    "\n",
    "    if os.path.exists(f\"./{scriptname}.out\"):\n",
    "        payoffvalue = read_payoff(f\"{scriptname}.out\")\n",
    "        write_log(f\"Payoff for {scriptname} is {payoffvalue}, calibration complete!\", logfile)\n",
    "        return payoffvalue # For optimisation runs, return payoff\n",
    "    return 0 # Set default payoff value for simtypes that don't generate one\n",
    "\n",
    "\n",
    "def run_vensim_script(scriptname, vensim7path, maxattempts, logfile):\n",
    "    \"\"\"Call Vensim (not Vengine) with command script using subprocess; \n",
    "    for instances when Vengine not necessary; try up to `maxattempts` \n",
    "    times; return payoff if applicable\"\"\"\n",
    "    attempts = 0\n",
    "    while attempts < maxattempts:\n",
    "        attempts += 1 # Track & update number of attempts to prevent infinite loop\n",
    "        if os.path.exists(f\"./{scriptname}.tab\"):\n",
    "            os.remove(f\"./{scriptname}.tab\") # Delete old output tabfile if needed\n",
    "        try:\n",
    "            subprocess.run(f\"{vensim7path} \\\"./{scriptname}.cmd\\\"\", check=True)\n",
    "            pass\n",
    "        except subprocess.CalledProcessError:\n",
    "            print(\"Vensim! Trying again...\")\n",
    "            continue\n",
    "        if os.path.exists(f\"./{scriptname}.tab\"): # Check for output tabfile\n",
    "            break\n",
    "        else:\n",
    "            write_log(f\"Help! {scriptname} is being repressed!\", logfile)\n",
    "            continue\n",
    "    \n",
    "    if os.path.exists(f\"./{scriptname}.out\"):\n",
    "        payoffvalue = read_payoff(f\"{scriptname}.out\")\n",
    "        write_log(f\"Payoff for {scriptname} is {payoffvalue}, calibration complete!\", logfile)\n",
    "        return payoffvalue # For optimisation runs, return payoff\n",
    "    return 0 # Set default payoff value for simtypes that don't generate one\n",
    "\n",
    "\n",
    "def create_syn_mdl(modelname, newmodelname):\n",
    "    \"\"\"Open .mdl as text and modify as needed for syndata analysis\"\"\"\n",
    "    with open(modelname,'r') as f:\n",
    "        filedata = f.read()\n",
    "        \n",
    "    synvarregex = regex.compile(r\"SynVar\\[Elm\\]=[\\s\\S]*?(\\n\\t~)\")\n",
    "    datavarregex = regex.compile(r\"MaxDataTime,\\s*DataVarBase\\[Elm\\]\")\n",
    "    \n",
    "    # Switch SynVar to read data from file (from previously generated synthetic dataset)\n",
    "    newdata = synvarregex.sub(\"SynVar[Elm]\\n\\t~\", filedata)\n",
    "    \n",
    "    # Switch DataVar to read data from SynVar\n",
    "    # DataVar default equation is `IF THEN ELSE(Time <= MaxDataTime, DataVarBase[Elm], :NA:)`\n",
    "    newdata = datavarregex.sub(\"MaxDataTime, SynVar[Elm]\", newdata)\n",
    "\n",
    "    with open(newmodelname,'w') as f:\n",
    "        f.write(newdata) # Save new .mdl for synthetic data calibration\n",
    "\n",
    "                       \n",
    "def split_voc(vocname, fractolfactor, linekey, mcsettings):\n",
    "    \"\"\"Splits .VOC file into multiple versions, for main, overdose, \n",
    "    initial values, and full model calibration\"\"\"\n",
    "    with open(vocname,'r') as f0:\n",
    "        filedata = f0.readlines()\n",
    "    \n",
    "    ### LINEKEY should include all OD-related parameter keywords e.g. fentanyl, overdose, OD\n",
    "    \n",
    "    vocmain = [line for line in filedata if line[0] == ':' or \n",
    "               ('Initial' not in line and not any(k in line for k in linekey))]\n",
    "    vocod = [line for line in filedata if line[0] == ':' or any(k in line for k in linekey)]\n",
    "    vocinit = [line for line in filedata if line[0] == ':' or 'Initial' in line]\n",
    "    vocfull = filedata.copy()\n",
    "    vocsens = vocfull.copy()\n",
    "    vocmc = ''.join(vocfull)\n",
    "    \n",
    "    # Identify and multiply fractional tolerance by fractolfactor for initial & sens runs\n",
    "    for l, line in enumerate(vocsens):\n",
    "        if ':FRACTIONAL_TOLERANCE' in line:\n",
    "            fractol = float(line.split('=')[1])\n",
    "            vocsens[l] = f\":FRACTIONAL_TOLERANCE={min(fractol*fractolfactor,0.1)}\\n\"\n",
    "\n",
    "    # Make necessary substitutions for MCMC settings\n",
    "    for k,v in mcsettings.items():\n",
    "        vocmc = regex.sub(f\":{regex.escape(k)}=.*\", f\":{k}={v}\", vocmc)\n",
    "        \n",
    "    # Write various voc versions to separate .voc files\n",
    "    for fname, suffix in zip([vocmain, vocod, vocinit, vocfull, vocsens, vocmc], \n",
    "                             ['m', 'o', 'i', 'f', 's', 'mc']):\n",
    "        with open(f\"{vocname[:-4]}_{suffix}.voc\", 'w') as f:\n",
    "            f.writelines(fname)\n",
    "\n",
    "\n",
    "def split_lk_vocs(vocname, lkdict):\n",
    "    \"\"\"Splits .VOC file into multiple versions for loop knockout tests \n",
    "    based on `lkdict`; one version for each dict key, setting all vars \n",
    "    in list matching key to 1e-06; also create matching .cin file\"\"\"\n",
    "    for key, varlist in lkdict.items():\n",
    "        with open(vocname,'r') as f0:\n",
    "            voclk = f0.read()\n",
    "\n",
    "        for var in varlist:\n",
    "            # Select and replace entire line containing var name for each var\n",
    "            lineregex = r'.*(?<=[^\\w ]|\\n)\\s?' + regex.escape(var) + r'\\s?(?=[^\\w ]).*'\n",
    "            voclk = regex.sub(lineregex, f\"1e-06<={var}<=1e-06\", voclk) # Fix var at 1e-06\n",
    "\n",
    "        with open(f\"{vocname[:-4]}_lk_{key}.voc\", 'w') as f1:\n",
    "            f1.writelines(voclk)\n",
    "        \n",
    "        # Also create .cin file for use in normal runs\n",
    "        cinlines = [f'{var} = 1e-06\\n' for var in varlist]\n",
    "        with open(f\"LK_{key}.cin\", 'w') as f2:\n",
    "            f2.writelines(cinlines)\n",
    "            \n",
    "\n",
    "def split_vpd(vpdname):\n",
    "    \"\"\"Splits .vpd file into main and overdose PMC versions\"\"\"\n",
    "    with open(vpdname,'r') as f0:\n",
    "        filedata = f0.read()\n",
    "    \n",
    "    vpdod = filedata.replace('Elm', 'ODElm')\n",
    "        \n",
    "    # Write various vpd versions to separate .vpd files\n",
    "    for fname, suffix in zip([vpdod], ['o']):\n",
    "        with open(f\"{vpdname[:-4]}_{suffix}.vpd\", 'w') as f:\n",
    "            f.write(fname)\n",
    "\n",
    "\n",
    "def check_zeroes(scriptname):\n",
    "    \"\"\"Check if an .out file has any parameters set to zero (indicates \n",
    "    Vengine error), return True if any parameters zeroed OR if # runs = \n",
    "    # restarts, and False otherwise\"\"\"\n",
    "    filename = f\"{scriptname}.out\"\n",
    "    with open(filename,'r') as f0:\n",
    "        filedata = f0.readlines()\n",
    "    \n",
    "    checklist = []\n",
    "    for line in filedata:\n",
    "        if line[0] != ':': # Include only parameter lines\n",
    "            if ' = 0 ' in line:\n",
    "                checklist.append(True)\n",
    "            else:\n",
    "                checklist.append(False)\n",
    "        elif ':RESTART_MAX' in line:\n",
    "            restarts = regex.findall(r'\\d+', line)[0]\n",
    "    \n",
    "    # Ensure number of simulations != number of restarts\n",
    "    if f\"After {restarts} simulations\" in filedata[0]:\n",
    "        checklist.append(True)\n",
    "    \n",
    "    return any(checklist)\n",
    "\n",
    "\n",
    "def read_payoff(outfile, line=1):\n",
    "    \"\"\"Identifies payoff value from .OUT or .REP file - \n",
    "    use line 1 (default) for .OUT, or use line 0 for .REP\"\"\"\n",
    "    with open(outfile, 'r') as f:\n",
    "        payoffline = f.readlines()[line]\n",
    "    payoffvalue = [float(s) for s in \n",
    "                   regex.findall(r'-?(?:0|[1-9]\\d*)(?:\\.\\d*)?(?:[eE][+\\-]?\\d+)?', payoffline)][0]\n",
    "    return payoffvalue\n",
    "\n",
    "\n",
    "def compare_payoff(scriptname):\n",
    "    \"\"\"Returns the difference in payoffs between .OUT and .REP file, \n",
    "    which should be zero in most cases except when MCMC bugs out\"\"\"\n",
    "    difference = read_payoff(f\"{scriptname}.out\") - read_payoff(f\"{scriptname}.rep\", 0)\n",
    "    print(\".OUT and .REP payoff difference is\", difference)\n",
    "    return difference\n",
    "\n",
    "\n",
    "def increment_seed(vocfile, logfile):\n",
    "    \"\"\"Increments random number seed in a .VOC file by 1\"\"\"\n",
    "    with open(vocfile, 'r') as f:\n",
    "        vocdata = f.read()\n",
    "    seedregex = regex.compile(r':SEED=\\d+')\n",
    "    try:\n",
    "        i = int(regex.search(r'\\d+', regex.search(seedregex, vocdata).group()).group())\n",
    "        newdata = seedregex.sub(f\":SEED={i+1}\", vocdata)\n",
    "        with open(vocfile, 'w') as f:\n",
    "            f.write(newdata)\n",
    "    except:\n",
    "        write_log(\"No seed found, skipping incrementing.\", logfile)\n",
    "\n",
    "\n",
    "def get_value(file, varname):\n",
    "    \"\"\"General purpose function for reading values from .mdl, .out, etc. \n",
    "    files; returns value matching `varname` in a 'var = val' syntax\"\"\"\n",
    "    varregex = regex.compile(r'(?<=([^\\w ]|\\n)\\s?' + regex.escape(varname)\n",
    "                             + r'\\s*=)\\s*-?(?:\\d*)(\\.\\d*)?([eE][+\\-]?\\d+)?')\n",
    "\n",
    "    with open(file, 'r') as f:\n",
    "        filetext = f.read()\n",
    "        value = float((regex.search(varregex, filetext))[0])\n",
    "\n",
    "    return value\n",
    "\n",
    "\n",
    "def read_outvals(filename):\n",
    "    \"\"\"Read .out file and return matching lists of varnames and values; \n",
    "    NOTE need to include file extension in `filename`; also only works \n",
    "    for .out files with lower bounds specified\"\"\"\n",
    "    with open(filename, 'r') as f:\n",
    "        output = [line for line in f.readlines() if line[0] != ':']\n",
    "        \n",
    "        varnames = [line.split(' <= ')[1].split(' = ')[0].strip() for line in output]\n",
    "        basevals = [float(line.split(' <= ')[1].split(' = ')[1]) for line in output]\n",
    "        \n",
    "    return varnames, basevals\n",
    "\n",
    "\n",
    "def compile_sensoutvals(baserunname, sensdict, name, sfxs=['']):\n",
    "    \"\"\"Compile dataframe / tabfile of estimated parameter values read \n",
    "    from .out files for each of multiple sensitivity runs as specified \n",
    "    in `sensdict`, reading from runnames specified with `name` and any \n",
    "    `sfxs` (remember to include underscores if needed)\"\"\"\n",
    "    # Read base param values and set up dataframe\n",
    "    varnames, basevals = read_outvals(f\"{baserunname}_main_full.out\")\n",
    "    paramdf = pd.DataFrame({'Value': basevals}, index=varnames)\n",
    "\n",
    "    # Read iterator from sensdict\n",
    "    try:\n",
    "        keys = sensdict.keys()\n",
    "    except AttributeError:\n",
    "        keys = sensdict\n",
    "    \n",
    "    # Loop through scenarios / runs, with sfxs if needed\n",
    "    for key in keys:\n",
    "        for sfx in sfxs:\n",
    "            varnames, basevals = read_outvals(f\"{baserunname}_{name}_{key}{sfx}.out\")\n",
    "            sensdf = pd.DataFrame({f'{key}{sfx}': basevals}, index=varnames)\n",
    "            \n",
    "            # Add values from current scenario / run to dataframe\n",
    "            paramdf = paramdf.merge(sensdf, how='left', left_index=True, right_index=True)\n",
    "\n",
    "    paramdf.to_csv(f\"./{baserunname}_{name}_params.tab\", sep='\\t')\n",
    "        \n",
    "    \n",
    "def export_fits(baserunname, fitlist, base='main_full', name='main'):\n",
    "    \"\"\"Write .cmd file to cleanly export only fit-to-data variables \n",
    "    (SimVar, DataVar, NormErr) for elements specified in `fitlist`\"\"\"\n",
    "    savelisttext = [f'SimVar[{var[0]}]\\nDataVar[{var[0]}]\\nNormErr[{var[0]}]\\n'\n",
    "                    for var in fitlist]\n",
    "\n",
    "    with open(\"FitList.lst\", 'w') as f:\n",
    "        f.writelines(savelisttext)\n",
    "    exporttext = []\n",
    "    exporttext.append(\"SPECIAL>NOINTERACTION\\n\")\n",
    "    exporttext.append(f\"MENU>VDF2TAB|{baserunname}_{base}.vdf|\"\n",
    "                      f\"{baserunname}_{name}_fits_raw.tab|FitList.lst||0.25|\\n\")\n",
    "    exporttext.append(\"MENU>EXIT\\n\")\n",
    "\n",
    "    with open(\"ExportFits.cmd\", 'w') as f:\n",
    "        f.writelines(exporttext)\n",
    "\n",
    "\n",
    "def clean_fits(baserunname, name='main'):\n",
    "    \"\"\"Clean fit-to-data tabfile generated using `export_fits` .cmd\"\"\"\n",
    "    cols = pd.read_csv(f\"{baserunname}_{name}_fits_raw.tab\", sep='\\t', index_col=0, \n",
    "                       nrows=1, error_bad_lines=False).columns\n",
    "    res = pd.read_csv(f\"{baserunname}_{name}_fits_raw.tab\", sep='\\t', names=cols, \n",
    "                      index_col=0, error_bad_lines=False)\n",
    "    res.drop('Time', inplace=True)\n",
    "    res.to_csv(f\"./{baserunname}_{name}_fits.tab\", sep='\\t')\n",
    "        \n",
    "        \n",
    "def downsample(scriptname, samplefrac, remove=True):\n",
    "    \"\"\"Downsamples an MCMC _sample tab file according to `samplefrac`, \n",
    "    then deletes MCMC _sample and _points files to free up disk space\"\"\"\n",
    "    rawdf = pd.read_csv(f\"{scriptname}_MCMC_sample.tab\", sep='\\t')\n",
    "    newdf = rawdf.sample(frac=samplefrac)\n",
    "    newdf.to_csv(f\"{scriptname}_MCMC_sample_frac.tab\", sep='\\t', index=False)\n",
    "    if remove:\n",
    "        os.remove(f\"{scriptname}_MCMC_sample.tab\")\n",
    "        os.remove(f\"{scriptname}_MCMC_points.tab\")\n",
    "\n",
    "    \n",
    "def merge_samples(baserunname, namelist, name='main', sfx=''):\n",
    "    \"\"\"Combines multiple downsampled MCMC outputs into a single \n",
    "    sensitivity input tabfile, and writes corresponding .vsc file\"\"\"\n",
    "    filelist = [f\"{baserunname}_{n}_MC{sfx}_MCMC_sample_frac.tab\" for n in namelist]\n",
    "    dflist = []\n",
    "    \n",
    "    for f in filelist:\n",
    "        filedf = pd.read_csv(f, sep='\\t')\n",
    "        dflist.append(filedf)\n",
    "    \n",
    "    sensdf = pd.concat(dflist, axis=1)\n",
    "    sensdf.dropna(axis=1, how='all', inplace=True)\n",
    "    sensdf.dropna().to_csv(f\"{baserunname}_{name}_sample_frac{sfx}.tab\", sep='\\t', index=False)\n",
    "\n",
    "    # Generate file input method .vsc file, reading from sample\n",
    "    with open(f\"{baserunname}_full{sfx}.vsc\", 'w') as f:\n",
    "        f.write(f\",F,,{baserunname}_{name}_sample_frac{sfx}.tab,0\")\n",
    "\n",
    "\n",
    "def recalc_weights(scriptname, fitlist, weights, wtsname='CalWtsAdj.cin', dto=False):\n",
    "    \"\"\"Recalculate calibration weights based on residuals from previous \n",
    "    calibration round, for all elms in `fitlist`; read existing weights \n",
    "    from `weights` and output to `wtsname`; if `dto` is specified, de-\n",
    "    trend residuals before resampling (use with caution)\"\"\"\n",
    "    t = pd.read_csv(f'{scriptname}.tab', sep='\\t', index_col=0)\n",
    "    t = t[t.columns[::16]]\n",
    "\n",
    "    # Calculate residuals for each variable\n",
    "    for var in fitlist:\n",
    "        t.loc[f'StDev[{var[0]}]'] = t.loc[f'SimVar[{var[0]}]'] - t.loc[f'DataVar[{var[0]}]']\n",
    "        t.replace({0: np.nan}, inplace=True)\n",
    "\n",
    "    if dto != False: # Detrend residuals using fitted polynomial of order `dto` (usually 1 or 2)\n",
    "        for var in fitlist:\n",
    "            t.loc[f'StDev[{var[0]}]'] = detrend(t.loc[f'StDev[{var[0]}]'].dropna(), order=dto)\n",
    "        \n",
    "    t = t.loc[[f'StDev[{var[0]}]' for var in fitlist]]\n",
    "    stdevs = t.std(axis=1) # Take standard deviations of residuals\n",
    "\n",
    "    with open(weights, 'r') as f0: # Read in any additional weight modifiers from original file\n",
    "        tail = [line for line in f0.readlines() \n",
    "                if all(f'StDev[{var[0]}]' not in line for var in fitlist)]\n",
    "\n",
    "    newwts = [f'{k}={v}\\n' for k, v in stdevs.items()] + tail\n",
    "\n",
    "    with open(wtsname, 'w') as f1:\n",
    "        f1.writelines(newwts)\n",
    "        \n",
    "        \n",
    "def generate_noise(baserunname, fitlist, yearsubs=range(2,33), name='main', sfx='', dto=False):\n",
    "    \"\"\"Calculate and resample normalised residuals for vars in `fitlist` \n",
    "    based on a fitted multivariate normal distribution, over yearsubs as \n",
    "    specified numerically in `yearsubs`, and generate sensitivity input \n",
    "    tabfile containing year x var normalised noise values with matching \n",
    "    .vsc file; note the tabfile may be large and generating may take a \n",
    "    while; if `dto` is specified, detrend normalised residuals before \n",
    "    resampling to remove systematic bias / error (use with caution)\"\"\"\n",
    "    table = pd.read_csv(f'{baserunname}_{name}_fits{sfx}.tab', sep='\\t', index_col=0)\n",
    "    table = table[table.columns[::4]]\n",
    "    \n",
    "    if dto == False:\n",
    "        table = table.loc[[f'NormErr[{var[0]}]' for var in fitlist]]\n",
    "\n",
    "    # Detrend residuals using fitted polynomial of order `dto` (usually 1 or 2)\n",
    "    else:\n",
    "        for var in fitlist:\n",
    "            table.loc[f'NormErrDT[{var[0]}]'] = detrend(\n",
    "                table.loc[f'NormErr[{var[0]}]'].dropna(), order=dto)\n",
    "        table = table.loc[[f'NormErrDT[{var[0]}]' for var in fitlist]]\n",
    "    \n",
    "    sample = pd.read_csv(f'{baserunname}_{name}_sample_frac{sfx}.tab', sep='\\t')\n",
    "\n",
    "    # Calculate means and covariance matrix for normalised residuals\n",
    "    means = table.mean(axis=1)\n",
    "    covs = table.T.cov().astype('float64')\n",
    "\n",
    "    # Generate flattened column list of yearsubs x fitlist vars\n",
    "    cols = [[f'RepErrRaw[{var[0]},Y{year:02d}]' for var in fitlist] for year in yearsubs]\n",
    "    cols_flat = [i for sub in cols for i in sub]\n",
    "\n",
    "    # Draw and flatten new residual sets from multivariate normal distribution\n",
    "    rng = np.random.default_rng()\n",
    "    noise_raw = rng.multivariate_normal(means, covs, (len(sample.index), len(yearsubs)))\n",
    "    noise_flat = [[i for sub in draw for i in sub] for draw in noise_raw]\n",
    "    noise = pd.DataFrame(noise_flat, columns=cols_flat)\n",
    "    noise.clip(-0.99, inplace=True) # Truncate noise values less than -1\n",
    "\n",
    "    # Merge new residuals with sensitivity input tabfile and save\n",
    "    merged = pd.concat([sample, noise], axis=1)\n",
    "\n",
    "    merged.to_csv(f'{baserunname}_{name}_sample_frac_noise{sfx}.tab', sep='\\t', index=False)\n",
    "\n",
    "    # Write corresponding new .vsc file\n",
    "    with open(f\"{baserunname}_full{sfx}.vsc\", 'w') as f:\n",
    "        f.write(f\",F,,{baserunname}_{name}_sample_frac_noise{sfx}.tab,0\")\n",
    "        \n",
    "\n",
    "def generate_param_intervals(baserunname, perc_list, base='main_full', name='main', sfx=''):\n",
    "    \"\"\"Compile dataframe / tabfile of estimated parameter values from \n",
    "    .out and quantiles from MCMC subsample specified in `perc_list`\"\"\"\n",
    "    samdf = pd.read_csv(f\"{baserunname}_{name}_sample_frac{sfx}.tab\", sep='\\t')\n",
    "\n",
    "    quants = samdf.quantile(perc_list).T\n",
    "    quants.index = quants.index.str.strip() # Necessary to strip any trailing spaces\n",
    "\n",
    "    varnames, basevals = read_outvals(f\"{baserunname}_{base}.out\")\n",
    "    \n",
    "    paramdf = pd.DataFrame({'Value': basevals}, index=varnames)\n",
    "    paramdf = paramdf.merge(quants, how='left', left_index=True, right_index=True)\n",
    "    \n",
    "    paramdf.to_csv(f\"{baserunname}_params{sfx}.tab\", sep='\\t')\n",
    "\n",
    "    \n",
    "def clean_sens_raw(senstab, dest=None, remove=True):\n",
    "    \"\"\"Do basic reformat and save output tabfile from sensitivity run\"\"\"\n",
    "    sens = pd.read_csv(senstab, sep='\\t')\n",
    "\n",
    "    senstable = pd.pivot_table(sens, index='Time', columns='Percentile')\n",
    "    senstable = senstable.T\n",
    "    senstable.index.names = ['Var', 'Perc']\n",
    "    senstable = senstable[senstable.columns[::4]] # Switch from 1/16 time step to 1/4\n",
    "\n",
    "    senstable.to_csv(f'{senstab[:-4]}_clean.tab', sep='\\t')\n",
    "    \n",
    "    if dest: # If specified, copy clean file to destination directory\n",
    "        copy(f'./{senstab[:-4]}_clean.tab', dest)\n",
    "        \n",
    "    if remove: # Unless turned off, delete sensitivity VDF to free up space\n",
    "        os.remove(f\"{senstab[:-4]}.vdf\")\n",
    "    del sens\n",
    "\n",
    "\n",
    "def clean_full_sens(senstab, polvars, dest=None, remove=True):\n",
    "    \"\"\"Do basic reformat and save output tabfile from sensitivity run \n",
    "    with full simulation sample saved, saving only `polvars`\"\"\"\n",
    "    colnames = ['Simulation', 'Time'] + polvars # List only relevant columns to subset\n",
    "    sens = pd.read_csv(senstab, sep='\\t', usecols=colnames)\n",
    "\n",
    "    senstable = pd.pivot_table(sens, index='Time', columns='Simulation')\n",
    "    senstable = senstable.T\n",
    "    senstable.index.names = ['Var', 'Run']\n",
    "    senstable = senstable[senstable.columns[::4]] # Switch from 1/16 time step to 1/4\n",
    "\n",
    "    senstable.to_csv(f'{senstab[:-4]}_clean.tab', sep='\\t')\n",
    "\n",
    "    if dest: # If specified, copy clean file to destination directory\n",
    "        copy(f'./{senstab[:-4]}_clean.tab', dest)\n",
    "\n",
    "    if remove: # Unless turned off, delete main sensitivity outputs to free up space\n",
    "        os.remove(senstab)\n",
    "        os.remove(f\"{senstab[:-4]}.vdf\")\n",
    "    del sens        \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Read specified controlfile and unpack into variables\n",
    "controlfilename = input(\"Enter control file name (with extension):\")\n",
    "cf = json.load(open(controlfilename, 'r'))\n",
    "\n",
    "for k,v in cf.items():\n",
    "    exec(k + '=v')\n",
    "\n",
    "for setting in [analysissettings]:\n",
    "    for k, v in setting.items():\n",
    "        exec(k + '=v')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "##### MAIN ITERATIVE CALIBRATION PROCESS AND MCMC #####\n",
    "\n",
    "# Set up files in run directory and initialise logfile\n",
    "master = Script(cf)\n",
    "master.changes.extend(basescens + scenariolist + policylist + testlist)\n",
    "master.copy_model_files(f\"{baserunname}_IterCal\")\n",
    "copy(f\"../{controlfilename}\", \"./\")\n",
    "copy(f\"../{graphs}\", \"./\")\n",
    "copy(f\"../{polsavelist}\", \"./\")\n",
    "basedir = os.getcwd()\n",
    "logfile = f\"{basedir}/{baserunname}.log\"\n",
    "write_log(f\"-----\\nStarting new log at {time.ctime()}\\nReady to work!\", logfile)\n",
    "\n",
    "# Initialise necessary .voc and .vpd files\n",
    "split_voc(simsettings['optparm'], fractolfactor, odlinekey, mcsettings)\n",
    "split_vpd(simsettings['payoff'])\n",
    "\n",
    "# Define encapsulation function for core calibration process\n",
    "def iterative_calibration():\n",
    "    # First do initial calibration round\n",
    "    compile_script(cf, SubdirScript, 'od', 0, {'optparm': '_o', 'payoff': '_o'}, \n",
    "                   logfile, subdir='OD')\n",
    "    payoff_list = [compile_script(cf, SubdirScript, 'main', 0, {'optparm': '_m'}, \n",
    "                                  logfile, chglist=[('od', 0)], subdir='Main')]\n",
    "    payoff_delta = abs(payoff_list[0])\n",
    "    compile_script(cf, SubdirScript, 'init', 0, {'optparm': '_i'}, logfile, \n",
    "                   chglist=[(['od', 'main'], 0)], subdir='Initials')\n",
    "    i = 1\n",
    "\n",
    "    # Then iterate until convergence or until limit is reached\n",
    "    while payoff_delta > threshold:\n",
    "        write_log(f\"More work? Okay! Starting iteration {i}\", logfile)\n",
    "        compile_script(cf, SubdirScript, 'od', i, {'optparm': '_o', 'payoff': '_o'}, \n",
    "                       logfile, chglist=[(['od','main','init'], i-1)], subdir='OD')\n",
    "        payoff_list.append(\n",
    "            compile_script(cf, SubdirScript, 'main', i, {'optparm': '_m'}, logfile, \n",
    "                           chglist=[(['main','init'], i-1), ('od', i)], subdir='Main'))\n",
    "        payoff_delta = abs(payoff_list[-1] - payoff_list[-2])\n",
    "        compile_script(cf, SubdirScript, 'init', i, {'optparm': '_i'}, logfile, \n",
    "                       chglist=[('init', i-1), (['od', 'main'], i)], subdir='Initials')\n",
    "        i += 1\n",
    "\n",
    "        # Increment random number seeds for VOC files        \n",
    "        for sfx in ['_o', '_m', '_i']:\n",
    "            increment_seed(f\"{simsettings['optparm'][:-4]}{sfx}.voc\", logfile)\n",
    "        write_log(f\"Payoff list thus far is {payoff_list}\", logfile)\n",
    "        write_log(f\"Payoff delta is {payoff_delta}\", logfile)\n",
    "\n",
    "        if i > cf['iterlimit']:\n",
    "            write_log(\"Iteration limit reached!\", logfile)\n",
    "            break\n",
    "    else:\n",
    "        write_log(\"Payoff delta is less than threshold. Moving on!\", logfile)\n",
    "\n",
    "    # Run one more full calibration with all parameters\n",
    "    compile_script(cf, LongScript, 'main', 'full', {'optparm': '_f'}, logfile, \n",
    "                   chglist=[(['od', 'main', 'init'], i-1)])\n",
    "\n",
    "# If iterlimit set to 0 (bypass), go straight to all-params Powell optimization\n",
    "if iterlimit == 0:\n",
    "    write_log(\"Iteration is no basis for a system of estimation. Bypassing!\", logfile)\n",
    "    # Skip all-params if previously already done\n",
    "    if os.path.exists(f\"./{baserunname}_main_full.out\"):\n",
    "        if twostep != 0:\n",
    "            write_log(\"Starting with the second time!\", logfile)\n",
    "            cf['simsettings']['changes'].append('CalWtsAdj.cin')\n",
    "        write_log(\"Hang on to outdated imperialist dogma! Using previous output...\", logfile)\n",
    "    else:\n",
    "        compile_script(cf, LongScript, 'main', 'full', {'optparm': '_f'}, logfile)\n",
    "        \n",
    "        # If two-step calibration needed, recalculate weights and calibrate again\n",
    "        if twostep != 0:\n",
    "            write_log(\"I shall do it a second time!\", logfile)\n",
    "            recalc_weights(f'{baserunname}_main_full', fitlist, \n",
    "                           simsettings['changes'][0], dto=noise_dto)\n",
    "            cf['simsettings']['changes'].append('CalWtsAdj.cin')\n",
    "            compile_script(cf, LongScript, 'main', 'full', {'optparm': '_f'}, logfile)\n",
    "            time.sleep(2)\n",
    "\n",
    "else: # Otherwise run iterative calibration process as normal\n",
    "    iterative_calibration() # Run encapsulation function defined above\n",
    "    \n",
    "    # If two-step calibration needed, recalculate weights and calibrate again\n",
    "    if twostep != 0:\n",
    "        write_log(\"I shall do it a second time!\", logfile)\n",
    "        recalc_weights(f'{baserunname}_main_full', fitlist, \n",
    "                       simsettings['changes'][0], dto=noise_dto)\n",
    "        cf['simsettings']['changes'].append('CalWtsAdj.cin')\n",
    "        iterative_calibration() # Run iterative calibration process again with new weights\n",
    "        time.sleep(2)\n",
    "        \n",
    "    # Export fit-to-data results only\n",
    "    export_fits(baserunname, fitlist)\n",
    "    subprocess.run(f\"{vensim7path} \\\"./ExportFits.cmd\\\"\", check=True)\n",
    "    clean_fits(baserunname)\n",
    "\n",
    "# If MCMC option is on, initialise MCMC\n",
    "if mccores != 0:\n",
    "    write_log(\"We're an anarcho-syndicalist commune!\\n\"\n",
    "              f\"Initiating MCMC at {time.ctime()}!\", logfile)\n",
    "    compile_script(cf, MCScript, 'main', 'MC', {'optparm': '_mc'}, \n",
    "                   logfile, chglist=[('main', 'full')])\n",
    "    write_log(f\"MCMC completed at {time.ctime()}!\", logfile)\n",
    "    if twostep != 0:\n",
    "        cf['simsettings']['changes'].pop() # Remember to remove new weights CIN file!\n",
    "\n",
    "    # Compile MCMC subsample & calculate quantiles on parameter estimates\n",
    "    merge_samples(baserunname, ['main']) # Necessary to drop empty column & fix name\n",
    "    generate_param_intervals(baserunname, param_percs)\n",
    "    \n",
    "\n",
    "    # If noise setting is on, generate resampled noise from residuals for MCMC subsample\n",
    "    if noise != 0:\n",
    "        yearsubs = [99] + list(range(0, 33)) # NOTE change this as needed to match years included\n",
    "        generate_noise(baserunname, fitlist, yearsubs=yearsubs, dto=noise_dto)\n",
    "        write_log(\"Ready for something completely different!\", logfile)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "##### BASELINE PROJECTIONS AND PROJECTION ASSUMPTIONS SENSITIVITY #####\n",
    "\n",
    "# Switch to Scenarios subfolder & copy necessary files\n",
    "os.chdir(basedir)\n",
    "os.makedirs('./Results', exist_ok=True)\n",
    "master.copy_model_files(\"Scenarios\")\n",
    "for file in [f\"{controlfilename}\", f\"{baserunname}_main_full.out\", \n",
    "             f\"{baserunname}_full.vsc\", f\"{baserunname}_main_sample_frac.tab\", \n",
    "             f\"{baserunname}_main_sample_frac_noise.tab\"]:\n",
    "    copy(f\"../{file}\", \"./\")\n",
    "\n",
    "# Run basic projections & scenarios with basic & sensitivity runs\n",
    "for cin in (basescens + scenariolist):\n",
    "    copy(f\"../{cin}\", \"./\")\n",
    "    chglist = [('main', 'full'), cin]\n",
    "    write_log(f\"Running scenario {cin}!\", logfile)\n",
    "    compile_script(cf, ScenRunScript, 'final', cin[:-4], {}, logfile, chglist=chglist)\n",
    "    compile_script(cf, ScenSensScript, 'sens', cin[:-4], {}, logfile, chglist=chglist)\n",
    "    time.sleep(2)\n",
    "    clean_sens_raw(f'{baserunname}_sens_{cin[:-4]}.tab', '../Results')\n",
    "    copy(f'./{baserunname}_final_{cin[:-4]}.tab', '../Results')\n",
    "\n",
    "# Run element-wise tests of sensitivity to projection assumptions\n",
    "for cin in basescens[0:2]: # Limit to just main two base scenarios\n",
    "    for proj in proj_subs:\n",
    "        # Create _chg cin files to test single proj element variations\n",
    "        with open(cin, 'r') as f0:\n",
    "            lines = f0.readlines()\n",
    "\n",
    "        # Read existing value of switch and reverse it for proj subscript element\n",
    "        switchval = [line.split('=')[1] for line in lines \n",
    "                     if \"Switch for constant projections[Proj]\" in line]\n",
    "        lines.append(f'\\nSwitch for constant projections[{proj}]={abs(int(switchval[0])-1)}\\n')\n",
    "        with open(f'{cin[:-4]}_chg.cin', 'w') as f1:\n",
    "            f1.writelines(lines)\n",
    "        \n",
    "        # Run basic and sensitivity runs with assumption test .cin\n",
    "        scen = f'{cin[:-4]}{proj}'\n",
    "        chglist = [('main', 'full'), f'{cin[:-4]}_chg.cin']\n",
    "        write_log(f\"Running scenario {scen}!\", logfile)\n",
    "        compile_script(cf, ScenRunScript, 'final', scen, {}, logfile, chglist=chglist)\n",
    "        compile_script(cf, ScenSensScript, 'sens', scen, {}, logfile, chglist=chglist)\n",
    "        time.sleep(2)\n",
    "        clean_sens_raw(f'{baserunname}_sens_{scen}.tab', '../Results')\n",
    "        copy(f'./{baserunname}_final_{scen}.tab', '../Results')\n",
    "    \n",
    "os.chdir(\"..\")\n",
    "\n",
    "print(\"Ready and waiting!\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "##### POLICY SCENARIOS & SENSITIVITY #####\n",
    "\n",
    "# Switch to Policies subfolder & copy necessary files\n",
    "os.chdir(basedir)\n",
    "os.makedirs('./Results', exist_ok=True)\n",
    "master.copy_model_files(\"Policies\")\n",
    "for file in [f\"{controlfilename}\", f\"{baserunname}_main_full.out\", \n",
    "             f\"{baserunname}_full.vsc\", f\"{baserunname}_main_sample_frac.tab\", \n",
    "             f\"{baserunname}_main_sample_frac_noise.tab\", polsavelist] + policylist:\n",
    "    copy(f\"../{file}\", \"./\")\n",
    "\n",
    "# Switch to policy savelist\n",
    "senssavelist = cf['simsettings']['senssavelist'] # Record regular savelist to switch back later\n",
    "cf['simsettings']['senssavelist'] = polsavelist\n",
    "\n",
    "# Run each baseline scenario in combination with each policy \n",
    "for cin in basescens:\n",
    "    for pol in policylist:\n",
    "        # Run basic and sensitivity runs for each combination\n",
    "        scen = f'{cin[:-4]}{pol[:-4]}'\n",
    "        chglist = [('main', 'full'), cin, pol]\n",
    "        write_log(f\"Running scenario {scen}!\", logfile)\n",
    "        compile_script(cf, ScenRunScript, 'final', scen, {}, logfile, chglist=chglist)\n",
    "        compile_script(cf, FullSensScript, 'sens', scen,\n",
    "                       {}, logfile, chglist=chglist) # NOTE: saves full sensitivity sample\n",
    "        time.sleep(2)\n",
    "        clean_full_sens(f'{baserunname}_sens_{scen}.tab', polvars, '../Results')\n",
    "        copy(f'./{baserunname}_final_{scen}.tab', '../Results')\n",
    "\n",
    "cf['simsettings']['senssavelist'] = senssavelist # Switch back to regular savelist\n",
    "os.chdir(\"..\")\n",
    "\n",
    "print(\"Ready and waiting!\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "##### TEMPORARY ######\n",
    "##### PAIRWISE POLICY SCENARIO COMBINATIONS ######\n",
    "\n",
    "# Switch to Policies subfolder & copy necessary files\n",
    "os.chdir(basedir)\n",
    "os.makedirs('./Results', exist_ok=True)\n",
    "master.copy_model_files(\"PolicyCombos\")\n",
    "for file in [f\"{controlfilename}\", f\"{baserunname}_main_full.out\", \n",
    "             f\"{baserunname}_full.vsc\", f\"{baserunname}_main_sample_frac.tab\", \n",
    "             f\"{baserunname}_main_sample_frac_noise.tab\", polsavelist] + policylist:\n",
    "    copy(f\"../{file}\", \"./\")\n",
    "\n",
    "# Switch to policy savelist\n",
    "senssavelist = cf['simsettings']['senssavelist'] # Record regular savelist to switch back later\n",
    "cf['simsettings']['senssavelist'] = polsavelist\n",
    "\n",
    "pollist = [\"BupProvBarriers.cin\", \"BupProv.cin\", \"DevOUD.cin\", \"DivRxInit.cin\", \n",
    "           \"FentODRisk.cin\", \"HInit.cin\", \"NxKits.cin\", \"OwnRxInit.cin\", \n",
    "           \"PeerRecovery.cin\",  \"RxRate.cin\", \"ReturntoOUD.cin\"]\n",
    "pollist10 = [\"BupProvBarriers10.cin\", \"BupProv10.cin\", \"DevOUD10.cin\", \"DivRxInit10.cin\", \n",
    "             \"FentODRisk10.cin\", \"HInit10.cin\", \"NxKits10.cin\", \"OwnRxInit10.cin\", \n",
    "             \"PeerRecovery10.cin\", \"RxRate10.cin\", \"ReturntoOUD10.cin\"]\n",
    "pollist50 = [\"BupProvBarriers50.cin\", \"BupProv50.cin\", \"DevOUD50.cin\", \"DivRxInit50.cin\", \n",
    "             \"FentODRisk50.cin\", \"HInit50.cin\", \"NxKits50.cin\", \"OwnRxInit50.cin\", \n",
    "             \"PeerRecovery50.cin\", \"RxRate50.cin\", \"ReturntoOUD50.cin\"]\n",
    "\n",
    "for cin in basescens[0:1]:\n",
    "    for pols in [pollist, pollist10, pollist50]:\n",
    "        for i, pol in enumerate(pols[:-1]):\n",
    "            for pol2 in pols[i+1:]:\n",
    "                scen = f'{cin[:-4]}{pol[:-4]}{pol2[:-4]}'\n",
    "                chglist = [('main', 'full'), cin, pol, pol2]\n",
    "                write_log(f\"Running scenario combo {pol} + {pol2}!\", logfile)\n",
    "                compile_script(cf, ScenRunScript, 'final', scen, {}, logfile, chglist=chglist)\n",
    "#                 compile_script(cf, ScenSensScript, 'sens', scen, {}, logfile, chglist=chglist)\n",
    "                time.sleep(2)\n",
    "#                 clean_full_sens(f'{baserunname}_sens_{scen}.tab', polvars, '../Results')\n",
    "                copy(f'./{baserunname}_final_{scen}.tab', '../Results')\n",
    "\n",
    "cf['simsettings']['senssavelist'] = senssavelist # Switch back to regular savelist\n",
    "os.chdir(\"..\")\n",
    "\n",
    "print(\"Ready and waiting!\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "##### LOOP KNOCKOUT ANALYSES #####\n",
    "\n",
    "# Switch to Sensitivity subfolder & copy necessary files\n",
    "os.chdir(basedir)\n",
    "master.copy_model_files(\"Sensitivity\")\n",
    "copy(f\"../{controlfilename}\", \"./\")\n",
    "copy(f\"../{baserunname}_main_full.out\", \"./\")\n",
    "\n",
    "split_lk_vocs(simsettings['optparm'], lkdict) # Create necessary .vocs for loop knockout tests\n",
    "\n",
    "# Iterate through sets of loops / vars to inactivate as specified in lkdict\n",
    "for key in lkdict.keys():\n",
    "    write_log(f\"It's {time.ctime()}! Time for a new... idiom!\", logfile)\n",
    "    \n",
    "    # First run crude (no re-estimation) loop knockout projection using .cin\n",
    "    compile_script(cf, ScenRunScript, 'lk_run', f'{key}_{basescens[0][:-4]}', {}, logfile,\n",
    "                   chglist=[('main', 'full'), basescens[0], f'LK_{key}.cin'])\n",
    "    time.sleep(2)\n",
    "    \n",
    "    # Then re-estimate with loops deactivated, and re-run projections\n",
    "    compile_script(cf, LongScript, 'lk', key, {'optparm': f'_lk_{key}'}, logfile)\n",
    "    time.sleep(2)\n",
    "    \n",
    "    # If two-step calibration needed, recalculate weights and calibrate again\n",
    "    if twostep != 0:\n",
    "        recalc_weights(f'{baserunname}_lk_{key}', fitlist, \n",
    "                       simsettings['changes'][0], wtsname=f'CalWtsAdj_{key}.cin')\n",
    "        cf['simsettings']['changes'].append(f'CalWtsAdj_{key}.cin')\n",
    "        compile_script(cf, LongScript, 'lk', key, {'optparm': f'_lk_{key}'}, logfile)\n",
    "        time.sleep(2)\n",
    "        cf['simsettings']['changes'].pop() # Remember to remove new weights CIN file!\n",
    "\n",
    "    compile_script(cf, ScenRunScript, 'lk', f'{key}_{basescens[0][:-4]}', \n",
    "                   {}, logfile, chglist=[('lk', key), basescens[0]]) # Re-run baseline projections\n",
    "    time.sleep(2)\n",
    "\n",
    "compile_sensoutvals(baserunname, lkdict, 'lk') # Compile summary sensitivity table\n",
    "write_log(f\"That's not right for my idiom! Leaving at {time.ctime()}!\", logfile)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "##### PARAMETRIC ASSUMPTION ELASTICITY ANALYSES #####\n",
    "\n",
    "# Switch to Sensitivity subfolder & copy necessary files\n",
    "os.chdir(basedir)\n",
    "master.copy_model_files(\"Sensitivity\")\n",
    "copy(f\"../{controlfilename}\", \"./\")\n",
    "copy(f\"../{baserunname}_main_full.out\", \"./\")\n",
    "copy(f\"../{simsettings['optparm'][:-4]}_f.voc\", \"./\")\n",
    "\n",
    "# Compile dict of short form names from first letters of sensvars\n",
    "sensdict = dict([[''.join([w[0] for w in regex.findall(r\"[\\w']+\", var)]), var] \n",
    "                 for var in sensvars])\n",
    "\n",
    "for key, var in sensdict.items():\n",
    "    baseval = get_value(simsettings['model'], var) # Get base value directly from .mdl\n",
    "    \n",
    "    # Create low and high value assumption .cin files\n",
    "    with open(f'{key}_L.cin', 'w') as f:\n",
    "        f.write(f'{var}={baseval*(1 - sensrange)}')\n",
    "    with open(f'{key}_H.cin', 'w') as f:\n",
    "        f.write(f'{var}={baseval*(1 + sensrange)}')\n",
    "    \n",
    "    write_log(f\"It's {time.ctime()} - how do we know if {var} is made of wood?\", logfile)\n",
    "    \n",
    "    # Run each scenario with re-estimation first then projections\n",
    "    for scen in [f'{key}_L', f'{key}_H']:\n",
    "        cf['simsettings']['changes'].append(f'{scen}.cin')\n",
    "        compile_script(cf, LongScript, 'assm', scen, {'optparm': '_f'}, \n",
    "                       logfile, chglist=[('main', 'full')])\n",
    "        time.sleep(2)\n",
    "        \n",
    "        # If two-step calibration needed, recalculate weights and calibrate again\n",
    "        if twostep != 0:\n",
    "            recalc_weights(f'{baserunname}_assm_{scen}', fitlist, \n",
    "                           simsettings['changes'][0], wtsname=f'CalWtsAdj_{scen}.cin')\n",
    "            cf['simsettings']['changes'].append(f'CalWtsAdj_{scen}.cin')\n",
    "            compile_script(cf, LongScript, 'assm', scen, {'optparm': '_f'}, \n",
    "                           logfile, chglist=[('main', 'full')])\n",
    "            time.sleep(2)\n",
    "            cf['simsettings']['changes'].pop() # Remember to remove new weights CIN file!\n",
    "        \n",
    "        compile_script(cf, ScenRunScript, 'assm', f'{scen}_{basescens[0][:-4]}', {}, \n",
    "                       logfile, chglist=[('assm', scen), basescens[0]]) # Run baseline projection\n",
    "        time.sleep(2)\n",
    "        cf['simsettings']['changes'].pop() # Remember to remove assumption CIN file!\n",
    "\n",
    "# Compile summary sensitivity table\n",
    "compile_sensoutvals(baserunname, sensdict, 'assm', sfxs=['_L', '_H'])\n",
    "write_log(f\"It's {time.ctime()}! Burn them!\", logfile)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "##### ALTERNATIVE DATA CONDITION ANALYSES #####\n",
    "\n",
    "# Switch to Sensitivity subfolder & copy necessary files\n",
    "os.chdir(basedir)\n",
    "master.copy_model_files(\"Sensitivity\")\n",
    "copy(f\"../{controlfilename}\", \"./\")\n",
    "copy(f\"../{baserunname}_main_full.out\", \"./\")\n",
    "copy(f\"../{simsettings['optparm'][:-4]}_f.voc\", \"./\")\n",
    "for file in testlist:\n",
    "    copy(f\"../{file}\", \"./\")\n",
    "\n",
    "# Run calibration for each alternative data condition scenario\n",
    "for cin in testlist:\n",
    "    cf['simsettings']['changes'].append(cin)\n",
    "    compile_script(cf, LongScript, 'test', cin[:-4], {'optparm': '_f'}, \n",
    "                   logfile, chglist=[('main', 'full')])\n",
    "    time.sleep(2)\n",
    "    \n",
    "    # If two-step calibration needed, recalculate weights and calibrate again\n",
    "    if twostep != 0:\n",
    "        recalc_weights(f'{baserunname}_test_{cin[:-4]}', fitlist, \n",
    "                       simsettings['changes'][0], wtsname=f'CalWtsAdj_{cin[:-4]}.cin')\n",
    "        cf['simsettings']['changes'].append(f'CalWtsAdj_{cin[:-4]}.cin')\n",
    "        compile_script(cf, LongScript, 'test', cin[:-4], {'optparm': '_f'}, \n",
    "                       logfile, chglist=[('main', 'full')])\n",
    "        time.sleep(2)\n",
    "        cf['simsettings']['changes'].pop() # Remember to remove new weights CIN file!    \n",
    "    \n",
    "    compile_script(cf, ScenRunScript, 'test', f'{cin[:-4]}_{basescens[0][:-4]}', {}, \n",
    "                   logfile, chglist=[('test', cin[:-4]), basescens[0]]) # Run baseline projection\n",
    "\n",
    "    time.sleep(2)\n",
    "    cf['simsettings']['changes'].pop() # Remember to remove assumption CIN file!\n",
    "\n",
    "# Compile summary sensitivity table\n",
    "compile_sensoutvals(baserunname, [cin[:-4] for cin in testlist], 'test')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "##### HOLDOUT DATA ANALYSIS #####\n",
    "\n",
    "# Switch to Sensitivity subfolder & copy necessary files\n",
    "os.chdir(basedir)\n",
    "master.copy_model_files(\"Sensitivity\")\n",
    "copy(f\"../{controlfilename}\", \"./\")\n",
    "copy(f\"../{simsettings['optparm'][:-4]}_f.voc\", \"./\")\n",
    "copy(f\"../{simsettings['optparm'][:-4]}_mc.voc\", \"./\")\n",
    "\n",
    "# Run initial calibration using data up to `holdoutyear`\n",
    "write_log(f\"At {time.ctime()}, none shall pass after {holdoutyear}!\", logfile)\n",
    "compile_script(cf, LongScript, 'hold', 'full', {'optparm': '_f'}, \n",
    "               logfile, setvals=[('MaxDataTime', holdoutyear)])\n",
    "time.sleep(2)\n",
    "    \n",
    "# If two-step calibration needed, recalculate weights and calibrate again\n",
    "if twostep != 0:\n",
    "    recalc_weights(f'{baserunname}_hold_full', fitlist, \n",
    "                   simsettings['changes'][0], wtsname='CalWtsAdj_hold.cin')\n",
    "    cf['simsettings']['changes'].append('CalWtsAdj_hold.cin')\n",
    "    compile_script(cf, LongScript, 'hold', 'full', {'optparm': '_f'}, \n",
    "                   logfile, setvals=[('MaxDataTime', holdoutyear)])\n",
    "    time.sleep(2)\n",
    "\n",
    "# Export fit-to-data results only\n",
    "export_fits(baserunname, fitlist, base='hold_full', name='hold')\n",
    "subprocess.run(f\"{vensim7path} \\\"./ExportFits.cmd\\\"\", check=True)\n",
    "clean_fits(baserunname, name='hold')\n",
    "\n",
    "# Run MCMC following initial calibration\n",
    "compile_script(cf, MCScript, 'hold', 'MC', {'optparm': '_mc'}, \n",
    "               logfile, chglist=[('hold', 'full')], setvals=[('MaxDataTime', holdoutyear)])\n",
    "write_log(f\"Just a flesh wound. MCMC completed at {time.ctime()}!\", logfile)\n",
    "if twostep != 0:\n",
    "    cf['simsettings']['changes'].pop() # Remember to remove new weights CIN file!\n",
    "\n",
    "# Compile MCMC subsample & calculate quantiles on parameter estimates\n",
    "merge_samples(baserunname, ['hold'], name='hold') # Necessary to drop empty column & fix name\n",
    "\n",
    "# Generate resampled noise from residuals for MCMC subsample\n",
    "yearsubs = [99] + list(range(0, 33)) # NOTE change this as needed to match years included\n",
    "generate_noise(baserunname, fitlist, yearsubs=yearsubs, name='hold', dto=1)\n",
    "\n",
    "# Run baseline projection and sensitivity\n",
    "compile_script(cf, ScenRunScript, 'hold', f'final_{basescens[0][:-4]}', {}, \n",
    "               logfile, chglist=[('hold', 'full'), basescens[0]])\n",
    "compile_script(cf, ScenSensScript, 'hold', f'sens_{basescens[0][:-4]}', {}, logfile, \n",
    "               chglist=[('hold', 'full'), basescens[0]], setvals=[('NoiseStartTime', holdoutyear), ('Switch for historical noise', 1)])\n",
    "time.sleep(2)\n",
    "\n",
    "clean_sens_raw(f'{baserunname}_hold_sens_{basescens[0][:-4]}.tab', '../Results', remove=False)\n",
    "copy(f'./{baserunname}_hold_final_{basescens[0][:-4]}.tab', '../Results')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "##### SYNTHETIC DATA GENERATION AND ANALYSES #####\n",
    "\n",
    "os.chdir(basedir)\n",
    "\n",
    "# Create syndata generation sample from main noise dataframe\n",
    "rawdf = pd.read_csv(f\"{baserunname}_main_sample_frac_noise.tab\", sep='\\t')\n",
    "syndatdf = rawdf.sample(synsample) # Randomly subsample `synsample` param & noise sets\n",
    "del rawdf\n",
    "syndatdf.index = range(synsample)\n",
    "syndatdf.columns = syndatdf.columns.str.strip() # Necessary to strip any trailing spaces\n",
    "syndatdf.to_csv(f\"{baserunname}_syndata_sample.tab\", sep='\\t', index=False)\n",
    "\n",
    "# Switch to SynData subfolder & copy necessary files\n",
    "master.copy_model_files(\"SynData\")\n",
    "\n",
    "for file in [f\"{controlfilename}\", f\"{simsettings['optparm'][:-4]}_f.voc\", \n",
    "             f\"{simsettings['optparm'][:-4]}_mc.voc\"]:\n",
    "    copy(f\"../{file}\", \"./\")\n",
    "\n",
    "# Create model version to read syndata as data\n",
    "synmodel = simsettings['model'][:-4] + '_Syn.mdl'\n",
    "create_syn_mdl(simsettings['model'], synmodel)\n",
    "syncf = json.load(open(controlfilename, 'r'))\n",
    "\n",
    "syn_res_list = [] # Initialise container for results of each syndata run\n",
    "\n",
    "for i in range(synsample):\n",
    "    # Create .cin files for syndata generation from sample\n",
    "    values = [f'NoiseStartTime = {styear}\\n', 'Switch for historical noise = 1\\n']\n",
    "    values.extend([f'{k} = {v}\\n' for k,v in syndatdf.loc[i].to_dict().items()])\n",
    "    with open(f'SynData{i}.cin', 'w') as f:\n",
    "        f.writelines(values)\n",
    "    \n",
    "    # Run base model with 'true' syndata values from .cin to generate syndata\n",
    "    write_log(f\"Building a large wooden rabbit, round {i}!\", logfile)\n",
    "    compile_script(cf, ScenRunScript, 'syndata', i, {}, logfile, chglist=[f'SynData{i}.cin'])\n",
    "    time.sleep(2)\n",
    "    \n",
    "    # Run estimation and MCMC on synthetic data\n",
    "    write_log(\"Now we leap out of the rabbit...\", logfile)\n",
    "    syncf['simsettings']['data'] = [f'{baserunname}_syndata_{i}.vdf', 'yearsubs.vdf']\n",
    "    compile_script(syncf, LongScript, 'syn', i, {'model': '_Syn', 'optparm': '_f'}, logfile)\n",
    "    time.sleep(2)\n",
    "\n",
    "    # If two-step calibration needed, recalculate weights and calibrate again\n",
    "    if twostep != 0:\n",
    "        recalc_weights(f'{baserunname}_syn_{i}', fitlist, \n",
    "                       simsettings['changes'][0], wtsname=f'CalWtsAdj_{i}.cin')\n",
    "        cf['simsettings']['changes'].append(f'CalWtsAdj_{i}.cin')\n",
    "        compile_script(syncf, LongScript, 'syn', i, {'model': '_Syn', 'optparm': '_f'}, logfile)\n",
    "        time.sleep(2)\n",
    "    \n",
    "    write_log(f\"...and at {time.ctime()}, take them by surprise!\", logfile)\n",
    "    compile_script(syncf, MCScript, 'syn', f'MC_{i}', {'model': '_Syn', 'optparm': '_mc'}, \n",
    "                   logfile, chglist=[('syn', i)])\n",
    "    write_log(f\"Surprise completed at {time.ctime()}!\", logfile)\n",
    "    if twostep != 0:\n",
    "        cf['simsettings']['changes'].pop() # Remember to remove new weights CIN file!\n",
    "    \n",
    "    # Generate CrIs for estimated values\n",
    "    merge_samples(baserunname, ['syn'], sfx=f'_{i}') # Necessary to drop empty column & fix name\n",
    "    generate_param_intervals(baserunname, syn_percs, base=f'syn_{i}', name='main', sfx=f'_{i}')\n",
    "    \n",
    "    # Combine estimated params and CrIs with 'true' syndata values from original sample\n",
    "    paramdf = pd.read_csv(f'{baserunname}_params_{i}.tab', sep='\\t', index_col=0)\n",
    "    resdf = pd.concat([syndatdf.loc[i], paramdf], axis=1)\n",
    "    resdf.dropna(inplace=True)\n",
    "    resdf.rename(columns={i: 'True'}, inplace=True) # Label 'true' param values\n",
    "    resdf = resdf.T\n",
    "    \n",
    "    syn_res_list.append(resdf)\n",
    "\n",
    "# Compile and export combined results of all syndata estimations w/ CIs\n",
    "synresdf = pd.concat(syn_res_list, keys=range(synsample), names=['Run', 'Perc'])\n",
    "synresdf.to_csv(f'{baserunname}_syndata_results.tab', sep='\\t')\n",
    "copy(f'./{baserunname}_syndata_results.tab', '../Results')\n",
    "display(synresdf)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "##### CALCULATE PSRF STATISTICS #####\n",
    "\n",
    "os.chdir(basedir)\n",
    "\n",
    "name = 'main' # Change if main run name is different\n",
    "\n",
    "# Compile .cmd script to convert stats .dat to tabfile\n",
    "cmdtext = [\n",
    "    \"SPECIAL>NOINTERACTION\\n\",\n",
    "    f\"MENU>DAT2VDF|{baserunname}_{name}_MC_MCMC_stats.dat\\n\",\n",
    "    f\"SIMULATE>RUNNAME|{baserunname}_{name}_MC\\n\",\n",
    "    f\"MENU>VDF2TAB|{baserunname}_{name}_MC_MCMC_stats|{baserunname}_{name}_MC_MCMC_stats|\\n\",\n",
    "    \"SPECIAL>CLEARRUNS\\n\",\n",
    "    \"MENU>EXIT\\n\"\n",
    "]\n",
    "    \n",
    "with open(f\"{baserunname}_{name}_PSRF.cmd\", 'w') as scriptfile:\n",
    "    scriptfile.writelines(cmdtext)\n",
    "\n",
    "# Run .dat conversion with Vensim\n",
    "while True:\n",
    "    subprocess.run(f\"{vensim7path} \\\"./{baserunname}_{name}_PSRF.cmd\\\"\")\n",
    "    time.sleep(1)\n",
    "    try:\n",
    "        copy(f\"./{baserunname}_{name}_MC_MCMC_stats.tab\", \"./Results\")\n",
    "        break\n",
    "    except FileNotFoundError:\n",
    "        print(f\"Help! {name} is being repressed!\")\n",
    "        continue\n",
    "\n",
    "os.chdir(\"..\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "000000000000000000000000000000000000000000000000000000000000000000000000\n",
    "000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}