<jupyter_start><jupyter_text>Import dependencies<jupyter_code>## Get dependencies ##

import numpy as np
import string
import math
import sys
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib
sys.path.append('..')
from GIR import *
import scipy as sp
import pickle
import time
import scipy as sp
from scipy import signal
from scipy.io.idl import readsav
import os
from mpl_toolkits.axes_grid1.inset_locator import inset_axes
import glob
import requests
import ftplib
import io
import cmocean
from bs4 import BeautifulSoup
import urllib.request
from io import StringIO, BytesIO
from zipfile import ZipFile
import tqdm
import seaborn as sn
import re
import palettable
import xarray as xr
import ipywidgets as widgets
from IPython.display import display

## get SSP emissions
from tools.RCMIP import *

GCP_historical_emissions=pd.read_csv('https://raw.githubusercontent.com/openclimatedata/global-carbon-budget/master/data/historical-budget.csv',index_col=0,usecols=[0,1,2])
# set 1750 value to 0
GCP_historical_emissions.loc[1750,'Fossil-Fuel-Industry'] = 0
# linearly interpolate between 0 and first land-use-change datapoint
GCP_historical_emissions.loc[1750:1849,'Land-Use-Change'] = np.linspace(0,0.99,100) * GCP_historical_emissions.loc[1850,'Land-Use-Change']

header = {
  "User-Agent": "Mozilla/5.0 (X11; Linux x86_64) AppleWebKit/537.36 (KHTML, like Gecko) Chrome/50.0.2661.75 Safari/537.36",
  "X-Requested-With": "XMLHttpRequest"
}
    
## Pandas options setup

pd.set_option('display.max_columns', 100)
    
## Matplotlib rcparams setup:

matplotlib.rcParams['font.family']='Helvetica'
matplotlib.rcParams['font.size']=11
# matplotlib.rcParams['font.weight']=400

matplotlib.rcParams['image.cmap']='cmo.ice'

matplotlib.rcParams['axes.prop_cycle']=matplotlib.cycler('color',['011936','FF7D00','225560','BFACAA','D72638','788C9B','A33F00','7CAAB0','685655','EB767C'])
matplotlib.rcParams['axes.formatter.limits']=-3,3
# matplotlib.rcParams['axes.labelweight']=300

matplotlib.rcParams['legend.frameon']=False

matplotlib.rcParams['boxplot.whiskers']=(5,95)
matplotlib.rcParams['boxplot.showfliers']=False
matplotlib.rcParams['boxplot.showfliers']=False
matplotlib.rcParams['boxplot.medianprops.color']='black'

matplotlib.rcParams['errorbar.capsize']=5

matplotlib.rcParams['hist.bins']='auto'

plt.rcParams['pdf.fonttype'] = 42

%matplotlib inline

%load_ext line_profiler
%load_ext memory_profiler
# erf_ar6 = pd.read_csv('https://raw.githubusercontent.com/Priestley-Centre/ssp_erf/master/SSPs/ERF_ssp245_1750-2500.csv',index_col=0,dtype=float)
# erf_ar6['solar'].loc[2000:2020].plot()<jupyter_output><empty_output><jupyter_text>NROY Ensemble Generation  Part I: the Global Warming Index  required OLSE script:<jupyter_code># definition to compute the regressions:
def OLSE_NORM(X,Y,add_intercept=True):
    
    ## computes a multiple OLS regression over a field against several indices. First dimension is time, second is features (X), or targets (Y)
    
    if add_intercept:
    
        X_1 = np.concatenate((np.ones(X.shape[0])[:,np.newaxis],X),axis=1)
        
    else:
        
        X_1 = X.copy()
    
    B = np.dot( np.linalg.inv( np.dot( X_1.T , X_1 ) ) , np.dot( X_1.T , Y ) )
    
    e = Y - np.dot(X_1,B)
    
    SSE = np.sum(e**2,axis=0)

    SST = np.sum((Y - np.mean(Y,axis=0))**2,axis=0)

    R2 = 1 - SSE / SST

    MSE_var = SSE / (X_1.shape[0] - X_1.shape[-1])

    SE_B = np.sqrt( np.diag( np.linalg.inv( np.dot( X_1.T , X_1 ) ) )[:,np.newaxis] * MSE_var[np.newaxis,:] )
    
    return {'coefs':B[1:],'coef_err':SE_B[1:],'res':e,'intercept':B[0],'intercept_err':SE_B[0],'R2':R2}<jupyter_output><empty_output><jupyter_text>Retrieve observational GMST data<jupyter_code># ## Import data
# baseline = ['1861','1880']

# url = "https://www.metoffice.gov.uk/hadobs/hadcrut4/data/current/time_series/HadCRUT.4.6.0.0.monthly_ns_avg.txt"
# req = requests.get(url, headers=header)
# data = StringIO(req.text)
# HadCRUT4 = pd.read_csv(data,delim_whitespace=True,index_col=0,names=['date','anom','lower_bias','upper_bias','lower_measure','upper_measure','lower_cover','upper_cover','lower_bias_measure','upper_bias_measure','lower_full','upper_full']).apply(pd.to_numeric)
# HadCRUT4.index = pd.to_datetime(['-'.join(x.split('/'))+'-01' for x in HadCRUT4.index])
# HadCRUT4 = HadCRUT4['anom'].rename('HadCRUT4')

# GISS = pd.read_csv("https://data.giss.nasa.gov/gistemp/tabledata_v4/GLB.Ts+dSST.txt",skiprows=7,skipfooter=7,delim_whitespace=True,index_col=0).iloc[:,:12]
# GISS = GISS.loc[GISS.index!='Year'].replace('****',np.nan).apply(pd.to_numeric)/100
# month_map = dict(zip(GISS.columns,["{:02d}".format(x) for x in np.arange(1,13)]))
# GISS = pd.DataFrame(GISS.values.flatten(),index = pd.to_datetime([x+'-'+month_map[y]+'-01' for x in GISS.index for y in GISS.columns]),columns=['anom'])
# GISS = GISS['anom'].rename('GISTEMP')

# CW_Had4_krig = pd.read_csv("https://www-users.york.ac.uk/~kdc3/papers/coverage2013/had4_krig_v2_0_0.txt",delim_whitespace=True,index_col=0,names=['date','anom','unc_full','unc_cover','unc_ens']).apply(pd.to_numeric)
# CW_Had4_krig.index = pd.to_datetime([str(int(x//1))+'-'+"{:02d}".format(int((x % 1)*12)+1)+'-01' for x in CW_Had4_krig.index])
# CW_Had4_krig = CW_Had4_krig['anom'].rename('CW')

# Berkeley = pd.read_csv("http://berkeleyearth.lbl.gov/auto/Global/Land_and_Ocean_complete.txt",skiprows=76,delim_whitespace=True,usecols=[0,1,2,3],names=['date','month','anom','unc_full'])
# Berkeley = Berkeley.iloc[:Berkeley.loc[Berkeley.date=='%'].index.min()].apply(pd.to_numeric)
# Berkeley.index = pd.to_datetime(Berkeley.date.astype(str) + Berkeley.month.astype(str),format='%Y%m')
# Berkeley = Berkeley['anom'].rename('Berkeley')

# parser = 'html.parser'  # or 'lxml' (preferred) or 'html5lib', if installed
# resp = urllib.request.urlopen("https://www.ncei.noaa.gov/data/noaa-global-surface-temperature/v5/access/timeseries/")
# soup = BeautifulSoup(resp, parser, from_encoding=resp.info().get_param('charset'))
# NOAA_fname = [x['href'] for x in soup.find_all('a', href=True) if x['href'].split('.')[1:5]==['mon','land_ocean','90S','90N']][0]
# NOAA = pd.read_csv("https://www.ncei.noaa.gov/data/noaa-global-surface-temperature/v5/access/timeseries/"+NOAA_fname,delim_whitespace=True,names=['date','month','anom','unc_full'],usecols=[0,1,2,3]).apply(pd.to_numeric)
# NOAA.index = pd.to_datetime(NOAA.date.astype(str) + NOAA.month.astype(str),format='%Y%m')
# NOAA = NOAA['anom'].rename('NOAA')

# ## Combine into a single df

# GMST_datasets = pd.concat([HadCRUT4,CW_Had4_krig,GISS,NOAA,Berkeley],axis=1)

# ## Baseline to 1861:1880 OR to HadCRUT4 20C mean if data begins in 1880 (GISTEMP / NOAA)

# GMST_datasets.loc[:,['HadCRUT4','CW','Berkeley']] -= GMST_datasets.loc[baseline[0]:baseline[1],['HadCRUT4','CW','Berkeley']].mean()
# GMST_datasets.loc[:,['GISTEMP','NOAA']] -= GMST_datasets.loc['1900':'2000',['GISTEMP','NOAA']].mean() -  GMST_datasets.loc['1900':'2000','HadCRUT4'].mean()
# ERA5 = xr.open_dataset('/home/leachl/Documents/Datasets/GMST/ERA5.nc').sel(expver=1).t2m#.squeeze().to_pandas()
# weights = np.cos(np.deg2rad(ERA5.latitude))
# ERA5_fldmean = ERA5.weighted(weights).mean(('longitude','latitude'))
# ERA5_glob = (ERA5_fldmean.groupby("time.month") - ERA5_fldmean.groupby("time.month").mean('time')).to_pandas()

# GSAT_to_GMST = OLSE_NORM(GMST_datasets.loc[ERA5_glob.index[0]:ERA5_glob.dropna().index[-1]].mean(axis=1).values[:,np.newaxis],ERA5_glob.dropna().values[:,np.newaxis])
# print('GMST to GSAT scaling = ',round(GSAT_to_GMST['coefs'].flatten()[0],3),'+-',round(GSAT_to_GMST['coef_err'].flatten()[0],4))

# GMST_datasets.rolling(12).mean().plot()
# ## Import the observational uncertainty from HadCRUT

# HadCRUT_ensmem_zip = ZipFile(BytesIO(requests.get('https://www.metoffice.gov.uk/hadobs/hadcrut4/data/current/time_series/HadCRUT.4.6.0.0.annual_ns_avg_realisations.zip').content))
# HadCRUT_ensmem_fnames = [x.filename for x in HadCRUT_ensmem_zip.infolist() if x.filename[-4:] == '.txt']
# HadCRUT_ensmem_index = pd.read_csv(HadCRUT_ensmem_zip.open(HadCRUT_ensmem_fnames[0]),header=None,delim_whitespace=True,usecols=[0],index_col=0).index
# # HadCRUT_ensmem_index = pd.to_datetime(['-'.join(x.split('/'))+'-01' for x in HadCRUT_ensmem_index])
# HadCRUT_ensmem = pd.concat([pd.read_csv(HadCRUT_ensmem_zip.open(x),header=None,delim_whitespace=True,usecols=[1]).iloc[:,0].rename(str(i)) for i,x in enumerate(HadCRUT_ensmem_fnames)],axis=1)
# HadCRUT_ensmem.index = HadCRUT_ensmem_index
## import forcing ensemble: for now, use the most up-to-date copy I have (.idlsave)

# forcing_dataset = readsav('/home/leachl/Documents/Datasets/RF/1750_Oct2017_forcings.idlsave')
# rf_ant = pd.DataFrame(forcing_dataset['rf_anthro'].byteswap().newbyteorder(),index = pd.to_datetime([str(int(x))+'-01-01' for x in forcing_dataset['year']])).apply(pd.to_numeric)
# rf_nat = pd.DataFrame(forcing_dataset['rf_nat'].byteswap().newbyteorder(),index = pd.to_datetime([str(int(x))+'-01-01' for x in forcing_dataset['year']])).apply(pd.to_numeric)<jupyter_output><empty_output><jupyter_text>Generate forcing profiles<jupyter_code>## creating a forcing response dataset

def generate_forcing(N):

    erf_ar6 = pd.read_csv('https://raw.githubusercontent.com/Priestley-Centre/ssp_erf/master/SSPs/ERF_ssp245_1750-2500.csv',index_col=0,dtype=float)
    erf_ar6 -= erf_ar6.loc[1750]
    
    ## generating many rf shapes for the GWI

    rf_factors = {}

    rf_factors['co2'] = sp.stats.norm(1,0.20/1.645).rvs(N)#ppf((erf_ar6.mean()>0).co2*random_percentiles + ~(erf_ar6.mean()>0).co2*(1-random_percentiles))
    rf_factors['ch4'] = sp.stats.norm(1,0.28/1.645).rvs(N)#ppf((erf_ar6.mean()>0).ch4*random_percentiles + ~(erf_ar6.mean()>0).ch4*(1-random_percentiles))
    rf_factors['n2o'] = sp.stats.norm(1,0.2/1.645).rvs(N)#ppf((erf_ar6.mean()>0).n2o*random_percentiles + ~(erf_ar6.mean()>0).n2o*(1-random_percentiles))
    rf_factors['other_wmghg'] = sp.stats.norm(1,0.20/1.645).rvs(N)#ppf((erf_ar6.mean()>0).other_wmghg*random_percentiles + ~(erf_ar6.mean()>0).other_wmghg*(1-random_percentiles))
    rf_factors['o3_tropospheric'] = sp.stats.norm(1,0.5/1.645).rvs(N)#ppf((erf_ar6.mean()>0).o3_tropospheric*random_percentiles + ~(erf_ar6.mean()>0).o3_tropospheric*(1-random_percentiles))
    rf_factors['o3_stratospheric'] = sp.stats.norm(1,2/1.645).rvs(N)#ppf((erf_ar6.mean()>0).o3_stratospheric*random_percentiles + ~(erf_ar6.mean()>0).o3_stratospheric*(1-random_percentiles))
    rf_factors['h2o_stratospheric'] = sp.stats.norm(1,0.72/1.645).rvs(N)#ppf((erf_ar6.mean()>0).h2o_stratospheric*random_percentiles + ~(erf_ar6.mean()>0).h2o_stratospheric*(1-random_percentiles))
    rf_factors['contrails'] = sp.stats.norm(1,0.75/1.645).rvs(N)#ppf((erf_ar6.mean()>0).contrails*random_percentiles + ~(erf_ar6.mean()>0).contrails*(1-random_percentiles))
    rf_factors['land_use'] = sp.stats.norm(1,0.75/1.645).rvs(N)#ppf((erf_ar6.mean()>0).land_use*random_percentiles + ~(erf_ar6.mean()>0).land_use*(1-random_percentiles))
    rf_factors['volcanic'] = sp.stats.norm(1,0.5/1.645).rvs(N)
    rf_factors['solar'] = sp.stats.norm(1,1/1.645).rvs(N)

    pct_5 = 0.04 / 0.08
    pct_95 = 0.18 / 0.08
    sigma = (np.log(pct_95) - np.log(pct_5)) / (sp.stats.norm().ppf(0.95)-sp.stats.norm().ppf(0.05))
    mu = np.log(pct_5) - sigma * sp.stats.norm().ppf(0.05)
    rf_factors['bc_on_snow'] = sp.stats.lognorm(s=sigma,scale=np.exp(mu)).rvs(N)

    ant_rf_ensemble = np.zeros((270,N))
    nat_rf_ensemble = np.zeros((270,N))

    for agent in ['co2','ch4','n2o','other_wmghg','o3_tropospheric','o3_stratospheric','h2o_stratospheric','contrails','land_use','bc_on_snow']:
        ant_rf_ensemble += rf_factors[agent]*erf_ar6.loc[1750:2019,agent].values[:,None]

    for agent in ['volcanic','solar']:
        nat_rf_ensemble += rf_factors[agent]*erf_ar6.loc[1750:2019,agent].values[:,None]

    ## aerosol distribution 

    ### generate the mean shape:

    RCMIP_so2_emms = RCMIP_emms.loc[('World','ssp245','Emissions|Sulfur')].iloc[4:].apply(pd.to_numeric).interpolate().loc[:'2019'].values[None]
    RCMIP_oc_emms = RCMIP_emms.loc[('World','ssp245','Emissions|OC')].iloc[4:].apply(pd.to_numeric).interpolate().loc[:'2019'].values[None]
    RCMIP_bc_emms = RCMIP_emms.loc[('World','ssp245','Emissions|BC')].iloc[4:].apply(pd.to_numeric).interpolate().loc[:'2019'].values[None]

    RCMIP_so2_emms -= RCMIP_so2_emms[0,0]
    RCMIP_oc_emms -= RCMIP_oc_emms[0,0]
    RCMIP_bc_emms -= RCMIP_bc_emms[0,0]

    ERFari_coefs = pd.read_csv('../Parameter_Sets/CMIP6_ERFari.csv',index_col=0)
    ERFaci_coefs = pd.read_csv('../Parameter_Sets/CMIP6_ERFaci.csv',index_col=0)

    beta = -1*np.exp(ERFaci_coefs['log_beta'].mean())
    s = np.exp(ERFaci_coefs['log_s'].mean())
    alpha = ERFaci_coefs['alpha'].mean()

    # ERFari_shape = ERFari_coefs.so2.mean() * RCMIP_so2_emms + ERFari_coefs.bc.mean() * RCMIP_bc_emms + ERFari_coefs.oc.mean() * RCMIP_oc_emms
    # ERFaci_shape = beta * np.log(1 + RCMIP_so2_emms / s) + alpha * ( RCMIP_oc_emms + RCMIP_bc_emms )

    ERFari_shape = erf_ar6.loc[1750:2019,['aerosol-radiation_interactions']].values.T
    ERFaci_shape = erf_ar6.loc[1750:2019,['aerosol-cloud_interactions']].values.T

    ## construct the Smith distributions 
    def fit_skewnorm(x,X,percentiles):
        distr = sp.stats.skewnorm(*x)
        return sum(abs(distr.ppf(percentiles) - X))

    ERFaci_smith_sknorm_params = sp.optimize.minimize(fit_skewnorm,x0=[1,1,1],args=([0.13,0.59,1.17],[0.05,0.5,0.95]),method='nelder-mead').x
    ERFari_smith_sknorm_params = sp.optimize.minimize(fit_skewnorm,x0=[1,1,1],args=([0.07,0.27,0.6],[0.05,0.5,0.95]),method='nelder-mead').x

    ERFaci_smith_sknorm_distr = -1*sp.stats.skewnorm(*ERFaci_smith_sknorm_params).rvs(N)[:,None]
    ERFari_smith_sknorm_distr = -1*sp.stats.skewnorm(*ERFari_smith_sknorm_params).rvs(N)[:,None]

    ERFari_samples = ERFari_smith_sknorm_distr / ERFari_shape[0,-1] * ERFari_shape
    ERFaci_samples = ERFaci_smith_sknorm_distr / ERFaci_shape[0,-1] * ERFaci_shape
    
    tot_aer = (ERFari_samples+ERFaci_samples).T

    ## Combine all anthro forcings
    ant_rf_ensemble += tot_aer
    
    ## 
    rf_nat = pd.DataFrame(nat_rf_ensemble,index=np.arange(1750,2020),columns=pd.MultiIndex.from_product([['forcing_'+str(x) for x in np.arange(N)],['forcing']])).loc[:2019]

    rf_ant = pd.DataFrame(ant_rf_ensemble,index=np.arange(1750,2020),columns=pd.MultiIndex.from_product([['forcing_'+str(x) for x in np.arange(N)],['forcing']])).loc[:2019]
    
    return rf_ant,rf_nat
N_forc = 10000

rf_ant,rf_nat = generate_forcing(N_forc)

rf_ant.plot(legend=None,lw=0.02,c='k')

rf_nat.plot(legend=None,lw=0.02,c='k')<jupyter_output><empty_output><jupyter_text>Combine HadCRUT uncertainty with 5-dataset mean signal<jupyter_code># ## generate the uncertain temperature array:

# temp_obs_med = GMST_datasets.mean(axis=1).loc[:str(int(rf_ant.index[-1]))].resample('Y').mean().values

# temp_obs_unc = HadCRUT_ensmem.loc[:rf_ant.index[-1]].values - HadCRUT_ensmem.loc[:rf_ant.index[-1]].median(axis=1)[:,np.newaxis]

# temp_obs_full = temp_obs_med[:,np.newaxis] + temp_obs_unc<jupyter_output><empty_output><jupyter_text>Alternative aerosol generation with varying shapes (work-in-progress):<jupyter_code># ## New aerosol generation methodology, following Smith et al (2020) ; multivariate gaussian distribution, sample parameters, generate shapes from emission data.
# ERFari_params = pd.DataFrame(data=np.array([[ 3.25934432e-02, -3.46653116e-04, -2.48761587e-03],
#        [ 2.48231885e-02, -1.26321908e-02, -9.42090973e-04],
#        [ 2.68560170e-02, -2.08950151e-03, -2.60579679e-03],
#        [ 1.01700746e-01, -3.04192033e-02, -2.63868473e-03],
#        [ 1.46132214e-01, -4.41485999e-02, -6.67827786e-03],
#        [ 1.95732170e-03,  4.15481407e-03, -2.90822562e-03],
#        [-5.61261154e-02,  8.85201623e-03, -7.47602910e-04],
#        [ 3.87211606e-02, -1.42408227e-02, -1.78155606e-03],
#        [ 3.02173214e-03, -3.40228173e-03, -1.25702781e-03],
#        [ 2.54832726e-03,  6.32322444e-05, -2.38863636e-03]]),index=['CanESM5', 'E3SM', 'GFDL-CM4', 'GFDL-ESM4', 'GISS-E2-1-G',
#        'HadGEM3-GC31-LL', 'IPSL-CM6A-LR', 'MIROC6', 'NorESM2-LM',
#        'UKESM1-0-LL'],columns=['bc', 'oc', 'so2'])

# ERFari_sfs = sp.stats.multivariate_normal(ERFari_params.mean(),ERFari_params.cov()).rvs(N)
# ERFari_samples = (ERFari_sfs[...,None]*RCMIP_emms.loc[('World','ssp245',['Emissions|BC','Emissions|OC','Emissions|Sulfur'])].droplevel((0,1)).iloc[:,4:].apply(pd.to_numeric).interpolate(axis=1).loc[:,:'2019'].values[None,...]).sum(axis=1)

# ### LINEAR ERFaci params
# # ERFaci_params = pd.DataFrame(data=np.array([[-0.23977316,  0.05033984, -0.0026841 ],
# #        [-0.03596736, -0.0061301 , -0.00999338],
# #        [-0.02306089, -0.00397075, -0.00440759],
# #        [-0.04100239, -0.00815999, -0.00296047],
# #        [-0.11952911,  0.01185454, -0.00136009],
# #        [-0.0751146 ,  0.0115965 , -0.00595064],
# #        [-0.00958696,  0.00027465, -0.00252801],
# #        [-0.10616509,  0.01464634, -0.00292868],
# #        [ 0.02490414, -0.02261304, -0.00721553],
# #        [ 0.06822656, -0.01548075, -0.01065098]])
# # ,index=['CanESM5', 'E3SM', 'GFDL-CM4', 'GFDL-ESM4', 'GISS-E2-1-G',
# #        'HadGEM3-GC31-LL', 'IPSL-CM6A-LR', 'MIROC6', 'NorESM2-LM',
# #        'UKESM1-0-LL'],columns=['bc', 'oc', 'so2'])

# ### LOG so2 ERFaci params
# ERFaci_params = pd.DataFrame(data=np.array([[-8.04936152e-01, -4.33765114e+00, -1.76516101e-02],
#        [ 4.20348167e-01, -4.19165818e+00, -1.57827154e-02],
#        [ 1.67034191e+00, -5.34562320e+00, -8.80674737e-03],
#        [ 2.25088976e+00, -5.57251763e+00, -1.19361347e-02],
#        [-1.80621064e+00, -4.82004944e+00, -1.88411062e-02],
#        [-1.48244252e-01, -4.48450777e+00, -6.70740113e-03],
#        [ 1.03705925e+00, -5.88364385e+00, -1.91422284e-03],
#        [-1.00837014e+00, -4.73118298e+00, -1.30437040e-02],
#        [ 2.07731610e+00, -5.03084334e+00, -1.39245898e-02],
#        [-2.52178641e-01, -3.93423009e+00,  3.35367830e-03]])

# ,index=['CanESM5', 'E3SM', 'GFDL-CM4', 'GFDL-ESM4', 'GISS-E2-1-G',
#        'HadGEM3-GC31-LL', 'IPSL-CM6A-LR', 'MIROC6', 'NorESM2-LM',
#        'UKESM1-0-LL'],columns=['log_beta', 'log_beta_s', 'alpha'])

# # ERFaci_params.loc[:,'so2_logcoef'] = np.log(-1*ERFaci_params.loc[:,'so2_logcoef'])

# ### LINEAR ERFaci params
# # ERFaci_sfs = sp.stats.multivariate_normal(ERFaci_params.mean(),ERFaci_params.cov()).rvs(N)
# # ERFaci_samples = (ERFaci_sfs[...,None]*RCMIP_emms.loc[('World','ssp245',['Emissions|BC','Emissions|OC','Emissions|Sulfur'])].droplevel((0,1)).iloc[:,4:].apply(pd.to_numeric).interpolate(axis=1).loc[:,:'2019'].values[None,...]).sum(axis=1)

# ### LOG so2 ERFaci params
# beta_samples = -1 * np.exp( sp.stats.norm(ERFaci_params['log_beta'].mean(),ERFaci_params['log_beta'].std()).rvs(N) )[...,None]
# beta_s_samples = np.exp( sp.stats.norm(ERFaci_params['log_beta_s'].mean(),ERFaci_params['log_beta_s'].std()).rvs(N) )[...,None]
# s_samples = -1 * beta_samples / beta_s_samples
# alpha_samples = sp.stats.norm(ERFaci_params['alpha'].mean(),ERFaci_params['alpha'].std()).rvs(N)[...,None]

# RCMIP_so2_emms = RCMIP_emms.loc[('World','ssp245','Emissions|Sulfur')].iloc[4:].apply(pd.to_numeric).interpolate().loc[:'2019'].values[None]
# RCMIP_oc_emms = RCMIP_emms.loc[('World','ssp245','Emissions|OC')].iloc[4:].apply(pd.to_numeric).interpolate().loc[:'2019'].values[None]
# RCMIP_bc_emms = RCMIP_emms.loc[('World','ssp245','Emissions|BC')].iloc[4:].apply(pd.to_numeric).interpolate().loc[:'2019'].values[None]

# ERFaci_samples = beta_samples * np.log(1+RCMIP_so2_emms/s_samples) + alpha_samples * ( RCMIP_oc_emms+RCMIP_bc_emms )<jupyter_output><empty_output><jupyter_text>Generate response parameters<jupyter_code>def generate_response_params():

    ## create the response parameter ranges:

    ### based on inferred ranges from CMIP6 - 18 combinations
    response_names = ['response_'+str(x) for x in np.arange(24)]
    response_params = pd.DataFrame(index=['d','q'],columns = pd.MultiIndex.from_product([response_names,[1,2,3]])).apply(pd.to_numeric)
    response_params.loc[:] = 0

    d1_range = [0.2,0.8,1.4,2]
    d2_range = [4,8,12,16]
    d3_range = [100,200,400,800]
    q1_range = [0.04,0.16,0.28,0.4]

    RWF_range = [0.3 , 0.4 , 0.5, 0.6, 0.7, 0.8]
    ECS=3

    i=0
    for d_num in np.arange(4):
        for RWF in RWF_range:
            response_params.loc['d',(response_names[i],1)] = d1_range[d_num]
            response_params.loc['d',(response_names[i],2)] = d2_range[d_num]
            response_params.loc['d',(response_names[i],3)] = d3_range[d_num]
            response_params.loc['q',(response_names[i],1)] = q1_range[d_num]

            q2,q3 = sp.optimize.root(lambda x: q_to_tcr_ecs(response_params.reindex([response_names[i]],axis=1,level=0) + np.array([[0,0],[0,x[0]],[0,x[1]]]).T,F_2x=3.76).values.flatten() - np.array([ECS,ECS*RWF]),x0=np.array([0.32,0.47]),method='lm').x

            response_params.loc['q',(response_names[i],2)] = q2
            response_params.loc['q',(response_names[i],3)] = q3

            i+=1

    ## remove any that are unphysical (negative parameters)

    response_params = response_params.reindex([i for i,x in ((response_params<0).sum().unstack().sum(axis=1)==1).iteritems() if not x],axis=1,level=0)
    
    return response_params
response_params = generate_response_params()<jupyter_output><empty_output><jupyter_text>Generate temperature profiles:<jupyter_code>## generate the temperature responses:

null_gas_params = pd.read_csv('../Parameter_Sets/Complete_gas_cycle_params.csv',header=[0,1],index_col=0).reindex(['carbon_dioxide'],axis=1,level=1)
##
temp_ant = run_GIR(emissions_in=return_empty_emissions(rf_ant,gases_in=['carbon_dioxide']),gas_parameters=null_gas_params,forcing_in=rf_ant,thermal_parameters=response_params)['T'].droplevel(1,axis=1)
temp_nat = run_GIR(emissions_in=return_empty_emissions(rf_nat,gases_in=['carbon_dioxide']),gas_parameters=null_gas_params,forcing_in=rf_nat,thermal_parameters=response_params)['T'].droplevel(1,axis=1)<jupyter_output>Integrating 10000 scenarios, 1 gas cycle parameter sets, 18 thermal response parameter sets, over ['carbon_dioxide'] forcing agents, between 1750 and 2019...<jupyter_text>save temperature series to netcdfs for combination with OLSE coefficients later<jupyter_code>xr.DataArray(temp_ant.values,
             dims=['time','index'],
             coords=dict(time=pd.to_datetime(temp_ant.index.astype('str')),index=temp_ant.columns)).unstack('index').rename({'Scenario':'forcing_mem','Thermal set':'response_mem'}).to_netcdf('../Parameter_Sets/NROY_GWI_informed/GWI_components/ant_temperature.nc')

xr.DataArray(temp_nat.values,
             dims=['time','index'],
             coords=dict(time=pd.to_datetime(temp_ant.index.astype('str')),index=temp_ant.columns)).unstack('index').rename({'Scenario':'forcing_mem','Thermal set':'response_mem'}).to_netcdf('../Parameter_Sets/NROY_GWI_informed/GWI_components/nat_temperature.nc')<jupyter_output><empty_output><jupyter_text>Retrieve CMIP6 internal variability & subsample<jupyter_code>## computing the uncertainty from internal variability

def generate_IV(series_length=170):
    ### get internal variability timeseries:

    piControl_data = pd.read_csv('../tools/cmip6_data/piControl.csv',index_col=0)
    for i,x in enumerate(['source','member','variable']):
        piControl_data.loc[x] = [x.split('_')[i+1] for x in piControl_data.columns]

    CMIP6_int_var=piControl_data.T.set_index(['source','member','variable']).T.xs('tas',axis=1,level=-1).apply(pd.to_numeric)

    ### subsample 100 * 170 year slices
    ### check if drift, discard if drift
    ### randomly pick one non-drifting sample

    CMIP6_int_var_samples = pd.DataFrame(index=np.arange(series_length))
    for model,data in CMIP6_int_var.iteritems():
        arr = data.dropna().values
        chunksize = arr.size
        if chunksize<series_length:
            continue
        starting_points = np.random.choice(np.arange(chunksize-series_length),100)
        nodrift_points = []
        for s in starting_points:
            subarr = arr[s:s+series_length]
            drift = sp.stats.linregress(np.arange(series_length),subarr).slope
            if abs(drift)*100>0.15:
                continue
            else:
                nodrift_points+=[s]
        chosen_start = np.random.choice(nodrift_points)
        index_name = '_'.join(list(model))+'_'+str(chosen_start)
        CMIP6_int_var_samples.loc[:,index_name] = arr[chosen_start:chosen_start+series_length] - arr[chosen_start:chosen_start+series_length].mean()

    ### remove model degeneracies (take first member of each model):
    models = []
    chosen_members = []
    for model,data in CMIP6_int_var_samples.dropna(axis=1).iteritems():
        if model.split('_')[0] in models:
            continue
        else:
            chosen_members+=[model]
            models+=[model.split('_')[0]]

    return CMIP6_int_var_samples[chosen_members]
CMIP6_int_var_samples_nondeg = pd.concat([generate_IV(),generate_IV()],axis=1)
CMIP6_int_var_samples_nondeg.index = np.arange(1850,2020)
CMIP6_int_var_samples_nondeg.plot(legend=None,lw=0.1,c='k')<jupyter_output><empty_output><jupyter_text>Choose and import the observational temperature dataset (updated to start of 2021)Allowed values:[GISTEMP](https://data.giss.nasa.gov/gistemp/) - single series[NOAA](https://www.ncei.noaa.gov/data/noaa-global-surface-temperature/v5/access/timeseries/) - single series[BERKELEY](http://berkeleyearth.org/data/) - singles series[CW](https://www-users.york.ac.uk/~kdc3/papers/coverage2013/series.html) - 100 member ensemble[HadCRUT5](http://data.ceda.ac.uk/badc/ukmo-hadobs/data/insitu/MOHC/HadOBS/HadCRUT/HadCRUT5/analysis/diagnostics) - 200 member ensemble[HadCRUT4](https://www.metoffice.gov.uk/hadobs/hadcrut4/data/current/download.html) - 100 member ensemble<jupyter_code>## button selector...
# select_dataset = widgets.ToggleButtons(options=['HadCRUT5','HadCRUT4','NOAA','GISTEMP','CW','BERKELEY'],value='HadCRUT5',description='Choose observational timeseries:')
# select_dataset
def get_GMST(select_dataset):

    ## retrieve chosen dataset:
    if select_dataset == 'HadCRUT5':
        GMST = xr.open_dataset('/home/leachl/Documents/Datasets/GMST/HadCRUT.5.0.0.0_analysis_ensemble-series_global_monthly.nc').tas.to_pandas().T
        GMST.columns.name = 'HadCRUT5_obs_mem'

    if select_dataset == 'HadCRUT4':
        GMST = pd.concat([pd.read_csv(x,header=None,delim_whitespace=True,index_col=0,usecols=[0,1]).iloc[:,0] for x in sorted(glob.glob('/home/leachl/Documents/Datasets/GMST/HadCRUT.4.6.0.0.monthly_ns_avg_realisations/*.txt'))],axis=1,keys=np.arange(100))
        GMST.index = pd.to_datetime(GMST.index)
        GMST.columns.name = 'HadCRUT4_obs_mem'

    if select_dataset == 'NOAA':
        GMST = pd.read_csv('/home/leachl/Documents/Datasets/GMST/aravg.mon.land_ocean.90S.90N.v5.0.0.202012.asc',delim_whitespace=True,names=['date','month','anom','unc_full'],usecols=[0,1,2,3]).apply(pd.to_numeric)
        GMST.index = pd.to_datetime(GMST.date.astype(str)+'-'+GMST.month.astype(str))
        GMST = pd.DataFrame(GMST['anom'].rename(0))
        GMST.columns.name = 'NOAA_obs_mem'

    if select_dataset == 'GISTEMP':
        GMST = pd.read_csv('/home/leachl/Documents/Datasets/GMST/GLB.Ts+dSST.csv',skiprows=2,names=['year']+list(range(1,13)),usecols=range(13),index_col=0).unstack().reset_index()
        GMST.index = pd.to_datetime(GMST.year.astype(str)+'-'+GMST.level_0.astype(str))
        GMST = pd.DataFrame(GMST.sort_index().iloc[:,-1].rename(0))
        GMST.columns.name = 'GISTEMP_obs_mem'

    if select_dataset == 'CW':
        GMST=pd.read_csv('/home/leachl/Documents/Datasets/GMST/had4_krig_ensemble_v2_0_0.txt',delim_whitespace=True,index_col=0,names=['date']+list(range(100)))
        GMST.index = pd.to_datetime(GMST.index.astype(int).astype(str)+'-'+np.ceil((GMST.index%1)*12).astype(int).astype(str))
        # current version of the dataset has a null member at 32
        GMST = GMST.drop(32,axis=1)
        GMST.columns.name = 'CW_obs_mem'

    if select_dataset == 'BERKELEY':
        # NB this uses air temps above sea ice, not water temps
        GMST = pd.read_csv('/home/leachl/Documents/Datasets/GMST/Land_and_Ocean_complete.txt',skiprows=76,skipfooter=2057,delim_whitespace=True,usecols=[0,1,2,3],names=['date','month','anom','unc_full'])
        GMST.index = pd.to_datetime(GMST.date.astype(str)+'-'+GMST.month.astype(str))
        GMST = pd.DataFrame(GMST['anom'].rename(0))
        GMST.columns.name = 'BERKELEY_obs_mem'
        
    return GMST
def run_GWI(select_dataset,temp_ant,temp_nat):

    GMST = get_GMST(select_dataset)

    ## resample GMST to annual:

    GMST = GMST.resample('y').mean()
    GMST.index = GMST.index.year.astype(int)
    GMST = GMST.loc[:2019]
    
    # regress temps vs observed temps
    
    ant_coefs = np.empty((GMST.shape[1],temp_ant.shape[1]))
    nat_coefs = np.empty((GMST.shape[1],temp_nat.shape[1]))
    R2_vals = np.empty((GMST.shape[1],temp_ant.shape[1]))

    temp_vstack = np.array([temp_ant.loc[GMST.index].values,temp_nat.loc[GMST.index].values]).T

    for i in tqdm.tqdm(np.arange(ant_coefs.shape[1])):
        _lreg = OLSE_NORM(temp_vstack[i],GMST.values)
        ant_coefs[:,i] = _lreg['coefs'][0]
        nat_coefs[:,i] = _lreg['coefs'][1]
        R2_vals[:,i] = _lreg['R2']
        
    ## save results
    
    xr.DataArray(ant_coefs,
                 dims=[GMST.columns.name,'index'],
                 coords={GMST.columns.name:GMST.columns,'index':temp_ant.columns}).unstack('index').rename({'Scenario':'forcing_mem','Thermal set':'response_mem'}).to_netcdf('../Parameter_Sets/NROY_GWI_informed/GWI_components/ant_coefs_forc_'+select_dataset+'.nc')

    xr.DataArray(nat_coefs,
                 dims=[GMST.columns.name,'index'],
                 coords={GMST.columns.name:GMST.columns,'index':temp_nat.columns}).unstack('index').rename({'Scenario':'forcing_mem','Thermal set':'response_mem'}).to_netcdf('../Parameter_Sets/NROY_GWI_informed/GWI_components/nat_coefs_forc_'+select_dataset+'.nc')
    
    xr.DataArray(R2_vals,
                 dims=[GMST.columns.name,'index'],
                 coords={GMST.columns.name:GMST.columns,'index':temp_nat.columns}).unstack('index').rename({'Scenario':'forcing_mem','Thermal set':'response_mem'}).to_netcdf('../Parameter_Sets/NROY_GWI_informed/GWI_components/R2_forc_'+select_dataset+'.nc')
    
    
    ## regress temps vs IV
    
    ant_coefs_iv = np.empty((CMIP6_int_var_samples_nondeg.shape[1],temp_ant.shape[1]))
    nat_coefs_iv = np.empty((CMIP6_int_var_samples_nondeg.shape[1],temp_nat.shape[1]))
    R2_vals_iv = np.empty((CMIP6_int_var_samples_nondeg.shape[1],temp_nat.shape[1]))
    
    IV_regress = CMIP6_int_var_samples_nondeg.loc[GMST.index].values

    for i in tqdm.tqdm(np.arange(ant_coefs_iv.shape[1])):
        _lreg = OLSE_NORM(temp_vstack[i],IV_regress)
        ant_coefs_iv[:,i] = _lreg['coefs'][0]
        nat_coefs_iv[:,i] = _lreg['coefs'][1]
        R2_vals_iv[:,i] = _lreg['R2']
        
    
    xr.DataArray(ant_coefs_iv,
                 dims=['int_var_mem','index'],
                 coords={'int_var_mem':CMIP6_int_var_samples_nondeg.columns,'index':temp_ant.columns}).unstack('index').rename({'Scenario':'forcing_mem','Thermal set':'response_mem'}).to_netcdf('../Parameter_Sets/NROY_GWI_informed/GWI_components/ant_coefs_IV_'+select_dataset+'.nc')

    xr.DataArray(nat_coefs_iv,
                 dims=['int_var_mem','index'],
                 coords={'int_var_mem':CMIP6_int_var_samples_nondeg.columns,'index':temp_nat.columns}).unstack('index').rename({'Scenario':'forcing_mem','Thermal set':'response_mem'}).to_netcdf('../Parameter_Sets/NROY_GWI_informed/GWI_components/nat_coefs_IV_'+select_dataset+'.nc')
    
    xr.DataArray(R2_vals_iv,
                 dims=['int_var_mem','index'],
                 coords={'int_var_mem':CMIP6_int_var_samples_nondeg.columns,'index':temp_nat.columns}).unstack('index').rename({'Scenario':'forcing_mem','Thermal set':'response_mem'}).to_netcdf('../Parameter_Sets/NROY_GWI_informed/GWI_components/R2_IV_'+select_dataset+'.nc')
for select_dataset in ['HadCRUT5','HadCRUT4','NOAA','GISTEMP','CW','BERKELEY']:
    print('calculating GWI for '+select_dataset)
    run_GWI(select_dataset,temp_ant,temp_nat)<jupyter_output>calculating GWI for HadCRUT5<jupyter_text>Regress temp profiles vs observed temps<jupyter_code># ant_coefs = np.empty((GMST.shape[1],temp_ant.shape[1]))
# nat_coefs = np.empty((GMST.shape[1],temp_nat.shape[1]))
# # R2_vals = np.zeros(1000)

# temp_vstack = np.array([temp_ant.loc[GMST.index].values,temp_nat.loc[GMST.index].values]).T

# for i in tqdm.tqdm(np.arange(ant_coefs.shape[1])):
#     _lreg_coefs = OLSE_NORM(temp_vstack[i],GMST.values)
#     ant_coefs[:,i] = _lreg_coefs['coefs'][0]
#     nat_coefs[:,i] = _lreg_coefs['coefs'][1]
# #     R2_vals = _lreg_coefs['R2']<jupyter_output>100%|██████████| 180000/180000 [00:21<00:00, 8566.76it/s]<jupyter_text>save OLSE coefficients to netcdfs<jupyter_code># xr.DataArray(ant_coefs,
#              dims=[GMST.columns.name,'index'],
#              coords={GMST.columns.name:GMST.columns,'index':temp_ant.columns}).unstack('index').rename({'Scenario':'forcing_mem','Thermal set':'response_mem'}).to_netcdf('../Parameter_Sets/NROY_GWI_informed/GWI_components/ant_coefs_forc_'+select_dataset.value+'.nc')

# xr.DataArray(nat_coefs,
#              dims=[GMST.columns.name,'index'],
#              coords={GMST.columns.name:GMST.columns,'index':temp_nat.columns}).unstack('index').rename({'Scenario':'forcing_mem','Thermal set':'response_mem'}).to_netcdf('../Parameter_Sets/NROY_GWI_informed/GWI_components/nat_coefs_forc_'+select_dataset.value+'.nc')<jupyter_output><empty_output><jupyter_text>Regress temp profiles vs internal variability<jupyter_code># ant_coefs_iv = np.empty((CMIP6_int_var_samples_nondeg.shape[1],temp_ant.shape[1]))
# nat_coefs_iv = np.empty((CMIP6_int_var_samples_nondeg.shape[1],temp_nat.shape[1]))

# for i in tqdm.tqdm(np.arange(ant_coefs_iv.shape[1])):
#     _lreg_coefs = OLSE_NORM(temp_vstack[i],CMIP6_int_var_samples_nondeg.loc[GMST.index].values)['coefs']
#     ant_coefs_iv[:,i] = _lreg_coefs[0]
#     nat_coefs_iv[:,i] = _lreg_coefs[1]
# xr.DataArray(ant_coefs_iv,
#              dims=['int_var_mem','index'],
#              coords={'int_var_mem':CMIP6_int_var_samples_nondeg.columns,'index':temp_ant.columns}).unstack('index').rename({'Scenario':'forcing_mem','Thermal set':'response_mem'}).to_netcdf('../Parameter_Sets/NROY_GWI_informed/GWI_components/ant_coefs_IV_'+select_dataset.value+'.nc')

# xr.DataArray(nat_coefs_iv,
#              dims=['int_var_mem','index'],
#              coords={'int_var_mem':CMIP6_int_var_samples_nondeg.columns,'index':temp_nat.columns}).unstack('index').rename({'Scenario':'forcing_mem','Thermal set':'response_mem'}).to_netcdf('../Parameter_Sets/NROY_GWI_informed/GWI_components/nat_coefs_IV_'+select_dataset.value+'.nc')<jupyter_output><empty_output><jupyter_text>Create subsamples<jupyter_code>ant_temps = xr.open_dataarray('../Parameter_Sets/NROY_GWI_informed/GWI_components/ant_temperature.nc',chunks={'forcing_mem':100}).astype(np.single)
nat_temps = xr.open_dataarray('../Parameter_Sets/NROY_GWI_informed/GWI_components/nat_temperature.nc').astype(np.single)
select_dataset = widgets.ToggleButtons(options=['HadCRUT5','HadCRUT4','NOAA','GISTEMP','CW','BERKELEY'],value='HadCRUT5',description='Choose observational timeseries:')
select_dataset
for select_dataset in ['HadCRUT5','HadCRUT4','NOAA','GISTEMP','CW','BERKELEY']:
    
    ant_coefs = xr.open_dataarray('../Parameter_Sets/NROY_GWI_informed/GWI_components/ant_coefs_forc_'+select_dataset+'.nc',chunks={'forcing_mem':100}).astype(np.single)
    nat_coefs = xr.open_dataarray('../Parameter_Sets/NROY_GWI_informed/GWI_components/nat_coefs_forc_'+select_dataset+'.nc',chunks={'forcing_mem':100}).astype(np.single)

    ant_coefs_iv = xr.open_dataarray('../Parameter_Sets/NROY_GWI_informed/GWI_components/ant_coefs_IV_'+select_dataset+'.nc',chunks={'forcing_mem':100}).astype(np.single)
    nat_coefs_iv = xr.open_dataarray('../Parameter_Sets/NROY_GWI_informed/GWI_components/nat_coefs_IV_'+select_dataset+'.nc',chunks={'forcing_mem':100}).astype(np.single)
    
    ## need to subsample otherwise RAM goes boo
    N_subsample = int(1e8)
    random_subsample = np.random.choice((ant_coefs*ant_coefs_iv).size,N_subsample)
    
    np.save(
        '../Parameter_Sets/NROY_GWI_informed/GWI_components/subsamples/AWI_'+select_dataset,
        ((ant_coefs_iv+ant_coefs)*(ant_temps.sel(time=slice('2010','2019')).mean('time')-ant_temps.sel(time=slice('1850','1900')).mean('time'))).values.flatten()[random_subsample]
    )
    
    np.save(
        '../Parameter_Sets/NROY_GWI_informed/GWI_components/subsamples/NWI_'+select_dataset,
        ((nat_coefs_iv+nat_coefs)*(nat_temps.sel(time=slice('2010','2019')).mean('time')-nat_temps.sel(time=slice('1850','1900')).mean('time'))).values.flatten()[random_subsample]
    )
    
    ## create rate arrays
    ant_rate=(ant_temps.sel(time=slice('2010','2019')).assign_coords(time=np.arange(10)).polyfit(dim='time',deg=1).polyfit_coefficients.sel(degree=1)).astype(np.single).compute()
    nat_rate=(nat_temps.sel(time=slice('2010','2019')).assign_coords(time=np.arange(10)).polyfit(dim='time',deg=1).polyfit_coefficients.sel(degree=1)).astype(np.single).compute()
    
    np.save(
        '../Parameter_Sets/NROY_GWI_informed/GWI_components/subsamples/AWI_r_'+select_dataset,
        ((ant_coefs_iv+ant_coefs)*ant_rate).values.flatten()[random_subsample]
    )
    
    np.save(
        '../Parameter_Sets/NROY_GWI_informed/GWI_components/subsamples/NWI_r_'+select_dataset,
        ((nat_coefs_iv+nat_coefs)*nat_rate).values.flatten()[random_subsample]
    )<jupyter_output><empty_output><jupyter_text>analyse subsamples<jupyter_code>## subsubsample for kde analysis
subsubsample = np.random.choice(int(1e8),int(1e6))

GWI_results = pd.concat([pd.concat([pd.Series(np.load('../Parameter_Sets/NROY_GWI_informed/GWI_components/subsamples/'+x+'_'+y+'.npy')[subsubsample]) for x in ['AWI','AWI_r','NWI','NWI_r']],axis=1,keys=['ant','ant_r','nat','nat_r']) for y in ['HadCRUT5','HadCRUT4','NOAA','GISTEMP','CW','BERKELEY']],axis=0,keys=['HadCRUT5','HadCRUT4','NOAA','GISTEMP','CW','BERKELEY']).reset_index(level=0).rename(dict(level_0='dataset'),axis=1)
def kdeplot_norm(y,norm=100,clip=False,vertical=False,quantiles=[[0.05,0.95]],display='mean',**kwargs):
    if clip:
        y = y[(y>clip[0])&(y<clip[1])]
    plot_range = np.linspace(y.min(),y.max(),500)
    y_kde = sp.stats.gaussian_kde(y)
    y_pdf = y_kde(plot_range)/norm
    
    if vertical:
        plt.fill_betweenx(plot_range,y_pdf,0,**kwargs,alpha=0.1,lw=0)
    else:
        plt.fill_between(plot_range,y_pdf,0,**kwargs,alpha=0.1,lw=0)
    
    for quantile_band in quantiles:
        y_quantiles = np.quantile(y,quantile_band)
        quantile_range = np.linspace(*y_quantiles,500)
        quantile_pdf = y_kde(quantile_range)/norm
        if vertical:
            plt.fill_betweenx(quantile_range,quantile_pdf,0,alpha=0.1,lw=0,**kwargs)
        else:
            plt.fill_between(quantile_range,quantile_pdf,0,alpha=0.1,lw=0,**kwargs)
        
    if display == 'mean':
        if vertical:
            plt.hlines(y=y.mean(),xmin=0,xmax=y_kde(y.mean())/norm,ls='--',lw=1,**kwargs)
        else:
            plt.vlines(x=y.mean(),ymin=0,ymax=y_kde(y.mean())/norm,ls='--',lw=1,**kwargs)
        
    elif type(display)==float:
        if vertical:
            plt.hlines(y=np.quantile(y,display),xmin=0,xmax=y_kde(np.quantile(y,display))/norm,ls='--',lw=1,**kwargs)
        else:
            plt.vlines(x=np.quantile(y,display),ymin=0,ymax=y_kde(np.quantile(y,display))/norm,ls='--',lw=1,**kwargs)
            
    if vertical:
        plt.plot(y_pdf,plot_range,zorder=2,lw=1.5,**kwargs)
    else:
        plt.plot(plot_range,y_pdf,zorder=2,lw=1.5,**kwargs)
def plot_crosstime(x,y,**kwargs):
    
    ## plot data lines:
    
    time_to_15 = (1.5-x)/y
    
    quantiles = np.quantile(time_to_15,[0.17,0.5,0.83])
    
#     plt.plot([1.5,0],[0,1.5/quantiles[1]],**kwargs)
    
    plt.fill_between([1.5,0],[0,1.5/quantiles[0]],[0,1.5/quantiles[2]],alpha=0.2,lw=0,**kwargs)
    
    ## add on labels:
    
    X = np.linspace(0,1.6,100)
    Y = np.linspace(0,0.1,100)
    X,Y=np.meshgrid(X,Y)
    Z = (1.5-X)/Y
    contours = plt.contour(X,Y,Z,levels=quantiles,linewidths=[0,1.5,0],colors=kwargs['color'])
    labels = [str(round(x,1))+' years'  for x in quantiles]
    y0=0
    x0=1.5
    rx=2
    ry=0.06
    def find_loc(y,Z):
        x=1.5-Z*y
        return 1 - ((y-y0)/ry)**2 - ((x-x0)/rx)**2
    y_vals_levels = np.array([sp.optimize.root(lambda x:find_loc(x,quantiles[i]),x0=0.5).x[0] for i in np.arange(3)])
    x_vals_levels = 1.5 - y_vals_levels*quantiles
    plt.clabel(contours,fmt=dict(zip(contours.levels,labels)),  manual = list(zip(x_vals_levels,y_vals_levels)),colors='k')
def hexbin(x, y, **kwargs):
    plt.hexbin(x, y, **kwargs)
g=sn.FacetGrid(data=GWI_results,col='dataset',col_wrap=3,height=5)
g.map(hexbin,'ant','ant_r',cmap='binary',gridsize=50,extent=[0,1.6,0,0.07],linewidths=0)
g.map(kdeplot_norm,'ant',norm=350,clip=[0,1.6],quantiles=[[0.05,0.95],[0.17,0.83]],display=0.5,color='k')
g.map(kdeplot_norm,'ant_r',norm=200,clip=[0,0.1],quantiles=[[0.05,0.95],[0.17,0.83]],display=0.5,vertical=True,color='k')

[a.set_xlim(0,1.6) for a in g.axes.flatten()]
[a.set_ylim(0,0.07) for a in g.axes.flatten()]

g.map(plot_crosstime,'ant','ant_r',color='cornflowerblue')

[a.set_ylabel('Anthropogenic contribution to 2010-2019 rate / K yr$^{-1}$') for a in g.axes[[0,3]]]
[a.set_xlabel('Anthropogenic contribution to 2010-2019 level / K') for a in g.axes[-3:]]
''
((1.5-GWI_results.set_index('dataset').ant) / GWI_results.set_index('dataset').ant_r).groupby(level=0).quantile([0.05,0.17,0.5,0.83,0.95]).unstack(level=0)<jupyter_output><empty_output><jupyter_text>R$^{2}$ for each dataset<jupyter_code>fig,ax = plt.subplots(1,figsize=(12,8))

[(xr.open_dataset('../Parameter_Sets/NROY_GWI_informed/GWI_components/R2_forc_'+x+'.nc').__xarray_dataarray_variable__).plot.hist(ax=ax,histtype='step',cumulative=True,density=True,label=x) for x in ['HadCRUT5','HadCRUT4','NOAA','GISTEMP','CW','BERKELEY']]
def forward(x):
    return np.log(1-x)
def inverse(x):
    return 1-np.exp(x)
plt.xscale('function',functions=(forward,inverse))
plt.xlim(0.95,0.05)
plt.ylim(0,1)
plt.axhline(0.5,lw=1,ls=':',alpha=0.5)
plt.title('forced regression R$^2$')
plt.legend(loc=2)
''
### OLD ANALYSIS<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1864093567251462, 'fraction_numerical': 0.05742690058479532, 'mean_word_length': 5.255633596722271, 'pattern_counts': {'":': 2, '<': 70, '<?xml version=': 0, '>': 86, 'https://': 13, 'lorem ipsum': 0, 'www.': 8, 'xml': 1}, 'pii_count': 5, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '17', 'dir': 'jupyter-structured-clean-dedup', 'id': '195643', 'license': 'CC-BY-4.0', 'n_tokens_mistral': 18904, 'n_tokens_neox': 17393, 'n_words': 1660, 'path': 'GIR/notebooks/GWI_generation.ipynb', 'repo_name': 'njleach/GIR', 'size': '42750'}