# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python [conda env:miniconda3-lens-py]
#     language: python
#     name: conda-env-miniconda3-lens-py-py
# ---

import numpy as np
import pandas as pd
import xarray as xr
import glob
import matplotlib
import matplotlib.pyplot as plt
import cftime
import dask
import statsmodels.api as sm

import os
from ppe_tools import Ensemble,Member,ParamInfo
from ppe_tools.utils import get_default, parse_val


# Updated for PBS scheduler
# this could go into utils.
# By default gets 1 core w/ 25 GB memory
def get_ClusterClient(ncores=1, nmem='25GB'):
    import dask
    from dask_jobqueue import PBSCluster
    from dask.distributed import Client
    ncores=ncores
    nmem = nmem

    cluster = PBSCluster(
        cores=ncores, # The number of cores you want
        memory=nmem, # Amount of memory
        processes=ncores, # How many processes
        queue='casper', # The type of queue to utilize (/glade/u/apps/dav/opt/usr/bin/execcasper)
        resource_spec='select=1:ncpus='+str(ncores)+':mem='+nmem, # Specify resources
        project='P93300641', # Input your project ID here
        walltime='2:30:00', # Amount of wall time
        interface='ib0', # Interface to use
    )

    client = Client(cluster)
    return cluster, client


cluster, client = get_ClusterClient()
cluster.scale(20)

client


def get_params(keys,paramkey):
    params=[]
    minmaxs=[]
    for key in keys:
        ix     = paramkey.key==key
        params.append(paramkey.param[ix].values[0])
        m      = paramkey.minmax[ix].values[0]
        if 'min' in m:
            minmaxs.append('min')
        elif 'max' in m:
            minmaxs.append('max')
        else:
            minmaxs.append(paramkey.minmax[ix].values[0])
    return params,minmaxs


def get_ensemble(name,data_vars,htape,keys,paramkey):
    '''
    Returns a dataset comprising the given ensemble
    name  = 'CTL2020' or 'C285' or 'AF1905'
    data_vars, e.g. ['GPP','HR','AR']
    htape, e.g. 'h0' 0/1/2/3/4/5/7 available
    '''
    
    #define the directory structure and find files
    topdir     = '/glade/scratch/djk2120/PPEn11/hist/' 
    thisdir    = topdir+name+'/'
    files      = [glob.glob(thisdir+'*'+key+'*'+htape+'*.nc')[0] for key in keys]

    def preprocess(ds):
        return ds[data_vars]

    #read in the dataset
    ds = xr.open_mfdataset(files,combine='nested',concat_dim='ens',
                           parallel=True,preprocess=preprocess)

    #fix up time dimension
    if htape=='h0' or htape=='h1':
        ds['time'] = xr.cftime_range(str(2005),periods=120,freq='MS') #fix time bug

    #specify extra variables    
    if htape=='h0':
        extras     = ['grid1d_lat','grid1d_lon']
    elif htape=='h1':
        extras     = ['pfts1d_lat','pfts1d_lon','pfts1d_wtgcell','pfts1d_itype_veg']
    
    #add in some extra variables
    ds0 = xr.open_dataset(files[0])
    for extra in extras:
        ds[extra]=ds0[extra]

    #append some info about key/param/minmax
    params,minmaxs = get_params(keys,paramkey) 
    ds['key']    = xr.DataArray(keys,dims='ens')
    ds['param']  = xr.DataArray(params,dims='ens')
    ds['minmax'] = xr.DataArray(minmaxs,dims='ens')
    
    return ds


# ### round 3

# +
#fetch the paraminfo
csv = '/glade/scratch/djk2120/PPEn11/pass2.csv' 
paramkey = pd.read_csv(csv)

#fetch the sparsegrid landarea
la_file = '/glade/scratch/djk2120/PPEn08/sparsegrid_landarea.nc'
la = xr.open_dataset(la_file).landarea  #km2

# +
envs = ['CTL2010','AF1855','C285']
topdir = '/glade/scratch/djk2120/PPEn11/hist/'

for env in envs:
    files = sorted(glob.glob(topdir+env+'/*h0*'))
    thiskeys  = [file.split('_')[2].split('.')[0] for file in files]
    if env=='CTL2010':
        keys = set(thiskeys)
    keys = keys&set(thiskeys)
    
keys = sorted(list(keys)) 
# -

dsets = []
data_vars=['TLAI']
for env in envs:
    ds = get_ensemble(env,data_vars,'h1',keys,paramkey)
    ds['pft']=ds.pfts1d_itype_veg
    dsets.append(ds)

lai_thresh  = 0.1
surv_thresh = 0.7
survs = []
for i in range(3):
    npix = (dsets[i].TLAI.max(dim='time')>lai_thresh).groupby('pft').sum().compute()
    surv = (npix/npix.isel(ens=0)).compute().isel(pft=np.arange(17)>0)>surv_thresh
    survs.append(surv)
    if i==0:
        surv3 = surv.values
    else:
        surv3 = surv3*surv.values
surv  = xr.DataArray(surv3,dims=['ens','pft'],name='survived')
sflag = surv.sum(dim='pft')==16

plt.figure(figsize=[11,3])
surv.T.plot()
plt.title('CTL2010 & AF1855 & C285');

uniques = []
for p,m in zip(ds.param.values,ds.minmax.values):
    ustr = p+'-'+m
    if ustr not in uniques:
        uniques.append(ustr)
#which params still need resubmission?
for u in uniques:
    p = u.split('-')[0]
    m = u.split('-')[1]
    ix = (ds.param==p) & (ds.minmax==m)
    if True not in sflag.isel(ens=ix):
        print(u)

# +
#create a paramkey with only the first good surviving key
goodkeys=[]
goodps  =[]
goodms  =[]
for u in uniques:
    p = u.split('-')[0]
    m = u.split('-')[1]
    ix = (ds.param==p) & (ds.minmax==m)
    s = sflag[ix]
    if True in s:
        firstk = str(ds.key[ix][s][0].values)  #take the first key where sflag is true
        goodkeys.append(firstk)
        goodps.append(p)
        goodms.append(m)
        
ranks = np.argsort(goodkeys)
goodkeys = [goodkeys[rank] for rank in ranks]
goodps   = [goodps[rank] for rank in ranks]
goodms   = [goodms[rank] for rank in ranks]
survkey = pd.DataFrame({'key':goodkeys,'param':goodps,'minmax':goodms})
# -

# ### examine medlynintercept-max and sucsat-both

data_vars = ['FCTR','FGEV','FCEV','QHR','SOILWATER_10CM','EFLX_LH_TOT']

ix1 = (paramkey.param=='medlynintercept')
ix2 = np.array(['max' in minmax for minmax in paramkey.minmax])
ix3 = (paramkey.param=='sucsat_sf')
ix4 = paramkey.param=='default'
ix = (ix1&ix2)|ix3|ix4
mkeys = paramkey.key[ix]

ds = get_ensemble('CTL2010',data_vars,'h0',mkeys,paramkey)


def month_wts(nyears):
    days_pm  = [31,28,31,30,31,30,31,31,30,31,30,31]
    return xr.DataArray(np.tile(days_pm,nyears),dims='time')


# how often is FCEV<0, as compared to the default case??
x=(la*(ds.FCEV<0)).sum(dim=['gridcell','time']).compute()
for i in range(20):
    if 'medlyn' in str(ds.param.isel(ens=i).values):
        k=str(ds.key.isel(ens=i).values)
        narea = x.isel(ens=i)/x.isel(ens=0)
        print(k,narea.values)

#OAAT0413 looks best
medkey ='OAAT0413'

qhr = 1/la.sum()/365*(month_wts(10)*(la*ds.QHR).sum(dim='gridcell')).groupby('time.year').sum().compute()
plt.figure(figsize=[12,6])
i=0
for ee in range(19):
    p=str(ds.param.isel(ens=ee).values)
    m=str(ds.minmax.isel(ens=ee).values)
    k=str(ds.key.isel(ens=ee).values)
    if ('sucsat' in p) and ('max' in m):
        i+=1
        plt.subplot(2,3,i)
        qhr.isel(ens=ee).plot()
        qhr.isel(ens=0).plot()
        if (i<4):
            plt.xlabel('')
        plt.title(k)
        plt.ylim([-5e-5,0])
plt.subplots_adjust(hspace=0.3)

qhr = 1/la.sum()/365*(month_wts(10)*(la*ds.EFLX_LH_TOT).sum(dim='gridcell')).groupby('time.year').sum().compute()
plt.figure(figsize=[12,6])
i=0
for ee in range(19):
    p=str(ds.param.isel(ens=ee).values)
    m=str(ds.minmax.isel(ens=ee).values)
    k=str(ds.key.isel(ens=ee).values)
    if ('sucsat' in p) and ('max' in m):
        i+=1
        plt.subplot(2,3,i)
        qhr.isel(ens=ee).plot()
        qhr.isel(ens=0).plot()
        if (i<4):
            plt.xlabel('')
        plt.title(k)
        plt.ylim([36,44])
plt.subplots_adjust(hspace=0.3)

qhr = 1/la.sum()/365*(month_wts(10)*(la*ds.QHR).sum(dim='gridcell')).groupby('time.year').sum().compute()
plt.figure(figsize=[12,6])
i=0
for ee in range(19):
    p=str(ds.param.isel(ens=ee).values)
    m=str(ds.minmax.isel(ens=ee).values)
    k=str(ds.key.isel(ens=ee).values)
    if ('sucsat' in p) and ('min' in m):
        i+=1
        plt.subplot(2,3,i)
        qhr.isel(ens=ee).plot()
        qhr.isel(ens=0).plot()
        if (i<4):
            plt.xlabel('')
        plt.title(k)
        plt.ylim([-5e-6,0])
plt.subplots_adjust(hspace=0.3)

qhr = 1/la.sum()/365*(month_wts(10)*(la*ds.EFLX_LH_TOT).sum(dim='gridcell')).groupby('time.year').sum().compute()
plt.figure(figsize=[12,6])
i=0
for ee in range(19):
    p=str(ds.param.isel(ens=ee).values)
    m=str(ds.minmax.isel(ens=ee).values)
    k=str(ds.key.isel(ens=ee).values)
    if ('sucsat' in p) and ('min' in m):
        i+=1
        plt.subplot(2,3,i)
        qhr.isel(ens=ee).plot()
        qhr.isel(ens=0).plot()
        if (i<4):
            plt.xlabel('')
        plt.title(k)
        plt.ylim([36,42])
plt.subplots_adjust(hspace=0.3)

skey_max = 'OAAT0449'
skey_min = 'OAAT0454'

ix = (survkey.param=='medlynintercept')&(survkey.minmax=='max')
survkey.key[ix]=medkey
ix = (survkey.param=='sucsat_sf')&(survkey.minmax=='min')
survkey.key[ix]=skey_min
ix = (survkey.param=='sucsat_sf')&(survkey.minmax=='max')
survkey.key[ix]=skey_max

survkey.to_csv('/glade/scratch/djk2120/PPEn11/surv.csv',index=False)

# + [markdown] tags=[]
# ## look at the original CTL,AF1855, and C285 ensembles

# +
envs = ['CTL2010','AF1855','C285']
topdir = '/glade/scratch/djk2120/PPEn11/hist/'

for env in envs:
    files = sorted(glob.glob(topdir+env+'/*h0*'))
    thiskeys  = [file.split('_')[2].split('.')[0] for file in files]
    if env=='CTL2010':
        keys = set(thiskeys)
    keys = keys&set(thiskeys)
    
keys = sorted(list(keys))       
# -

dsets = []
for env in envs:
    ds = get_ensemble(env,data_vars,'h1',keys,paramkey)
    ds['pft']=ds.pfts1d_itype_veg
    dsets.append(ds)

lai_thresh  = 0.1
surv_thresh = 0.7
survs = []
for i in range(3):
    npix = (dsets[i].TLAI.max(dim='time')>lai_thresh).groupby('pft').sum().compute()
    surv = (npix/npix.isel(ens=0)).compute().isel(pft=np.arange(17)>0)>surv_thresh
    survs.append(surv)
    if i==0:
        surv3 = surv.values
    else:
        surv3 = surv3*surv.values
surv  = xr.DataArray(surv3,dims=['ens','pft'],name='survived')
sflag = surv.sum(dim='pft')==16

plt.figure(figsize=[11,3])
surv.T.plot()
plt.title('CTL2010 & AF1855 & C285');

# ### create new paramfiles to contrain ranges
#  - submitting three more per dead member (to save time)
#  - it happened to be that no namelist params died
#      - so this will break if trying to constrain namelist params

# +
#instantiate the Ensemble object
basefile = '/glade/p/cesm/cseg/inputdata/lnd/clm2/paramdata/ctsm51_params.c210507.nc'
lndin    = '/glade/work/oleson/lmbirch_wkattge.n01_ctsm5.1.dev006/cime/scripts/clm51_lmbirchwkattgen01ctsm51d006_2deg_GSWP3V1_PPE_1850pAD/CaseDocs/lnd_in'
pdir     = '/glade/scratch/djk2120/PPEn11/paramfiles/'
ndir     = '/glade/scratch/djk2120/PPEn11/namelist_mods/'
surv2    = Ensemble(basefile,pdir,ndir)


defp = xr.open_dataset(pdir+'OAAT0000.nc')
ot = 1+np.arange(16)
ds = dsets[0]
nextk = 400
kvals=[];pvals=[];mvals=[]

for w1,howmany in zip([0.5,0.25,0.125],[2,3,4]):
    w0 = 1-w1
    for i,k in zip(ds.ens.values, ds.key.values):
        if not sflag.isel(ens=i):
            p=ds.param.isel(ens=i).values.ravel()[0]
            m=ds.minmax.isel(ens=i).values.ravel()[0]

            p0 = defp[p]
            p1 = xr.open_dataset(pdir+k+'.nc')[p]

            if p1.shape: #param has pft dimension                    
                ix = ot[~surv.isel(ens=i).values] #which pfts died?
                #constrain dead pfts
                newp = p1.copy(deep=True)
                newp[ix] = w0*p0[ix]+w1*p1[ix] 
            else: #param is scalar
                newp = w0*p0+w1*p1
                
            loc    = 'P'
            defval = get_default(p, loc, defp, lndin)
            value  = parse_val(loc,defval,newp)
            minmax = m+str(howmany)
            key    = 'OAAT'+str(nextk).zfill(4)
            
            paraminfo = ParamInfo(p, loc, defval, value)
            paramdict = {p:paraminfo}
            member    = Member(key,paramdict,basefile,minmax)   
            surv2.add_member(member)
            
            kvals.append(key)
            pvals.append(p)
            mvals.append(minmax)
            
            nextk+=1
            
surv2.write()
newcsv=pd.DataFrame({'key':kvals,'param':pvals,'minmax':mvals})
newcsv.to_csv('/glade/scratch/djk2120/PPEn11/surv_resubmit.csv',index=False)

# +
#spot check a pft variable
pk = pd.read_csv('/glade/scratch/djk2120/PPEn11/surv_resubmit.csv')
thisp = 'medlynslope'
plt.plot(defp[thisp][:17])
p1 = xr.open_dataset(pdir+'OAAT0117.nc')[thisp]

plt.plot(p1[:17])
for p,k in zip(pk.param.values,pk.key.values):
    if p==thisp:
        newp = xr.open_dataset(pdir+k+'.nc')[p]
        plt.plot(newp[:17])
plt.legend(range(5));
# -

#spot check the non-pft param
print('          original perturbation -----> default value')
parms = list(np.unique(pk.param.values))
parms.remove('medlynslope')
parms.remove('medlynintercept')
for parm in parms:
    ix = pk.param==parm
    mm = pk.minmax[ix].values[0][:3]
    ix = np.logical_and(ds.param==parm,ds.minmax==mm)
    k1 = ds.key.isel(ens=ix).values[0]
    p1 = xr.open_dataset(pdir+k1+'.nc')[parm]
    
    thisstr = parm+': '+str(p1.values)
    ix = pk.param==parm
    for key in pk.key[ix]:        
        newp = xr.open_dataset(pdir+key+'.nc')[parm]
        thisstr = thisstr+', '+str(newp.values)
    thisstr = thisstr+', '+str(defp[parm].values)    
    
    print(thisstr)
    

x = [75325.0, 112987.5, 131818.75, 141234.375, 150650.0]
np.array(x)/150650

# ### second round of resubmissions
# - I handled medlynintercept separately, in the next section

# +
#fetch the paraminfo
csv = '/glade/scratch/djk2120/PPEn11/full.csv' 
paramkey = pd.read_csv(csv)

#fetch the sparsegrid landarea
la_file = '/glade/scratch/djk2120/PPEn08/sparsegrid_landarea.nc'
la = xr.open_dataset(la_file).landarea  #km2
# -

ix = paramkey.param=='medlynintercept'
paramkey.minmax[ix]

envs = ['CTL2010','C285','AF1855']
datavars=['TLAI']
dsets = []
for env in envs:
    ds = get_ensemble(env,datavars,'h1',paramkey.key,paramkey)
    ds['pft']=ds.pfts1d_itype_veg
    dsets.append(ds)

lai_thresh  = 0.1
surv_thresh = 0.7
survs = []
for i in range(3):
    npix = (dsets[i].TLAI.max(dim='time')>lai_thresh).groupby('pft').sum().compute()
    surv = (npix/npix.isel(ens=0)).compute().isel(pft=np.arange(17)>0)>surv_thresh
    survs.append(surv)
    if i==0:
        surv3 = surv.values
    else:
        surv3 = surv3*surv.values
surv  = xr.DataArray(surv3,dims=['ens','pft'],name='survived')
sflag = surv.sum(dim='pft')==16

uniques = []
for p,m in zip(ds.param.values,ds.minmax.values):
    ustr = p+'-'+m
    if ustr not in uniques:
        uniques.append(ustr)
#which params still need resubmission?
for u in uniques:
    p = u.split('-')[0]
    m = u.split('-')[1]
    ix = (ds.param==p) & (ds.minmax==m)
    if True not in sflag.isel(ens=ix):
        print(u)

#which params still need resubmission?
for u in uniques:
    p = u.split('-')[0]
    m = u.split('-')[1]
    ix = (ds.param==p) & (ds.minmax==m)
    if True not in sflag.isel(ens=ix):
        print(u)

# +
#create a paramkey with only the first good surviving key
goodkeys=[]
goodps  =[]
goodms  =[]
for u in uniques:
    p = u.split('-')[0]
    m = u.split('-')[1]
    ix = (ds.param==p) & (ds.minmax==m)
    s = sflag[ix]
    if True in s:
        firstk = str(ds.key[ix][s][0].values)  #take the first key where sflag is true
        goodkeys.append(firstk)
        goodps.append(p)
        goodms.append(m)
        
ranks = np.argsort(goodkeys)
goodkeys = [goodkeys[rank] for rank in ranks]
goodps   = [goodps[rank] for rank in ranks]
goodms   = [goodms[rank] for rank in ranks]

survkey = pd.DataFrame({'key':goodkeys,'param':goodps,'minmax':goodms})
survkey.to_csv('/glade/scratch/djk2120/PPEn11/surv.csv',index=False)


# +

nextk

# +
constrains=['vcmaxhd-min','lmrhd-min','vcmaxse_sf-max']
nextk = int(sorted(glob.glob('/glade/scratch/djk2120/PPEn11/paramfiles/OAAT*.nc'))[-1].split('/')[-1].split('.')[0][4:])+1

#instantiate the Ensemble object
basefile = '/glade/p/cesm/cseg/inputdata/lnd/clm2/paramdata/ctsm51_params.c210507.nc'
lndin    = '/glade/work/oleson/lmbirch_wkattge.n01_ctsm5.1.dev006/cime/scripts/clm51_lmbirchwkattgen01ctsm51d006_2deg_GSWP3V1_PPE_1850pAD/CaseDocs/lnd_in'
pdir     = '/glade/scratch/djk2120/PPEn11/paramfiles/'
ndir     = '/glade/scratch/djk2120/PPEn11/namelist_mods/'
surv3    = Ensemble(basefile,pdir,ndir)

f = open('/glade/scratch/djk2120/PPEn11/surv_resubmit2.csv', "a")
f.write('key,param,minmax\n')


defp = xr.open_dataset(pdir+'OAAT0000.nc')

for w1,howmany in zip([0.5,0.25,0.125],[5,6,7]):
    w0 = 1-w1

    for u in constrains:
        p = u.split('-')[0]
        m = u.split('-')[1]
        ix = (ds.param==p) & (ds.minmax==m)
        k = str(ds.key[ix][-1].values)
        
        p0 = defp[p]
        p1 = xr.open_dataset(pdir+k+'.nc')[p]
        newp = w0*p0+w1*p1
        
        
        loc    = 'P'
        defval = get_default(p, loc, defp, lndin)
        value  = parse_val(loc,defval,newp)
        minmax = m+str(howmany)
        key    = 'OAAT'+str(nextk).zfill(4)

        paraminfo = ParamInfo(p, loc, defval, value)
        paramdict = {p:paraminfo}
        member    = Member(key,paramdict,basefile,minmax)   
        
        
        surv3.add_member(member)
        
        f.write(key+','+p+','+minmax+'\n')
        nextk+=1
      
f.close()
# -

surv3.write()


# + [markdown] tags=[]
# ### investigate medlynintercept
#  - can lead to global negative FCEV
# -

datavars = ['FCEV']
envs = ['CTL2010','C285','AF1855']
fcevs=[]
for env in envs:
    ds = get_ensemble(env,datavars,'h0',paramkey.key,paramkey)
    fcevs.append(ds)

ds = fcevs[0]
p = 'medlynintercept'
m = 'max'
ix = np.logical_and(ds.param==p,ds.minmax==mm)
ix = np.logical_or(ix,ds.param=='default')

plt.figure(figsize=[10,6])
for i in range(5):
    plt.subplot(2,3,i+1)
    fcev_neg = ds.FCEV.isel(ens=ix).isel(ens=i)<-5
    fcev_neg.plot()


def month_wts(nyears):
    days_pm  = [31,28,31,30,31,30,31,31,30,31,30,31]
    return xr.DataArray(np.tile(days_pm,nyears),dims='time')


cf = 1/la.sum()/365
fcev_glob = cf*(month_wts(10)*(la*ds.FCEV.isel(ens=ix)).sum(dim='gridcell')).groupby('time.year').sum().mean(dim='year').compute()
fcev_glob

medkeys = ds.key[ix].values
for key in medkeys:
    pdir = '/glade/scratch/djk2120/PPEn11/paramfiles/'
    f = pdir+key+'.nc'
    p = xr.open_dataset(f)['medlynintercept']
    m = np.mean(p[1:17]).values
    print(key,m)

# ### I think we should redo medlynintercept constraining all pfts at once

# +
kvals = [401,413,425,436,437,438]
mvals = 100000/(2**np.arange(6))
iters = [2,3,4,5,6,7]

#instantiate the Ensemble object
basefile = '/glade/p/cesm/cseg/inputdata/lnd/clm2/paramdata/ctsm51_params.c210507.nc'
lndin    = '/glade/work/oleson/lmbirch_wkattge.n01_ctsm5.1.dev006/cime/scripts/clm51_lmbirchwkattgen01ctsm51d006_2deg_GSWP3V1_PPE_1850pAD/CaseDocs/lnd_in'
pdir     = '/glade/scratch/djk2120/PPEn11/paramfiles/'
ndir     = '/glade/scratch/djk2120/PPEn11/namelist_mods/'
medg0    = Ensemble(basefile,pdir,ndir)

for mval,kval,howmany in zip(mvals,kvals,iters):

    defp = xr.open_dataset(pdir+'OAAT0000.nc')
    p = 'medlynintercept'
    p0 = defp[p]
    newp = p0.copy(deep=True)
    newp[:] = mval


    loc    = 'P'
    defval = get_default(p, loc, defp, lndin)
    value  = parse_val(loc,defval,newp)
    minmax = 'max'+str(howmany)
    key    = 'OAAT'+str(kval).zfill(4)

    paraminfo = ParamInfo(p, loc, defval, value)
    paramdict = {p:paraminfo}
    member    = Member(key,paramdict,basefile,minmax)   
    medg0.add_member(member)
# -

medg0.write()

pp = xr.open_dataset('/glade/scratch/djk2120/PPEn11/paramfiles/OAAT0401.nc')['medlynintercept']
pp

u  = 'medlynintercept'
mm = 'max'
ix = np.logical_and(ds.param==u,ds.minmax==mm)
ix = np.logical_or(ix,ds.param=='default')
npixs = []
survs = []
for ds in dsets:
    npix = (ds.TLAI.isel(ens=ix).max(dim='time')>lai_thresh).groupby('pft').sum().compute()
    npixs.append(npix)
    surv = (npix/npix.isel(ens=0)).compute().isel(pft=np.arange(17)>0)>surv_thresh
    survs.append(surv)

surv

plt.figure(figsize=[10,8])
for i in 1+np.arange(16):
    plt.subplot(4,4,i)
    for j in range(3):
        plt.plot(survs[j].isel(pft=i-1),'-x')

plt.figure(figsize=[10,8])
for i in 1+np.arange(16):
    plt.subplot(4,4,i)
    for j in range(3):
        plt.plot(npixs[j].isel(pft=i),'-x')

# ### investigate sucsat_sf

ds.key[ds.param=='sucsat_sf']

datavars = ['FCTR','FGEV','QHR']
keys = ['OAAT0000','OAAT0066','OAAT0067']
ds = get_ensemble('CTL2010',datavars,'h0',keys,paramkey)
ds

cf = 24*60*60/la.sum()/365
qhr = cf*(month_wts(10)*(la*ds.QHR).sum(dim='gridcell')).groupby('time.year').sum()
plt.plot(1+np.arange(10),-qhr,'-x')
plt.ylim([0,5])
plt.ylabel('HR (mm/d)')
plt.xlabel('year');

# +
sucsat_high = 1+8/2**np.arange(6)[1:]
sucsat_low = 1/sucsat_high
svals = [*sucsat_high,*sucsat_low]
iters = 2+np.arange(5)
iters = [*iters,*iters]

nextk = int(sorted(glob.glob('/glade/scratch/djk2120/PPEn11/paramfiles/OAAT*.nc'))[-1].split('/')[-1].split('.')[0][4:])+1
#instantiate the Ensemble object
basefile = '/glade/p/cesm/cseg/inputdata/lnd/clm2/paramdata/ctsm51_params.c210507.nc'
lndin    = '/glade/work/oleson/lmbirch_wkattge.n01_ctsm5.1.dev006/cime/scripts/clm51_lmbirchwkattgen01ctsm51d006_2deg_GSWP3V1_PPE_1850pAD/CaseDocs/lnd_in'
pdir     = '/glade/scratch/djk2120/PPEn11/paramfiles/'
ndir     = '/glade/scratch/djk2120/PPEn11/namelist_mods/'
scsat    = Ensemble(basefile,pdir,ndir)

for sval,howmany in zip(svals,iters):

    if sval>1:
        m='max'
    else:
        m='min'
    
    defp = xr.open_dataset(pdir+'OAAT0000.nc')
    p = 'sucsat_sf'
    p0 = defp[p]
    newp = p0.copy(deep=True)
    newp = sval


    loc    = 'P'
    defval = get_default(p, loc, defp, lndin)
    value  = parse_val(loc,defval,newp)
    minmax = m+str(howmany)
    key    = 'OAAT'+str(nextk).zfill(4)

    paraminfo = ParamInfo(p, loc, defval, value)
    paramdict = {p:paraminfo}
    member    = Member(key,paramdict,basefile,minmax)   
    scsat.add_member(member)
    
    nextk+=1
# -

scsat.write(csvfile='/glade/scratch/djk2120/PPEn11/sucsat.csv')

cf = 1/la.sum()/365
fctr = cf*(month_wts(10)*(la*ds.FCTR).sum(dim='gridcell')).groupby('time.year').sum()
plt.plot(1+np.arange(10),fctr,'-x')
#plt.ylim([0,5])
plt.ylabel('FCTR (W/m2)')
plt.xlabel('year');

cf = 1/la.sum()/365
fgev = cf*(month_wts(10)*(la*ds.FGEV).sum(dim='gridcell')).groupby('time.year').sum()
plt.plot(1+np.arange(10),fgev,'-x')
#plt.ylim([0,5])
plt.ylabel('FCTR (W/m2)')
plt.xlabel('year');

# ## deprecated

env = 'CTL2010'
thedir = '/glade/scratch/djk2120/PPEn11/hist/'+env+'/'
laithresh = 0.1
files = sorted(glob.glob(thedir+'*h1*nc'))
keys  = [file.split('_')[2].split('.')[0] for file in files]

ctl    = get_ensemble('CTL2010',['TLAI','GPP'],'h1',paramkey)
c285   = get_ensemble('C285'   ,['TLAI','GPP'],'h1',paramkey)
af1905 = get_ensemble('AF1905' ,['TLAI','GPP'],'h1',paramkey)
dsets  = [ctl,c285,af1905]
envs   = ['CTL','C285','AF1905']

survs = []
for ds in dsets:
    npix = (ds.TLAI.max(dim='time')>laithresh).groupby('pft').sum().compute()
    surv = (npix/npix.isel(ens=0)).compute().isel(pft=np.arange(17)>0)>0.7
    survs.append(surv)

ot = 1+np.arange(16)
resubmit = dict()
for ds,surv,env in zip(dsets,survs,envs):
    dead = surv.sum(dim='pft')<16
    ivals = ds.ens.isel(ens=dead)
    for ival in ivals:
        p  = str(ds.param.isel(ens=ival).values)
        m  = ds.minmax.isel(ens=ival).values
        k  = ds.key.isel(ens=ival).values
        kn = int(str(k)[4:])
        if kn<400:
            ix = np.logical_and(ds.param==p,ds.minmax==m) 
            if False not in dead.isel(ens=ix):
                i = np.max(ds.ens.isel(ens=ix))
                deadpfts = surv.isel(ens=i).values<0.7
                thisp = p+'-'+str(m)
                print(env+':',p,m,'unresolved')
                if thisp not in resubmit:
                    resubmit[thisp] = list(ot[deadpfts])
                else:
                    x0 = resubmit[thisp]
                    x  = x0+list(ot[deadpfts])
                    resubmit[thisp] = list(np.unique(x))

topdir = '/glade/scratch/djk2120/PPEn11/hist/'
envs = ['CTL2010','C285', 'AF1905']
for env in envs:
    tmp = [file.split('_')[2].split('.')[0] for file in glob.glob(topdir+env+'/*h1*')]
    if env=='CTL2010':
        keys=set(tmp)
    else:
        keys=keys&set(tmp)
keys = sorted(list(keys))

dsets = []
for env in envs:
    ds=get_ensemble(env,keys,data_vars,'h1',paramkey)
    dsets.append(ds)

survs = []
for ds in dsets:
    npix = (ds.TLAI.max(dim='time')>laithresh).groupby('pft').sum().compute()
    surv = (npix/npix.isel(ens=0)).compute().isel(pft=np.arange(17)>0)>0.7
    survs.append(surv)

alive = survs[0].sum(dim='pft')>-1
for surv in survs:
    alive = np.logical_and(alive,surv.sum(dim='pft')==16)

goods = dict()
bads  = []
goodkeys=[]
goodmins=[]
goodps=[]
for i in alive.ens:
    if alive.sel(ens=i):
        p=str(ds.param.sel(ens=i).values)+'-'+str(ds.minmax.sel(ens=i).values)
        if p not in goods:
            goods[p]=str(ds.key.sel(ens=i).values)
            goodps.append(  str(ds.param.sel(ens=i).values))
            goodkeys.append(str(ds.key.sel(ens=i).values))
            goodmins.append(str(ds.minmax.sel(ens=i).values))
        if p in bads:
            bads.remove(p)
    elif p not in goods:
        bads.append(p)

survkey = pd.DataFrame({'key':goodkeys,'param':goodps,'minmax':goodmins})
survkey

survkey.to_csv('/glade/scratch/djk2120/PPEn11/surv.csv',index=False)

npix0 = npix.isel(ens=0).values[1:]
print(npix0)
for i in ds.ens:
    if ds.param.isel(ens=i)=='vcmaxhd':
        if ds.minmax.isel(ens=i)=='min':
            print(ds.key.sel(ens=i).values)
            print(np.round(npix.isel(ens=i).values[1:]/npix0,2))

# +
for ival in [0,363,416,427,434,441]:
    k = 'OAAT'+str(ival).zfill(4)
    f = '/glade/scratch/djk2120/PPEn11/paramfiles/'+k+'.nc'
    p = xr.open_dataset(f)['vcmaxhd']
    print(p.values)
    
    
    
    

# +
#instantiate the Ensemble object
basefile = '/glade/p/cesm/cseg/inputdata/lnd/clm2/paramdata/ctsm51_params.c210507.nc'
lndin    = '/glade/work/oleson/lmbirch_wkattge.n01_ctsm5.1.dev006/cime/scripts/clm51_lmbirchwkattgen01ctsm51d006_2deg_GSWP3V1_PPE_1850pAD/CaseDocs/lnd_in'
pdir     = '/glade/scratch/djk2120/PPEn11/paramfiles/'
ndir     = '/glade/scratch/djk2120/PPEn11/namelist_mods/'
pfiles = glob.glob(pdir+'*.nc')
keys = [int(p.split('/')[-1].split('.')[0][4:]) for p in pfiles]
surv2    = Ensemble(basefile,pdir,ndir)
i  = max(keys)




ds = dsets[0]
pdef = xr.open_dataset(pdir+'OAAT0000.nc')
wds = [0.5,0.75,0.875]
nxs = ['3','4','5']
kvals=[];pvals=[];mvals=[]
for wd,nx in zip(wds,nxs):
    w0 = 1-wd
    for r in resubmit:
        i     += 1
        newk   = 'OAAT'+str(i).zfill(4)

        r      = r.split('-')
        param  = r[0]
        minmax = r[1]
        ix     = np.logical_and(ds.param==param,ds.minmax==minmax)
        k      = str(ds.key[ix][1].values)
        xd     = pdef[param].values

        p0 = xr.open_dataset(pdir+k+'.nc')
        x0 = p0[param].values



        loc    = 'P'
        defval = get_default(param, loc, pdef, lndin)


        newx   = wd*xd+w0*x0
        value  = parse_val(loc,defval,newx)
        paraminfo = ParamInfo(param, loc, defval, value)
        paramdict = {param:paraminfo}
        member = Member(newk,paramdict,basefile,minmax)   
        surv2.add_member(member)

        kvals.append(newk)
        pvals.append(param)
        mvals.append(minmax+nx)
        print(newk,param,value)
# -

surv2.write()
newcsv=pd.DataFrame({'key':kvals,'param':pvals,'minmax':mvals})
newcsv.to_csv('/glade/scratch/djk2120/PPEn11/surv_resubmission2.csv',index=False)



# ### older stuff











# +
env    = 'C285'
thedir = '/glade/scratch/djk2120/PPEn11/hist/'+env+'/'
k      = 'OAAT0000'
f      = glob.glob(thedir+'PPEn*_'+env+'_'+k+'.clm2.h1.2005-02-01-00000.nc')[0]
tlai   = xr.open_dataset(f)['TLAI']

tlai['pft']  = dstmp['pfts1d_itype_veg']

laithresh = 0.1
npix0     = (tlai.max(dim='time')>laithresh).groupby('pft').sum().values

# +
k2  = 'OAAT0401'
ix2 = paramkey.key==k2
p   = paramkey.param[ix2].values[0]
m   = paramkey.minmax[ix2].values[0][:3]
ix1 = np.logical_and(paramkey.param==p,paramkey.minmax==m)
k1  = paramkey.key[ix1].values[0]

f            = glob.glob(thedir+'PPEn*_'+env+'_'+k1+'.clm2.h1.2005-02-01-00000.nc')[0]
tlai         = xr.open_dataset(f)['TLAI']
tlai['pft']  = dstmp['pfts1d_itype_veg']
npix1        = (tlai.max(dim='time')>laithresh).groupby('pft').sum().values

f            = glob.glob(thedir+'PPEn*_'+env+'_'+k2+'.clm2.h1.2005-02-01-00000.nc')[0]
tlai         = xr.open_dataset(f)['TLAI']
tlai['pft']  = dstmp['pfts1d_itype_veg']
npix2        = (tlai.max(dim='time')>laithresh).groupby('pft').sum().values
# -

resubmits=[]
for i in 400+np.arange(21):
    k2 = 'OAAT'+str(i).zfill(4)
    ix2 = paramkey.key==k2
    p   = paramkey.param[ix2].values[0]
    m   = paramkey.minmax[ix2].values[0][:3]
    ix1 = np.logical_and(paramkey.param==p,paramkey.minmax==m)
    k1  = paramkey.key[ix1].values[0]
    resubmits.append(k1)

resubmits



npix0

npix1

npix2







pp = xr.open_dataset('/glade/scratch/djk2120/PPEn11/paramfiles/OAAT0000.nc')
pftname = pp['pftname']

for i in 1+np.arange(16):
    ix = dstmp['pfts1d_itype_veg']==i
    print(i,pftname[i].values,ix.sum().values)

file    = '/glade/scratch/djk2120/PPEn11/hist/CTL2010/PPEn08_CTL2010_OAAT0000.clm2.h1.2005-02-01-00000.nc'
tlai = xr.open_dataset(file)['TLAI']
tlai['time'] = xr.cftime_range(str(2005),periods=120,freq='MS') #fix time bug
tlai['pft']  = dstmp['pfts1d_itype_veg']
laithresh    = 0.1
npix0        = (tlai.max(dim='time')>laithresh).groupby('pft').sum().compute()

alives = []
dsets  = []
for env in ['CTL2010','C285','AF1905']:
    print(env)
    ds        = get_ensemble(env,['TLAI','GPP'],'h1',paramkey)
    ds['pft'] = dstmp['pfts1d_itype_veg'] 
    npix0     = (ds['TLAI'].isel(ens=0).max(dim='time')>laithresh).groupby('pft').sum().compute()
    npixels   = (ds['TLAI'].max(dim='time')>laithresh).groupby('pft').sum().compute()
    surv      = 100*npixels/npix0.compute().isel(pft=1+np.arange(15))
    thresh    = 70
    alive     = surv.T>thresh
    alives.append(alive)
    dsets.append(ds)

alives[0]

fig,ax = plt.subplots(figsize=[12,3])
plt.pcolormesh(surv.T,vmin=0,vmax=100)
plt.colorbar()
plt.title('Survival % relative to default \n maxLAI>'+str(thresh)+' threshold')
plt.yticks([0,3,8,11,15])
ax.set_yticks([1.5,5.5,9.5,13], minor=True)
ax.yaxis.set_major_formatter(plt.NullFormatter())
ax.set_yticklabels(['NT','BT','S','G'],minor=True)
for tick in ax.yaxis.get_minor_ticks():
    tick.tick1line.set_markersize(0)
    tick.tick2line.set_markersize(0)
plt.ylabel('PFT')
plt.xlabel('Ensemble number');

# +
vals = np.ones([2,4])
vals[0,:3]=np.array([223,194,125])/256
vals[1,:3]=np.array([50,130,50])/256
newcmp = matplotlib.colors.ListedColormap(vals)


fig,ax = plt.subplots(figsize=[12,3])
plt.pcolormesh(alive,cmap=newcmp)

cbar=plt.colorbar(ticks=[0, 1])
cbar.ax.set_yticklabels(['dead', 'alive']) 
cbar.ax.tick_params(labelsize=14) 
plt.yticks([0,3,8,11,15])
ax.set_yticks([1.5,5.5,9.5,13], minor=True)
for y in [3,8,11]:
    plt.plot([0,len(ds.ens)],[y,y],':',color=[0.9,0.9,0.9],linewidth=0.5)

ax.yaxis.set_major_formatter(plt.NullFormatter())
ax.set_yticklabels(['NT','BT','S','G'],minor=True)
for tick in ax.yaxis.get_minor_ticks():
    tick.tick1line.set_markersize(0)
    tick.tick2line.set_markersize(0)
plt.title('C285: '+str(thresh)+'% threshold (relative to default)')
plt.ylabel('PFT')
plt.xlabel('Ensemble number');
# -

deads = dict()
for alive,ds in zip(alives,dsets):
    for i in alive.ens:
        if alive.isel(ens=i).sum()<15:
            p = ds.param[i].values
            m = ds.minmax[i].values
            k = str(ds.key[i].values)
            for f in alive.pft.values:
                if not alive.isel(ens=i).sel(pft=f):
                    if k in deads:
                        if f not in deads[k]:
                            deads[k].append(f)
                    else:
                        deads[k]=[f]


deads

for r in resubmits:
    if r in deads:
        print(r,'in deads')
    else:
        print(r,'not in')

for d in deads:
    if d in resubmits:
        print(d,'resubmitted')
    else:
        print(d,'not')

laithresh=0.1

ds = dsets[2]
npix0=(ds.TLAI.isel(ens=0).max(dim='time')>laithresh).groupby('pft').sum().compute()
for i in 400+np.arange(1):
    k   = 'OAAT'+str(i).zfill(4)
    ix  = ds.key==k
    p   = ds.param[ix].values[0]
    m   = ds.minmax[ix].values[0][:3]
    ix2 = np.logical_and(ds.param==p,ds.minmax==m)
    k2  = ds.key[ix2].values[0]
    
    npix1=(ds.TLAI.isel(ens=ix).max(dim='time')>laithresh).groupby('pft').sum().values.ravel()
    npix2=(ds.TLAI.isel(ens=ix2).max(dim='time')>laithresh).groupby('pft').sum().values.ravel()
    
    
    print(k2)

p

npix0

npix2

plt.plot(npix0,'-x')
plt.plot(npix1,'-x')
plt.plot(npix2,'-x')

deads

for k in deads:
    ix    = paramkey.key==k
    param = paramkey.param[ix].values[0]
    print(k,param)

p0 = xr.open_dataset('/glade/scratch/djk2120/PPEn11/paramfiles/OAAT0000.nc')



# +
#instantiate the Ensemble object
basefile = '/glade/p/cesm/cseg/inputdata/lnd/clm2/paramdata/ctsm51_params.c210507.nc'
lndin    = '/glade/work/oleson/lmbirch_wkattge.n01_ctsm5.1.dev006/cime/scripts/clm51_lmbirchwkattgen01ctsm51d006_2deg_GSWP3V1_PPE_1850pAD/CaseDocs/lnd_in'
pdir     = '/glade/scratch/djk2120/PPEn11/paramfiles/'
ndir     = '/glade/scratch/djk2120/PPEn11/namelist_mods/'
nextnum  = 399
x        = Ensemble(basefile,pdir,ndir)


for k in deads:
    nextnum +=1
    newk  = 'OAAT'+str(nextnum).zfill(4)
    
    p     = xr.open_dataset('/glade/scratch/djk2120/PPEn11/paramfiles/'+k+'.nc')
    ix    = paramkey.key==k
    param = paramkey.param[ix].values[0]
    minmax = paramkey.minmax[ix].values[0]+'2'
    
    if param != 'KCN':
        if param in p.data_vars:
            loc    = 'P'
            defval = get_default(param, loc, p0, lndin)    
            nx     = p[param].shape

            if not nx:
                newx   = 0.5*(p0[param]+p[param]).values
            else:
                newx = p[param].copy(deep=True).values
                for f in deads[k]:
                    f       = int(f)
                    newx[f] = 0.5*(p[param][f]+p0[param][f])

            value  = parse_val(loc,defval,newx)
        else:
            loc    = 'N'
            nfile  = '/glade/scratch/djk2120/PPEn11/namelist_mods/'+k+'.txt'
            cmd = 'grep '+param+' '+nfile
            tmp = os.popen(cmd).read()
            defval = get_default(param, loc, p0, lndin)
            x0     = float(tmp.split('=')[1].split('\n')[0])
            newx   = 0.5*(defval+x0) 
            value  = parse_val(loc,defval,newx)

        print(newk,param,loc,defval,value)
        paraminfo = ParamInfo(param, loc, defval, value)
        paramdict = {param:paraminfo}
        member = Member(newk,paramdict,basefile,minmax)   
        x.add_member(member)
# -

csvfile='/glade/scratch/djk2120/PPEn11/resubmission_071221.csv'
x.write(csvfile=csvfile)


newx

x0

cmd = 'grep '+param+' '+lndin

deads['a'].append(2)

deads['a']

alive.pft.values

resubmits = alive.isel(ens=ix)
    

# ### AF1905

ds = get_ensemble('AF1905',['TLAI','GPP'],'h1',paramkey)
#reassign pft variable for easier analysis
afile = '/glade/scratch/djk2120/PPEn11/hist/C285/PPEn11_C285_OAAT0398.clm2.h1.2005-02-01-00000.nc'
tmp = xr.open_dataset(afile)
ds['pft']=tmp['pfts1d_itype_veg'] 

# +
thresh  = 0.1
npixels = (ds['TLAI'].max(dim='time')>thresh).groupby('pft').sum().compute()
surv    = 100*npixels/npixels.isel(ens=0).compute().isel(pft=1+np.arange(15))

fig,ax = plt.subplots(figsize=[12,3])
plt.pcolormesh(surv.T,vmin=0,vmax=100)
plt.colorbar()
plt.title('Survival % relative to default \n maxLAI>'+str(thresh)+' threshold')
plt.yticks([0,3,8,11,15])
ax.set_yticks([1.5,5.5,9.5,13], minor=True)
ax.yaxis.set_major_formatter(plt.NullFormatter())
ax.set_yticklabels(['NT','BT','S','G'],minor=True)
for tick in ax.yaxis.get_minor_ticks():
    tick.tick1line.set_markersize(0)
    tick.tick2line.set_markersize(0)
plt.ylabel('PFT')
plt.xlabel('Ensemble number');

# +
vals = np.ones([2,4])
vals[0,:3]=np.array([223,194,125])/256
vals[1,:3]=np.array([50,130,50])/256
newcmp = matplotlib.colors.ListedColormap(vals)

thresh = 70
alive = surv.T>thresh
fig,ax = plt.subplots(figsize=[12,3])
plt.pcolormesh(alive,cmap=newcmp)

cbar=plt.colorbar(ticks=[0, 1])
cbar.ax.set_yticklabels(['dead', 'alive']) 
cbar.ax.tick_params(labelsize=14) 
plt.yticks([0,3,8,11,15])
ax.set_yticks([1.5,5.5,9.5,13], minor=True)
for y in [3,8,11]:
    plt.plot([0,len(ds.ens)],[y,y],':',color=[0.9,0.9,0.9],linewidth=0.5)

ax.yaxis.set_major_formatter(plt.NullFormatter())
ax.set_yticklabels(['NT','BT','S','G'],minor=True)
for tick in ax.yaxis.get_minor_ticks():
    tick.tick1line.set_markersize(0)
    tick.tick2line.set_markersize(0)
plt.title('AF1905: '+str(thresh)+'% threshold (relative to default)')
plt.ylabel('PFT')
plt.xlabel('Ensemble number');
# -

ix     = alive.sum(dim='pft')<15
ndeads = 15-alive.sum(dim='pft').values[ix]
for p,m,nd in zip(ds.param[ix].values,ds.minmax[ix].values,ndeads):
    print(p+'-'+m+': '+str(nd)+' PFTs dead')

# ### CTL2010

ds = get_ensemble('CTL2010',['TLAI','GPP'],'h1',paramkey)
#reassign pft variable for easier analysis
afile = '/glade/scratch/djk2120/PPEn11/hist/C285/PPEn11_C285_OAAT0398.clm2.h1.2005-02-01-00000.nc'
tmp = xr.open_dataset(afile)
ds['pft']=tmp['pfts1d_itype_veg'] 

# +
thresh  = 0.1
npixels = (ds['TLAI'].max(dim='time')>thresh).groupby('pft').sum().compute()
surv    = 100*npixels/npixels.isel(ens=0).compute().isel(pft=1+np.arange(15))

fig,ax = plt.subplots(figsize=[12,3])
plt.pcolormesh(surv.T,vmin=0,vmax=100)
plt.colorbar()
plt.title('Survival % relative to default \n maxLAI>'+str(thresh)+' threshold')
plt.yticks([0,3,8,11,15])
ax.set_yticks([1.5,5.5,9.5,13], minor=True)
ax.yaxis.set_major_formatter(plt.NullFormatter())
ax.set_yticklabels(['NT','BT','S','G'],minor=True)
for tick in ax.yaxis.get_minor_ticks():
    tick.tick1line.set_markersize(0)
    tick.tick2line.set_markersize(0)
plt.ylabel('PFT')
plt.xlabel('Ensemble number');

# +
vals = np.ones([2,4])
vals[0,:3]=np.array([223,194,125])/256
vals[1,:3]=np.array([50,130,50])/256
newcmp = matplotlib.colors.ListedColormap(vals)

thresh = 70
alive = surv.T>thresh
fig,ax = plt.subplots(figsize=[12,3])
plt.pcolormesh(alive,cmap=newcmp)

cbar=plt.colorbar(ticks=[0, 1])
cbar.ax.set_yticklabels(['dead', 'alive']) 
cbar.ax.tick_params(labelsize=14) 
plt.yticks([0,3,8,11,15])
ax.set_yticks([1.5,5.5,9.5,13], minor=True)
for y in [3,8,11]:
    plt.plot([0,len(ds.ens)],[y,y],':',color=[0.9,0.9,0.9],linewidth=0.5)

ax.yaxis.set_major_formatter(plt.NullFormatter())
ax.set_yticklabels(['NT','BT','S','G'],minor=True)
for tick in ax.yaxis.get_minor_ticks():
    tick.tick1line.set_markersize(0)
    tick.tick2line.set_markersize(0)
plt.title('CTL2010: '+str(thresh)+'% threshold (relative to default)')
plt.ylabel('PFT')
plt.xlabel('Ensemble number');
# -

ix     = alive.sum(dim='pft')<15
ndeads = 15-alive.sum(dim='pft').values[ix]
for p,m,nd in zip(ds.param[ix].values,ds.minmax[ix].values,ndeads):
    print(p+'-'+m+': '+str(nd)+' PFTs dead')



StarcoderdataJupyterScriptsDedupFiltered