optimisation-unified-canada.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "\n",
    "import sys\n",
    "\n",
    "from smrt import make_snowpack, make_model, sensor_list, make_soil, SMRTError\n",
    "from smrt.inputs.make_medium import make_medium\n",
    "\n",
    "%load_ext autoreload\n",
    "%autoreload 2\n",
    "\n",
    "%matplotlib widget\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sp_xls = pd.ExcelFile(\"data/CanadaData/SP.xlsx\")\n",
    "\n",
    "obs = pd.read_csv(\"data/CanadaData/TB.csv\")\n",
    "\n",
    "# each sheet is one snowpit\n",
    "sps = [sp for sp in sp_xls.sheet_names]\n",
    "len(sps)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "def prepare_snowpit(sp):\n",
    "\n",
    "    profile = sp_xls.parse(sp)\n",
    "\n",
    "    temp_substrate = min(profile['TempSol'].iloc[0], 273.15)\n",
    "\n",
    "    profile['z'] = profile['Hneige'] * 0.01\n",
    "\n",
    "    profile = profile.drop(columns=['TempSol', 'TEMPERATURE', 'IRIS', 'Hneige', 'Hneige_cm', 'Ropt', 'Hneige (m)']).dropna()\n",
    "\n",
    "    profile = profile.rename(columns={'Temp_Kel': 'temperature',\n",
    "                                      'Density': 'density'})\n",
    "\n",
    "    return temp_substrate, profile\n",
    "        \n",
    "snowpits = {sp: prepare_snowpit(sp) for sp in sps}\n",
    "\n",
    "#database_with_result['depth'] = [snowpits[sp][1]['z'].abs().max() for sp in database_with_result.index]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "freqs = ['19', '37']\n",
    "\n",
    "\n",
    "def run_model(temp_substrate, profile, microstructure, coef):\n",
    "    density_ice = 917.\n",
    "    porod_length = 4 * (1 - profile.density / density_ice) / (density_ice * profile.SSA)\n",
    "    \n",
    "    K_TS = np.where(profile.SSA > 15, 0.6, coef)\n",
    "    K_SEXP = np.where(profile.SSA > 15, 0.625, coef)\n",
    "    K_SHS = np.where(profile.SSA > 15, 0.64, coef)\n",
    "        \n",
    "    microstructure_args = {\n",
    "        \"USEXP\": dict(microstructure_model=\"unified_scaled_exponential\", porod_length=porod_length, polydispersity=K_SEXP),\n",
    "        \"UTS\": dict(microstructure_model=\"unified_teubner_strey\", porod_length=porod_length, polydispersity=K_TS),\n",
    "        \"USHS\": dict(microstructure_model=\"unified_sticky_hard_spheres\", porod_length=porod_length, polydispersity=K_SHS)\n",
    "    }\n",
    "\n",
    "    substrate = make_soil(\"soil_wegmuller\", permittivity_model=\"montpetit2008\",\n",
    "                          temperature=temp_substrate, roughness_rms=1.93e-3)\n",
    "\n",
    "    sp = {m: make_medium(profile, **microstructure_args[m], substrate=substrate) for m in microstructure}\n",
    "\n",
    "    model = make_model(\"iba\", \"dort\", rtsolver_options=dict(prune_deep_snowpack=8, error_handling='nan'))\n",
    "    sensor = sensor_list.amsr2(freqs)\n",
    "\n",
    "    results = model.run(sensor, sp, parallel_computation=True)\n",
    "\n",
    "    results = results.to_dataframe()\n",
    "    results = results.unstack()\n",
    "    results.index = [\"_\".join(ind) for ind in results.index]  # collapse multiindex by joining names with _\n",
    "\n",
    "    return results\n",
    "\n",
    "\n",
    "def run_all_site(microstructure, coef):\n",
    "    results = []\n",
    "    for sp in snowpits:\n",
    "        if snowpits[sp] is None:\n",
    "            continue\n",
    "        res = run_model(*snowpits[sp], microstructure, coef)\n",
    "        if res is None:\n",
    "            continue\n",
    "        results.append(res)\n",
    "\n",
    "    results = pd.DataFrame(results, index=snowpits.keys())\n",
    "    database_with_result = results.join(obs.set_index(\"CorresSP\"))\n",
    "\n",
    "    return database_with_result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "def compute_error(microstructure, coef_list):\n",
    "\n",
    "    errors = []\n",
    "    for coef in coef_list:\n",
    "        database_with_result = run_all_site(microstructure=[microstructure], coef=coef)\n",
    "        #print(database_with_result)\n",
    "        d = {'polydispersity': coef}\n",
    "\n",
    "        for freq in freqs:\n",
    "            d[f'rmse_{freq}'] = np.sqrt(((database_with_result[f'{freq}V MES'] - database_with_result[f'{freq}V_{microstructure}'])**2).mean())\n",
    "            d[f'abse_{freq}'] = (database_with_result[f'{freq}V MES'] - database_with_result[f'{freq}V_{microstructure}']).abs().mean()\n",
    "            d[f'bias_{freq}'] = (database_with_result[f'{freq}V MES'] - database_with_result[f'{freq}V_{microstructure}']).mean()\n",
    "            d[f'count_{freq}'] = database_with_result[f'{freq}V_{microstructure}'].dropna().count()\n",
    "        errors.append(d)\n",
    "    return pd.DataFrame(errors).set_index('polydispersity')\n",
    "\n",
    "polydispersity = np.arange(1, 3, 0.05)\n",
    "# polydispersity = np.arange(1.3, 1.4, 0.01)  # refine\n",
    "\n",
    "lazzy = True\n",
    "update = False\n",
    "\n",
    "errors = {}\n",
    "for m in ['UTS', 'USEXP', 'USHS']:\n",
    "    \n",
    "    filename = f\"results/simulations-errors-canada-{m.lower()}.csv\"\n",
    "    file_exists = os.path.exists(filename)\n",
    "    \n",
    "    if (lazzy or update) and file_exists:\n",
    "        errors[m] = pd.read_csv(filename).set_index('polydispersity').drop_duplicates()\n",
    "    \n",
    "    if update or not file_exists:\n",
    "        err = compute_error(m, polydispersity)\n",
    "        if update and file_exists:\n",
    "            err = pd.concat((err, errors[m])).sort_index().drop_duplicates()\n",
    "        errors[m] = err\n",
    "\n",
    "        errors[m].to_csv(filename)\n",
    "        \n",
    "    errors[m]['rmse'] = np.sqrt((errors[m]['rmse_19']**2 + errors[m]['rmse_37']**2) / 2)\n",
    "    errors[m]['bias'] = (errors[m]['bias_19'] + errors[m]['bias_37']) / 2\n",
    "    \n",
    "    imin = errors[m]['rmse'].argmin()\n",
    "    #print(m, errors[m].index[imin])\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for m in errors:\n",
    "    Koptimal = errors[m]['rmse'].idxmin()\n",
    "    print(m, Koptimal, errors[m].loc[Koptimal]['rmse'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f, axs = plt.subplots(1, 3, figsize=(6, 3), sharey=True, sharex=True)\n",
    "\n",
    "axs[2].plot(errors['UTS'].index, errors['UTS']['rmse'], '-',)\n",
    "axs[2].plot(errors['UTS'].index, abs(errors['UTS']['bias']), '-',)\n",
    "#axs[2].set_xlabel('Repeat coefficient $q$')  # = d / 2\\\\pi \\\\xi$')\n",
    "axs[2].set_ylabel('RMSE 19 and 37 GHz, V-pol (K)')\n",
    "#axs[2].set_ylabel('Error (K)')\n",
    "#axs[2].set_ylim((0, 35))\n",
    "axs[2].set_title('Teubner-Strey')\n",
    "\n",
    "axs[0].plot(errors['USEXP'].index, errors['USEXP']['rmse'], '-', label='RMSE')\n",
    "axs[0].plot(errors['USEXP'].index, abs(errors['USEXP']['bias']), '-', label='|bias|')\n",
    "axs[0].set_title('Scaled exponential')\n",
    "#axs[0].set_xlabel('Scaling coefficient $\\\\phi$')\n",
    "axs[0].legend()\n",
    "\n",
    "axs[1].plot(errors['USHS'].index, errors['USHS']['rmse'], '-')\n",
    "axs[1].plot(errors['USHS'].index, abs(errors['USHS']['bias']), '-')\n",
    "\n",
    "#axs[2].set_xlabel('Stickiness $\\\\tau$')\n",
    "axs[1].set_title('Sticky Hard Sphere')\n",
    "\n",
    "for i in [0, 1, 2]:\n",
    "    axs[i].grid(alpha=0.2)\n",
    "    axs[i].set_xlabel('Polydispersity $K$')\n",
    "\n",
    "axs[2].set_xlim((1, 2.3))\n",
    "f.tight_layout()\n",
    "plt.savefig(\"fig-global-optimisation-coefs-arctic.pdf\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python (smrt)",
   "language": "python",
   "name": "smrt"
  },
  "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.9.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}