DeepDA_allMC.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "The config (yml) file saved : \n",
      "\n",
      "/Users/mingsongli/Dropbox/PETM/posterior/petmproxy3slices_v0.1.5shift.csv_petm38_v0.1.5shift_20240728_All.bays_MC100_pHcor_frac0.99_Ca75_biogem_SatuBen_setLim.yml\n",
      "\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/mingsongli/miniconda3/envs/deepda/lib/python3.6/site-packages/ipykernel_launcher.py:374: DeprecationWarning: tostring() is deprecated. Use tobytes() instead.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "##########  Monte Carlo 1 / 100 => okay   ##########\n",
      "\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/mingsongli/miniconda3/envs/deepda/lib/python3.6/site-packages/pandas/core/generic.py:2449: PerformanceWarning: \n",
      "your performance may suffer as PyTables will pickle object types that it cannot\n",
      "map directly to c-types [inferred_type->mixed,key->block1_values] [items->Index(['Site', 'Type', 'DepthOri', 'SourcePaleodepth', 'Unnamed: 12', 'Mask',\n",
      "       'assessment', 'File', 'doi', 'Proxy', 'Glassy', 'mg_bcp_red',\n",
      "       'JESS NOTES'],\n",
      "      dtype='object')]\n",
      "\n",
      "  encoding=encoding,\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "\n",
      "##########  Monte Carlo 2 / 100 => okay   ##########\n",
      "\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "\n",
      "##########  Monte Carlo 3 / 100 => okay   ##########\n",
      "\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "\n",
      "##########  Monte Carlo 4 / 100 => okay   ##########\n",
      "\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "\n",
      "##########  Monte Carlo 5 / 100 => okay   ##########\n",
      "\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "\n",
      "##########  Monte Carlo 6 / 100 => okay   ##########\n",
      "\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "\n",
      "##########  Monte Carlo 7 / 100 => okay   ##########\n",
      "\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "\n",
      "##########  Monte Carlo 8 / 100 => okay   ##########\n",
      "\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "\n",
      "##########  Monte Carlo 9 / 100 => okay   ##########\n",
      "\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "\n",
      "##########  Monte Carlo 10 / 100 => okay   ##########\n",
      "\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "\n",
      "##########  Monte Carlo 11 / 100 => okay   ##########\n",
      "\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "\n",
      "##########  Monte Carlo 12 / 100 => okay   ##########\n",
      "\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "\n",
      "##########  Monte Carlo 13 / 100 => okay   ##########\n",
      "\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "\n",
      "##########  Monte Carlo 14 / 100 => okay   ##########\n",
      "\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "\n",
      "##########  Monte Carlo 15 / 100 => okay   ##########\n",
      "\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "\n",
      "##########  Monte Carlo 16 / 100 => okay   ##########\n",
      "\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "\n",
      "##########  Monte Carlo 17 / 100 => okay   ##########\n",
      "\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "\n",
      "##########  Monte Carlo 18 / 100 => okay   ##########\n",
      "\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "\n",
      "##########  Monte Carlo 19 / 100 => okay   ##########\n",
      "\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "\n",
      "##########  Monte Carlo 20 / 100 => okay   ##########\n",
      "\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "\n",
      "##########  Monte Carlo 21 / 100 => okay   ##########\n",
      "\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "\n",
      "##########  Monte Carlo 22 / 100 => okay   ##########\n",
      "\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "\n",
      "##########  Monte Carlo 23 / 100 => okay   ##########\n",
      "\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "\n",
      "##########  Monte Carlo 24 / 100 => okay   ##########\n",
      "\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "\n",
      "##########  Monte Carlo 25 / 100 => okay   ##########\n",
      "\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "\n",
      "##########  Monte Carlo 26 / 100 => okay   ##########\n",
      "\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "\n",
      "##########  Monte Carlo 27 / 100 => okay   ##########\n",
      "\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "\n",
      "##########  Monte Carlo 28 / 100 => okay   ##########\n",
      "\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "\n",
      "##########  Monte Carlo 29 / 100 => okay   ##########\n",
      "\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "\n",
      "##########  Monte Carlo 30 / 100 => okay   ##########\n",
      "\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "\n",
      "##########  Monte Carlo 31 / 100 => okay   ##########\n",
      "\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "\n",
      "##########  Monte Carlo 32 / 100 => okay   ##########\n",
      "\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "\n",
      "##########  Monte Carlo 33 / 100 => okay   ##########\n",
      "\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "\n",
      "##########  Monte Carlo 34 / 100 => okay   ##########\n",
      "\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "\n",
      "##########  Monte Carlo 35 / 100 => okay   ##########\n",
      "\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "\n",
      "##########  Monte Carlo 36 / 100 => okay   ##########\n",
      "\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "\n",
      "##########  Monte Carlo 37 / 100 => okay   ##########\n",
      "\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "\n",
      "##########  Monte Carlo 38 / 100 => okay   ##########\n",
      "\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "\n",
      "##########  Monte Carlo 39 / 100 => okay   ##########\n",
      "\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "\n",
      "##########  Monte Carlo 40 / 100 => okay   ##########\n",
      "\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "\n",
      "##########  Monte Carlo 41 / 100 => okay   ##########\n",
      "\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "\n",
      "##########  Monte Carlo 42 / 100 => okay   ##########\n",
      "\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "\n",
      "##########  Monte Carlo 43 / 100 => okay   ##########\n",
      "\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "\n",
      "##########  Monte Carlo 44 / 100 => okay   ##########\n",
      "\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "\n",
      "##########  Monte Carlo 45 / 100 => okay   ##########\n",
      "\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "\n",
      "##########  Monte Carlo 46 / 100 => okay   ##########\n",
      "\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "\n",
      "##########  Monte Carlo 47 / 100 => okay   ##########\n",
      "\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "\n",
      "##########  Monte Carlo 48 / 100 => okay   ##########\n",
      "\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "\n",
      "##########  Monte Carlo 49 / 100 => okay   ##########\n",
      "\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "\n",
      "##########  Monte Carlo 50 / 100 => okay   ##########\n",
      "\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "\n",
      "##########  Monte Carlo 51 / 100 => okay   ##########\n",
      "\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "\n",
      "##########  Monte Carlo 52 / 100 => okay   ##########\n",
      "\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "\n",
      "##########  Monte Carlo 53 / 100 => okay   ##########\n",
      "\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "\n",
      "##########  Monte Carlo 54 / 100 => okay   ##########\n",
      "\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "\n",
      "##########  Monte Carlo 55 / 100 => okay   ##########\n",
      "\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "\n",
      "##########  Monte Carlo 56 / 100 => okay   ##########\n",
      "\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "\n",
      "##########  Monte Carlo 57 / 100 => okay   ##########\n",
      "\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "\n",
      "##########  Monte Carlo 58 / 100 => okay   ##########\n",
      "\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "\n",
      "##########  Monte Carlo 59 / 100 => okay   ##########\n",
      "\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "\n",
      "##########  Monte Carlo 60 / 100 => okay   ##########\n",
      "\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "\n",
      "##########  Monte Carlo 61 / 100 => okay   ##########\n",
      "\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "\n",
      "##########  Monte Carlo 62 / 100 => okay   ##########\n",
      "\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "\n",
      "##########  Monte Carlo 63 / 100 => okay   ##########\n",
      "\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "\n",
      "##########  Monte Carlo 64 / 100 => okay   ##########\n",
      "\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "\n",
      "##########  Monte Carlo 65 / 100 => okay   ##########\n",
      "\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "\n",
      "##########  Monte Carlo 66 / 100 => okay   ##########\n",
      "\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "\n",
      "##########  Monte Carlo 67 / 100 => okay   ##########\n",
      "\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "\n",
      "##########  Monte Carlo 68 / 100 => okay   ##########\n",
      "\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "\n",
      "##########  Monte Carlo 69 / 100 => okay   ##########\n",
      "\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "\n",
      "##########  Monte Carlo 70 / 100 => okay   ##########\n",
      "\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "\n",
      "##########  Monte Carlo 71 / 100 => okay   ##########\n",
      "\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "\n",
      "##########  Monte Carlo 72 / 100 => okay   ##########\n",
      "\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "\n",
      "##########  Monte Carlo 73 / 100 => okay   ##########\n",
      "\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "\n",
      "##########  Monte Carlo 74 / 100 => okay   ##########\n",
      "\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "\n",
      "##########  Monte Carlo 75 / 100 => okay   ##########\n",
      "\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "\n",
      "##########  Monte Carlo 76 / 100 => okay   ##########\n",
      "\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "\n",
      "##########  Monte Carlo 77 / 100 => okay   ##########\n",
      "\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "\n",
      "##########  Monte Carlo 78 / 100 => okay   ##########\n",
      "\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "\n",
      "##########  Monte Carlo 79 / 100 => okay   ##########\n",
      "\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "\n",
      "##########  Monte Carlo 80 / 100 => okay   ##########\n",
      "\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "\n",
      "##########  Monte Carlo 81 / 100 => okay   ##########\n",
      "\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "\n",
      "##########  Monte Carlo 82 / 100 => okay   ##########\n",
      "\n",
      "\n",
      "##########  Monte Carlo 83 / 100 => okay   ##########\n",
      "\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "\n",
      "##########  Monte Carlo 84 / 100 => okay   ##########\n",
      "\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "\n",
      "##########  Monte Carlo 85 / 100 => okay   ##########\n",
      "\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "\n",
      "##########  Monte Carlo 86 / 100 => okay   ##########\n",
      "\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "\n",
      "##########  Monte Carlo 87 / 100 => okay   ##########\n",
      "\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "\n",
      "##########  Monte Carlo 88 / 100 => okay   ##########\n",
      "\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "\n",
      "##########  Monte Carlo 89 / 100 => okay   ##########\n",
      "\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "\n",
      "##########  Monte Carlo 90 / 100 => okay   ##########\n",
      "\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "\n",
      "##########  Monte Carlo 91 / 100 => okay   ##########\n",
      "\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "\n",
      "##########  Monte Carlo 92 / 100 => okay   ##########\n",
      "\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "\n",
      "##########  Monte Carlo 93 / 100 => okay   ##########\n",
      "\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "\n",
      "##########  Monte Carlo 94 / 100 => okay   ##########\n",
      "\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "\n",
      "##########  Monte Carlo 95 / 100 => okay   ##########\n",
      "\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "\n",
      "##########  Monte Carlo 96 / 100 => okay   ##########\n",
      "\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "\n",
      "##########  Monte Carlo 97 / 100 => okay   ##########\n",
      "\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "\n",
      "##########  Monte Carlo 98 / 100 => okay   ##########\n",
      "\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "\n",
      "##########  Monte Carlo 99 / 100 => okay   ##########\n",
      "\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "\n",
      "##########  Monte Carlo 100 / 100 => okay   ##########\n",
      "\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "  bayspar Warning. search_tol may be too small. try a larger number + 10\n",
      "\n",
      ">>    ID from 0 to 1296: field is ocn_sur_temp\n",
      ">>               Xb mean is 34.14978790, std is 2.70977869, var is 7.34290056\n",
      ">>      Recon 0. Xa mean is 34.34038162, std is 0.97226858, var is 0.94530618\n",
      "\n",
      ">>    ID from 1296 to 2592: field is atm_temp\n",
      ">>               Xb mean is 31.44051170, std is 3.14972261, var is 9.92075253\n",
      ">>      Recon 0. Xa mean is 31.66066933, std is 1.12927341, var is 1.27525842\n",
      "\n",
      ">>    ID from 2592 to 3888: field is atm_pCO2\n",
      ">>               Xb mean is 1945.77429199, std is 480.85505512, var is 231221.58403803\n",
      ">>      Recon 0. Xa mean is 1972.54321289, std is 172.66026306, var is 29811.56835938\n",
      "\n",
      ">>    ID from 3888 to 5184: field is ocn_sur_sal\n",
      ">>               Xb mean is 33.69539261, std is 0.26982926, var is 0.07280783\n",
      ">>      Recon 0. Xa mean is 33.68947220, std is 0.10282134, var is 0.01057223\n",
      "\n",
      ">>    ID from 5184 to 6480: field is misc_pH\n",
      ">>               Xb mean is 7.37344885, std is 0.13349122, var is 0.01781990\n",
      ">>      Recon 0. Xa mean is 7.45184326, std is 0.05731646, var is 0.00328518\n",
      "\n",
      ">>    ID from 6480 to 7776: field is carb_sur_ohm_cal\n",
      ">>               Xb mean is 3.19574952, std is 1.21246111, var is 1.47006194\n",
      ">>      Recon 0. Xa mean is 3.82881570, std is 0.63210016, var is 0.39955065\n",
      "\n",
      ">>    ID from 7776 to 9072: field is carb_sur_ohm_arg\n",
      ">>               Xb mean is 1.59029484, std is 0.60065628, var is 0.36078796\n",
      ">>      Recon 0. Xa mean is 1.90873921, std is 0.31451771, var is 0.09892138\n",
      "\n",
      ">>    ID from 9072 to 10368: field is ocn_ben_temp\n",
      ">>               Xb mean is 15.15946770, std is 2.06973846, var is 4.28381728\n",
      ">>      Recon 0. Xa mean is 15.46140194, std is 1.22251832, var is 1.49455094\n",
      "\n",
      ">>    ID from 10368 to 11664: field is ocn_sur_ALK\n",
      ">>               Xb mean is 0.00179770, std is 0.00029902, var is 0.00000009\n",
      ">>      Recon 0. Xa mean is 0.00209920, std is 0.00018412, var is 0.00000003\n",
      "\n",
      ">>    ID from 11664 to 12960: field is sed_CaCO3\n",
      ">>               Xb mean is 18.09581566, std is 28.28526131, var is 800.05600730\n",
      ">>      Recon 0. Xa mean is 15.35647488, std is 5.74948978, var is 33.05663300\n",
      "\n",
      ">>    ID from 12960 to 14256: field is carb_ben_ohm_cal\n",
      ">>               Xb mean is 0.63249749, std is 0.28290039, var is 0.08003263\n",
      ">>      Recon 0. Xa mean is 0.69795787, std is 0.10257214, var is 0.01052104\n",
      "\n",
      ">>    ID from 14256 to 15552: field is carb_ben_ohm_arg\n",
      ">>               Xb mean is 0.35208413, std is 0.14998174, var is 0.02249452\n",
      ">>      Recon 0. Xa mean is 0.38852242, std is 0.05302934, var is 0.00281211\n",
      "\n",
      "########## All Done ##########\n",
      "\n",
      "This web page saved as DeepDA_allMC.html in the working directory : /Users/mingsongli/Dropbox/PETM/posterior/petmproxy3slices_v0.1.5shift.csv_petm38_v0.1.5shift_20240728_All.bays_MC100_pHcor_frac0.99_Ca75_biogem_SatuBen_setLim/\n"
     ]
    }
   ],
   "source": [
    "'''\n",
    "DeepDA - Data assimilation framework for deep time paleoclimate projects\n",
    "\n",
    "Prior: cgenie simulation ensemble\n",
    "Proxy: deepda\n",
    "PSM: baysian psm\n",
    "\n",
    "By :  Mingsong Li\n",
    "      Peking University @ Penn State\n",
    "      msli@pku.edu.cn\n",
    "      limingsonglms@gmail.com\n",
    "     \n",
    "Date: Feb 26, 2020\n",
    "\n",
    "Updated\n",
    "\n",
    "    Mar. 3. 2020\n",
    "        partly clean the code; add MC for local_rad, withheld_rate, and scaled Rg\n",
    "    June 2020\n",
    "        include d18O of CESM by Zhu et al., 2019 Sci Adv\n",
    "    August 2020\n",
    "        Two options for the proxy order: all random & use the given list\n",
    "    Oct. 12, 2020\n",
    "        Add multi_seed for Monte Carlo simulations\n",
    "    Nov. 2, 2020\n",
    "        Add DeepMIP PSMs\n",
    "        Add log_level\n",
    "        Note:\n",
    "        if Mg/Ca proxy is included, may need to run\n",
    "            correct_cgenie_carb_ohm_cal.ipynb for the estimation of carb_ohm_cal_ben field\n",
    "        if d13C proxy is included, may need to run\n",
    "            correct_cgenie_sed_caco3_13c.ipynb for the correction of d13C because cgenie simulations may have nan values\n",
    "    May 6, 2021\n",
    "        d18O correction -1.42 per mil per pH (Zeebe, 2001)\n",
    "    June 29, 2021\n",
    "        pH correction using the bayfox 1.1\n",
    "    July 15, 2021\n",
    "        Rscale_style: tuple or use dict  \n",
    "        v0.10\n",
    "    Mar. 30, 2022\n",
    "        Add comments\n",
    "    Jul. 23, 2022\n",
    "        Update proxy to v0.1.2: correct location for Site U1443\n",
    "        Update prior to ML.petm031: correct LHS code\n",
    "    Dec. 26, 2023\n",
    "        Update proxy to v0.1.3: Add ZS site in Tibet\n",
    "        Update prior to ML.petm035, 36 for pre and peak PETM, \n",
    "        consider 3 uncertainties: pCO2, ALK, rainratio\n",
    "    July 4, 2024\n",
    "        Update prior reading, if CaCO3_d13C less than -100, revise to NaN\n",
    "        \n",
    "'''\n",
    "### ===================================\n",
    "### Import packages\n",
    "from DeepDA_lib import LMR_DA\n",
    "from DeepDA_lib import modules_nc\n",
    "from DeepDA_lib import DeepDA_tools\n",
    "from DeepDA_lib import DeepDA_psm\n",
    "import h5py\n",
    "import time\n",
    "import yaml\n",
    "import numpy as np\n",
    "import random\n",
    "import pandas\n",
    "import os\n",
    "import shutil\n",
    "from netCDF4 import Dataset\n",
    "import numpy.ma as ma\n",
    "import numpy.matlib as mat\n",
    "from sys import platform as sys_pf\n",
    "import matplotlib.pyplot as plt\n",
    "if sys_pf == 'darwin':\n",
    "    import matplotlib\n",
    "    matplotlib.use(\"TkAgg\")\n",
    "    import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "#from mpl_toolkits.basemap import Basemap, shiftgrid, cm\n",
    "\n",
    "try:\n",
    "    import bayspar\n",
    "except ImportError as e2:\n",
    "    print('Warning:', e2)\n",
    "try:\n",
    "    import bayfox\n",
    "except ImportError as e3:\n",
    "    print('Warning:', e3)\n",
    "try:\n",
    "    import baymag\n",
    "    # fetch modern pH and omega from baymag\n",
    "    import baymag.omgph\n",
    "except ImportError as e4:\n",
    "    print('Warning:', e4)\n",
    "### ===================================\n",
    "### Read config file\n",
    "\n",
    "config_name = \"DeepDA_config.yml\"\n",
    "\n",
    "f = open(config_name, 'r')\n",
    "yml_dict = yaml.load(f, Loader=yaml.FullLoader)\n",
    "f.close()\n",
    "log_level = yml_dict['log_level']\n",
    "if log_level > 1:\n",
    "    print('>>  Import packages...  => Okay')\n",
    "\n",
    "### ===================================\n",
    "# user-defined\n",
    "t = -1  # last time slice, for cGENIE\n",
    "#t = -2  # last time slice, for cGENIE  ML.petm041 only\n",
    "k = 0   # surface layer, for SST\n",
    "\n",
    "kcov_saving = 0 # save covariance??? 0=no, 1 = yes\n",
    "\n",
    "### ===================================\n",
    "### ===================================\n",
    "# read config.yml settings\n",
    "if log_level > 1:\n",
    "    print('')\n",
    "    print(' ########## Load yml config file ########## ')\n",
    "    print('')\n",
    "########## Proxy + PSM #########\n",
    "\n",
    "MCn= yml_dict['MonteCarlo']['number']\n",
    "multi_seed = yml_dict['MonteCarlo']['multi_seed']\n",
    "dir_proxy         = yml_dict['core']['proxy_dir']\n",
    "dir_proxy_data    = dir_proxy +'/'+ yml_dict['proxies'][yml_dict['proxies']['use_from'][0]]['dbversion']\n",
    "dir_proxy_save_dir= yml_dict['core']['wrkdir'] + '/'\n",
    "dir_proxy_save    = yml_dict['proxies'][yml_dict['proxies']['use_from'][0]]['dbversion']\n",
    "proxy_psm_type    = yml_dict['proxies'][yml_dict['proxies']['use_from'][0]]['proxy_psm_type']\n",
    "proxy_assim2      = yml_dict['proxies'][yml_dict['proxies']['use_from'][0]]['proxy_assim2']\n",
    "proxy_order       = yml_dict['proxies'][yml_dict['proxies']['use_from'][0]]['proxy_order']\n",
    "proxy_err_eval    = yml_dict['proxies'][yml_dict['proxies']['use_from'][0]]['proxy_err_eval']\n",
    "proxy_blacklist   = yml_dict['proxies'][yml_dict['proxies']['use_from'][0]]['proxy_blacklist']\n",
    "proxy_list = [item for item in proxy_order if item not in proxy_blacklist]\n",
    "psm_d18osw_adjust = yml_dict['psm']['bayesreg_d18o_pooled']['psm_d18osw_adjust']\n",
    "d18osw_local_choice = yml_dict['psm']['bayesreg_d18o_pooled']['d18osw_local_choice']\n",
    "d18osw_icesm_pco2 = yml_dict['psm']['bayesreg_d18o_pooled']['d18osw_icesm_pco2']\n",
    "\n",
    "try:\n",
    "    d18o_phcor_water = yml_dict['psm']['bayesreg_d18o_pooled']['d18o_phcor_water']\n",
    "except:\n",
    "    d18o_phcor_water = 'deeptime'\n",
    "    \n",
    "proxy_qc          = yml_dict['proxies']['proxy_qc']\n",
    "lon_label = yml_dict['proxies'][yml_dict['proxies']['use_from'][0]]['lon_label']\n",
    "lat_label = yml_dict['proxies'][yml_dict['proxies']['use_from'][0]]['lat_label']\n",
    "\n",
    "prior_source = yml_dict['prior']['prior_source'] #\n",
    "prior_state_variable = yml_dict['prior'][prior_source]['state_variable']  # note: ['2d': xxx; '3d': xxx]\n",
    "dum_lon_offset = yml_dict['prior'][prior_source]['dum_lon_offset'] # longitude offset\n",
    "dir_prior = yml_dict['core']['prior_dir']\n",
    "dir_prior_full = os.listdir(dir_prior)\n",
    "dir_prior_full = sorted(dir_prior_full)\n",
    "prior_len = len(dir_prior_full)\n",
    "\n",
    "nexp = yml_dict['core']['nexp']\n",
    "data_period_id    = yml_dict['proxies'][yml_dict['proxies']['use_from'][0]]['data_period_id']\n",
    "data_period_idstd = yml_dict['proxies'][yml_dict['proxies']['use_from'][0]]['data_period_idstd']\n",
    "recon_period = yml_dict['core']['recon_period']\n",
    "recon_timescale = yml_dict['core']['recon_timescale_interval']\n",
    "recon_period_full = np.arange(recon_period[0],recon_period[1]+1,recon_timescale)\n",
    "recon_period_len = recon_period_full.shape[0]\n",
    "geologic_age = yml_dict['core']['geologic_age']\n",
    "limit_hard_keys = list(yml_dict['prior'][prior_source]['limit_hard'].keys())\n",
    "\n",
    "# save config\n",
    "config_save_name = dir_proxy_save_dir + dir_proxy_save + nexp + '.yml'\n",
    "configos = 'cp ' + config_name + ' ' +  config_save_name\n",
    "os.system(configos)\n",
    "if log_level > 0:\n",
    "    print('')\n",
    "    print('The config (yml) file saved : ')\n",
    "    print('')\n",
    "    print(config_save_name)\n",
    "    print('')\n",
    "if log_level > 1:\n",
    "    print('Set limit for {}'.format(limit_hard_keys))\n",
    "    print(' directory of prior: {}'.format(dir_prior_full))\n",
    "\n",
    "nens = yml_dict['core']['nens']\n",
    "save_ens_full = yml_dict['core']['save_ens_full']\n",
    "save_mc_full = yml_dict['core']['save_mc_full']\n",
    "\n",
    "proxy_err_eval = yml_dict['proxies'][yml_dict['proxies']['use_from'][0]]['proxy_err_eval']\n",
    "# glassy d18O blacklist\n",
    "proxy_d18o_glassy  = yml_dict['proxies']['proxy_d18o_glassy']\n",
    "proxy_assim3 = yml_dict['proxies'][yml_dict['proxies']['use_from'][0]]['proxy_assim3']\n",
    "data_glassy_label_blacklist = proxy_assim3['Marine sediments_d18o_pooled_glassy']\n",
    "# bayspar\n",
    "search_tol_i = yml_dict['psm']['bayesreg_tex86']['search_tol']\n",
    "nens_i = yml_dict['psm']['bayesreg_tex86']['nens']\n",
    "\n",
    "# ========= dataset for plot =========\n",
    "cGENIEGrid = yml_dict['core']['proj_dir'] + '/data_misc/cGENIEGrid.csv'\n",
    "cGENIEGrid = pandas.read_csv(cGENIEGrid)\n",
    "cGENIEGridB_lat36 = cGENIEGrid['lat']\n",
    "cGENIEGridB_lon36 = cGENIEGrid['lon']\n",
    "cGENIEGrid = cGENIEGrid.to_numpy()\n",
    "if log_level > 2:\n",
    "    print('>>  Load dataset for plot => Okay')\n",
    "\n",
    "# ========= Monte Carlo =========\n",
    "local_rad_list = yml_dict['core']['local_rad_list'] #\n",
    "locRadn= len(local_rad_list)\n",
    "local_rad_list = np.asarray(local_rad_list)\n",
    "\n",
    "# proxy fraction\n",
    "proxy_frac_list   = yml_dict['proxies']['proxy_frac']\n",
    "proxy_fracn = len(proxy_frac_list)\n",
    "proxy_frac_list = np.asarray(proxy_frac_list)\n",
    "proxy_order_type = yml_dict['proxies'][yml_dict['proxies']['use_from'][0]]['proxy_order_type']\n",
    "\n",
    "# R scale\n",
    "Rscale_style = yml_dict['core']['Rscale_style']\n",
    "if  Rscale_style==1:\n",
    "    Rscale_list = yml_dict['core']['Rscale']\n",
    "    Rscalen = len(Rscale_list)\n",
    "    Rscale_list = np.asarray(Rscale_list)\n",
    "elif Rscale_style ==2:\n",
    "    Rscalen = 1\n",
    "    Rscale_d18o = yml_dict['core']['Rscale_dict']['Rscale_d18o']\n",
    "    Rscale_mgca = yml_dict['core']['Rscale_dict']['Rscale_mgca']\n",
    "    Rscale_tex = yml_dict['core']['Rscale_dict']['Rscale_tex']\n",
    "    Rscale_caco3 = yml_dict['core']['Rscale_dict']['Rscale_caco3']\n",
    "\n",
    "# output\n",
    "if log_level > 1:\n",
    "    print('>>  Prior member size: {}'.format(prior_len))\n",
    "    print('>>  Recon_period {} - {}. '.format(recon_period[0], recon_period[1]))\n",
    "    if log_level > 1:\n",
    "        print('      List: {}'.format(recon_period_full))\n",
    "    print('>>  Proxy error evaluation: {}'.format(proxy_err_eval))\n",
    "    if log_level > 2:\n",
    "        print('>>  Proxy full list:')\n",
    "        print('      {}'.format(proxy_order))\n",
    "        print('>>  Proxy blacklist:')\n",
    "        print('      {}'.format(proxy_blacklist))\n",
    "    print('>>  Proxy to be assimilated (some may not exist)')\n",
    "    print('      {}'.format(proxy_list))\n",
    "    print('>>  Proxy quality control selection: {}'.format(proxy_qc))\n",
    "\n",
    "data_psm_d18o_find = 0\n",
    "data_psm_mgca_find = 0\n",
    "if 'Marine sediments_mgca_pooled_bcp' in proxy_list or 'Marine sediments_mgca_pooled_red' in proxy_list:\n",
    "    data_psm_mgca_find = 1\n",
    "    if log_level > 2:\n",
    "        print('>>    Mg/Ca proxy found ')\n",
    "        \n",
    "# debug:   force the DA to calculate Xb_sal, Xb_ph and Xb_sal, etc.\n",
    "# data_psm_mgca_find = 1\n",
    "# debug:   end\n",
    "\n",
    "if 'Marine sediments_d18o_pooled' in proxy_list:\n",
    "    data_psm_d18o_find = 1\n",
    "    if log_level > 2:\n",
    "        print('>>    d18O proxy found ')\n",
    "    \n",
    "\n",
    "########## Prior #########\n",
    "\n",
    "if log_level > 1:\n",
    "    print('')\n",
    "    print('########## Read prior ######### ')\n",
    "    print('')\n",
    "\n",
    "# prior variable list, for saving outputs\n",
    "prior_variable_dict = []  # variable list\n",
    "prior_nc_file_list = []  # nc file list\n",
    "prior_variable_dict_3d = []  # variable list\n",
    "prior_nc_file_list_3d = []  # nc file list\n",
    "\n",
    "for key, value in prior_state_variable.items():\n",
    "    nc_keyvalue = prior_state_variable[key]['ncname']  # note: 2d or 3d dict\n",
    "    \n",
    "    for key1, value1 in nc_keyvalue.items():\n",
    "        for i in range(len(prior_state_variable[key][value1])):\n",
    "            if key in ['2d']:\n",
    "                prior_variable_dict.append(prior_state_variable[key][value1][i])\n",
    "                prior_nc_file_list.append(key1+'/'+value1+'.nc')\n",
    "            elif key in ['3d']:\n",
    "                prior_variable_dict_3d.append(prior_state_variable[key][value1][i])\n",
    "                prior_nc_file_list_3d.append(key1+'/'+value1+'.nc')\n",
    "\n",
    "# prepare variable list for Xb\n",
    "\n",
    "prior_variable2d_len = len(prior_variable_dict)\n",
    "prior_variable3d_len = len(prior_variable_dict_3d)\n",
    "if log_level > 2:\n",
    "    print('>>  Number of 2d prior variables is: {}.'.format(prior_variable2d_len))\n",
    "if prior_variable2d_len>0:\n",
    "    if log_level > 2:\n",
    "        print('      List:')\n",
    "        for i in range(prior_variable2d_len):\n",
    "            print('        {}/{}'.format(prior_nc_file_list[i], prior_variable_dict[i]))\n",
    "if log_level > 2:\n",
    "    print('>>  Number of 3d prior variables is: {}'.format(prior_variable3d_len))\n",
    "    if prior_variable3d_len>0:\n",
    "        print('      List:')\n",
    "        for i in range(prior_variable3d_len):\n",
    "            print('        {}/{}'.format(prior_nc_file_list_3d[i], prior_variable_dict_3d[i]))\n",
    "   \n",
    "\n",
    "# If there is no field in the model, convert model unit to proxy unit\n",
    "\n",
    "if log_level > 2:\n",
    "    print('>>  Reading prior state variables')\n",
    "    \n",
    "# read first variable data, first time slice, to get the shape of prior grid\n",
    "try:\n",
    "    x1 = Dataset(dir_prior+'/'+dir_prior_full[0]+'/'+ prior_nc_file_list_3d[0]).variables[prior_variable_dict_3d[0]][0,:,:,:]\n",
    "    dum_dmax = x1.shape[0] # depth\n",
    "    dum_imax = x1.shape[1]  # lon\n",
    "    dum_jmax = x1.shape[2]  # lat\n",
    "except:\n",
    "    try:\n",
    "        x0 = Dataset(dir_prior+'/'+dir_prior_full[0]+'/'+ prior_nc_file_list[0]).variables[prior_variable_dict[0]][0,:,:]\n",
    "        dum_imax = 36 #x1.shape[0]  # lon\n",
    "        dum_jmax = 36 #x1.shape[1]  # lat\n",
    "        dum_dmax = 16\n",
    "    except:\n",
    "        dum_dmax = 16\n",
    "        dum_imax = 36\n",
    "        dum_jmax = 36\n",
    "        \n",
    "# prepare 2d Xb for lon-lat state \n",
    "dum_ijmax = dum_imax*dum_jmax  # lonn * latn\n",
    "\n",
    "if log_level > 3:\n",
    "    print('>>  Shape of dum_dmax {}, dum_imax {}, dum_jmax {}, dum_ijmax {}'.format(dum_dmax,dum_imax,dum_jmax,dum_ijmax))\n",
    "\n",
    "# save units of each variable\n",
    "prior_variable_units = list()\n",
    "prior_variable_units_init = 0\n",
    "\n",
    "# define nan matrix for storing 2d and 3d variables\n",
    "# Xb for 2D variables\n",
    "if prior_variable2d_len>0:\n",
    "    Xb_shape = (prior_variable2d_len*dum_jmax*dum_imax, prior_len)  # lonn * latn * varn\n",
    "    Xb   = np.full(Xb_shape,np.nan)\n",
    "\n",
    "# prepare 3D version of Xb\n",
    "if prior_variable3d_len > 0:\n",
    "    Xb3d_shape = (prior_variable3d_len*dum_dmax*dum_jmax*dum_imax, prior_len)  # lonn * latn * varn\n",
    "    Xb3d = np.full(Xb3d_shape,np.nan)\n",
    "    # read prior and save Xb\n",
    "    #Xb = np.full((dum_ijmax, prior_len),np.nan)\n",
    "    \n",
    "if log_level > 2:\n",
    "    print('>>  Reading prior ...')\n",
    "    \n",
    "if data_psm_d18o_find == 1:\n",
    "    if log_level > 3:\n",
    "        print('>>  Prepare d18O related state variable ...')\n",
    "    Xb_ph        = np.full(Xb_shape,np.nan)\n",
    "    Xb_sal       = np.full(Xb_shape,np.nan)\n",
    "    \n",
    "if data_psm_mgca_find == 1:\n",
    "    if log_level > 3:\n",
    "        print('>>  Prepare Mg/Ca related state variable ...')\n",
    "    # for Mg/Ca SST proxy salinity, ph, omega\n",
    "    Xb_sal       = np.full(Xb_shape,np.nan)\n",
    "    Xb_ph        = np.full(Xb_shape,np.nan)\n",
    "    Xb_omega     = np.full(Xb_shape,np.nan)\n",
    "    spp = 'all'\n",
    "    # ``1`` for reductive, ``0`` for BCP (Barker).\n",
    "    cleaningr = np.tile(np.array([1]),prior_len).reshape((prior_len,1))\n",
    "    cleaningb = np.tile(np.array([0]),prior_len).reshape((prior_len,1))\n",
    "\n",
    "# read units of each variable from prior and save as prior_variable_units\n",
    "if prior_variable2d_len > 0:\n",
    "    for j in range(prior_variable2d_len):\n",
    "        name_nc_2d = dir_prior+'/'+dir_prior_full[0]+'/'+ prior_nc_file_list[j]\n",
    "        nc_field = prior_variable_dict[j]\n",
    "        try:\n",
    "            unit_j = Dataset(name_nc_2d).variables[nc_field].units\n",
    "        except:\n",
    "            unit_j ='unit'\n",
    "        prior_variable_units.append((unit_j))\n",
    "if prior_variable3d_len > 0:\n",
    "    for j in range(prior_variable3d_len):\n",
    "        name_nc_3d = dir_prior+'/'+dir_prior_full[0]+'/'+ prior_nc_file_list_3d[j]\n",
    "        nc_field = prior_variable_dict_3d[j]\n",
    "        try:\n",
    "            try:\n",
    "                unit_j = Dataset(name_nc_3d).variables[nc_field].units\n",
    "            except:\n",
    "                unit_j ='unit'\n",
    "            prior_variable_units.append((unit_j))\n",
    "        except:\n",
    "            prior_variable_units.append((''))\n",
    "\n",
    "data_psm_mgca_find_count = 0\n",
    "\n",
    "# loop for each member of an ensemble\n",
    "\n",
    "for i in range(prior_len):\n",
    "    \n",
    "    # loop for each variable of each member\n",
    "    # if Mg/Ca data is to be assimilated, then\n",
    "    if data_psm_mgca_find == 1:\n",
    "        water_saturation    = yml_dict['psm']['bayesreg_mgca_pooled_red']['water_saturation']        \n",
    "        psm_required_nc     = yml_dict['psm']['bayesreg_mgca_pooled_red']['psm_required_nc']\n",
    "        psm_required_nc_ohm = yml_dict['psm']['bayesreg_mgca_pooled_red']['psm_required_nc_ohm']\n",
    "        psm_required_nc_mg  = yml_dict['psm']['bayesreg_mgca_pooled_red']['psm_required_nc_mg']\n",
    "        psm_baymag_ln       =  yml_dict['psm']['bayesreg_mgca_pooled_red']['psm_baymag_ln']\n",
    "        \n",
    "        name_nc_2d = dir_prior+'/'+dir_prior_full[i]+psm_required_nc\n",
    "        name_nc_2d_mgca = dir_prior+'/'+dir_prior_full[i]+psm_required_nc_mg\n",
    "        psm_required_nc_ohm = dir_prior+'/'+dir_prior_full[i]+psm_required_nc_ohm\n",
    "        \n",
    "        try:\n",
    "            x00 = Dataset(name_nc_2d).variables['ocn_sur_sal'][t,:,:] # time-lat-lon\n",
    "        except:\n",
    "            print(' prior error: no ocn_sur_sal field; run /utils/correct_cgenie_biogem_2d.ipynb first')\n",
    "            \n",
    "        x01 = Dataset(name_nc_2d).variables['misc_pH'][t,:,:] # time-lat-lon | core top pH\n",
    "        \n",
    "        if water_saturation in ['bottom']:\n",
    "            if log_level > 2:\n",
    "                if data_psm_mgca_find_count < 1:\n",
    "                    print('>>  Water_saturation_field :  bottom')\n",
    "            water_saturation_field = yml_dict['psm']['bayesreg_mgca_pooled_red']['water_saturation_field']\n",
    "            x02 = Dataset(psm_required_nc_ohm).variables[water_saturation_field][-1,:,:]   # only 1 time slice \n",
    "             \n",
    "        if water_saturation in ['fixed']:\n",
    "            water_saturation_value = yml_dict['psm']['bayesreg_mgca_pooled_red']['water_saturation_value']\n",
    "            if log_level > 2:\n",
    "                if data_psm_mgca_find_count < 1:\n",
    "                    print('>>  Water_saturation_field :  fixed {}'.format(water_saturation_value))\n",
    "                    \n",
    "            x02 = np.full((dum_jmax,dum_imax),water_saturation_value)\n",
    "            \n",
    "        if water_saturation in ['modern']:\n",
    "            if log_level > 2:\n",
    "                print('>>  Water_saturation_field :  modern from baymag')\n",
    "                \n",
    "        data_psm_mgca_find_count += 1\n",
    "        \n",
    "    # no Mg/Ca data is to be assimilated\n",
    "    else:\n",
    "        if data_psm_d18o_find == 1:\n",
    "            psm_required_nc = yml_dict['psm']['bayesreg_d18o_pooled']['psm_required_nc']\n",
    "            name_nc_2d = dir_prior+'/'+dir_prior_full[i]+psm_required_nc\n",
    "            try:\n",
    "                x00 = Dataset(name_nc_2d).variables['ocn_sur_sal'][t,:,:] # time-lat-lon\n",
    "            except:\n",
    "                print(' prior error: no ocn_sur_sal field; run /utils/correct_cgenie_biogem_2d.ipynb first')\n",
    "                \n",
    "            x01 = Dataset(name_nc_2d).variables['misc_pH'][t,:,:] # time-lat-lon | core top pH\n",
    "            \n",
    "    if log_level > 2:\n",
    "        print(' reading prior: member {}. {}'.format(i, name_nc_2d))\n",
    "        \n",
    "    # if more than 1 variable is to be assimilated, then\n",
    "    if prior_variable2d_len > 0:\n",
    "        for j in range(prior_variable2d_len):\n",
    "            # full directory of netcdf file\n",
    "            name_nc_2d = dir_prior+'/'+dir_prior_full[i]+'/'+ prior_nc_file_list[j]\n",
    "            j0 = dum_ijmax * j\n",
    "            j1 = dum_ijmax * (j+1)\n",
    "            nc_field = prior_variable_dict[j]\n",
    "            if log_level > 4:\n",
    "                print('>>    name nc 2d is {}, nc field is {}'.format(name_nc_2d, nc_field))\n",
    "                \n",
    "            x = Dataset(name_nc_2d).variables[nc_field][t,:,:]  # time-lat-lon\n",
    "            \n",
    "            # pCO2 from 1 to ppm\n",
    "            if nc_field in ['atm_pCO2']:\n",
    "                x = x * 1.0e+06\n",
    "\n",
    "            # CaCO3_d13C ensure no bad values (-999999999 etc.) \n",
    "            \n",
    "            if nc_field in ['sed_CaCO3_13C']:\n",
    "                x[x<-100] = np.nan\n",
    "                if log_level > 1 and i == 1:\n",
    "                    print('>>    CaCO3_d13C: ensure (any values < -100) = np.nan ')\n",
    "                \n",
    "            # fburial_CaCO3: ensure no negative values\n",
    "            if nc_field in ['fburial_CaCO3']:\n",
    "                if log_level > 4:\n",
    "                    print('>>      min of fburial_CaCO3 raw')\n",
    "                    print(np.min(x))\n",
    "                x[x<0] = 0\n",
    "                if log_level > 4:\n",
    "                    print('>>      min of fburial_CaCO3 new')\n",
    "                    print(np.min(x))\n",
    "            Xb[j0:j1,i] = np.copy(x.reshape(dum_ijmax))  # var-lat-lon: Nx x 1\n",
    "            if data_psm_d18o_find == 1:\n",
    "                try:\n",
    "                    Xb_sal[j0:j1,i] = np.copy(x00.reshape(dum_ijmax)) # var-lat-lon: Nx x 1  | surface water salinity\n",
    "                    Xb_ph[j0:j1,i]  = np.copy(x01.reshape(dum_ijmax)) # var-lat-lon: Nx x 1\n",
    "                except:\n",
    "                    if i == 0:\n",
    "                        # warning one time\n",
    "                        if log_level > 1:\n",
    "                            print('>>  Warning: reading state variable error. ocn_sur_sal or misc_pH')\n",
    "            if data_psm_mgca_find == 1:\n",
    "                try:\n",
    "                    Xb_sal[j0:j1,i] = np.copy(x00.reshape(dum_ijmax)) # var-lat-lon: Nx x 1  | surface water salinity\n",
    "                    Xb_ph[j0:j1,i] = np.copy(x01.reshape(dum_ijmax)) # var-lat-lon: Nx x 1\n",
    "                    if water_saturation in ['bottom','fixed']:\n",
    "                        Xb_omega[j0:j1,i] = np.copy(x02.reshape(dum_ijmax)) # var-lat-lon: Nx x 1\n",
    "                except:\n",
    "                    if i == 0:\n",
    "                        # warning one time\n",
    "                        if log_level > 1:\n",
    "                            print('>>  Warning: reading state variable error. ocn_sur_sal, misc_pH, carb_ohm_cal')\n",
    "            # print the last one data\n",
    "            if log_level > 2:\n",
    "                if i > prior_len-2:\n",
    "                    print('    Last member: {}: {}: {}'.format(i, dir_prior_full[i], prior_variable_dict[j]))\n",
    "        \n",
    "        # mask for the nan values\n",
    "        Xb = np.ma.MaskedArray(Xb, Xb >= 9.9692e+36)\n",
    "        \n",
    "    # if 3d variables are used\n",
    "    if prior_variable3d_len > 0:\n",
    "        for k in range(prior_variable3d_len):\n",
    "            name_nc_3d = dir_prior+'/'+dir_prior_full[i]+'/'+ prior_nc_file_list_3d[k]\n",
    "            nc_field = prior_variable_dict_3d[k]\n",
    "            k0 = dum_ijmax*dum_dmax * k\n",
    "            k1 = dum_ijmax*dum_dmax * (k+1)\n",
    "            x = Dataset(name_nc_3d).variables[nc_field][t,:,:,:]  # time-depth-lat-lon\n",
    "            Xb3d[k0:k1,i] = np.copy(x.reshape(dum_dmax*dum_ijmax)) # var-depth-lat-lon\n",
    "        # mask for the nan values\n",
    "        Xb3d = np.ma.MaskedArray(Xb3d, Xb3d >= 9.9692e+36)\n",
    "\n",
    "if log_level > 1:\n",
    "    print('>>  Units of state variables {}: {}'.format(prior_variable_dict+prior_variable_dict_3d,prior_variable_units))\n",
    "\n",
    "# hard copy of 2D variable\n",
    "Xb_prior = np.copy(Xb)\n",
    "# hard copy of 3D variable\n",
    "if prior_variable3d_len > 0:\n",
    "    Xb3d_prior = np.copy(Xb3d)\n",
    "    \n",
    "# hard copy of variables that are useful for the d18O\n",
    "if data_psm_d18o_find == 1:\n",
    "    if log_level > 3:\n",
    "        print('>>  Prepare d18O related state variable ...')\n",
    "    Xb_sal_prior       = np.copy(Xb_sal)\n",
    "    Xb_ph_prior        = np.copy(Xb_ph)\n",
    "# hard copy of variables that are useful for the Mg/Ca\n",
    "if data_psm_mgca_find == 1:\n",
    "    if log_level > 3:\n",
    "        print('>>  Prepare Mg/Ca related state variable ...')\n",
    "    # for Mg/Ca SST proxy salinity, ph, omega\n",
    "    Xb_sal_prior       = np.copy(Xb_sal)\n",
    "    Xb_ph_prior        = np.copy(Xb_ph)\n",
    "    Xb_omega_prior     = np.copy(Xb_omega)\n",
    "    # ``1`` for reductive, ``0`` for BCP (Barker).\n",
    "    cleaningr_prior = np.copy(cleaningr)\n",
    "    cleaningb_prior = np.copy(cleaningb)\n",
    "if log_level > 1:\n",
    "    print('>>  Reading Prior => Okay')\n",
    "    print('')\n",
    "    print(' ########## Read proxies database ########## ')\n",
    "    print('')\n",
    "    \n",
    "\n",
    "### ===================================\n",
    "### read proxies database ###\n",
    "### ===================================\n",
    "\n",
    "proxies = pandas.read_csv(dir_proxy_data)\n",
    "if log_level > 1:\n",
    "    print('>>  All proxy: shape')\n",
    "    print(proxies.shape)\n",
    "\n",
    "proxies_len0 = len(proxies)\n",
    "\n",
    "if log_level > 3:\n",
    "    print('>>  All proxy: '.format(proxies))\n",
    "\n",
    "proxy_select_0 = 0\n",
    "\n",
    "### remove proxy data that is in the blacklist ###\n",
    "for j in range(proxies_len0):\n",
    "    # Read proxy type from the database\n",
    "    data_psm_type = proxies['Proxy'][j]\n",
    "    # initial default 0 : this proxy is not included\n",
    "    data_assimilate_i = 0\n",
    "    for jlist in range(len(proxy_list)):\n",
    "        if data_psm_type in proxy_assim2[proxy_list[jlist]]:\n",
    "            # find and save this proxy\n",
    "            data_assimilate_i = 1\n",
    "            \n",
    "    if data_assimilate_i == 1:\n",
    "        if log_level > 3:\n",
    "            print('>>    File {}, {} included'.format(proxies.loc[j,'File'], data_psm_type))\n",
    "        if proxy_select_0 == 0:\n",
    "            proxy_select0 = proxies.iloc[[j]]\n",
    "            proxy_select0 = proxy_select0.reset_index(drop=True) # reset_index, avoid index error\n",
    "            proxy_select_0 = 1\n",
    "        else:\n",
    "            proxy_select0 = proxy_select0.append(proxies.iloc[[j]], ignore_index=True)\n",
    "            \n",
    "proxies_select_len0 = len(proxy_select0)\n",
    "\n",
    "if log_level > 1:\n",
    "    print('>>  Proxy: selected proxy dataset number {}: those in blacklist removed!'.format(proxies_select_len0))\n",
    "\n",
    "### remove data of non-glassy fora if proxy_d18o_glassy is TRUE\n",
    "\n",
    "proxy_select_0 = 0\n",
    "if proxy_d18o_glassy:\n",
    "    for jj in range(proxies_select_len0):\n",
    "        data_glassy_label = proxy_select0['Glassy'][jj]\n",
    "        if data_glassy_label not in data_glassy_label_blacklist:\n",
    "            if proxy_select_0 == 0:\n",
    "                proxy_select = proxy_select0.iloc[[jj]]\n",
    "                proxy_select = proxy_select.reset_index(drop=True) # reset_index, avoid index error\n",
    "                proxy_select_0 = 1\n",
    "            else:\n",
    "                proxy_select = proxy_select.append(proxy_select0.iloc[[jj]], ignore_index=True)\n",
    "\n",
    "    if log_level > 3:\n",
    "        print(proxy_select)\n",
    "    proxies_select_len0 = len(proxy_select)\n",
    "    if log_level > 1:\n",
    "        print('>>  Proxy: selected proxy dataset number {}: those unknown/frosty removed!'.format(proxies_select_len0))\n",
    "else:\n",
    "    proxy_select = proxy_select0.copy()\n",
    "    \n",
    "#######     ########     #######     ########     #######     ########     #######     ######## \n",
    "#######                  OKAY, setting read, now let's DA                              ######## \n",
    "#######     ########     #######     ########     #######     ########     #######     ######## \n",
    "\n",
    "# DA for each loop using a user-defined localization radius\n",
    "for locRadi in range(locRadn):\n",
    "    locRad = local_rad_list[locRadi]\n",
    "    if log_level > 2:\n",
    "        print('')\n",
    "        print('>>  Localization id {} radius distance {} km'.format(locRadi, locRad))\n",
    "    if locRad is None:\n",
    "        locRadv = 0 # for filename only\n",
    "    else:\n",
    "        locRadv = locRad\n",
    "    \n",
    "    # DA for each loop using a user-defined proxy fraction\n",
    "    for proxy_fraci in range(proxy_fracn):\n",
    "        proxy_frac = proxy_frac_list[proxy_fraci]\n",
    "        \n",
    "        # DA for each loop using a user-defined R scale\n",
    "        for Rscalei in range(Rscalen):\n",
    "            if Rscale_style == 1:\n",
    "                Rscale_use = Rscale_list[Rscalei]\n",
    "                \n",
    "            #######     ########     #######     ########     #######     ########     #######     ########     \n",
    "            if log_level > 1:\n",
    "                print('')\n",
    "                print('>>  Starting Monte Carlo ... ')\n",
    "            #######     ########     #######     ########     #######     ########     #######     ########\n",
    "            \n",
    "            \n",
    "            ### TEST SAVE FULL\n",
    "            # Prepare empty matrix for saving all DA results\n",
    "            if prior_variable2d_len>0:\n",
    "                Xa_shape_MC = (prior_variable2d_len*dum_jmax*dum_imax, prior_len, recon_period_len, MCn)  # lonn * latn * varn\n",
    "                Xa_MC   = np.full(Xa_shape_MC,np.nan)\n",
    "\n",
    "            # prepare 3D version of Xb\n",
    "            if prior_variable3d_len > 0:\n",
    "                Xa3d_shape_MC = (prior_variable3d_len*dum_dmax*dum_jmax*dum_imax, prior_len, recon_period_len, MCn)  # lonn * latn * varn\n",
    "                Xa3d_MC = np.full(Xa3d_shape_MC,np.nan)\n",
    "            \n",
    "            ### TEST END\n",
    "            # Monte Carlo simulations\n",
    "            for MCi in range(MCn):\n",
    "            # debug\n",
    "            #for MCj in range(100-59):\n",
    "            #    MCi = MCj + 59\n",
    "            \n",
    "                # copy back:\n",
    "                Xb = np.copy(Xb_prior)\n",
    "                if prior_variable3d_len > 0:\n",
    "                    Xb3d = np.copy(Xb3d_prior)\n",
    "                if data_psm_d18o_find == 1:\n",
    "                    Xb_ph        = np.copy(Xb_ph_prior)\n",
    "                    Xb_ph1 = np.copy(Xb_ph)\n",
    "                    Xb_ph1[Xb_ph1> 3.0e+36] = np.nan\n",
    "                    Xb_ph_mean = np.nanmean(Xb_ph1)\n",
    "                    \n",
    "                    Xb_sal       = np.copy(Xb_sal_prior)\n",
    "                    Xb_sal1 = np.copy(Xb_sal)\n",
    "                    Xb_sal1[Xb_sal1> 3.0e+36] = np.nan\n",
    "                    Xb_sal_mean = np.nanmean(Xb_sal1)\n",
    "                    \n",
    "                    if log_level > 3:\n",
    "                        print('')\n",
    "                        print('>>    mean of Xb_sal {}, Xb_ph {}'.format(Xb_sal_mean, Xb_ph_mean))\n",
    "                \n",
    "                if data_psm_mgca_find == 1:\n",
    "                    Xb_sal       = np.copy(Xb_sal_prior)\n",
    "                    Xb_ph        = np.copy(Xb_ph_prior)\n",
    "                    Xb_omega     = np.copy(Xb_omega_prior)\n",
    "                    # ``1`` for reductive, ``0`` for BCP (Barker).\n",
    "                    cleaningr = np.copy(cleaningr_prior)\n",
    "                    cleaningb = np.copy(cleaningb_prior)\n",
    "                    # maybe used in Mg/Ca PSM\n",
    "                    Xb_sal1 = np.copy(Xb_sal)\n",
    "                    Xb_sal1[Xb_sal1> 3.0e+36] = np.nan\n",
    "                    Xb_sal_mean = np.nanmean(Xb_sal1)\n",
    "                    Xb_ph1 = np.copy(Xb_ph)\n",
    "                    Xb_ph1[Xb_ph1> 3.0e+36] = np.nan\n",
    "                    Xb_ph_mean = np.nanmean(Xb_ph1)\n",
    "                    Xb_omega1 = np.copy(Xb_omega)\n",
    "                    Xb_omega1[Xb_omega1> 3.0e+36] = np.nan\n",
    "                    Xb_omega_mean = np.nanmean(Xb_omega1)\n",
    "                    if log_level > 3:\n",
    "                        print('')\n",
    "                        print('>>    mean of Xb_sal {}, Xb_ph {}, Xb_omega {}'.format(Xb_sal_mean, Xb_ph_mean, Xb_omega_mean))\n",
    "                \n",
    "                ### Select a fraction of proxy sites ###\n",
    "                if proxy_frac <= 1.0:\n",
    "                    if log_level > 2:\n",
    "                        print('')\n",
    "                        print('>>  Proxy fraction is {}'.format(proxy_frac))\n",
    "                    ## Seed\n",
    "                    \n",
    "                    curr_seed = multi_seed[MCi]\n",
    "                    if log_level > 2:\n",
    "                        print('Setting current prior iteration seed: {}'.format(curr_seed))\n",
    "                    random.seed(curr_seed)\n",
    "                    sites_assim, sites_eval = DeepDA_psm.proxy_frac_4da_eval(proxy_select,proxy_frac,log_level)\n",
    "                else:\n",
    "                    sites_assim = proxy_select.copy()\n",
    "                    sites_eval = []\n",
    "                \n",
    "                ### log\n",
    "                if log_level > 3:\n",
    "                    print('>>  Selected proxy sties: ')\n",
    "                    print(sites_assim)\n",
    "                    print('>>  Un-selected proxy sties: ')\n",
    "                    print(sites_eval)                \n",
    "                \n",
    "                ### sort proxy data using the given order ###\n",
    "                proxies_frac_len = len(sites_assim)\n",
    "                proxy_select_1 = 0\n",
    "                if proxy_order_type in ['use_list']:\n",
    "                    for i in range(len(proxy_order)):\n",
    "                        proxy_order_i = proxy_assim2[proxy_order[i]]\n",
    "                        for j in range(proxies_frac_len):\n",
    "                            # Read proxy type from the database\n",
    "                            data_psm_type = sites_assim['Proxy'][j]\n",
    "                            if data_psm_type in proxy_order_i:\n",
    "                                if proxy_select_1 == 0:\n",
    "                                    # first element\n",
    "                                    #proxy_select_sort = proxy_select.iloc[[j]]\n",
    "                                    proxy_select_sort = sites_assim.iloc[[j]]\n",
    "                                    proxy_select_1 = 1  # proxy included\n",
    "                                else:\n",
    "                                    # rest elements\n",
    "                                    proxy_select_sort = proxy_select_sort.append(sites_assim.iloc[[j]], ignore_index=True)\n",
    "                    if log_level > 2:\n",
    "                        print('>>  Proxy order: use user-defined list.')\n",
    "                else:\n",
    "                    proxy_select_sort = sites_assim.sample(frac=1).reset_index(drop=True);\n",
    "                    if log_level > 2:\n",
    "                        print('>>  Proxy order: use random list.')\n",
    "                \n",
    "                ### update proxies using sorted proxy order ###\n",
    "                proxies =   proxy_select_sort.copy()\n",
    "                proxies_len = len(proxies)\n",
    "\n",
    "                if proxies_len0 > proxies_len:\n",
    "                    if log_level > 1:\n",
    "                        print('>>  Selected proxy data length {}'.format(proxies_len))\n",
    "\n",
    "                ######## Prepare Ye   ########\n",
    "                # for saving proxy unit data Ye\n",
    "                Ye       = np.full((proxies_len,prior_len),np.nan)\n",
    "                obvalue  = np.full((proxies_len,recon_period_len),np.nan)\n",
    "                ob_err   = np.full((proxies_len,recon_period_len),np.nan) # data obs error\n",
    "                ob_err0  = np.full((proxies_len,recon_period_len),np.nan) # PSM obs error\n",
    "                ob_err_comb  = np.full((proxies_len,recon_period_len),np.nan) # PSM obs error\n",
    "                yo_all = np.full((proxies_len,2),np.nan) # PSM obs error\n",
    "                if log_level > 2:\n",
    "                    print('>>  OKAY.')\n",
    "                    print('')\n",
    "                # check the consistency of the config.yml file and proxy database\n",
    "                # AND get obs R\n",
    "                if log_level > 2:\n",
    "                    print('########## Check the consistency of the config.yml file and proxy database ##########')\n",
    "                    print('')\n",
    "                \n",
    "                proxy_psm_type_dict = {}\n",
    "                for j in range(proxies_len):\n",
    "                    # Read proxy type from the database\n",
    "                    data_psm_type = proxies['Proxy'][j]\n",
    "                    # Read allowed proxy from the DTDA-config.yml\n",
    "                    data_psm_type_find = 0\n",
    "                    for key, value in proxy_assim2.items():\n",
    "                        if log_level > 5:\n",
    "                            print(key,value)\n",
    "                        # check this proxy type exist or not, and how many times it occurrs\n",
    "                        if data_psm_type in proxy_assim2[key]:\n",
    "                            data_psm_type_find = data_psm_type_find + 1\n",
    "                    if data_psm_type_find == 1:\n",
    "                        for key, value in proxy_psm_type.items():\n",
    "                            if data_psm_type in proxy_assim2[key]:\n",
    "                                data_psm_key = key\n",
    "                        proxy_psm_type_i = proxy_psm_type[data_psm_key]\n",
    "\n",
    "                        proxy_psm_type_dict[j] =proxy_psm_type_i\n",
    "\n",
    "                        if log_level > 3:\n",
    "                            print('>>  {}. PSM for {} is {}'.format(j, data_psm_type,proxy_psm_type_i))\n",
    "\n",
    "                    elif data_psm_type_find == 0:\n",
    "                        if log_level > 3:\n",
    "                            print('>>  Warning, {} in database is not find in config.yml dictionary'.format(data_psm_type))\n",
    "                    else:\n",
    "                        if log_level > 3:\n",
    "                            print('>>  Warning, {} in database appears more than one time in config.yml dictionary'.format(data_psm_type))\n",
    "\n",
    "                    # Now PSM type has been found. Let's precal Ye\n",
    "\n",
    "                    if proxy_psm_type_i in ['bayesreg_mgca_pooled_red','bayesreg_mgca_pooled_bcp','deepmip_mgca']:\n",
    "                        data_psm_mgca_find = 1\n",
    "\n",
    "                if log_level > 3:\n",
    "                    print('>>  Proxy_psm_type_dict: ')\n",
    "                    print(proxy_psm_type_dict)\n",
    "                if log_level > 2:\n",
    "                    print('')\n",
    "                    print('>>  All looks good.')\n",
    "                    print('')\n",
    "\n",
    "                    ##### Ye calculation ####\n",
    "\n",
    "                    print('##########  Ye calculation  ##########')\n",
    "                    print('')\n",
    "                    \n",
    "                    \n",
    "                # precal_Ye\n",
    "                proi = 0\n",
    "                for j in range(proxies_len):\n",
    "                    # Read proxy type from the database\n",
    "                    data_psm_type = proxies['Proxy'][j]\n",
    "                    proxy_psm_type_i = proxy_psm_type_dict[j]\n",
    "                    psm_required_variable_key = list(yml_dict['psm'][proxy_psm_type_i]['psm_required_variables'].keys())[0]\n",
    "                    if log_level > 3:\n",
    "                        print(psm_required_variable_key)\n",
    "                    # ID-ID match: proxy type matches with the prior type. This allows assimilate multiple proxy types for multiple state variables\n",
    "                    if psm_required_variable_key in prior_variable_dict:\n",
    "                        psm_required_variable_key_index = prior_variable_dict.index(psm_required_variable_key)\n",
    "                    elif psm_required_variable_key in prior_variable_dict_3d:\n",
    "                        psm_required_variable_key_index = prior_variable_dict_3d.index(psm_required_variable_key)\n",
    "\n",
    "                ######################## FOR 2D field ONLY TO DO: adjusted to include 3d proxies ##############\n",
    "                \n",
    "                    # read lon lat for each line of proxy\n",
    "                    dum_lat = proxies[lat_label][j]  # (paleo)latitude of this site\n",
    "                    dum_lon = proxies[lon_label][j]  # (paleo)longitude of this site\n",
    "                    yo_all[proi,:] = np.array([dum_lon, dum_lat])  # save location of this site\n",
    "\n",
    "                    lonlat = modules_nc.cal_find_ij(dum_lon,dum_lat,dum_lon_offset,dum_imax,dum_jmax)\n",
    "                    # output [lon, lat], \n",
    "                    # lon ranges from 0 (-180) to 35 (180), lat ranges from 0 (-90) to 35 (90)\n",
    "\n",
    "                    Filei = proxies['File'][j]\n",
    "                    Sitei = proxies['Site'][j]\n",
    "                    # find 1d grid location\n",
    "                    lonlati = lonlat[1] * dum_jmax + lonlat[0] + psm_required_variable_key_index * dum_ijmax\n",
    "                    if log_level > 3:\n",
    "                        print('>>  lonlat id is {}'.format(lonlati))\n",
    "                    # read prior\n",
    "                    prior_1grid = np.copy(Xb[lonlati,:])   # prior\n",
    "                    #print(prior_1grid.shape)\n",
    "                    if log_level > 3:\n",
    "                        print(prior_1grid)\n",
    "                    \n",
    "                ######################## FOR 2D field ONLY. TO DO: adjusted to include 3d priors ##############\n",
    "                    if log_level > 2:\n",
    "                        print('')\n",
    "                        print('>>  {}. grid [lon lat] {}, index {}, PSM for {} is {}, File: {}, Site: {}'.format(j,lonlat,lonlati,data_psm_type,proxy_psm_type_i,Filei,Sitei))\n",
    "                    if log_level > 3:\n",
    "                        if psm_required_variable_key in prior_variable_dict:                     \n",
    "                            print('>>    Key Found: {} in prior_variable_dict 2d list, index = {}'.format(psm_required_variable_key, psm_required_variable_key_index))\n",
    "                        elif psm_required_variable_key in prior_variable_dict_3d:\n",
    "                            print('>>    Key Found: {} in prior_variable_dict_3d list, index = {}'.format(psm_required_variable_key, psm_required_variable_key_index))\n",
    "                    if log_level > 3:\n",
    "                        print('>>      Mean of Prior is {:.6f}, variance is {:.6f}'.format(np.mean(prior_1grid), np.var(prior_1grid)))\n",
    "\n",
    "                    # Now PSM type has been found. Let's precal Ye\n",
    "\n",
    "                    if proxy_psm_type_i in ['bayesreg_d18o_pooled']:\n",
    "                        if Rscale_style == 1:\n",
    "                            Rscale = Rscale_use\n",
    "                        elif Rscale_style == 2:\n",
    "                            Rscale = Rscale_d18o\n",
    "                        type_i = 1\n",
    "                        salinity_i = np.copy(Xb_sal1[lonlati,:])\n",
    "                        # water_saturation option\n",
    "                        if d18o_phcor_water in ['modern']:\n",
    "                            latlon_j = (dum_lat, dum_lon)\n",
    "                            ph_i = baymag.omgph.fetch_ph(latlon=latlon_j)\n",
    "                            #ph_i = np.full((prior_len,1),ph_j)\n",
    "                            if log_level > 3:\n",
    "                                print('>>    {} use modern pH is {}'.format(proxy_psm_type_i, ph_i))\n",
    "                        else:\n",
    "                            ph_i       = np.mean(np.copy(Xb_ph1[lonlati,:]))\n",
    "                            \n",
    "                        if d18osw_local_choice in ['zachos94']:\n",
    "                            # d18o_localsw using method by Zachos et al., 1994 PALEOCEANOGRAPHY\n",
    "                            d18o_localsw = DeepDA_psm.d18o_localsw(abs(dum_lat))\n",
    "                            \n",
    "                        else:\n",
    "                            if d18osw_icesm_pco2 == 1.0:\n",
    "                                proxy_col_d18osw = 'd18osw_1x'\n",
    "                            elif d18osw_icesm_pco2 == 3.0:\n",
    "                                proxy_col_d18osw = 'd18osw_3x'\n",
    "                            elif d18osw_icesm_pco2 == 6.0:\n",
    "                                proxy_col_d18osw = 'd18osw_6x'\n",
    "                            elif d18osw_icesm_pco2 == 9.0:\n",
    "                                proxy_col_d18osw = 'd18osw_9x'\n",
    "                            else:\n",
    "                                proxy_col_d18osw = 'd18osw_3x'\n",
    "                            d18o_localsw = proxies[proxy_col_d18osw][j]\n",
    "                            \n",
    "                        # total d18osw = d18o_localsw + d18o_adj + psm_d18osw_adjust\n",
    "                        # d18o_adj has been included in the bayfox model\n",
    "                        #print('>>  Prior is {}'.format(prior_1grid))\n",
    "                        if d18osw_local_choice in ['zachos94']:\n",
    "                            prediction_d18O_raw = bayfox.predict_d18oc(prior_1grid,d18o_localsw + psm_d18osw_adjust) # pool model for bayfox\n",
    "                            prediction_d18O     = bayfox.pHCorrect(prediction_d18O_raw.ensemble,prior_1grid,salinity_i,ph_i,type_i)\n",
    "                            if log_level > 3:\n",
    "                                print('>>        Sea water d18O is {:.6f}, d18osw_adjust is {:.6f} '.format(d18o_localsw, psm_d18osw_adjust))\n",
    "                        else:\n",
    "                            prediction_d18O_raw = bayfox.predict_d18oc(prior_1grid,d18o_localsw) # pool model for bayfox\n",
    "                            prediction_d18O     = bayfox.pHCorrect(prediction_d18O_raw.ensemble,prior_1grid,salinity_i,ph_i,type_i)\n",
    "                            if log_level > 3:\n",
    "                                print('>>        Sea water d18O is {:.6f}'.format(d18o_localsw))\n",
    "                        if log_level > 4:\n",
    "                            print('>>  prediction_d18O.ensemble shape {}'.format(prediction_d18O.ensemble.shape))\n",
    "                        \n",
    "                        Ye[proi,:] = np.mean(prediction_d18O.ensemble, axis = 1)\n",
    "                        if log_level > 4:\n",
    "                            print('>>  Ye is {}'.format(Ye[proi,:]))\n",
    "                        if log_level > 3:\n",
    "                            print('>>      Mean of  Ye  is {:.6f}, variance is {:.6f} '.format(np.mean(Ye[proi,:]), np.var(Ye[proi,:],ddof=1)))\n",
    "                        for reconi in range(recon_period_len):\n",
    "                            reconid = recon_period_full[reconi]\n",
    "                            obvalue[proi,reconi] = proxies[data_period_id[reconid]][j]\n",
    "                            ob_err[proi,reconi] = proxies[data_period_idstd[reconid]][j] ** 2\n",
    "                            if proxy_err_eval in ['proxy_err_psm']:\n",
    "                                if d18osw_local_choice in ['zachos94']:\n",
    "                                    ob_err0[proi,reconi]= DeepDA_psm.obs_estimate_r_d18o(obvalue[proi,reconi], d18o_localsw+psm_d18osw_adjust) * Rscale\n",
    "                                else:\n",
    "                                    ob_err0[proi,reconi]= DeepDA_psm.obs_estimate_r_d18o(obvalue[proi,reconi], d18o_localsw) * Rscale\n",
    "                            else:\n",
    "                                ob_err0[proi,reconi]= DeepDA_psm.obs_estimate_r_fixed_d18o(15) * Rscale\n",
    "                            ob_err_comb[proi,reconi] = np.nansum([ob_err[proi,reconi], ob_err0[proi,reconi]])\n",
    "                            if ob_err_comb[proi,reconi] == 0: ob_err_comb[proi,reconi] = np.nan\n",
    "                            if log_level > 3:\n",
    "                                print('>>   {}. Proxy variance from PSM is {:.6f} vs. from PSM + time variance is {:.6f} '.format(reconi,ob_err0[proi,reconi], ob_err_comb[proi,reconi]))\n",
    "\n",
    "                            # Quality control\n",
    "                            if log_level > 1:\n",
    "                                if proxy_qc is not None:\n",
    "                                    print('>>   Quality Control (QC) ...')\n",
    "                            if proxy_err_eval in ['proxy_err_psm']:\n",
    "                                qc_i = DeepDA_psm.obs_qc(Ye[proi,:], obvalue[proi,reconi], ob_err_comb[proi,reconi], proxy_qc)\n",
    "                            else:\n",
    "                                qc_i = DeepDA_psm.obs_qc(Ye[proi,:], obvalue[proi,reconi], ob_err0[proi,reconi], proxy_qc)\n",
    "                            #print(qc_i)\n",
    "                            if qc_i:\n",
    "                                if proxy_qc is not None:\n",
    "                                    if log_level > 1:\n",
    "                                        print('    Pass QC. ye {}, obs {}, obs_var {}, qc {}'.format(np.mean(Ye[proi,:]), obvalue[proi,reconi], ob_err_comb[proi,reconi], proxy_qc))\n",
    "                            else:\n",
    "                                ob_err_comb[proi,reconi] = np.nan\n",
    "                                if proxy_qc is not None:   \n",
    "                                    if log_level > 1:\n",
    "                                        print('    Failed QC. ye {}, obs {}, obs_var {}, qc {}'.format(np.mean(Ye[proi,:]), obvalue[proi,reconi], ob_err_comb[proi,reconi], proxy_qc))\n",
    "                        proi = proi + 1  # increasement\n",
    "                        \n",
    "                        \n",
    "                    elif proxy_psm_type_i in ['deepmip_d18o']:\n",
    "                        if Rscale_style == 1:\n",
    "                            Rscale = Rscale_use\n",
    "                        elif Rscale_style == 2:\n",
    "                            Rscale = Rscale_d18o\n",
    "                        if d18osw_local_choice in ['zachos94']:\n",
    "                            # d18o_localsw using method by Zachos et al., 1994 PALEOCEANOGRAPHY\n",
    "                            d18o_localsw = DeepDA_psm.d18o_localsw(abs(dum_lat))\n",
    "                        else:\n",
    "                            if d18osw_icesm_pco2 == 1.0:\n",
    "                                proxy_col_d18osw = 'd18osw_1x'\n",
    "                            elif d18osw_icesm_pco2 == 3.0:\n",
    "                                proxy_col_d18osw = 'd18osw_3x'\n",
    "                            elif d18osw_icesm_pco2 == 6.0:\n",
    "                                proxy_col_d18osw = 'd18osw_6x'\n",
    "                            elif d18osw_icesm_pco2 == 9.0:\n",
    "                                proxy_col_d18osw = 'd18osw_9x'\n",
    "                            else:\n",
    "                                proxy_col_d18osw = 'd18osw_3x'\n",
    "                            d18o_localsw = proxies[proxy_col_d18osw][j]\n",
    "                            \n",
    "                        # total d18osw = d18o_localsw + d18o_adj + psm_d18osw_adjust\n",
    "                        # d18o_adj has been included in the bayfox model\n",
    "                        #print('>>  Prior is {}'.format(prior_1grid))\n",
    "                        \n",
    "                        if d18osw_local_choice in ['zachos94']:\n",
    "                            Ye[proi,:] = DeepDA_psm.d18oc_linear_forward(prior_1grid,d18o_localsw + psm_d18osw_adjust)\n",
    "                            if log_level > 4:\n",
    "                                print('>>        Sea water d18O is {:.6f}, d18osw_adjust is {:.6f} '.format(d18o_localsw, psm_d18osw_adjust))\n",
    "                        else:\n",
    "                            Ye[proi,:] = DeepDA_psm.d18oc_linear_forward(prior_1grid,d18o_localsw)\n",
    "                            if log_level > 4:\n",
    "                                print('>>        Sea water d18O is {:.6f}'.format(d18o_localsw))\n",
    "                        if log_level > 4:\n",
    "                            print('>>  prediction_d18O.ensemble shape {}'.format(prediction_d18O.ensemble.shape))\n",
    "                        \n",
    "                        if log_level > 4:\n",
    "                            print('>>  Ye is {}'.format(Ye[proi,:]))\n",
    "                        if log_level > 3:\n",
    "                            print('>>      Mean of  Ye  is {:.6f}, variance is {:.6f} '.format(np.mean(Ye[proi,:]), np.var(Ye[proi,:],ddof=1)))\n",
    "                        for reconi in range(recon_period_len):\n",
    "                            reconid = recon_period_full[reconi]\n",
    "                            obvalue[proi,reconi] = proxies[data_period_id[reconid]][j]\n",
    "                            ob_err[proi,reconi] = proxies[data_period_idstd[reconid]][j] ** 2\n",
    "                            if proxy_err_eval in ['proxy_err_psm']:\n",
    "                                if d18osw_local_choice in ['zachos94']:\n",
    "                                    ob_err0[proi,reconi]= DeepDA_psm.obs_estimate_r_d18o(obvalue[proi,reconi], d18o_localsw+psm_d18osw_adjust) * Rscale\n",
    "                                else:\n",
    "                                    ob_err0[proi,reconi]= DeepDA_psm.obs_estimate_r_d18o(obvalue[proi,reconi], d18o_localsw) * Rscale\n",
    "                                    \n",
    "                            elif proxy_err_eval in ['proxy_err_psm_fixed']:\n",
    "                                ob_err0[proi,reconi]= DeepDA_psm.obs_estimate_r_fixed_d18o(15) * Rscale\n",
    "                            else:\n",
    "                                ob_err0[proi,reconi] = yml_dict['psm'][proxy_psm_type_i]['psm_error'] * Rscale\n",
    "                                \n",
    "                            ob_err_comb[proi,reconi] = np.nansum([ob_err[proi,reconi], ob_err0[proi,reconi]])\n",
    "                            if ob_err_comb[proi,reconi] == 0: ob_err_comb[proi,reconi] = np.nan\n",
    "                            if log_level > 3:\n",
    "                                print('>>   {}. Proxy variance from PSM is {:.6f} vs. from PSM + time variance is {:.6f} '.format(reconi,ob_err0[proi,reconi], ob_err_comb[proi,reconi]))\n",
    "\n",
    "                            # Quality control\n",
    "                            if log_level > 1:\n",
    "                                if proxy_qc is not None:\n",
    "                                    print('>>   Quality Control (QC) ...')\n",
    "                            if proxy_err_eval in ['proxy_err_psm']:\n",
    "                                qc_i = DeepDA_psm.obs_qc(Ye[proi,:], obvalue[proi,reconi], ob_err_comb[proi,reconi], proxy_qc)\n",
    "                            elif proxy_err_eval in ['proxy_err_psm_fixed']:\n",
    "                                qc_i = DeepDA_psm.obs_qc(Ye[proi,:], obvalue[proi,reconi], ob_err0[proi,reconi], proxy_qc)\n",
    "                            else:\n",
    "                                qc_i = DeepDA_psm.obs_qc(Ye[proi,:], obvalue[proi,reconi], ob_err0[proi,reconi], proxy_qc)\n",
    "                            #print(qc_i)\n",
    "                            if qc_i:\n",
    "                                if proxy_qc is not None:\n",
    "                                    if log_level > 1:\n",
    "                                        print('    Pass QC. ye {}, obs {}, obs_var {}, qc {}'.format(np.mean(Ye[proi,:]), obvalue[proi,reconi], ob_err_comb[proi,reconi], proxy_qc))\n",
    "                            else:\n",
    "                                ob_err_comb[proi,reconi] = np.nan\n",
    "                                if proxy_qc is not None:      \n",
    "                                    if log_level > 1:\n",
    "                                        print('    Failed QC. ye {}, obs {}, obs_var {}, qc {}'.format(np.mean(Ye[proi,:]), obvalue[proi,reconi], ob_err_comb[proi,reconi], proxy_qc))\n",
    "                        proi = proi + 1  # increasement\n",
    "                        \n",
    "                    elif proxy_psm_type_i in ['cgenie_caco3', 'cgenie_caco3_13c']:\n",
    "                        if Rscale_style == 1:\n",
    "                            Rscale = Rscale_use\n",
    "                        elif Rscale_style == 2:\n",
    "                            Rscale = Rscale_caco3\n",
    "                            \n",
    "                        Ye[proi,:] = np.nanmean(prior_1grid)\n",
    "                        for reconi in range(recon_period_len):\n",
    "                            reconid = recon_period_full[reconi]\n",
    "                            obvalue[proi,reconi] = proxies[data_period_id[reconid]][j]\n",
    "                            ob_err[proi,reconi] = proxies[data_period_idstd[reconid]][j] ** 2\n",
    "                            ob_err0[proi,reconi] = yml_dict['psm'][proxy_psm_type_i]['psm_error'] * Rscale\n",
    "                            ob_err_comb[proi,reconi] = np.nansum([ob_err[proi,reconi], ob_err0[proi,reconi]])\n",
    "                            # Quality control\n",
    "                            if proxy_err_eval in ['proxy_err_psm']:\n",
    "                                qc_i = DeepDA_psm.obs_qc(Ye[proi,:], obvalue[proi,reconi], ob_err_comb[proi,reconi], proxy_qc)\n",
    "                            else:\n",
    "                                qc_i = DeepDA_psm.obs_qc(Ye[proi,:], obvalue[proi,reconi], ob_err0[proi,reconi], proxy_qc)\n",
    "                            if qc_i:\n",
    "                                if proxy_qc is not None:\n",
    "                                    if log_level > 1:\n",
    "                                        print('    Pass QC. ye {}, obs {}, obs_var {}, qc {}'.format(np.nanmean(Ye[proi,:]), obvalue[proi,reconi], ob_err_comb[proi,reconi], proxy_qc))\n",
    "                            else:\n",
    "                                ob_err_comb[proi,reconi] = np.nan\n",
    "                                if proxy_qc is not None: \n",
    "                                    if log_level > 1:\n",
    "                                        print('    Failed QC. ye {}, obs {}, obs_var {}, qc {}'.format(np.nanmean(Ye[proi,:]), obvalue[proi,reconi], ob_err_comb[proi,reconi], proxy_qc))\n",
    "                        proi = proi + 1  # increasement\n",
    "                        \n",
    "                    elif proxy_psm_type_i in ['bayesreg_tex86']:\n",
    "                        if Rscale_style == 1:\n",
    "                            Rscale = Rscale_use\n",
    "                        elif Rscale_style == 2:\n",
    "                            Rscale = Rscale_tex\n",
    "                        # bayspar\n",
    "                        try:\n",
    "                            prediction = bayspar.predict_tex_analog(prior_1grid, temptype = 'sst', search_tol = search_tol_i, nens=nens_i)\n",
    "                        except:\n",
    "                            prediction = bayspar.predict_tex_analog(prior_1grid, temptype = 'sst', search_tol = search_tol_i+10, nens=nens_i)\n",
    "                            \n",
    "                        Ye[proi,:] = np.mean(prediction.ensemble, axis = 1)\n",
    "                        if log_level > 3:\n",
    "                            print('>>      Mean of  Ye   is {:.6f}, variance is {:.6f} '.format(np.mean(Ye[proi,:]), np.var(Ye[proi,:],ddof=1)))\n",
    "                        for reconi in range(recon_period_len):\n",
    "                            reconid = recon_period_full[reconi]\n",
    "                            obvalue[proi,reconi] = proxies[data_period_id[reconid]][j]\n",
    "                            ob_err[proi,reconi] = proxies[data_period_idstd[reconid]][j] ** 2\n",
    "                            if proxy_err_eval in ['proxy_err_psm']:\n",
    "                                ob_err0[proi,reconi]= DeepDA_psm.obs_estimate_r_tex86(np.array([31]), 'sst', 15)  * Rscale\n",
    "                            else:\n",
    "                                ob_err0[proi,reconi]= DeepDA_psm.obs_estimate_r_fixed_tex86(31)  * Rscale\n",
    "                            #obvalue[proi,] = proxies['Lat'][j]\n",
    "                            ob_err_comb[proi,reconi] = np.nansum([ob_err[proi,reconi], ob_err0[proi,reconi]])\n",
    "                            if ob_err_comb[proi,reconi] == 0: ob_err_comb[proi,reconi] = np.nan\n",
    "                            if log_level > 3:\n",
    "                                print('>>   {}. Proxy variance from PSM is {:.6f}, from PSM and selected interval is {:.6f} '.format(reconi,ob_err0[proi,reconi], ob_err_comb[proi,reconi]))\n",
    "                            # Quality control\n",
    "                            if proxy_err_eval in ['proxy_err_psm']:\n",
    "                                qc_i = DeepDA_psm.obs_qc(Ye[proi,:], obvalue[proi,reconi], ob_err_comb[proi,reconi], proxy_qc)\n",
    "                            else:\n",
    "                                qc_i = DeepDA_psm.obs_qc(Ye[proi,:], obvalue[proi,reconi], ob_err0[proi,reconi], proxy_qc)\n",
    "                            if qc_i:\n",
    "                                if proxy_qc is not None:\n",
    "                                    if log_level > 1:\n",
    "                                        print('    Pass QC. ye {}, obs {}, obs_var {}, qc {}'.format(np.mean(Ye[proi,:]), obvalue[proi,reconi], ob_err_comb[proi,reconi], proxy_qc))\n",
    "                            else:\n",
    "                                ob_err_comb[proi,reconi] = np.nan\n",
    "                                if proxy_qc is not None: \n",
    "                                    if log_level > 1:\n",
    "                                        print('    Failed QC. ye {}, obs {}, obs_var {}, qc {}'.format(np.mean(Ye[proi,:]), obvalue[proi,reconi], ob_err_comb[proi,reconi], proxy_qc))\n",
    "                        proi = proi + 1  # increasement\n",
    "                        #except:\n",
    "                        #    if log_level > 2:\n",
    "                        #        print('>>  Warning {}'.format(proxy_psm_type_i))\n",
    "                        #        print('>>  search_tol too small for {}: mean sst is {}'.format(j, np.mean(prior_1grid)))\n",
    "                        \n",
    "                    elif proxy_psm_type_i in ['tex86h_forward']:\n",
    "                        if Rscale_style == 1:\n",
    "                            Rscale = Rscale_use\n",
    "                        elif Rscale_style == 2:\n",
    "                            Rscale = Rscale_tex\n",
    "                        Ye[proi,:] = DeepDA_psm.tex86h_forward(prior_1grid)\n",
    "                        if log_level > 3:\n",
    "                            print('>>      Mean of  Ye   is {:.6f}, variance is {:.6f} '.format(np.mean(Ye[proi,:]), np.var(Ye[proi,:],ddof=1)))\n",
    "                        for reconi in range(recon_period_len):\n",
    "                            reconid = recon_period_full[reconi]\n",
    "                            obvalue[proi,reconi] = proxies[data_period_id[reconid]][j]\n",
    "                            ob_err[proi,reconi] = proxies[data_period_idstd[reconid]][j] ** 2\n",
    "                            if proxy_err_eval in ['proxy_err_psm']:\n",
    "                                ob_err0[proi,reconi]= DeepDA_psm.obs_estimate_r_tex86(np.array([31]), 'sst', 15)  * Rscale\n",
    "                            elif proxy_err_eval in ['proxy_err_psm_fixed']:\n",
    "                                ob_err0[proi,reconi]= DeepDA_psm.obs_estimate_r_fixed_tex86(31)  * Rscale\n",
    "                            else:\n",
    "                                ob_err0[proi,reconi] = yml_dict['psm'][proxy_psm_type_i]['psm_error'] * Rscale\n",
    "                            #obvalue[proi,] = proxies['Lat'][j]\n",
    "                            ob_err_comb[proi,reconi] = np.nansum([ob_err[proi,reconi], ob_err0[proi,reconi]])\n",
    "                            if ob_err_comb[proi,reconi] == 0: ob_err_comb[proi,reconi] = np.nan\n",
    "                            if log_level > 3:\n",
    "                                print('>>   {}. Proxy variance from PSM is {:.6f}, from PSM and selected interval is {:.6f} '.format(reconi,ob_err0[proi,reconi], ob_err_comb[proi,reconi]))\n",
    "                            # Quality control\n",
    "                            if proxy_err_eval in ['proxy_err_psm']:\n",
    "                                qc_i = DeepDA_psm.obs_qc(Ye[proi,:], obvalue[proi,reconi], ob_err_comb[proi,reconi], proxy_qc)\n",
    "                            else:\n",
    "                                qc_i = DeepDA_psm.obs_qc(Ye[proi,:], obvalue[proi,reconi], ob_err0[proi,reconi], proxy_qc)\n",
    "                            if qc_i:\n",
    "                                if proxy_qc is not None:\n",
    "                                    if log_level > 1:\n",
    "                                        print('    Pass QC. ye {}, obs {}, obs_var {}, qc {}'.format(np.mean(Ye[proi,:]), obvalue[proi,reconi], ob_err_comb[proi,reconi], proxy_qc))\n",
    "                            else:\n",
    "                                ob_err_comb[proi,reconi] = np.nan\n",
    "                                if proxy_qc is not None:      \n",
    "                                    if log_level > 1:\n",
    "                                        print('    Failed QC. ye {}, obs {}, obs_var {}, qc {}'.format(np.mean(Ye[proi,:]), obvalue[proi,reconi], ob_err_comb[proi,reconi], proxy_qc))\n",
    "                        proi = proi + 1  # increasement\n",
    "\n",
    "                    #elif proxy_psm_type_i in ['bayesreg_uk37']:\n",
    "                        # \n",
    "                        #print('... bayesreg_uk37: To be done ...')\n",
    "                        \n",
    "                    elif proxy_psm_type_i in ['deepmip_mgca']:\n",
    "                        if Rscale_style == 1:\n",
    "                            Rscale = Rscale_use\n",
    "                        elif Rscale_style == 2:\n",
    "                            Rscale = Rscale_mgca\n",
    "                        salinity =  np.copy(Xb_sal1[lonlati,:])\n",
    "                        ph       =  np.copy(Xb_ph1[lonlati,:])\n",
    "                        mgcasw = yml_dict['psm'][proxy_psm_type_i]['mgcasw']\n",
    "                        mgcacorr = DeepDA_psm.mgca_evans18_forward(prior_1grid,ph,mgcasw)\n",
    "                        Ye[proi,:] = DeepDA_psm.mgca_sal_corr_forward(mgcacorr,salinity)\n",
    "                        \n",
    "                        for reconi in range(recon_period_len):\n",
    "                            reconid = recon_period_full[reconi]\n",
    "                            obvalue[proi,reconi] = proxies[data_period_id[reconid]][j]\n",
    "                            ob_err[proi,reconi]  = proxies[data_period_idstd[reconid]][j] ** 2\n",
    "                            if proxy_err_eval in ['proxy_err_psm', 'proxy_err_psm_fixed']:\n",
    "                                clearning_one = cleaningb  # use barker cleaning model\n",
    "                                if log_level > 3:\n",
    "                                    print('>>    mean of Xb_sal {}, Xb_ph {}, Xb_omega {}, cleaning {}'.format(Xb_sal_mean, Xb_ph_mean, Xb_omega_mean, clearning_one[0]))\n",
    "                            if proxy_err_eval in ['proxy_err_psm']:\n",
    "                                ob_err0[proi,reconi] = DeepDA_psm.obs_estimate_r_mgca_pooled(obvalue[proi,reconi], clearning_one[0], np.mean(salinity), np.mean(ph), np.mean(omega), spp, geologic_age) * Rscale\n",
    "                            elif proxy_err_eval in ['proxy_err_psm_fixed']:\n",
    "                                ob_err0[proi,reconi] = DeepDA_psm.obs_estimate_r_fixed_mgca_pooled((15, 16), clearning_one[0], Xb_sal_mean, Xb_ph_mean, Xb_omega_mean, spp, geologic_age) * Rscale\n",
    "                            else:\n",
    "                                ob_err0[proi,reconi] = yml_dict['psm'][proxy_psm_type_i]['psm_error'] * Rscale\n",
    "                                \n",
    "                            ob_err_comb[proi,reconi] = np.nansum([ob_err[proi,reconi], ob_err0[proi,reconi]])\n",
    "                            if ob_err_comb[proi,reconi] == 0: ob_err_comb[proi,reconi] = np.nan\n",
    "                            if log_level > 3:\n",
    "                                print('>>   {}. Proxy variance from PSM is {:.6f}, from PSM and selected interval is {:.6f} '.format(reconi,ob_err0[proi,reconi], ob_err_comb[proi,reconi]))\n",
    "                            # Quality control\n",
    "                            if proxy_err_eval in ['proxy_err_psm']:\n",
    "                                qc_i = DeepDA_psm.obs_qc(Ye[proi,:], obvalue[proi,reconi], ob_err_comb[proi,reconi], proxy_qc)\n",
    "                            else:\n",
    "                                qc_i = DeepDA_psm.obs_qc(Ye[proi,:], obvalue[proi,reconi], ob_err0[proi,reconi], proxy_qc)\n",
    "                            if qc_i:\n",
    "                                if proxy_qc is not None:\n",
    "                                    if log_level > 1:\n",
    "                                        print('      Pass QC. ye {}, obs {}, obs_var {}, qc {}'.format(np.mean(Ye[proi,:]), obvalue[proi,reconi], ob_err_comb[proi,reconi], proxy_qc))\n",
    "                            else:\n",
    "                                ob_err_comb[proi,reconi] = np.nan\n",
    "                                if proxy_qc is not None:  \n",
    "                                    if log_level > 1:\n",
    "                                        print('      Failed QC. ye {}, obs {}, obs_var {}, qc {}'.format(np.mean(Ye[proi,:]), obvalue[proi,reconi], ob_err_comb[proi,reconi], proxy_qc))\n",
    "                        proi = proi + 1  # increasement\n",
    "\n",
    "                    elif proxy_psm_type_i in ['bayesreg_mgca_pooled_red', 'bayesreg_mgca_pooled_bcp']:\n",
    "                        if Rscale_style == 1:\n",
    "                            Rscale = Rscale_use\n",
    "                        elif Rscale_style == 2:\n",
    "                            Rscale = Rscale_mgca\n",
    "                        if proxy_psm_type_i in ['bayesreg_mgca_pooled_red']:\n",
    "                            clearning_one = cleaningr\n",
    "                            proxy_explain = 'reductive'\n",
    "                        elif proxy_psm_type_i in ['bayesreg_mgca_pooled_bcp']:\n",
    "                            clearning_one = cleaningb\n",
    "                            proxy_explain = 'barker'\n",
    "                        #try:\n",
    "                        # prior_1grid = np.copy(Xb[lonlati,:])   # prior\n",
    "                        salinity =  np.copy(Xb_sal1[lonlati,:])#.reshape((prior_len,1))\n",
    "                        \n",
    "                        if water_saturation in ['modern']:\n",
    "                            \n",
    "                            latlon_i = (dum_lat, dum_lon)\n",
    "                            ph_i = baymag.omgph.fetch_ph(latlon=latlon_i)\n",
    "                            depth_i = proxies['Depth'][j]\n",
    "                            omg_i = baymag.omgph.fetch_omega(latlon=latlon_i, depth=depth_i)\n",
    "        \n",
    "                            ph       = np.full((prior_len,),ph_i)\n",
    "                            omega    = np.full((prior_len,),omg_i)\n",
    "                            if log_level > 3:\n",
    "                                print('>>  {} use modern pH and omega'.format(proxy_psm_type_i))\n",
    "                        else:\n",
    "                            ph       =  np.copy(Xb_ph1[lonlati,:])#.reshape((prior_len,1))\n",
    "                            omega    =  np.copy(Xb_omega1[lonlati,:])#.reshape((prior_len,1))\n",
    "                            \n",
    "                        prior_1grid = prior_1grid\n",
    "                        \n",
    "                        if log_level > 3:\n",
    "                            print('>>    mean of Xb_sal {}, Xb_ph {}, Xb_omega {}, cleaning {}'.format(Xb_sal_mean, Xb_ph_mean, Xb_omega_mean, clearning_one[0]))\n",
    "                            \n",
    "                            if log_level > 4:\n",
    "                                if MCi == 0:\n",
    "                                    print('shape of prior_1grid {}, clearning_one {}, salinity {}, ph {}, omega {}'\n",
    "                                          .format(prior_1grid.shape, clearning_one.shape, salinity.shape, ph.shape,omega.shape))\n",
    "                                    \n",
    "                        if psm_baymag_ln in ['no']:\n",
    "                            prediction_mgca = baymag.predict_mgca(prior_1grid, clearning_one[0], salinity, ph, omega, spp) # pool model for baymag reductive\n",
    "                            pred_mgca_adj = baymag.sw_correction(prediction_mgca, np.array([geologic_age]))\n",
    "                        \n",
    "                        if psm_baymag_ln in ['yes']:\n",
    "                            if log_level > 4:\n",
    "                                if MCi == 0:\n",
    "                                    print('baymag ln model for deep time ...')\n",
    "                            pred_mgca_adj = baymag.predict_mgca_ln_dt(geologic_age,prior_1grid, clearning_one[0], salinity, ph, omega, spp)\n",
    "                            \n",
    "                        if log_level > 4:\n",
    "                            if MCi == 0: \n",
    "                                #print(prediction_mgca.ensemble.shape)\n",
    "                                print(pred_mgca_adj.ensemble.shape)\n",
    "                        Ye[proi,:] = np.mean(pred_mgca_adj.ensemble, axis = 1)\n",
    "                        if log_level > 3:\n",
    "                            print('>>      Mean of  Ye   is {:.6f}, variance is {:.6f} '.format(np.mean(Ye[proi,:]), np.var(Ye[proi,:],ddof=1)))\n",
    "\n",
    "                        for reconi in range(recon_period_len):\n",
    "                            reconid = recon_period_full[reconi]\n",
    "                            \n",
    "                            obvalue[proi,reconi] = proxies[data_period_id[reconid]][j]\n",
    "                            \n",
    "                            ob_err[proi,reconi]  = proxies[data_period_idstd[reconid]][j] ** 2\n",
    "                            #obs_estimate_r_mgca_pooled(obs, cleaning, salinity, ph, omega, spp, age):\n",
    "                            if proxy_err_eval in ['proxy_err_psm']:\n",
    "                                ob_err0[proi,reconi] = DeepDA_psm.obs_estimate_r_mgca_pooled(obvalue[proi,reconi], clearning_one[0], np.mean(salinity), np.mean(ph), np.mean(omega), spp, geologic_age) * Rscale\n",
    "                            else:\n",
    "                                #ob_err0[proi,reconi] = DeepDA_psm.obs_estimate_r_fixed_mgca_pooled((15, 16), clearning_one[0], np.mean(salinity), np.mean(ph), np.mean(omega), spp, geologic_age)\n",
    "                                ob_err0[proi,reconi] = DeepDA_psm.obs_estimate_r_fixed_mgca_pooled((15, 16), clearning_one[0], Xb_sal_mean, Xb_ph_mean, Xb_omega_mean, spp, geologic_age) * Rscale\n",
    "                            ob_err_comb[proi,reconi] = np.nansum([ob_err[proi,reconi], ob_err0[proi,reconi]])\n",
    "                            \n",
    "                            if ob_err_comb[proi,reconi] == 0: ob_err_comb[proi,reconi] = np.nan\n",
    "                            \n",
    "                            #if psm_baymag_ln in ['yes']:\n",
    "                            #    obvalue[proi,reconi] = np.log(obvalue[proi,reconi])\n",
    "                            #    # need adjustment?\n",
    "                            #    ob_err0[proi,reconi] = np.log(ob_err0[proi,reconi])\n",
    "                            #    ob_err_comb[proi,reconi] = np.log(ob_err_comb[proi,reconi])\n",
    "                            \n",
    "                            if log_level > 3:\n",
    "                                print('>>   {}. Proxy variance from PSM is {:.6f}, from PSM and selected interval is {:.6f} '.format(reconi,ob_err0[proi,reconi], ob_err_comb[proi,reconi]))\n",
    "                            # Quality control\n",
    "                            if proxy_err_eval in ['proxy_err_psm']:\n",
    "                                qc_i = DeepDA_psm.obs_qc(Ye[proi,:], obvalue[proi,reconi], ob_err_comb[proi,reconi], proxy_qc)\n",
    "                            else:\n",
    "                                qc_i = DeepDA_psm.obs_qc(Ye[proi,:], obvalue[proi,reconi], ob_err0[proi,reconi], proxy_qc)\n",
    "                            if qc_i:\n",
    "                                if proxy_qc is not None:\n",
    "                                    if log_level > 1:\n",
    "                                        print('      Pass QC. ye {}, obs {}, obs_var {}, qc {}'.format(np.mean(Ye[proi,:]), obvalue[proi,reconi], ob_err_comb[proi,reconi], proxy_qc))\n",
    "                            else:\n",
    "                                ob_err_comb[proi,reconi] = np.nan\n",
    "                                if proxy_qc is not None:                    \n",
    "                                    if log_level > 1:\n",
    "                                        print('      Failed QC. ye {}, obs {}, obs_var {}, qc {}'.format(np.mean(Ye[proi,:]), obvalue[proi,reconi], ob_err_comb[proi,reconi], proxy_qc))\n",
    "                        if log_level > 3:\n",
    "                            print('        {}: mean salinity {}, ph {}, omega {}'.format(proxy_explain,np.mean(salinity), np.mean(ph), np.mean(omega)))\n",
    "                        proi = proi + 1  # increasement\n",
    "\n",
    "                    else:\n",
    "                        a = 1\n",
    "                \n",
    "                Ye_mean_print = (np.mean(Ye,axis=1))[np.newaxis]\n",
    "                if log_level > 1:\n",
    "                    print('')\n",
    "                    print('>>  Summary of this Monte Carlo simulation')\n",
    "                    if log_level > 2:\n",
    "                        print('')\n",
    "                        print('>>  Ye mean')\n",
    "                    \n",
    "                        print('>>  {}'.format(Ye_mean_print.T))\n",
    "                        print('>>  Observation value')\n",
    "                        print('>>  {}'.format(obvalue))\n",
    "                        print('>>  Observation error (ob_err0)')\n",
    "                        print('>>  {}'.format(ob_err0))\n",
    "                        print('>>  Observation error (ob_err_comb, from PSM and time-variance)')\n",
    "                        print('>>  {}'.format( ob_err_comb))\n",
    "                        \n",
    "                    print('[Ye, Ob_value]')\n",
    "                Yeobvalue = np.concatenate((Ye_mean_print.T, obvalue), axis=1)\n",
    "                if log_level > 1:\n",
    "                    print('>>  {}'.format(Yeobvalue))\n",
    "                if log_level > 0:\n",
    "                    print('')\n",
    "                    print('##########  Monte Carlo {} / {} => okay   ##########'.format(MCi+1, MCn))\n",
    "                    print('')\n",
    "\n",
    "                MC_dir = dir_proxy_save_dir + dir_proxy_save + nexp +'/'\n",
    "                \n",
    "                if not os.path.exists(MC_dir):\n",
    "                    os.makedirs(MC_dir)\n",
    "                # NetCDF file name\n",
    "                if Rscale_style == 1:\n",
    "                    filename_short = '_loc_', str(locRadv),'_proxy_frac_', str(proxy_frac),'_Rscale_',str(Rscale),'_MC_'\n",
    "                elif Rscale_style == 2:\n",
    "                    filename_short = '_loc_', str(locRadv),'_proxy_frac_', str(proxy_frac),'_Rscale_',str(0),'_MC_' \n",
    "                nc_filename = MC_dir + ''.join(filename_short) + str(MCi) + '.nc'\n",
    "                hdf5name    = MC_dir + ''.join(filename_short) + str(MCi) + '.hdf5'\n",
    "                proxy_psm_type_dict_df = pandas.DataFrame.from_dict(proxy_psm_type_dict, orient='index')\n",
    "\n",
    "                with h5py.File(hdf5name, 'w') as f:\n",
    "                    f.create_dataset('obvalue', data=obvalue)\n",
    "                    f.create_dataset('Ye', data=np.transpose(Ye))\n",
    "                    f.create_dataset('ob_err', data=ob_err)\n",
    "                    f.create_dataset('ob_err0', data=ob_err0)\n",
    "                    f.create_dataset('ob_err_comb', data=ob_err_comb)\n",
    "                    f.create_dataset('yo_all', data=yo_all)\n",
    "                    if MCi == 0:\n",
    "                        # if any 2d field selected\n",
    "                        if prior_variable2d_len > 0:\n",
    "                            f.create_dataset('Xb', data=Xb)\n",
    "                        # If any 3d field saved\n",
    "                        if prior_variable3d_len>0:\n",
    "                            f.create_dataset('Xb3d', data=Xb3d)\n",
    "                        if data_psm_d18o_find == 1:\n",
    "                            f.create_dataset('Xb_ph', data=Xb_ph1)\n",
    "                            f.create_dataset('Xb_sal', data=Xb_sal1)\n",
    "                        # if Mg/Ca proxy are used\n",
    "                        if data_psm_mgca_find == 1:\n",
    "                            if data_psm_d18o_find == 0:\n",
    "                                f.create_dataset('Xb_sal', data=Xb_sal1)\n",
    "                            if data_psm_d18o_find == 0:\n",
    "                                f.create_dataset('Xb_ph', data=Xb_ph1)\n",
    "                            f.create_dataset('Xb_omega', data=Xb_omega1)\n",
    "\n",
    "                    metadata = {'Date': time.time(),\n",
    "                                'proxy_dbversion':yml_dict['proxies'][yml_dict['proxies']['use_from'][0]]['dbversion'],\n",
    "                                'exp_dir':dir_prior,\n",
    "                                'Nens':str(prior_len)}\n",
    "\n",
    "                    f.attrs.update(metadata)\n",
    "\n",
    "                # append proxy to hdf5 file\n",
    "                proxies.to_hdf(hdf5name, key='proxies')\n",
    "                proxy_psm_type_dict_df.to_hdf(hdf5name, key='proxy_psm_type_dict_df')\n",
    "\n",
    "                if proxy_frac < 1.0:\n",
    "                    sites_eval.to_hdf(hdf5name, key='sites_eval')\n",
    "                pandas.DataFrame(prior_variable_dict).to_hdf(hdf5name, key='prior_variable_dict')\n",
    "                pandas.DataFrame(prior_variable_dict_3d).to_hdf(hdf5name, key='prior_variable_dict_3d')\n",
    "                if log_level > 1:\n",
    "                    print('>>  prior2proxyunit hdf5 file saved')\n",
    "                    print('      {}'.format(hdf5name))\n",
    "                    print('>>  Finished Step 1. Preparation')\n",
    "                    print('>>  Now Run: Step 2. Data Assimilation ...')\n",
    "                    print('')\n",
    "\n",
    "\n",
    "                #######     ########     #######     ########     #######     ########     #######     ########     \n",
    "                # STEP 2 Data Assimilation\n",
    "                #######     ########     #######     ########     #######     ########     #######     ########\n",
    "\n",
    "\n",
    "                ########## Prior #########\n",
    "                # prior variable list\n",
    "                prior_variable_dict = []  # variable list\n",
    "                prior_nc_file_list = []  # nc file list\n",
    "                prior_variable_dict_3d = []  # variable list\n",
    "                prior_nc_file_list_3d = []  # nc file list\n",
    "\n",
    "                for key, value in prior_state_variable.items():\n",
    "                    nc_keyvalue = prior_state_variable[key]['ncname']  # note: 2d dict\n",
    "                    if log_level > 3:\n",
    "                        print('>>  nc_keyvalue {}...'.format(nc_keyvalue))\n",
    "                    for key1, value1 in nc_keyvalue.items():\n",
    "                        if log_level > 3:\n",
    "                            print('>>  {}: {}'.format(key1,value1))\n",
    "\n",
    "                        for i in range(len(prior_state_variable[key][value1])):\n",
    "                            if key in ['2d']:\n",
    "                                prior_variable_dict.append(prior_state_variable[key][value1][i])\n",
    "                                prior_nc_file_list.append(key1+'/'+value1+'.nc')\n",
    "                            elif key in ['3d']:\n",
    "                                prior_variable_dict_3d.append(prior_state_variable[key][value1][i])\n",
    "                                prior_nc_file_list_3d.append(key1+'/'+value1+'.nc')\n",
    "\n",
    "                # variable list\n",
    "                prior_variable_len = len(prior_variable_dict)\n",
    "                prior_variable3d_len = len(prior_variable_dict_3d)\n",
    "                if log_level > 2:\n",
    "                    print('>>  Number of 2d prior variables is: {}. List:'.format(prior_variable_len))\n",
    "                    print('      {}'.format(prior_variable_dict))\n",
    "                    print('>>  Number of 3d prior variables is: {}. List:'.format(prior_variable3d_len))\n",
    "                    print('      {}'.format(prior_variable_dict_3d))\n",
    "\n",
    "                # for saving DA product Xa\n",
    "                if prior_variable_len > 0:\n",
    "                    Xa_output   = np.full((dum_ijmax * prior_variable_len, nens, recon_period_len),np.nan)\n",
    "                    Xa_output_all = Xa_output\n",
    "                    if prior_variable3d_len > 0:\n",
    "                        Xa3d_output   = np.full((dum_ijmax * dum_dmax * prior_variable_len, nens, recon_period_len),np.nan)\n",
    "                        Xa_output_all = np.concatenate((Xa_output, Xa3d_output), axis=0)\n",
    "                    else:\n",
    "                        if log_level > 1:\n",
    "                            print('>>  No 3d variable listed in {}'.format(config_name))\n",
    "                elif prior_variable_len == 0:\n",
    "                    if prior_variable3d_len > 0:\n",
    "                        Xa3d_output   = np.full((dum_ijmax * dum_dmax * prior_variable_len, nens, recon_period_len),np.nan)\n",
    "                        Xa_output_all = Xa3d_output\n",
    "                    if log_level > 1:\n",
    "                        print('>>  No 2d variable listed in {}'.format(config_name))\n",
    "                else:\n",
    "                    if log_level > 0:\n",
    "                        print('>>  Error! No 3d or 2d variables are listed in {}'.format(config_name))\n",
    "\n",
    "                # DA core script\n",
    "\n",
    "                proxies=pandas.read_hdf(hdf5name, 'proxies')\n",
    "                proxy_psm_type_dict_df = pandas.read_hdf(hdf5name, 'proxy_psm_type_dict_df')\n",
    "                proxy_psm_type_dict_list = proxy_psm_type_dict_df[0].values.tolist()\n",
    "                \n",
    "                # read the first hdf5 file for repeated xb and xb3d\n",
    "                hdf5name0    = MC_dir + ''.join(filename_short) + str(0) + '.hdf5'\n",
    "                with h5py.File(hdf5name0, 'r') as f:\n",
    "                    Xb = f.get('Xb')  # read Xb, background 2d field data\n",
    "                    Xb3d = f.get('Xb3d')  # read Xb, background 3d field data order: lon-lat-depth\n",
    "                    if Xb and Xb3d:\n",
    "                        Xball = np.concatenate((Xb, Xb3d), axis=0)\n",
    "                    elif Xb and Xb3d is None:\n",
    "                        Xball = Xb\n",
    "                    elif Xb is None and Xb3d:\n",
    "                        Xball = Xb3d\n",
    "                    else:\n",
    "                        print('>>  Error! No 3d or 2d variables are listed in {}'.format(config_name))\n",
    "                    \n",
    "                    Xb0 = np.copy(Xball)  # default Xb\n",
    "                    \n",
    "\n",
    "                    d18o_psm_list = ['bayesreg_d18o_pooled', 'deepmip_d18o']\n",
    "                    if any(item in d18o_psm_list for item in proxy_psm_type_dict_list):\n",
    "                        Xb_sal = f.get('Xb_sal')\n",
    "                        Xb_ph = f.get('Xb_ph')\n",
    "                        if log_level > 1:\n",
    "                            print('')\n",
    "                            print('>>  d18O proxy found. Loading salinity, pH and omega')\n",
    "                    mgca_psm_list = ['bayesreg_mgca_pooled_bcp', 'bayesreg_mgca_pooled_red', 'deepmip_mgca']\n",
    "                    if any(item in mgca_psm_list for item in proxy_psm_type_dict_list):\n",
    "                        #if 'bayesreg_mgca_pooled_bcp' in proxy_psm_type_dict_list or 'bayesreg_mgca_pooled_red' in proxy_psm_type_dict_list:\n",
    "                        Xb_sal   = f.get('Xb_sal')\n",
    "                        Xb_omega = f.get('Xb_omega')\n",
    "                        Xb_ph    = f.get('Xb_ph')\n",
    "                        if log_level > 1:\n",
    "                            print('')\n",
    "                            print('>>  Mg/Ca proxy found. Loading salinity, pH and omega')\n",
    "                            \n",
    "                            \n",
    "                    with h5py.File(hdf5name, 'r') as f:\n",
    "                        obvalue_full = f.get('obvalue')\n",
    "                        #Ye_full = f.get('Ye')   # not used???   Used for verification?\n",
    "                        #ob_err_full = f.get('ob_err')  # not used???   Used for verification?\n",
    "                        ob_err0_full = f.get('ob_err0')\n",
    "                        ob_err_comb = f.get('ob_err_comb')\n",
    "                        yo_all = f.get('yo_all')  # read location data\n",
    "\n",
    "                        Xa_output_all = np.full((Xball.shape[0], Xball.shape[1], recon_period_len),np.nan)\n",
    "                        ob_len = obvalue_full.shape[0]\n",
    "\n",
    "                        if log_level > 1:\n",
    "                            print('>>  Reconstruction time intervals: {}; Observation data set number: {}'.format(recon_period_len,ob_len))\n",
    "                        for reconi in range(recon_period_len):\n",
    "                            Xball = Xb0.copy()  # initialize Xball\n",
    "                            for obi in range(ob_len):\n",
    "\n",
    "                                yo_loc = yo_all[obi,:]  # read location\n",
    "                                obvalue  = obvalue_full[obi, reconi]  # read observation value\n",
    "                                if proxy_err_eval in ['proxy_err_psm_mp']:\n",
    "                                    ob_err = ob_err_comb[obi, reconi]  # read observation error, use PSM model + interval data uncertainty\n",
    "                                else:\n",
    "                                    ob_err = ob_err0_full[obi, reconi] # read observation error, use PSM model only\n",
    "\n",
    "                                # proxy type\n",
    "                                proxy_psm_type_i = proxy_psm_type_dict_df[0][obi]\n",
    "\n",
    "                                if log_level > 2:\n",
    "\n",
    "                                    print('>>ID Recon: {}, obser: {}, proxy psm type: {}'.format(reconi,obi, proxy_psm_type_i))\n",
    "\n",
    "                                if proxy_psm_type_i in ['bayesreg_tex86','tex86h_forward', 'cgenie_caco3', 'cgenie_caco3_13c']:\n",
    "                                    Ye = DeepDA_psm.cal_ye_cgenie(yml_dict,proxies,obi,Xball,proxy_assim2,proxy_psm_type,dum_lon_offset,dum_imax,dum_jmax)\n",
    "\n",
    "                                elif proxy_psm_type_i in ['bayesreg_d18o_pooled', 'deepmip_d18o']:\n",
    "                                    Ye = DeepDA_psm.cal_ye_cgenie_d18O(yml_dict,proxies,obi,Xball,Xb_sal,Xb_ph,proxy_assim2,proxy_psm_type,dum_lon_offset,dum_imax,dum_jmax)\n",
    "\n",
    "                                elif proxy_psm_type_i in ['bayesreg_mgca_pooled_bcp', 'bayesreg_mgca_pooled_red','deepmip_mgca']:\n",
    "                                    Ye = DeepDA_psm.cal_ye_cgenie_mgca(yml_dict,proxies,obi,Xball,proxy_psm_type_i,dum_lon_offset,dum_imax,dum_jmax,Xb_sal,Xb_ph,Xb_omega,geologic_age)\n",
    "\n",
    "                                if ~np.isnan(obvalue) and ~np.isnan(ob_err_comb[obi, reconi]):\n",
    "                                    if log_level > 1:\n",
    "                                        print('>>  ID Recon: {}, obser: {}. Loc: {}. Mean of Ye {:.6f}, var {:.6f}, obs {:.6f}, obs_err {:.6f}'.format(reconi,obi,yo_loc,np.mean(Ye),np.var(Ye,ddof=1), obvalue, ob_err))\n",
    "                                    if locRad:\n",
    "                                        covloc = modules_nc.covloc_eval(locRad, yo_loc, dum_jmax, dum_imax, cGENIEGrid)\n",
    "                                        covlocext = int(Xball.shape[0] / covloc.shape[0])\n",
    "                                        covloc = np.matlib.repmat(covloc, covlocext, 1).reshape((Xball.shape[0],))\n",
    "                                    else:\n",
    "                                        covloc = np.full((Xball.shape[0],),1)\n",
    "                                    if log_level > 4:\n",
    "                                        print('>>  Shape of Xball {}, ye {}, ob_err {}, covloc {}'.format(Xball.shape, Ye.shape, ob_err.shape, covloc.shape))\n",
    "\n",
    "                                    if proxy_psm_type_i in ['bayesreg_mgca_pooled_bcp', 'bayesreg_mgca_pooled_red']:\n",
    "                                        if psm_baymag_ln in ['yes']:\n",
    "                                            Xa = LMR_DA.enkf_update_array(Xball, np.log(obvalue), Ye, ob_err, loc = covloc)\n",
    "                                        if psm_baymag_ln in ['no']:\n",
    "                                            Xa = LMR_DA.enkf_update_array(Xball, obvalue, Ye, ob_err, loc = covloc)\n",
    "                                    else:\n",
    "                                        Xa = LMR_DA.enkf_update_array(Xball, obvalue, Ye, ob_err, loc = covloc)\n",
    "\n",
    "                                    # June 17, 2020 set hard limit for given variables\n",
    "                                    Xa = DeepDA_tools.deepda_hard_limit(Xa, yml_dict, prior_variable_dict, dum_ijmax,log_level)\n",
    "\n",
    "                                    if reconi == 0 and obi == 0:\n",
    "                                        kcov_saving = 1\n",
    "                                        ye = np.subtract(Ye, np.mean(Ye))\n",
    "                                        xbm = np.mean(Xball,axis=1)\n",
    "                                        Xbp = np.subtract(Xball,xbm[:,None])  # \"None\" means replicate in this dimension\n",
    "                                        kcov = np.dot(Xbp,np.transpose(ye)) / (nens-1)\n",
    "                                    # update Xb using Xa, to assimilate next observation\n",
    "                                    Xball = np.copy(Xa)\n",
    "                                else:\n",
    "                                    if log_level > 1:\n",
    "                                        print('>>  ID Recon: {}, obser: {}. Skip invalid obs.'.format(reconi,obi))\n",
    "                            if log_level > 2:\n",
    "                                print('>>  ... global mean is {}'.format(np.nanmean(Xa)))\n",
    "                            Xa_output_all[:,:,reconi] = np.copy(Xa) # for each reconi, all observations were assimilated. Save final result for this reconi\n",
    "\n",
    "\n",
    "                        if Xb is not None:\n",
    "                            #lenn1 = f.get('Xb').shape[0]\n",
    "                            lenn1 = Xb.shape[0]\n",
    "                            Xa_output_2d = Xa_output_all[0:lenn1,:,:]\n",
    "\n",
    "                            # save all MC\n",
    "                            # Xa_shape_MC = (prior_variable2d_len*dum_jmax*dum_imax, prior_len, recon_period_len, MCn)\n",
    "                            Xa_MC[:,:,:,MCi] = np.copy(Xa_output_all[0:lenn1,:,:])\n",
    "\n",
    "                            if Xb3d:\n",
    "                                lenn2 = f.get('Xb3d').shape[0]\n",
    "                                Xa_output_3d = Xa_output_all[lenn1:lenn2+lenn1,:,:]\n",
    "                                Xa3d_MC[:,:,:,MCi] = np.copy(Xa_output_all[lenn1:lenn2+lenn1,:,:])\n",
    "                        elif Xb is None:\n",
    "                            if Xb3d:\n",
    "                                #lenn2 = f.get('Xb3d').shape[0]\n",
    "                                lenn2 = Xb3d.shape[0]\n",
    "                                Xa_output_3d = Xa_output_all[0:lenn2,:,:]\n",
    "                                # save all MC\n",
    "                                Xa3d_MC[:,:,:,MCi] = np.copy(Xa_output_all[0:lenn2,:,:])\n",
    "                        else:\n",
    "                            print('>>  Error! No 3d or 2d variables are listed in {}'.format(config_name))\n",
    "                    if log_level > 1:\n",
    "                        print('>>')\n",
    "                        print('>>  Finished Step 2. Data Assimilation')\n",
    "                        print('>>  Now      Step 3. Save results')\n",
    "                        print('')\n",
    "                    if log_level > 4:\n",
    "                        print(Xa_output_all.shape)\n",
    "\n",
    "                    # DA save output in the netCDF file\n",
    "                    if save_mc_full:\n",
    "                        #with h5py.File(hdf5name, 'r') as f:\n",
    "                        if log_level > 1:\n",
    "                            print('>>  Start writing netCDF ...')\n",
    "                        # save netCDF file\n",
    "                        nf = Dataset(nc_filename, 'w', format='NETCDF4')\n",
    "                        nf.description = 'DeepDA' + nc_filename\n",
    "                        #Specifying dimensions\n",
    "                        nf.createDimension('lon', len(cGENIEGridB_lon36))\n",
    "                        nf.createDimension('lat', len(cGENIEGridB_lat36))\n",
    "                        z = np.arange(0,1,1) # level 2d\n",
    "                        nf.createDimension('z', len(z))  # level\n",
    "                        nf.createDimension('nens', nens)  # number of ens\n",
    "                        nf.createDimension('time', recon_period_len)\n",
    "                        # Building variables\n",
    "                        longitude = nf.createVariable('Longitude', 'f4', 'lon')\n",
    "                        # Passing data into variables\n",
    "                        longitude[:] = cGENIEGridB_lon36.values\n",
    "\n",
    "                        latitude = nf.createVariable('Latitude', 'f4', 'lat')\n",
    "                        latitude[:] = cGENIEGridB_lat36.values\n",
    "\n",
    "                        levels = nf.createVariable('Levels', 'i4', 'z')\n",
    "                        levels[:] = z  # 2d level\n",
    "                        if Xb3d is not None:\n",
    "                            nf.createDimension('zt', len(zt))\n",
    "                            levels = nf.createVariable('zt', 'f4', 'zt')\n",
    "                            levels[:] = zt\n",
    "\n",
    "                        if locRad:\n",
    "                            #nf.createDimension('prior_var', prior_variable_len)  # level\n",
    "                            covloc_nc = nf.createVariable('covloc', 'f4', ('lat', 'lon'))\n",
    "                            covloc_nc[:,:] = np.copy(covloc[0:dum_ijmax].reshape(dum_jmax,dum_imax))\n",
    "\n",
    "                        if Xb is not None:\n",
    "                            if log_level > 2:\n",
    "                                print('Writing 2d field.')\n",
    "                            for nc_var_i in range(prior_variable_len):\n",
    "                                nc_var_name = prior_variable_dict[nc_var_i]\n",
    "\n",
    "                                j0 = dum_ijmax * nc_var_i\n",
    "                                j1 = dum_ijmax * (nc_var_i+1)\n",
    "                                if log_level > 0:\n",
    "                                    print('')                            \n",
    "                                    print('>>    ID from {} to {}: field is {}'.format(j0, j1,nc_var_name))\n",
    "\n",
    "                                Xb0_i = np.copy(Xb[j0:j1,:])\n",
    "\n",
    "                                Xa_output_i = np.copy(Xa_output_2d[j0:j1,:,:])\n",
    "                                Xa_outputi = Xa_output_i.reshape(dum_imax,dum_jmax,nens,recon_period_len)\n",
    "\n",
    "                                XbNC_mean = nf.createVariable(nc_var_name+'_Xb_mean', 'f4', ('lat', 'lon','z'))\n",
    "                                xbm = np.mean(Xb0_i,axis=1)\n",
    "                                XbNC_mean[:,:,:] = np.copy(xbm.reshape(dum_jmax,dum_imax,1))\n",
    "\n",
    "                                XbNC_variance = nf.createVariable(nc_var_name+'_Xb_variance', 'f4', ('lat', 'lon','z'))\n",
    "                                Xb_temp = np.copy(np.var(Xb0_i,axis=1).reshape(dum_jmax,dum_imax,1))\n",
    "                                Xb_temp = np.ma.MaskedArray(Xb_temp, np.copy(xbm.reshape(dum_jmax,dum_imax,1)) >= 9.9692e+36)\n",
    "                                XbNC_variance[:,:,:] = Xb_temp\n",
    "                                if log_level > 0:\n",
    "                                    print('>>               Xb mean is {:.8f}, std is {:.8f}, var is {:.8f}'.format(np.nanmean(XbNC_mean), np.sqrt(np.nanmean(Xb_temp)), np.nanmean(Xb_temp)))\n",
    "\n",
    "                                XaNC_mean = nf.createVariable(nc_var_name+'_Xa_mean', 'f4', ('lat', 'lon','z','time'))\n",
    "                                Xam_temp = np.copy(np.nanmean(Xa_outputi,axis=2).reshape(dum_jmax,dum_imax,1,recon_period_len))\n",
    "                                XaNC_mean[:,:,:,:] = Xam_temp\n",
    "\n",
    "                                XaNC_variance = nf.createVariable(nc_var_name+'_Xa_variance', 'f4', ('lat', 'lon','z','time'))\n",
    "                                if log_level > 4:\n",
    "                                    print(Xa_outputi[0,0:36,0,0])\n",
    "                                Xa_temp = np.copy(np.ma.var(Xa_outputi,axis=2).reshape(dum_jmax,dum_imax,1,recon_period_len))\n",
    "                                Xa_temp = np.ma.MaskedArray(Xa_temp, Xam_temp >= 9.9692e+36)\n",
    "                                XaNC_variance[:,:,:,:] = Xa_temp\n",
    "\n",
    "                                for reconii in range(recon_period_len):\n",
    "                                    XaNC_mean_i = XaNC_mean[:,:,:,reconii]\n",
    "                                    XaNC_var_i = XaNC_variance[:,:,:,reconii]\n",
    "                                    if log_level > 0:\n",
    "                                        print('>>      Recon {}. Xa mean is {:.8f}, std is {:.8f}, var is {:.8f}'.format(reconii, np.nanmean(XaNC_mean_i),np.sqrt(np.nanmean(XaNC_var_i)), np.nanmean(XaNC_var_i)))\n",
    "\n",
    "                                if save_ens_full:\n",
    "                                    XaNC_full = nf.createVariable(nc_var_name+'_Xa_full', 'f4', ('lat', 'lon', 'nens', 'z','time'))\n",
    "                                    XaNC_full[:,:,:,:,:] = np.copy(Xa_outputi.reshape(dum_jmax,dum_imax,nens,1,recon_period_len))\n",
    "\n",
    "                                    XbNC_full = nf.createVariable(nc_var_name+'_Xb_full', 'f4', ('lat', 'lon', 'nens', 'z'))\n",
    "                                    XbNC_full[:,:,:,:] = np.copy(Xb0_i.reshape(dum_jmax,dum_imax,nens,1))\n",
    "\n",
    "                                if kcov_saving > 0:\n",
    "                                    kcov_i = np.copy(kcov[j0:j1]).reshape(dum_imax,dum_jmax,1)\n",
    "                                    kcov_i = np.ma.MaskedArray(kcov_i, np.copy(xbm.reshape(dum_jmax,dum_imax,1)) >= 9.9692e+36)\n",
    "                                    cov_ob0 = nf.createVariable(nc_var_name+'_obs0'+'_cov', 'f4', ('lat', 'lon','z'))\n",
    "                                    cov_ob0[:,:,:] = kcov_i\n",
    "\n",
    "                                #Add local attributes to variable instances\n",
    "                                longitude.units = '°'\n",
    "                                latitude.units = '°'\n",
    "                                levels.units = 'm'\n",
    "                                #XbNC_mean.units = '°C'\n",
    "                                #XbNC_variance.units = '°C^2'\n",
    "                                #if save_ens_full:\n",
    "                                #    XaNC_full.units = '°C'\n",
    "                                #    XbNC_full.units = '°C'\n",
    "\n",
    "                                #variance.warning = 'test ...'\n",
    "                        if Xb3d is not None:\n",
    "                            if log_level > 1:\n",
    "                                print('Writing 3d field.')\n",
    "                            for nc_var_i in range(prior_variable3d_len):\n",
    "                                nc_var_name = prior_variable_dict_3d[nc_var_i]\n",
    "\n",
    "                                j0 = dum_ijmax * dum_dmax * nc_var_i\n",
    "                                j1 = dum_ijmax * dum_dmax * (nc_var_i+1)\n",
    "                                if log_level > 0:\n",
    "                                    print('')\n",
    "                                    print('>>    ID from {} to {}: field is {}'.format(j0, j1,nc_var_name))\n",
    "\n",
    "                                Xb0_i = np.copy(Xb3d[j0:j1,:])\n",
    "                                Xa_output_i = np.copy(Xa_output_3d[j0:j1,:,:])\n",
    "                                Xa_outputi = Xa_output_i.reshape(dum_imax, dum_jmax,dum_dmax, nens,recon_period_len)\n",
    "\n",
    "                                XbNC_mean = nf.createVariable(nc_var_name+'_Xb_3d_mean', 'f4', ( 'zt', 'lat','lon'))\n",
    "                                xbm = np.mean(Xb0_i,axis=1)\n",
    "                                XbNC_mean[:,:,:] = np.copy(xbm.reshape(dum_dmax,dum_jmax,dum_imax))\n",
    "\n",
    "                                XbNC_variance = nf.createVariable(nc_var_name+'_Xb_3d_variance', 'f4', ( 'zt', 'lat','lon'))\n",
    "                                Xb_temp = np.copy(np.var(Xb0_i,axis= 1).reshape(dum_dmax,dum_jmax,dum_imax))\n",
    "                                Xb_temp = np.ma.MaskedArray(Xb_temp, np.copy(xbm.reshape(dum_dmax,dum_jmax,dum_imax)) >= 9.9692e+36)\n",
    "                                XbNC_variance[:,:,:] = Xb_temp\n",
    "                                if log_level > 0:\n",
    "                                    print('>>               Xb mean is {:.8f}, std is {:.8f}, var is {:.8f}'.format(np.nanmean(XbNC_mean),np.sqrt(np.nanmean(Xb_temp)), np.nanmean(Xb_temp)))\n",
    "\n",
    "                                XaNC_mean = nf.createVariable(nc_var_name+'_Xa_3d_mean', 'f4', ('zt','lat', 'lon','time'))\n",
    "                                Xam_temp = np.copy(np.nanmean(Xa_outputi,axis=3).reshape(dum_dmax,dum_jmax,dum_imax,recon_period_len))\n",
    "                                XaNC_mean[:,:,:,:] = Xam_temp\n",
    "\n",
    "                                XaNC_variance = nf.createVariable(nc_var_name+'_Xa_3d_variance', 'f4', ('zt','lat', 'lon','time'))\n",
    "                                Xa_temp = np.copy(np.ma.var(Xa_outputi,axis=3).reshape(dum_dmax,dum_jmax,dum_imax,recon_period_len))\n",
    "                                Xa_temp = np.ma.MaskedArray(Xa_temp, Xam_temp >= 9.9692e+36)\n",
    "                                XaNC_variance[:,:,:,:] = Xa_temp\n",
    "\n",
    "                                for reconii in range(recon_period_len):\n",
    "                                    XaNC_mean_i = XaNC_mean[:,:,:,reconii]\n",
    "                                    XaNC_var_i = XaNC_variance[:,:,:,reconii]\n",
    "                                    if log_level > 0:\n",
    "                                        print('>>      Recon {}. Xa mean is {:.8f}, std is {:.8f}, var is {:.8f}'.format(reconii, np.nanmean(XaNC_mean_i), np.sqrt(np.nanmean(XaNC_var_i)), np.nanmean(XaNC_var_i)))\n",
    "\n",
    "                                if save_ens_full:\n",
    "                                    XaNC_full = nf.createVariable(nc_var_name+'_Xa_3d_full', 'f4', ('zt','lat', 'lon', 'nens', 'time'))\n",
    "                                    XaNC_full[:,:,:,:,:] = np.copy(Xa_outputi.reshape(dum_dmax,dum_jmax,dum_imax,nens,recon_period_len))\n",
    "\n",
    "                                    XbNC_full = nf.createVariable(nc_var_name+'_Xb_3d_full', 'f4', ('zt','lat', 'lon', 'nens'))\n",
    "                                    XbNC_full[:,:,:,:] = np.copy(Xb0_i.reshape(dum_dmax,dum_jmax,dum_imax,nens))\n",
    "\n",
    "                                if kcov_saving > 0:\n",
    "                                    kcov_i = np.copy(kcov[lenn1:lenn1+dum_ijmax*dum_dmax]).reshape(dum_dmax,dum_jmax,dum_imax)\n",
    "                                    kcov_i = np.ma.MaskedArray(kcov_i, np.copy(xbm.reshape(dum_dmax,dum_jmax,dum_imax)) >= 9.9692e+36)\n",
    "                                    cov_ob0 = nf.createVariable(nc_var_name+'_3d_obs0'+'_cov', 'f4', ( 'zt', 'lat','lon'))\n",
    "                                    cov_ob0[:,:,:] = kcov_i\n",
    "\n",
    "                                #Add local attributes to variable instances\n",
    "                                longitude.units = '°'\n",
    "                                latitude.units = '°'\n",
    "                                levels.units = 'm'\n",
    "                                #XbNC_mean.units = '°C'\n",
    "                                #XbNC_variance.units = '°C^2'\n",
    "                                #if save_ens_full:\n",
    "                                #    XaNC_full.units = '°C'\n",
    "                                #    XbNC_full.units = '°C'\n",
    "                        # Closing the dataset\n",
    "                        nf.close()  # close the new file\n",
    "                        if log_level > 1:\n",
    "                            print('')\n",
    "                            print('>>  Data saved in netCDF file:')\n",
    "                        if log_level > 1:\n",
    "                            print('')\n",
    "                            print('netCDF file saved : ')\n",
    "                            print('')\n",
    "                            print(nc_filename)\n",
    "                    \n",
    "                    else:\n",
    "                        if MCi>0:\n",
    "                            file1 = hdf5name\n",
    "                            # Check and remove the files if they exist\n",
    "                            if os.path.exists(file1):\n",
    "                                os.remove(file1)\n",
    "                                if log_level > 2:\n",
    "                                    print(f\"{file1} has been removed.\")\n",
    "                            else:\n",
    "                                if log_level > 2:\n",
    "                                    print(f\"{file1} does not exist.\")\n",
    "\n",
    "                        \n",
    "                    if log_level > 1:\n",
    "                        print('')\n",
    "                        print('##########                    ##########')\n",
    "                        print('##########   This loop done   ##########')\n",
    "                        print('##########                    ##########')\n",
    "                        print('')\n",
    "                        print('')\n",
    "\n",
    "            \n",
    "            \n",
    "            # SAVE all MC result\n",
    "            # change index and reshape\n",
    "            with h5py.File(hdf5name0, 'r') as f:\n",
    "                Xb = f.get('Xb')  # read Xb, background 2d field data\n",
    "                Xb3d = f.get('Xb3d')  # read Xb, background 3d field data order: lon-lat-depth\n",
    "                \n",
    "                if log_level > 4:\n",
    "                    print(Xa_shape_MC.shape)\n",
    "                if Xb is not None:\n",
    "                    Xa_MC_s = np.swapaxes(Xa_MC,2,3)\n",
    "                    Xa_MC_s = Xa_MC_s.reshape((prior_variable2d_len*dum_jmax*dum_imax, prior_len*MCn, recon_period_len))\n",
    "                if Xb3d is not None:\n",
    "                    Xa3d_MC_s = np.swapaxes(Xa3d_MC,2,3)\n",
    "                    Xa3d_MC_s = Xa3d_MC_s.reshape((prior_variable3d_len*dum_dmax*dum_jmax*dum_imax, prior_len*MCn, recon_period_len))\n",
    "\n",
    "                nc_filename_all = MC_dir + ''.join(filename_short) + 'all.nc'\n",
    "\n",
    "                if log_level > 1:\n",
    "                    print('>>  Start writing netCDF for all Monte Carlo simulations...')\n",
    "                # save netCDF file\n",
    "                nf = Dataset(nc_filename_all, 'w', format='NETCDF4')\n",
    "                nf.description = 'DeepDA' + nc_filename_all\n",
    "                #Specifying dimensions\n",
    "                nf.createDimension('lon', len(cGENIEGridB_lon36))\n",
    "                nf.createDimension('lat', len(cGENIEGridB_lat36))\n",
    "                z = np.arange(0,1,1) # level 2d\n",
    "                nf.createDimension('z', len(z))  # level\n",
    "                nf.createDimension('nens', nens)  # number of ens   ###### revision\n",
    "                nf.createDimension('MC', nens*MCn)  # number of ens   ###### revision\n",
    "                nf.createDimension('time', recon_period_len)\n",
    "                # Building variables\n",
    "                longitude = nf.createVariable('Longitude', 'f4', 'lon')\n",
    "                # Passing data into variables\n",
    "                longitude[:] = cGENIEGridB_lon36.values\n",
    "\n",
    "                latitude = nf.createVariable('Latitude', 'f4', 'lat')\n",
    "                latitude[:] = cGENIEGridB_lat36.values\n",
    "\n",
    "                levels = nf.createVariable('Levels', 'i4', 'z')\n",
    "                levels[:] = z  # 2d level\n",
    "                if Xb3d is not None:\n",
    "                    nf.createDimension('zt', len(zt))\n",
    "                    levels = nf.createVariable('zt', 'f4', 'zt')\n",
    "                    levels[:] = zt\n",
    "\n",
    "                if locRad:\n",
    "                    #nf.createDimension('prior_var', prior_variable_len)  # level\n",
    "                    covloc_nc = nf.createVariable('covloc', 'f4', ('lat', 'lon'))\n",
    "                    covloc_nc[:,:] = np.copy(covloc[0:dum_ijmax].reshape(dum_jmax,dum_imax))\n",
    "\n",
    "                if Xb is not None:\n",
    "                    if log_level > 2:\n",
    "                        print('Writing 2d field.')\n",
    "                    for nc_var_i in range(prior_variable_len):\n",
    "                        nc_var_name = prior_variable_dict[nc_var_i]\n",
    "\n",
    "                        j0 = dum_ijmax * nc_var_i\n",
    "                        j1 = dum_ijmax * (nc_var_i+1)\n",
    "                        if log_level > 0:\n",
    "                            print('')                            \n",
    "                            print('>>    ID from {} to {}: field is {}'.format(j0, j1,nc_var_name))\n",
    "\n",
    "                        Xb0_i = np.copy(Xb[j0:j1,:])\n",
    "\n",
    "                        Xa_output_i = np.copy(Xa_MC_s[j0:j1,:,:])\n",
    "                        Xa_outputi = Xa_output_i.reshape(dum_imax,dum_jmax,nens*MCn,recon_period_len)\n",
    "\n",
    "                        XbNC_mean = nf.createVariable(nc_var_name+'_Xb_mean', 'f4', ('lat', 'lon','z'))\n",
    "                        xbm = np.mean(Xb0_i,axis=1)\n",
    "                        XbNC_mean[:,:,:] = np.copy(xbm.reshape(dum_jmax,dum_imax,1))\n",
    "\n",
    "                        XbNC_variance = nf.createVariable(nc_var_name+'_Xb_variance', 'f4', ('lat', 'lon','z'))\n",
    "                        Xb_temp = np.copy(np.var(Xb0_i,axis=1).reshape(dum_jmax,dum_imax,1))\n",
    "                        Xb_temp = np.ma.MaskedArray(Xb_temp, np.copy(xbm.reshape(dum_jmax,dum_imax,1)) >= 9.9692e+36)\n",
    "                        XbNC_variance[:,:,:] = Xb_temp\n",
    "                        if log_level > 0:\n",
    "                            print('>>               Xb mean is {:.8f}, std is {:.8f}, var is {:.8f}'.format(np.nanmean(XbNC_mean), np.sqrt(np.nanmean(Xb_temp)), np.nanmean(Xb_temp)))\n",
    "\n",
    "                        XaNC_mean = nf.createVariable(nc_var_name+'_Xa_mean', 'f4', ('lat', 'lon','z','time'))\n",
    "                        Xam_temp = np.copy(np.nanmean(Xa_outputi,axis=2).reshape(dum_jmax,dum_imax,1,recon_period_len))\n",
    "                        XaNC_mean[:,:,:,:] = Xam_temp\n",
    "\n",
    "                        XaNC_variance = nf.createVariable(nc_var_name+'_Xa_variance', 'f4', ('lat', 'lon','z','time'))\n",
    "                        if log_level > 4:\n",
    "                            print(Xa_outputi[0,0:36,0,0])\n",
    "                        Xa_temp = np.copy(np.ma.var(Xa_outputi,axis=2).reshape(dum_jmax,dum_imax,1,recon_period_len))\n",
    "                        Xa_temp = np.ma.MaskedArray(Xa_temp, Xam_temp >= 9.9692e+36)\n",
    "                        XaNC_variance[:,:,:,:] = Xa_temp\n",
    "\n",
    "                        for reconii in range(recon_period_len):\n",
    "                            XaNC_mean_i = XaNC_mean[:,:,:,reconii]\n",
    "                            XaNC_var_i = XaNC_variance[:,:,:,reconii]\n",
    "                            if log_level > 0:\n",
    "                                print('>>      Recon {}. Xa mean is {:.8f}, std is {:.8f}, var is {:.8f}'.format(reconii, np.nanmean(XaNC_mean_i),np.sqrt(np.nanmean(XaNC_var_i)), np.nanmean(XaNC_var_i)))\n",
    "\n",
    "                        if save_ens_full:\n",
    "                            XaNC_full = nf.createVariable(nc_var_name+'_Xa_full', 'f4', ('lat', 'lon', 'MC', 'z','time'))\n",
    "                            XaNC_full[:,:,:,:,:] = np.copy(Xa_outputi.reshape(dum_jmax,dum_imax,nens*MCn,1,recon_period_len))\n",
    "\n",
    "                            XbNC_full = nf.createVariable(nc_var_name+'_Xb_full', 'f4', ('lat', 'lon', 'nens', 'z'))\n",
    "                            XbNC_full[:,:,:,:] = np.copy(Xb0_i.reshape(dum_jmax,dum_imax,nens,1))\n",
    "\n",
    "                        if kcov_saving > 0:\n",
    "                            kcov_i = np.copy(kcov[j0:j1]).reshape(dum_imax,dum_jmax,1)\n",
    "                            kcov_i = np.ma.MaskedArray(kcov_i, np.copy(xbm.reshape(dum_jmax,dum_imax,1)) >= 9.9692e+36)\n",
    "                            cov_ob0 = nf.createVariable(nc_var_name+'_obs0'+'_cov', 'f4', ('lat', 'lon','z'))\n",
    "                            cov_ob0[:,:,:] = kcov_i\n",
    "\n",
    "                        #Add local attributes to variable instances\n",
    "                        longitude.units = '°'\n",
    "                        latitude.units = '°'\n",
    "                        levels.units = 'm'\n",
    "                        #XbNC_mean.units = '°C'\n",
    "                        #XbNC_variance.units = '°C^2'\n",
    "                        #if save_ens_full:\n",
    "                        #    XaNC_full.units = '°C'\n",
    "                        #    XbNC_full.units = '°C'\n",
    "\n",
    "                        #variance.warning = 'test ...'\n",
    "                if Xb3d is not None:\n",
    "                    if log_level > 1:\n",
    "                        print('Writing 3d field.')\n",
    "                    for nc_var_i in range(prior_variable3d_len):\n",
    "                        nc_var_name = prior_variable_dict_3d[nc_var_i]\n",
    "\n",
    "                        j0 = dum_ijmax * dum_dmax * nc_var_i\n",
    "                        j1 = dum_ijmax * dum_dmax * (nc_var_i+1)\n",
    "                        if log_level > 0:\n",
    "                            print('')\n",
    "                            print('>>    ID from {} to {}: field is {}'.format(j0, j1,nc_var_name))\n",
    "\n",
    "                        Xb0_i = np.copy(Xb3d[j0:j1,:])\n",
    "                        Xa_output_i = np.copy(Xa3d_MC_s[j0:j1,:,:])\n",
    "                        Xa_outputi = Xa_output_i.reshape(dum_imax, dum_jmax,dum_dmax, nens*MCn,recon_period_len)\n",
    "\n",
    "                        XbNC_mean = nf.createVariable(nc_var_name+'_Xb_3d_mean', 'f4', ( 'zt', 'lat','lon'))\n",
    "                        xbm = np.mean(Xb0_i,axis=1)\n",
    "                        XbNC_mean[:,:,:] = np.copy(xbm.reshape(dum_dmax,dum_jmax,dum_imax))\n",
    "\n",
    "                        XbNC_variance = nf.createVariable(nc_var_name+'_Xb_3d_variance', 'f4', ( 'zt', 'lat','lon'))\n",
    "                        Xb_temp = np.copy(np.var(Xb0_i,axis= 1).reshape(dum_dmax,dum_jmax,dum_imax))\n",
    "                        Xb_temp = np.ma.MaskedArray(Xb_temp, np.copy(xbm.reshape(dum_dmax,dum_jmax,dum_imax)) >= 9.9692e+36)\n",
    "                        XbNC_variance[:,:,:] = Xb_temp\n",
    "                        if log_level > 0:\n",
    "                            print('>>               Xb mean is {:.8f}, std is {:.8f}, var is {:.8f}'.format(np.nanmean(XbNC_mean),np.sqrt(np.nanmean(Xb_temp)), np.nanmean(Xb_temp)))\n",
    "\n",
    "                        XaNC_mean = nf.createVariable(nc_var_name+'_Xa_3d_mean', 'f4', ('zt', 'lat', 'lon','time'))\n",
    "                        Xam_temp = np.copy(np.nanmean(Xa_outputi,axis=3).reshape(dum_dmax,dum_jmax,dum_imax,recon_period_len))\n",
    "                        XaNC_mean[:,:,:,:] = Xam_temp\n",
    "\n",
    "                        XaNC_variance = nf.createVariable(nc_var_name+'_Xa_3d_variance', 'f4', ('zt','lat', 'lon','time'))\n",
    "                        Xa_temp = np.copy(np.ma.var(Xa_outputi,axis=3).reshape(dum_dmax,dum_jmax,dum_imax,recon_period_len))\n",
    "                        Xa_temp = np.ma.MaskedArray(Xa_temp, Xam_temp >= 9.9692e+36)\n",
    "                        XaNC_variance[:,:,:,:] = Xa_temp\n",
    "\n",
    "                        for reconii in range(recon_period_len):\n",
    "                            XaNC_mean_i = XaNC_mean[:,:,:,reconii]\n",
    "                            XaNC_var_i = XaNC_variance[:,:,:,reconii]\n",
    "                            if log_level > 0:\n",
    "                                print('>>      Recon {}. Xa mean is {:.8f}, std is {:.8f}, var is {:.8f}'.format(reconii, np.nanmean(XaNC_mean_i), np.sqrt(np.nanmean(XaNC_var_i)), np.nanmean(XaNC_var_i)))\n",
    "\n",
    "                        if save_ens_full:\n",
    "                            XaNC_full = nf.createVariable(nc_var_name+'_Xa_3d_full', 'f4', ('zt','lat', 'lon', 'MC', 'time'))\n",
    "                            XaNC_full[:,:,:,:,:] = np.copy(Xa_outputi.reshape(dum_dmax,dum_jmax,dum_imax,nens*MCn,recon_period_len))\n",
    "\n",
    "                            XbNC_full = nf.createVariable(nc_var_name+'_Xb_3d_full', 'f4', ('zt','lat', 'lon', 'nens'))\n",
    "                            XbNC_full[:,:,:,:] = np.copy(Xb0_i.reshape(dum_dmax,dum_jmax,dum_imax,nens))\n",
    "\n",
    "                        if kcov_saving > 0:\n",
    "                            kcov_i = np.copy(kcov[lenn1:lenn1+dum_ijmax*dum_dmax]).reshape(dum_dmax,dum_jmax,dum_imax)\n",
    "                            kcov_i = np.ma.MaskedArray(kcov_i, np.copy(xbm.reshape(dum_dmax,dum_jmax,dum_imax)) >= 9.9692e+36)\n",
    "                            cov_ob0 = nf.createVariable(nc_var_name+'_3d_obs0'+'_cov', 'f4', ( 'zt', 'lat','lon'))\n",
    "                            cov_ob0[:,:,:] = kcov_i\n",
    "\n",
    "                        #Add local attributes to variable instances\n",
    "                        longitude.units = '°'\n",
    "                        latitude.units = '°'\n",
    "                        levels.units = 'm'\n",
    "                        #XbNC_mean.units = '°C'\n",
    "                        #XbNC_variance.units = '°C^2'\n",
    "                        #if save_ens_full:\n",
    "                        #    XaNC_full.units = '°C'\n",
    "                        #    XbNC_full.units = '°C'\n",
    "                # Closing the dataset\n",
    "                nf.close()  # close the new file\n",
    "                if log_level > 1:\n",
    "                    print('')\n",
    "                    print('>>  Data saved in netCDF file:')\n",
    "                if log_level > 1:\n",
    "                    print('')\n",
    "                    print('netCDF file saved : ')\n",
    "                    print('')\n",
    "                    print(nc_filename_all)\n",
    "                    print('')\n",
    "                    print('##########                    ##########')\n",
    "                    print('##########   All loop done   ##########')\n",
    "                    print('##########                    ##########')\n",
    "                    print('')\n",
    "                    print('')\n",
    "# remove all hdf5 if required\n",
    "if not save_mc_full:\n",
    "    file1 = hdf5name0\n",
    "    # Check and remove the files if they exist\n",
    "    if os.path.exists(file1):\n",
    "        os.remove(file1)\n",
    "        if log_level > 2:\n",
    "            print(f\"{file1} has been removed.\")\n",
    "    else:\n",
    "        if log_level > 2:\n",
    "            print(f\"{file1} does not exist.\")\n",
    "\n",
    "# export jupyter notebook as html, for reference\n",
    "os.system('jupyter nbconvert --to html DeepDA_allMC.ipynb')\n",
    "shutil.move(\"DeepDA_allMC.html\", MC_dir+\"DeepDA_allMC.html\")\n",
    "if log_level > 0:\n",
    "    print('')\n",
    "    print('########## All Done ##########')\n",
    "    print('')\n",
    "    print('This web page saved as DeepDA_allMC.html in the working directory : {}'.format(MC_dir))\n",
    "########## Check the consistency of the config.yml file and proxy database ##########"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.6.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}