MWP-1A_sources_inversion.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 1. Basic Set Up"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "#load modules\n",
    "import numpy as np #version 1.18.1\n",
    "import pandas as pd #version 1.0.1\n",
    "from scipy.optimize import nnls #version 1.4.1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "#load RSL rise magnitude at each site\n",
    "#------------------Empirical Senario------------------\n",
    "Emp_mag = pd.read_excel('RSL magnitude/Lin et al., Sup data.xlsx',sheet_name=1) \n",
    "Tahiti_emp_dis = Emp_mag['Tahiti']\n",
    "Tahiti_emp_std = np.std(Tahiti_emp_dis)\n",
    "Barbados_emp_dis =  Emp_mag['Barbados']\n",
    "Barbados_emp_std = np.std(Barbados_emp_dis)\n",
    "Sunda_emp_dis =  Emp_mag['Sunda Shelf']\n",
    "Sunda_emp_std = np.std(Sunda_emp_dis)\n",
    "HYD_emp_dis =  Emp_mag['Hydrographer\\'s Passage (HYD)']\n",
    "HYD_emp_std = np.std(HYD_emp_dis)\n",
    "NOG_emp_dis =  Emp_mag['Noggin Pass (NOG)']\n",
    "NOG_emp_std = np.std(NOG_emp_dis)\n",
    "Scot_emp_dis =  Emp_mag['Northwest Scotland']\n",
    "Scot_emp_std = np.std(Scot_emp_dis)\n",
    "total_emp_std = np.array([Tahiti_emp_std,Barbados_emp_std,Sunda_emp_std,HYD_emp_std,NOG_emp_std,\n",
    "                        Scot_emp_std]) \n",
    "#------------------Uniform Senario------------------\n",
    "Uni_mag = pd.read_excel('RSL magnitude/Lin et al., Sup data.xlsx',sheet_name=2) #Uniform Senario\n",
    "Tahiti_uni_dis = Uni_mag['Tahiti']\n",
    "Tahiti_uni_std = np.std(Tahiti_uni_dis)\n",
    "Barbados_uni_dis =  Uni_mag['Barbados']\n",
    "Barbados_uni_std = np.std(Barbados_uni_dis)\n",
    "Sunda_uni_dis =  Uni_mag['Sunda Shelf']\n",
    "Sunda_uni_std = np.std(Sunda_uni_dis)\n",
    "HYD_uni_dis =  Uni_mag['Hydrographer\\'s Passage (HYD)']\n",
    "HYD_uni_std = np.std(HYD_uni_dis)\n",
    "NOG_uni_dis =  Uni_mag['Noggin Pass (NOG)']\n",
    "NOG_uni_std = np.std(NOG_uni_dis)\n",
    "Scot_uni_dis =  Uni_mag['Northwest Scotland']\n",
    "Scot_uni_std = np.std(Scot_uni_dis)\n",
    "total_uni_std = np.array([Tahiti_uni_std,Barbados_uni_std,Sunda_uni_std,HYD_uni_std,NOG_uni_std,\n",
    "                        Scot_uni_std])\n",
    "#load sea-level fingerprint matrix\n",
    "#each row corresponds to each sea-level site and each column corresponds to each ice sheet\n",
    "A =  np.array(pd.read_csv('RSL magnitude/fingerprint.csv',index_col=0))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def find_maximum_prob(values,percent):\n",
    "    '''This function is used to analytically find the smallest confidence interval  \n",
    "    range from a data ensemble \n",
    "    values:  data samples\n",
    "    percent: probablity range to find (e.g., 95 is 95% CI)'''\n",
    "    \n",
    "    v_max,v_min = np.max(values),np.min(values)\n",
    "    ranges = []\n",
    "    for i in np.linspace(0,100-percent,1000):\n",
    "        _b,_e = np.percentile(values,[i,i+percent])\n",
    "        ranges.append([_b,_e,np.abs(_b-_e),i,i+percent])\n",
    "    ranges = np.array(ranges)\n",
    "    min_index = np.argmin(ranges[:,2])\n",
    "    \n",
    "    return ranges[min_index]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 2. MWP-1A Sources Inversion\n",
    "\n",
    "All loaded 20,000 MWP-1A RSL rise magnitudes samples represent the elastic induced relative sea-level rise at each location, which can be used to invert the MWP-1A sources following the equation:\n",
    "\n",
    "\n",
    "\n",
    "$$ \\Delta {RSL}_{Elastic}(\\varphi)  = {ESL}_{NAIS}  \\times F_{NAIS}(\\varphi) + {ESL}_{AIS} \\times F_{AIS}(\\varphi) +{ESL}_{SIS}\\times F_{SIS}(\\varphi)$$\n",
    "\n",
    "where $\\Delta {RSL}_{Elastic}(\\varphi)$ represent the instantaneous response of the solid Earth and sea surface to an influx of meltwater at location $\\varphi$; three $ESL$ terms represent eustatic sea-level (ESL) contributions from the NAIS, AIS, and SIS ice sheets and the $F(\\varphi)$ terms denote ice-sheet-specific, site-specific sea-level fingerprint values.\n",
    "\n",
    "\n",
    "By optimising this equation above using the weighted non-negative least sqaure approch based on Lawson-Hanson algorithm (from a [scipy](https://docs.scipy.org/doc/scipy/reference/index.html) function [nnls](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.nnls.html)), three optimised $ESL$ paramters in the equation above can be obtained. We then repeated this process for 20,000 times (i.e., Monte Carlo sampling) to calculate the probability density functions for these ice sheets contribution to MWP-1A as well as the total MWP-1A magnitude. \n",
    "\n",
    "\n",
    "### Empirical senario for coral records\n",
    "\n",
    "#### ALL INVERSION RESULTS BELOW ARE EXPRESSED IN UNIT METER\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "----------------MWP1A Magnitude & Sources with Empirical Distribution -\n",
      "----------------Mean Result (95% confidence interval)----------------\n",
      "NAIS: 12.225 [0.000,18.573]\n",
      "AIS: 2.574 [0.000,10.764]\n",
      "SIS:   3.481 [0.000,6.582]\n",
      "Total MWP-1A Magnitude:   18.280 [14.377,22.068]\n"
     ]
    }
   ],
   "source": [
    "emp_result = []\n",
    "while len(emp_result)<20000:\n",
    "    i = np.random.randint(0,20000,1)[0]\n",
    "    obs = np.array([Tahiti_emp_dis.iloc[i],Barbados_emp_dis.iloc[i],Sunda_emp_dis.iloc[i],HYD_emp_dis.iloc[i],\n",
    "                   NOG_emp_dis.iloc[i],Scot_emp_dis.iloc[i]]) #compile all sample together\n",
    "    opti = list(nnls(np.sqrt(1/np.array(total_emp_std)**2)[:,None] * A,np.sqrt(1/np.array(total_emp_std)**2)*obs)[0]) #non-negative least square\n",
    "    emp_result.append(opti) #append results from this sample\n",
    "    \n",
    "emp_result= np.array(emp_result)\n",
    "print('----------------MWP1A Magnitude & Sources with Empirical Distribution -')\n",
    "print('----------------Mean Result (95% confidence interval)----------------')\n",
    "print('NAIS: {0:5.3f} [{1:5.3f},{2:5.3f}]'.format(np.mean(emp_result[:,0]),*find_maximum_prob(emp_result[:,0],95)[:2]))\n",
    "print('AIS: {0:5.3f} [{1:5.3f},{2:5.3f}]'.format(np.mean(emp_result[:,1]),*find_maximum_prob(emp_result[:,1],95)[:2]))\n",
    "print('SIS:   {0:5.3f} [{1:5.3f},{2:5.3f}]'.format(np.mean(emp_result[:,2]),*find_maximum_prob(emp_result[:,2],95)[:2]))\n",
    "print('Total MWP-1A Magnitude:   {0:5.3f} [{1:5.3f},{2:5.3f}]'.format(np.mean(np.sum(emp_result,axis=1)),*find_maximum_prob(np.sum(emp_result,axis=1),95)[:2]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Uniform senario coral records"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "----------------MWP1A Magnitude & Sources with Uniform Distribution -\n",
      "----------------Mean Result (95% confidence interval)----------------\n",
      "NAIS: 12.416 [2.956,19.109]\n",
      "AIS: 2.170 [0.000,8.775]\n",
      "SIS:   3.137 [0.000,6.007]\n",
      "Total MWP-1A Magnitude:   17.723 [14.851,20.490]\n"
     ]
    }
   ],
   "source": [
    "uni_result = []\n",
    "while len(uni_result)<20000:\n",
    "    i = np.random.randint(0,20000,1)[0]\n",
    "    obs = np.array([Tahiti_uni_dis.iloc[i],Barbados_uni_dis.iloc[i],Sunda_uni_dis.iloc[i],HYD_uni_dis.iloc[i],\n",
    "                   NOG_uni_dis.iloc[i],Scot_uni_dis.iloc[i]]) #compile all sample together\n",
    "    opti = list(nnls(np.sqrt(1/np.array(total_uni_std)**2)[:,None] * A,np.sqrt(1/np.array(total_uni_std)**2)*obs)[0]) #non-negative least square\n",
    "    uni_result.append(opti)\n",
    "    \n",
    "uni_result= np.array(uni_result)\n",
    "print('----------------MWP1A Magnitude & Sources with Uniform Distribution -')\n",
    "print('----------------Mean Result (95% confidence interval)----------------')\n",
    "print('NAIS: {0:5.3f} [{1:5.3f},{2:5.3f}]'.format(np.mean(uni_result[:,0]),*find_maximum_prob(uni_result[:,0],95)[:2]))\n",
    "print('AIS: {0:5.3f} [{1:5.3f},{2:5.3f}]'.format(np.mean(uni_result[:,1]),*find_maximum_prob(uni_result[:,1],95)[:2]))\n",
    "print('SIS:   {0:5.3f} [{1:5.3f},{2:5.3f}]'.format(np.mean(uni_result[:,2]),*find_maximum_prob(uni_result[:,2],95)[:2]))\n",
    "print('Total MWP-1A Magnitude:   {0:5.3f} [{1:5.3f},{2:5.3f}]'.format(np.mean(np.sum(uni_result,axis=1)),*find_maximum_prob(np.sum(uni_result,axis=1),95)[:2]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 3. Jackknife Resampling\n",
    "\n",
    "Jackknife resampling approach (also called leave one out cross validation) is used to quantify the overall stability of our inversion results as well as the potential bias contained in our original inversion results. \n",
    "\n",
    "### Empirical senario for coral records"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "----------------Without Tahiti -------------------------------------\n",
      "----------------MWP1A Magnitude & Sources with Empirical Distribution -----------\n",
      "----------------Mean Result [95% confidence interval] (1 std)-----------------\n",
      "North American Ice Sheet: 11.500 [0.000,19.039]  (6.201)\n",
      "West Antarctic Ice Sheet: 3.005 [0.000,11.753] (4.378)\n",
      "Scandinavian Ice Sheet:   3.472 [0.000,6.571]  (2.010)\n",
      "----------------------------------------------------------------------\n",
      "----------------Without Barbados -------------------------------------\n",
      "----------------MWP1A Magnitude & Sources with Empirical Distribution -----------\n",
      "----------------Mean Result [95% confidence interval] (1 std)-----------------\n",
      "North American Ice Sheet: 7.559 [0.000,18.042]  (7.282)\n",
      "West Antarctic Ice Sheet: 6.844 [0.000,15.840] (6.249)\n",
      "Scandinavian Ice Sheet:   5.190 [0.000,9.995]  (2.896)\n",
      "----------------------------------------------------------------------\n",
      "----------------Without Sunda Shelf -------------------------------------\n",
      "----------------MWP1A Magnitude & Sources with Empirical Distribution -----------\n",
      "----------------Mean Result [95% confidence interval] (1 std)-----------------\n",
      "North American Ice Sheet: 12.566 [0.000,18.604]  (4.983)\n",
      "West Antarctic Ice Sheet: 2.356 [0.000,10.458] (3.594)\n",
      "Scandinavian Ice Sheet:   3.480 [0.000,6.585]  (2.019)\n",
      "----------------------------------------------------------------------\n",
      "----------------Without Hydrographer -------------------------------------\n",
      "----------------MWP1A Magnitude & Sources with Empirical Distribution -----------\n",
      "----------------Mean Result [95% confidence interval] (1 std)-----------------\n",
      "North American Ice Sheet: 11.612 [0.000,18.318]  (5.346)\n",
      "West Antarctic Ice Sheet: 2.897 [0.000,11.067] (3.896)\n",
      "Scandinavian Ice Sheet:   3.429 [0.000,6.513]  (1.998)\n",
      "----------------------------------------------------------------------\n",
      "----------------Without Noggin Pass -------------------------------------\n",
      "----------------MWP1A Magnitude & Sources with Empirical Distribution -----------\n",
      "----------------Mean Result [95% confidence interval] (1 std)-----------------\n",
      "North American Ice Sheet: 11.961 [0.000,18.414]  (5.150)\n",
      "West Antarctic Ice Sheet: 2.711 [0.000,10.780] (3.759)\n",
      "Scandinavian Ice Sheet:   3.450 [0.000,6.480]  (1.991)\n",
      "----------------------------------------------------------------------\n",
      "----------------Without Scotland -------------------------------------\n",
      "----------------MWP1A Magnitude & Sources with Empirical Distribution -----------\n",
      "----------------Mean Result [95% confidence interval] (1 std)-----------------\n",
      "North American Ice Sheet: 13.962 [0.000,21.003]  (6.133)\n",
      "West Antarctic Ice Sheet: 2.911 [0.000,13.207] (4.477)\n",
      "Scandinavian Ice Sheet:   1.429 [0.000,10.709]  (3.665)\n",
      "----------------------------------------------------------------------\n",
      "----------------Overall Jackknife/Original Results --------------------\n",
      "NAIS: 11.527, 12.225  Bias: -0.698 Bias Corrected: 12.923 m\n",
      "AIS: 3.454, 2.574  Bias: 0.880 Bias Corrected: 1.694 m\n",
      "SIS:   3.408, 3.481  Bias: -0.073 Bias Corrected: 3.554 m\n",
      "----------------------------------------------------------------------\n"
     ]
    }
   ],
   "source": [
    "jack_knife_result =np.zeros([60000,3])\n",
    "site_name = ['Tahiti',\"Barbados\",'Sunda Shelf','Hydrographer','Noggin Pass',\n",
    "            'Scotland']\n",
    "for jack in range(6):\n",
    "    \n",
    "    all_emp_result = [Tahiti_emp_dis,Barbados_emp_dis,Sunda_emp_dis,NOG_emp_dis,\n",
    "                           HYD_emp_dis,Scot_emp_dis]\n",
    "    total_emp_std = [Tahiti_emp_std,Barbados_emp_std,Sunda_emp_std,HYD_emp_std,NOG_emp_std,\n",
    "                         Scot_emp_std]\n",
    "    A_index = ~(np.arange(0,6)==jack)\n",
    "    A_jack = A[A_index]\n",
    "\n",
    "    del all_emp_result[jack] #remove one site's observation\n",
    "    del total_emp_std[jack] #remove one site's standard deviation \n",
    "    emp_test_result = []\n",
    "    all_emp_result = np.array(all_emp_result)\n",
    "    total_emp_std = np.array(total_emp_std)\n",
    "    while len(emp_test_result)<10000:\n",
    "        \n",
    "        i = np.random.randint(0,20000,1)[0]\n",
    "        obs = np.array([Tahiti_emp_dis.iloc[i],Barbados_emp_dis.iloc[i],Sunda_emp_dis.iloc[i],HYD_emp_dis.iloc[i],\n",
    "                       NOG_emp_dis.iloc[i],Scot_emp_dis.iloc[i]])\n",
    "        obs = obs[A_index] \n",
    "        opti = list(nnls(np.sqrt(1/np.array(total_emp_std)**2)[:,None] * A_jack,np.sqrt(1/np.array(total_emp_std)**2)*obs)[0])\n",
    "        emp_test_result.append(opti)\n",
    "       \n",
    "    \n",
    "    emp_test_result= np.array(emp_test_result)\n",
    "    jack_knife_result[jack*10000:(jack+1)*10000,:]=emp_test_result\n",
    "    print('----------------Without {:} -------------------------------------'.format(site_name[jack]))\n",
    "    print('----------------MWP1A Magnitude & Sources with Empirical Distribution -----------')\n",
    "    print('----------------Mean Result [95% confidence interval] (1 std)-----------------')\n",
    "    print('North American Ice Sheet: {0:5.3f} [{1:5.3f},{2:5.3f}]  ({3:5.3f})'.format(np.mean(emp_test_result[:,0]),find_maximum_prob(emp_test_result[:,0],95)[0],find_maximum_prob(emp_test_result[:,0],95)[1],\n",
    "                                                                                   np.std(emp_test_result[:,0])))\n",
    "    print('West Antarctic Ice Sheet: {0:5.3f} [{1:5.3f},{2:5.3f}] ({3:5.3f})'.format(np.mean(emp_test_result[:,1]),find_maximum_prob(emp_test_result[:,1],95)[0],find_maximum_prob(emp_test_result[:,1],95)[1],\n",
    "                                                                                  np.std(emp_test_result[:,1])))\n",
    "    print('Scandinavian Ice Sheet:   {0:5.3f} [{1:5.3f},{2:5.3f}]  ({3:5.3f})'.format(np.mean(emp_test_result[:,2]),find_maximum_prob(emp_test_result[:,2],95)[0],find_maximum_prob(emp_test_result[:,2],95)[1],\n",
    "                                                                                  np.std(emp_test_result[:,2])))\n",
    "    print('----------------------------------------------------------------------')\n",
    "\n",
    "print('----------------Overall Jackknife/Original Results --------------------')\n",
    "print('NAIS: {0:5.3f}, {1:5.3f}  Bias: {2:5.3f} Bias Corrected: {3:5.3f} m'.format(np.mean(jack_knife_result[:,0]),np.mean(emp_result[:,0]),\n",
    "                                                                           np.mean(jack_knife_result[:,0])-np.mean(emp_result[:,0]),\n",
    "                                                                            np.mean(emp_result[:,0])-(np.mean(jack_knife_result[:,0])-np.mean(emp_result[:,0]))))\n",
    "print('AIS: {0:5.3f}, {1:5.3f}  Bias: {2:5.3f} Bias Corrected: {3:5.3f} m'.format(np.mean(jack_knife_result[:,1]),np.mean(emp_result[:,1]),\n",
    "                                                                           np.mean(jack_knife_result[:,1])-np.mean(emp_result[:,1]),\n",
    "                                                                           np.mean(emp_result[:,1])-(np.mean(jack_knife_result[:,1])-np.mean(emp_result[:,1]))))\n",
    "print('SIS:   {0:5.3f}, {1:5.3f}  Bias: {2:5.3f} Bias Corrected: {3:5.3f} m'.format(np.mean(jack_knife_result[:,2]),np.mean(emp_result[:,2]),\n",
    "                                                                          np.mean(jack_knife_result[:,2])-np.mean(emp_result[:,2]),\n",
    "                                                                        np.mean(emp_result[:,2])-( np.mean(jack_knife_result[:,2])-np.mean(emp_result[:,2]))))\n",
    "print('----------------------------------------------------------------------')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Uniform senario coral records"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "----------------Without Tahiti -------------------------------------\n",
      "----------------MWP1A Magnitude & Sources with Uniform Distribution -----------\n",
      "----------------Mean Result [95% confidence interval] (1 std)-----------------\n",
      "North American Ice Sheet: 10.683 [0.000,18.344]  (4.205)\n",
      "West Antarctic Ice Sheet: 3.269 [0.000,11.533] (4.283)\n",
      "Scandinavian Ice Sheet:   3.204 [0.000,6.002]  (1.891)\n",
      "----------------------------------------------------------------------\n",
      "----------------Without Barbados -------------------------------------\n",
      "----------------MWP1A Magnitude & Sources with Uniform Distribution -----------\n",
      "----------------Mean Result [95% confidence interval] (1 std)-----------------\n",
      "North American Ice Sheet: 8.817 [0.000,17.160]  (4.205)\n",
      "West Antarctic Ice Sheet: 5.460 [0.000,15.313] (5.941)\n",
      "Scandinavian Ice Sheet:   4.375 [0.000,8.759]  (2.607)\n",
      "----------------------------------------------------------------------\n",
      "----------------Without Sunda Shelf -------------------------------------\n",
      "----------------MWP1A Magnitude & Sources with Uniform Distribution -----------\n",
      "----------------Mean Result [95% confidence interval] (1 std)-----------------\n",
      "North American Ice Sheet: 12.527 [3.245,19.315]  (4.205)\n",
      "West Antarctic Ice Sheet: 2.107 [0.000,8.827] (2.996)\n",
      "Scandinavian Ice Sheet:   3.155 [0.000,5.996]  (1.881)\n",
      "----------------------------------------------------------------------\n",
      "----------------Without Hydrographer -------------------------------------\n",
      "----------------MWP1A Magnitude & Sources with Uniform Distribution -----------\n",
      "----------------Mean Result [95% confidence interval] (1 std)-----------------\n",
      "North American Ice Sheet: 12.322 [2.654,18.897]  (4.205)\n",
      "West Antarctic Ice Sheet: 2.190 [0.000,8.681] (3.049)\n",
      "Scandinavian Ice Sheet:   3.153 [0.000,5.946]  (1.863)\n",
      "----------------------------------------------------------------------\n",
      "----------------Without Noggin Pass -------------------------------------\n",
      "----------------MWP1A Magnitude & Sources with Uniform Distribution -----------\n",
      "----------------Mean Result [95% confidence interval] (1 std)-----------------\n",
      "North American Ice Sheet: 12.327 [2.702,19.412]  (4.205)\n",
      "West Antarctic Ice Sheet: 2.237 [0.000,9.139] (3.132)\n",
      "Scandinavian Ice Sheet:   3.089 [0.000,5.956]  (1.905)\n",
      "----------------------------------------------------------------------\n",
      "----------------Without Scotland -------------------------------------\n",
      "----------------MWP1A Magnitude & Sources with Uniform Distribution -----------\n",
      "----------------Mean Result [95% confidence interval] (1 std)-----------------\n",
      "North American Ice Sheet: 13.713 [3.050,20.777]  (4.205)\n",
      "West Antarctic Ice Sheet: 2.226 [0.000,9.989] (3.546)\n",
      "Scandinavian Ice Sheet:   1.825 [0.000,11.352]  (3.822)\n",
      "----------------------------------------------------------------------\n",
      "----------------Overall Jackknife/Original Results --------------------\n",
      "NAIS: 11.732, 12.416  Bias: -0.684 Bias Corrected: 13.100 m\n",
      "AIS: 2.915, 2.170  Bias: 0.745 Bias Corrected: 1.424 m\n",
      "SIS:   3.134, 3.137  Bias: -0.004 Bias Corrected: 3.141 m\n",
      "----------------------------------------------------------------------\n"
     ]
    }
   ],
   "source": [
    "jack_knife_result =np.zeros([60000,3])\n",
    "site_name = ['Tahiti',\"Barbados\",'Sunda Shelf','Hydrographer','Noggin Pass',\n",
    "            'Scotland']\n",
    "for jack in range(6):\n",
    "    \n",
    "    all_uni_result = [Tahiti_uni_dis,Barbados_uni_dis,Sunda_uni_dis,NOG_uni_dis,\n",
    "                           HYD_uni_dis,Scot_uni_dis]\n",
    "    total_uni_std = [Tahiti_uni_std,Barbados_uni_std,Sunda_uni_std,HYD_uni_std,NOG_uni_std,\n",
    "                         Scot_uni_std]\n",
    "    A_index = ~(np.arange(0,6)==jack)\n",
    "    A_jack = A[A_index]\n",
    "\n",
    "    del all_uni_result[jack] #remove one site's observation\n",
    "    del total_uni_std[jack] #remove one site's standard error \n",
    "    uni_test_result = []\n",
    "    all_uni_result = np.array(all_uni_result)\n",
    "    total_uni_std = np.array(total_uni_std)\n",
    "    while len(uni_test_result)<10000:\n",
    "        \n",
    "        i = np.random.randint(0,20000,1)[0]\n",
    "        obs = np.array([Tahiti_uni_dis.iloc[i],Barbados_uni_dis.iloc[i],Sunda_uni_dis.iloc[i],HYD_uni_dis.iloc[i],\n",
    "                       NOG_uni_dis.iloc[i],Scot_uni_dis.iloc[i]])\n",
    "        obs = obs[A_index] \n",
    "        opti = list(nnls(np.sqrt(1/np.array(total_uni_std)**2)[:,None] * A_jack,np.sqrt(1/np.array(total_uni_std)**2)*obs)[0])\n",
    "        uni_test_result.append(opti)\n",
    "       \n",
    "    \n",
    "    uni_test_result= np.array(uni_test_result)\n",
    "    jack_knife_result[jack*10000:(jack+1)*10000,:]=uni_test_result\n",
    "    print('----------------Without {:} -------------------------------------'.format(site_name[jack]))\n",
    "    print('----------------MWP1A Magnitude & Sources with Uniform Distribution -----------')\n",
    "    print('----------------Mean Result [95% confidence interval] (1 std)-----------------')\n",
    "    print('North American Ice Sheet: {0:5.3f} [{1:5.3f},{2:5.3f}]  ({3:5.3f})'.format(np.mean(uni_test_result[:,0]),find_maximum_prob(uni_test_result[:,0],95)[0],find_maximum_prob(uni_test_result[:,0],95)[1],\n",
    "                                                                                   np.std(uni_result[:,0])))\n",
    "    print('West Antarctic Ice Sheet: {0:5.3f} [{1:5.3f},{2:5.3f}] ({3:5.3f})'.format(np.mean(uni_test_result[:,1]),find_maximum_prob(uni_test_result[:,1],95)[0],find_maximum_prob(uni_test_result[:,1],95)[1],\n",
    "                                                                                  np.std(uni_test_result[:,1])))\n",
    "    print('Scandinavian Ice Sheet:   {0:5.3f} [{1:5.3f},{2:5.3f}]  ({3:5.3f})'.format(np.mean(uni_test_result[:,2]),find_maximum_prob(uni_test_result[:,2],95)[0],find_maximum_prob(uni_test_result[:,2],95)[1],\n",
    "                                                                                  np.std(uni_test_result[:,2])))\n",
    "    print('----------------------------------------------------------------------')\n",
    "\n",
    "print('----------------Overall Jackknife/Original Results --------------------')\n",
    "print('NAIS: {0:5.3f}, {1:5.3f}  Bias: {2:5.3f} Bias Corrected: {3:5.3f} m'.format(np.mean(jack_knife_result[:,0]),np.mean(uni_result[:,0]),\n",
    "                                                                           np.mean(jack_knife_result[:,0])-np.mean(uni_result[:,0]),\n",
    "                                                                            np.mean(uni_result[:,0])-(np.mean(jack_knife_result[:,0])-np.mean(uni_result[:,0]))))\n",
    "print('AIS: {0:5.3f}, {1:5.3f}  Bias: {2:5.3f} Bias Corrected: {3:5.3f} m'.format(np.mean(jack_knife_result[:,1]),np.mean(uni_result[:,1]),\n",
    "                                                                           np.mean(jack_knife_result[:,1])-np.mean(uni_result[:,1]),\n",
    "                                                                           np.mean(uni_result[:,1])-(np.mean(jack_knife_result[:,1])-np.mean(uni_result[:,1]))))\n",
    "print('SIS:   {0:5.3f}, {1:5.3f}  Bias: {2:5.3f} Bias Corrected: {3:5.3f} m'.format(np.mean(jack_knife_result[:,2]),np.mean(uni_result[:,2]),\n",
    "                                                                          np.mean(jack_knife_result[:,2])-np.mean(uni_result[:,2]),\n",
    "                                                                        np.mean(uni_result[:,2])-( np.mean(jack_knife_result[:,2])-np.mean(uni_result[:,2]))))\n",
    "print('----------------------------------------------------------------------')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 4. MWP-1A Sources Inversion with sea-level oscillation limit\n",
    "\n",
    "### Empirical senario for coral records"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "----------------MWP1A Magnitude & Sources with Empirical Distribution -----------\n",
      "----------------Mean Result (95% confidence interval)----------------\n",
      "NAIS: 11.315 [0.000,15.649]\n",
      "AIS:  2.041 [0.000,9.396]\n",
      "SIS:  4.927 [2.466,7.188]\n",
      "Total MWP-1A Magnitude:    18.284 [14.503,21.803]\n"
     ]
    }
   ],
   "source": [
    "emp_result = []\n",
    "total_emp_std = np.array([Tahiti_emp_std,Barbados_emp_std,Sunda_emp_std,HYD_emp_std,NOG_emp_std,\n",
    "                        Scot_emp_std]) \n",
    "while len(emp_result)<20000:\n",
    "    i = np.random.randint(0,20000,1)[0]\n",
    "    obs = np.array([Tahiti_emp_dis.iloc[i],Barbados_emp_dis.iloc[i],Sunda_emp_dis.iloc[i],HYD_emp_dis.iloc[i],\n",
    "                   NOG_emp_dis.iloc[i],Scot_emp_dis.iloc[i]]) #compile all sample together\n",
    "    opti = list(nnls(np.sqrt(1/np.array(total_emp_std)**2)[:,None] * A,np.sqrt(1/np.array(total_emp_std)**2)*obs)[0]) #non-negative least square\n",
    "    if (opti[0]*0.75+opti[1]*1.09-opti[2]*0.74)<9: # filter out the results producing a sea-level oscillation\n",
    "        emp_result.append(opti)\n",
    "    \n",
    "emp_result= np.array(emp_result)\n",
    "print('----------------MWP1A Magnitude & Sources with Empirical Distribution -----------')\n",
    "print('----------------Mean Result (95% confidence interval)----------------')\n",
    "print('NAIS: {0:5.3f} [{1:5.3f},{2:5.3f}]'.format(np.mean(emp_result[:,0]),*find_maximum_prob(emp_result[:,0],95)[:2]))\n",
    "print('AIS:  {0:5.3f} [{1:5.3f},{2:5.3f}]'.format(np.mean(emp_result[:,1]),*find_maximum_prob(emp_result[:,1],95)[:2]))\n",
    "print('SIS:  {0:5.3f} [{1:5.3f},{2:5.3f}]'.format(np.mean(emp_result[:,2]),*find_maximum_prob(emp_result[:,2],95)[:2]))\n",
    "print('Total MWP-1A Magnitude:    {0:5.3f} [{1:5.3f},{2:5.3f}]'.format(np.mean(np.sum(emp_result,axis=1)),*find_maximum_prob(np.sum(emp_result,axis=1),95)[:2]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Uniform senario coral records"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "----------------MWP1A Magnitude & Sources with Uniform Distribution -----------\n",
      "----------------Mean Result (95% confidence interval)----------------\n",
      "NAIS: 11.422 [3.294,16.114]\n",
      "AIS:  1.762 [0.000,7.519]\n",
      "SIS:  4.559 [2.595,6.606]\n",
      "Total MWP-1A Magnitude: 17.736 [15.100,20.419]\n"
     ]
    }
   ],
   "source": [
    "uni_result = []\n",
    "total_uni_std = np.array([Tahiti_uni_std,Barbados_uni_std,Sunda_uni_std,HYD_uni_std,NOG_uni_std,\n",
    "                        Scot_uni_std])\n",
    "while len(uni_result)<20000:\n",
    "    i = np.random.randint(0,20000,1)[0]\n",
    "    obs = np.array([Tahiti_uni_dis.iloc[i],Barbados_uni_dis.iloc[i],Sunda_uni_dis.iloc[i],HYD_uni_dis.iloc[i],\n",
    "                   NOG_uni_dis.iloc[i],Scot_uni_dis.iloc[i]]) #compile all sample together\n",
    "    opti = list(nnls(np.sqrt(1/np.array(total_uni_std)**2)[:,None] * A,np.sqrt(1/np.array(total_uni_std)**2)*obs)[0]) #non-negative least square\n",
    "    if (opti[0]*0.75+opti[1]*1.09-opti[2]*0.74)<9: # filter out the results producing a sea-level oscillation\n",
    "\n",
    "        uni_result.append(opti)\n",
    "    \n",
    "uni_result= np.array(uni_result)\n",
    "print('----------------MWP1A Magnitude & Sources with Uniform Distribution -----------')\n",
    "print('----------------Mean Result (95% confidence interval)----------------')\n",
    "print('NAIS: {0:5.3f} [{1:5.3f},{2:5.3f}]'.format(np.mean(uni_result[:,0]),*find_maximum_prob(uni_result[:,0],95)[:2]))\n",
    "print('AIS:  {0:5.3f} [{1:5.3f},{2:5.3f}]'.format(np.mean(uni_result[:,1]),*find_maximum_prob(uni_result[:,1],95)[:2]))\n",
    "print('SIS:  {0:5.3f} [{1:5.3f},{2:5.3f}]'.format(np.mean(uni_result[:,2]),*find_maximum_prob(uni_result[:,2],95)[:2]))\n",
    "print('Total MWP-1A Magnitude: {0:5.3f} [{1:5.3f},{2:5.3f}]'.format(np.median(np.sum(uni_result,axis=1)),*find_maximum_prob(np.sum(uni_result,axis=1),95)[:2]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Averaged results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "----------------Averaged MWP1A Magnitude & Sources -----------\n",
      "----------------Median Result (95% confidence interval)----------------\n",
      "NAIS: 12.097 [5.752,15.468]\n",
      "AIS:  1.300 [0.000,5.831]\n",
      "SIS:  4.712 [3.233,6.366]\n",
      "Total MWP-1A Magnitude: 18.027 [15.100,20.419]\n"
     ]
    }
   ],
   "source": [
    "print('----------------Averaged MWP1A Magnitude & Sources -----------')\n",
    "print('----------------Median Result (95% confidence interval)----------------')\n",
    "print('NAIS: {0:5.3f} [{1:5.3f},{2:5.3f}]'.format(np.median((uni_result[:,0]+emp_result[:,0])/2),*find_maximum_prob((uni_result[:,0]+emp_result[:,0])/2,95)[:2]))\n",
    "print('AIS:  {0:5.3f} [{1:5.3f},{2:5.3f}]'.format(np.median((uni_result[:,1]+emp_result[:,1])/2),*find_maximum_prob((uni_result[:,1]+emp_result[:,1])/2,95)[:2]))\n",
    "print('SIS:  {0:5.3f} [{1:5.3f},{2:5.3f}]'.format(np.median((uni_result[:,2]+emp_result[:,2])/2),*find_maximum_prob((uni_result[:,2]+emp_result[:,2])/2,95)[:2]))\n",
    "print('Total MWP-1A Magnitude: {0:5.3f} [{1:5.3f},{2:5.3f}]'.format(np.median((np.sum(uni_result,axis=1)+np.sum(emp_result,axis=1))/2),*find_maximum_prob(np.sum(uni_result,axis=1),95)[:2]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Citation\n",
    "\n",
    "Lin, Y., Hibbert, F.D., Whitehouse, P.L. *et al.* A reconciled solution of Meltwater Pulse 1A sources using sea-level fingerprinting. *Nat Commun* **12**, 2015 (2021).\n",
    "\n",
    "https://doi.org/10.1038/s41467-021-21990-y"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Final notes:\n",
    "\n",
    "- Due to the different ramdom seeds, this code can possibly result in a slightly different results comparing to the results shown in our manuscript\n",
    "- If you have any questions, feedbacks or comments, please feel free to contact [the corresponding author](yc-lin.com): yucheng.lin@durham.ac.uk\n"
   ]
  }
 ],
 "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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}