#!/usr/bin/env python3

print('\n###############################################################\n')
print('     A probabilistic framework for Windows of Opportunity:\n')
print('         The role of variance in critical transitions.\n')
print('         Coded by Jim van Belzen\n')
print('################################################################\n')

# Version 1.0 (19/10/2019) - Python implementation
# Version 1.1 (23/10/2019)
# Version 1.2 (05/11/2019)
# Version 1.3 (02/01/2020)
# Version 1.4 (11/01/2020)
# Version 1.5 (09/02/2020)
# Version 1.6 (03/07/2020)
# Version 1.7 (23/07/2020)
# Version 1.8 (13/08/2020)
# Version 1.9 (19/08/2020)
# Version 2.0 (10/08/2021)
# Version 2.1 (19/09/2021)
# Version 2.2 (21/09/2021)
# Version 2.3 (16/10/2021)
# Version 3.1 (20/10/2021)
# Version 3.2 (14/01/2022) - final version
# Coded and developed by Jim van Belzen
# published under the Creative Commons Attribution Non-Commercial license
# which allows users to read, copy, distribute and make derivative works
# for noncommercial purposes from the material, as long as the author of
# the original work is cited.
# This code comes with no warranties or support
# http://people.zeelandnet.nl/jbelzen/



###=== SETUP PYTHON ENVIRONMENT ===###
import numpy as np
import numpy.matlib
import matplotlib.pyplot as plt
#from scipy import stats
from scipy import ndimage
import pandas as pd



###=== DEFINE FUNCTIONS ===###
def fRDS(E0, E1, k0, k1, Tend, dt):
    """
    Function to create Regular Dichotomous Signal (RDS)
    """
    N = int(Tend/dt + 1)
    t = np.linspace(0, N-1, N)*dt
    E = np.linspace(0, N-1, N)*0

    #- initial state
    if np.random.rand() <= k0/(k0+k1):
        q = 0;
    else:
        q = 1;

    #- create signal
    c = 0
    for i in range(N):
        # report current state
        if q == 0:
            E[i] = E0
        elif q == 1:
            E[i] = E1
        
        # determine new state
        if q == 0:
            if c >= 1./k0:
                q = 1
                c = 0
        elif q == 1:
            if c >= 1./k1:
                q = 0
                c = 0
        
        c += dt

    return E, t


def fDMN(E0, E1, k0, k1, Tend, dt):
    """
    Function to create Dichotomous Markov Noise (DMN) signal
    """
    #- temporal rescaling of transistion probabilities
    k0 = 1-(1-k0)**dt
    k1 = 1-(1-k1)**dt
    
    N = int(Tend/dt + 1)
    t = np.linspace(0, N-1, N)*dt
    E = np.linspace(0, N-1, N)*0
    
    #- initial state
    if np.random.rand() <= k0/(k0+k1):
        q = 0;
    else:
        q = 1;

    #- create signal
    for i in range(N):
        # report current state
        if q == 0:
            E[i] = E0
        elif q == 1:
            E[i] = E1
        
        # determine new state
        if q == 0:
            if np.random.rand() <= k0:
                q = 1
        elif q == 1:
            if np.random.rand() <= k1:
               q = 0

    return E, t


def fOUP(mu, Theta, sigma, Tend, dt):
    """
    Function to create correlated noise signal based on Ornstein-Uhlenbeck process (OUP)
    mu      average of stationary signal
    Theta   rate at which signal restores to average
    sigma   magnitude of Gaussian noise
    """
    #- setting up output variable
    N = int(Tend/dt + 1)
    t = np.linspace(0, N-1, N)*dt
    E = np.linspace(0, N-1, N)*0

    #- initial state (based on 6sigma rule for 99.9%)
    xt = mu + 3*sigma*np.random.normal()

    #- create signal
    for i in range(N):
        E[i] = xt
        dxt = Theta*(mu - xt)*dt + sigma*np.random.normal()
        xt += dxt

    return E, t
    
    
def fHRN(mu, tau, sigma, Tend, dt):
    """
    Function to create correlated noise signal based on Hasselmann (1976) variation of the
    Ornstein-Uhlenbeck process (OUP)
    mu      average of stationary signal
    tau     approximate period of the signal
    sigma   magnitude of Gaussian noise
    """
    
    #- setting up output variable
    N = int(Tend/dt + 1)
    t = np.linspace(0, N-1, N)*dt
    E = np.linspace(0, N-1, N)*0

    #- initial state (based on 6sigma rule for 99.9%)
    xt = mu + 3*sigma*np.random.normal()

    #- create signal
    for i in range(N):
        E[i] = xt
        #dxt = (1 - (1/tau))*(mu - xt)*dt + sigma*np.random.normal()
        #dxt = (1 - (1/l))*(mu - xt)*dt + sigma*np.random.normal()
        #dxt = (1 - (1/l))*(xt - mu)*dt + sigma*np.random.normal()
        #dxt = (1 - (tau))*(xt - mu)*dt + sigma*np.random.normal()
        dxt = (1/tau)*(mu - xt)*dt + sigma*np.random.normal()
        xt += dxt

    return E, t


def fpotwin(E, T0, dt):
    """
    Analyze environmental signal for potential windows
    """
    lbl, nlbl = ndimage.label(E<=T0)
    w = np.zeros((nlbl+1,1))
    wlbl = range(1, nlbl+1)
    for i in wlbl:
        w[i] = np.sum(lbl == i)*dt

    return w


def fnumwoo(E, T0, dt, lamb):
    """
    Numerically derive Windows of Opportunity from any environmental signal
    """
    if T0>=np.min(E) and lamb>0:
        num = np.array((np.max(E) - T0)/(lamb*dt), dtype=int)
        
        #- run growth
        T = (E<=T0)*T0
        B = (T>0)
        for i in range(num):
            ii = i+1
            T += lamb*B*dt
            T[:-ii] = np.where(T[:-ii] > E[ii:], T[:-ii], 0)
            B = (T>0)
        
        #- get Windows of Opportunity
        lbl, nlbl = ndimage.label(B==1)
        w = np.zeros((nlbl+1,1))
        wlbl = range(1, nlbl+1)
        for i in wlbl:
            w[i] = np.sum(lbl == i)*dt
    else:
        w = 0
    
    return w



def fTt(T, lamb, t, n_model):
    """
    Tolerance after time t
    """
    if n_model == 0:
        """ linear model """
        T = T + lamb*t
    elif n_model == 1:
        """ exponential model """
        T = T*np.exp(lamb*t)

    return T



###=== FIGURE 1 ===###
#- set seed
np.random.seed(7)

# Figure panel properties
nrows   = 1
ncols   = 1
dpi     = 300
tmp = plt.rcParams['figure.figsize']

plt.rcParams['figure.figsize'] = [tmp[0] * 0.8, tmp[1] * .6]

##-- Build panel Figure 1
fig, ax1 = plt.subplots(nrows, ncols)

##-- Figure 1.a: Conceptual model
##- Parameters & settings
#- Environmental condition, E
E0      = 0.        # Emin
E1      = 1         # Emax
tau0    = 10.       # waiting time in Emin
tau1    = 2.5       # waiting time in Emax
k0      = 1./tau0   # transition rate from Emin to Emax
k1      = 1./tau1   # transition rate from Emin to Emax

tau     = tau0 + tau1   # total waiting time of the environmental condition
Eavg    = (tau0*E0 + tau1*E1)/tau

#- Tolerance of young organism, T
T0      = E0 + (Eavg-E0)*.5   # initial tolerance (for illustration it is in between Emin and Eavg)
delta   = E1 - T0

#- Temporal settings
Tend    = 300
dt      = .01

##-- Data/results to show
# build regular dichotomous sigal
[E, t] = fRDS(E0, E1, k0, k1, Tend, dt)

# arrow setting
hl = tau/25     # head_length
hw = (E1-E0)/20  # head_width

##- Building up panel Figure 1.a
# set-up environmental signal
ax1.plot(t, E, color='gray', lw=2)
ax1.plot([0, Tend],[Eavg, Eavg],'--', color='gray', lw=1.5)
ax1.plot(t[int(tau/dt):int(2*tau/dt)], E[int(tau/dt):int(2*tau/dt)], '-k', lw=2)

# add establishing organisms
ax1.plot(tau+tau0, T0+delta*.53, 'or', markersize=4)
ax1.arrow(tau, T0, tau0 - hl, delta*.5, head_width=hw, head_length=hl, lw=2, fc='g', ec='g')
ax1.text(tau+tau0, T0+delta*.51, r'${X}$', color='red', fontsize=16, \
    fontweight='bold', horizontalalignment = 'center', \
    verticalalignment = 'center')
ax1.arrow(tau, T0, tau, delta/tau0*tau, head_width=hw, head_length=hl, lw=2, fc='g', ec='g')
ax1.annotate("", xy=(tau+tau0, E0-0.1*(E1-E0)), xytext=(tau, E0-0.1*(E1-E0)), \
    arrowprops=dict(arrowstyle="<->"))
ax1.annotate("", xy=(2*tau, E0-0.1*(E1-E0)), xytext=(tau+tau0, E0-0.1*(E1-E0)), \
    arrowprops=dict(arrowstyle="<->"))
ax1.annotate("", xy=(tau+tau0+tau1/6, E1), xytext=(tau+tau0+tau1/6, T0), \
    arrowprops=dict(arrowstyle="<->"))
ax1.plot(tau, T0, 'og')

# add labels
ax1.text(tau*.6 + .025*tau, E0+.05*(E1-E0), r'$E_{low}$')
ax1.text(tau*.6 + .025*tau, Eavg+.05*(E1-E0), r'$E_{avg}$')
ax1.text(tau*.6 + .025*tau, E1+.05*(E1-E0), r'$E_{high}$')
ax1.text(tau+tau0/2, E0 + delta*.9, r'$\lambda_{crit}$', horizontalalignment = 'center')
ax1.text(tau+tau0/2, E0 + delta*.75, '\u2220', horizontalalignment = 'center')
ax1.text(tau+tau0/2, E0 - (E1-E0)*.25, r'$\tau_{low}$', horizontalalignment = 'center')
ax1.text(tau+tau0+tau1/2, E0 - (E1-E0)*.25, r'$\tau_{high}$', horizontalalignment = 'center')
ax1.text(tau+tau0+tau1/4, E0 + delta*.5, r'$\delta$')
ax1.text(tau-.09*tau, T0-.05*(E1-E0), r'T$_0$')

# figure axis
ax1.set_xlabel(r'Time, $t$')
ax1.set_ylabel(r'Environment, $E$', labelpad=20)
ax1.text(tau*.6 - 1.5, (E0-.4*(E1-E0) + 1.5*E1)/2, r'Tolerance, $T$', verticalalignment = 'center', rotation=90)
ax1.set_xlim((tau*.6,2.2*tau))
ax1.set_ylim((E0-.4*(E1-E0),1.5*E1))
# turn off ticks & labels
ax1.set_yticklabels([])
ax1.set_xticklabels([])
ax1.xaxis.set_ticks_position('none')
ax1.yaxis.set_ticks_position('none')


##- Finish & saving figure
plt.subplots_adjust(wspace = -0.5, hspace = .4)
plt.savefig('Figure1.png', dpi = dpi, bbox_inches = 'tight')
plt.show()
plt.close()
print('---### Figure 1 created and saved ###---')



###=== FIGURE 2 ===###
# Figure panel properties
nrows   = 2
ncols   = 1
dpi     = 300
plt.rcParams['figure.figsize'] = [tmp[0] * 0.8, tmp[1] * 1.2]

##-- Build panel Figure 2
fig, (ax1, ax2) = plt.subplots(nrows, ncols)

##-- Figure 2.a: Bifurcation diagram
##- Parameters & settings (additional)
#- Environmental condition, E
E0      = 1.2       # Emin
E1      = 1.2       # Emax
tau0    = 10        # waiting time in Emin
tau1    = 2.5       # waiting time in Emax
k0      = 1./tau0   # transition rate from Emin to Emax
k1      = 1./tau1   # transition rate from Emax to Emin

tau     = tau0 + tau1   # total waiting time of the environmental condition
Eavg    = (tau0*E0 + tau1*E1)/tau

#- parameters phase-plain
nb_bins     = 20
nb_grpnts   = 500
nb_arrows   = 7

n_model     = 0     # choose model

levels1 = np.arange(-100, 100, 100)

#- initialize
e0 = 1
e1 = 2
# lambda along x-axis
xmin = -0.01
xmax = 0.16
# T along x-axis
ymin = e0 - .25*(e1-e0)
ymax = e1 + .25*(e1-e0)

xplt = np.array([xmin, 0])
yplt = np.array([Eavg, Eavg])
xplt2 = np.array([xmin, xmax])
yplt2 = np.array([T0, T0])

##- Phase-plain
#- contour
x = np.linspace(xmin, xmax, nb_grpnts)
y = np.linspace(ymin, ymax, nb_grpnts)

X, Y = np.meshgrid(x, y)    # create x- and y-grids

T = fTt(Y, X, tau0, n_model = 0) * (Y>E0)
T[T<=E1] = 0
DT = T - Y

levels2 = np.linspace(np.min(DT), np.max(DT), nb_bins)

#- quiver
x2 = np.linspace(xmin-0.025, xmax+0.025, nb_arrows)
y2 = np.linspace(ymin-0.025, ymax+0.4, nb_arrows)

X2, Y2 = np.meshgrid(x2, y2) # create x- and y-grids

T2 = fTt(Y2, X2, tau0, n_model = 0) * (Y2>E0)
T2[T2<=E1] = 0
DT2 = T2 - Y2

ax1.plot(xplt, yplt, linestyle = '--', color='gray', lw=1.05)
ax1.quiver(X2, Y2, DT2*0, DT2*(DT2>0), scale_units='xy', scale=6, \
    pivot = 'tail', cmap = plt.cm.jet, color='green')
ax1.quiver(X2, Y2, DT2*0, DT2*(DT2<0), scale_units='xy', scale=6, \
    pivot = 'tail', cmap = plt.cm.jet, color='black')
ax1.contour(X, Y, DT, levels1, colors = 'black', linewidths = 2.5)

# figure axis
ax1.set_xlabel(r'Growth rate, $\lambda$', labelpad=-5)
ax1.set_ylabel(r'Initial tolerance, $T_{0}$')
ax1.set_xlim((xmin,xmax))
ax1.set_ylim((ymin,ymax))
ax1.set_xticks([0])
ax1.set_xticklabels(['0'])
ax1.set_yticks([Eavg])
ax1.set_yticklabels([r'$E_{avg}$'])

# add labels & arrow
ax1.text(xmin+(xmax-xmin)*.01, ymax+(ymax-ymin)*.03, \
    r'No variance, $\sigma^2(E)=0$')

# add explaining text
ax1.text(xmin+(xmax-xmin)*.32, ymin + (ymax-ymin)*0.02, \
    'No establishment')
ax1.text(xmin+(xmax-xmin)*.32, ymin + (ymax-ymin)*0.9, \
    'Establishment')
    
# panel label
ax1.text(xmin-(xmax-xmin)*.175, ymax+(ymax-ymin)*.01, '(a)', \
    fontsize=12, fontweight='bold')



##-- Figure 2.b: Bifurcation diagram
##- Parameters & settings (additional)
#- Environmental condition, E
E0      = 1         # Emin
E1      = 2         # Emax
tau0    = 10        # waiting time in Emin
tau1    = 5         # waiting time in Emax
k0      = 1./tau0   # transition rate from Emin to Emax
k1      = 1./tau1   # transition rate from Emin to Emax

tau     = tau0 + tau1   # total waiting time of the environmental condition
Eavg    = (tau0*E0 + tau1*E1)/tau

#- parameters phase-plain
nb_bins     = 20
nb_grpnts   = 500
nb_arrows   = 7

n_model     = 0     # choose model

levels1 = np.arange(-100, 100, 100)

#- initialize
# lambda along x-axis
xmin = -0.01
xmax = 0.16
# T along x-axis
ymin = E0 - .25*(E1-E0)
ymax = E1 + .25*(E1-E0)

xplt = np.array([xmin, xmax])
yplt = np.array([Eavg, Eavg])
xplt2 = np.array([xmin, xmax])
yplt2 = np.array([E0, E0])
xplt3 = np.array([xmin, xmax])
yplt3 = np.array([E1, E1])

##- Phase-plain
#- contour
x = np.linspace(xmin, xmax, nb_grpnts)
y = np.linspace(ymin, ymax, nb_grpnts)

X, Y = np.meshgrid(x, y)    # create x- and y-grids

T = fTt(Y, X, tau0, n_model = 0) * (Y>E0)
T[T<=E1] = 0
DT = T - Y

levels2 = np.linspace(np.min(DT), np.max(DT), nb_bins)

#- quiver
x2 = np.linspace(xmin-0.025, xmax+0.025, nb_arrows)
y2 = np.linspace(ymin-0.025, ymax+0.4, nb_arrows)

X2, Y2 = np.meshgrid(x2, y2) # create x- and y-grids

T2 = fTt(Y2, X2, tau0, n_model = 0) * (Y2>E0)
T2[T2<=E1] = 0
DT2 = T2 - Y2

ax2.plot(xplt, yplt, linestyle = '--', color='gray', lw=1.05)
ax2.plot(xplt2, yplt2, linestyle = '--', color='gray', lw=1.05)
ax2.plot(xplt3, yplt3, linestyle = '--', color='gray', lw=1.05)
ax2.quiver(X2, Y2, DT2*0, DT2*(DT2>0), scale_units='xy', scale=6, \
    pivot = 'tail', cmap = plt.cm.jet, color='green')
ax2.quiver(X2, Y2, DT2*0, DT2*(DT2<0), scale_units='xy', scale=6, \
    pivot = 'tail', cmap = plt.cm.jet, color='black')
ax2.contour(X, Y, DT, levels1, colors = 'black', linewidths = 2.5)


# figure axis
ax2.set_xlabel(r'Growth rate, $\lambda$')
ax2.set_ylabel(r'Initial tolerance, $T_{0}$', labelpad=-.1)
ax2.set_xlim((xmin,xmax))
ax2.set_ylim((ymin,ymax))
# turn off ticks & labels
ax2.set_xticks([0])
ax2.set_xticklabels(['0'])
ax2.set_yticks([E0, Eavg, E1])
ax2.set_yticklabels([r'$E_{low}$', r'$E_{avg}$', r'$E_{high}$'])

# add labels
ax2.text(xmin+(xmax-xmin)*.01, ymax+(ymax-ymin)*.03, \
    r'With variance, $\sigma^2(E)>0$')
    
# add explaining text
ax2.text(xmin+(xmax-xmin)*.32, ymin + (ymax-ymin)*0.02, \
    'No establishment')
ax2.text(xmin+(xmax-xmin)*.32, ymin + (ymax-ymin)*0.9, \
    'Establishment')
    
# Panel labels
ax2.text(xmin-(xmax-xmin)*.175, ymax+(ymax-ymin)*.01, '(b)', \
    fontsize=12, fontweight='bold')


##- Finish & saving figure
plt.subplots_adjust(wspace = -0.5, hspace = .4)
plt.savefig('Figure2.png', dpi = dpi, bbox_inches = 'tight')
plt.show()
plt.close()
print('---### Figure 2 created and saved ###---')



###=== FIGURE S1 ===###
#- set seed
np.random.seed(7)

# Figure panel properties
nrows   = 3
ncols   = 1
dpi     = 300
plt.rcParams['figure.figsize'] = [tmp[0] * .8, tmp[1] * 1.8]

##-- Build panel Figure 1
fig, (ax1, ax2, ax3) = plt.subplots(nrows, ncols)

##-- Figure 1.a: Conceptual model
##- Parameters & settings
#- Environmental condition, E
E0      = 0.        # E-
E1      = 1         # E+
tau0    = 10.       # waiting time in E-
tau1    = 2.5       # waiting time in E+
k0      = 1./tau0   # transition rate from E0 to E1
k1      = 1./tau1   # transition rate from E1 to E0

tau     = tau0 + tau1   # total waiting time of the environmental condition
Eavg    = (tau0*E0 + tau1*E1)/tau

#- Tolerance of young organism, T
T0      = E0 + (Eavg-E0)*.5   # initial tolerance (for illustration it is in between E0 and Eavg)
delta   = E1 - T0
print('T0 = ', T0, 'delta = ', delta)

#- Temporal settings
Tend    = 300
dt      = .01

##-- Data/results to show
# build regular dichotomous sigal
[E, t] = fRDS(E0, E1, k0, k1, Tend, dt)

# exponentially developing young organism
tt = np.linspace(0, int((tau-2)/dt + 1)-1, int((tau-2)/dt + 1))*dt
lcrit = (np.log(E1) - np.log(T0))/tau0
Tt = T0*np.exp(lcrit*tt)
print(tt)

# arrow setting
hl = tau/25      # head_length
hw = (E1-E0)/20  # head_width

##- Building up panel Figure S1.a
# set-up environmental signal
ax1.plot(t, E, color='gray', lw=2)
ax1.plot([0, Tend],[Eavg, Eavg],'--', color='gray', lw=1.5)
ax1.plot(t[int(tau/dt):int(2*tau/dt)], E[int(tau/dt):int(2*tau/dt)], '-k', lw=2)

# add establishing organisms
ax1.arrow(tau, T0, tau, delta/tau0*tau, head_width=hw, head_length=hl, lw=2, fc='g', ec='g')
ax1.annotate("", xy=(tau+tau0, E0-0.1*(E1-E0)), xytext=(tau, E0-0.1*(E1-E0)), \
    arrowprops=dict(arrowstyle="<->"))
ax1.annotate("", xy=(2*tau, E0-0.1*(E1-E0)), xytext=(tau+tau0, E0-0.1*(E1-E0)), \
    arrowprops=dict(arrowstyle="<->"))
ax1.annotate("", xy=(tau+tau0+tau1/6, E1), xytext=(tau+tau0+tau1/6, T0), \
    arrowprops=dict(arrowstyle="<->"))
ax1.plot(tt+tau,Tt, lw=2, color='darkblue')
ax1.arrow(tt[-2]+tau, Tt[-2], dt, Tt[-1]-Tt[-2], head_width=hw*1.75, head_length=hl, \
    lw=2, fc='darkblue', ec='darkblue')
ax1.plot(tau, T0, 'og')

# add labels
ax1.text(tau*.6 + .025*tau, E0+.05*(E1-E0), r'$E_{low}}$')
ax1.text(tau*.6 + .025*tau, Eavg+.05*(E1-E0), r'$E_{avg}$')
ax1.text(tau*.6 + .025*tau, E1+.05*(E1-E0), r'$E_{high}$')
ax1.text(tau+tau0/2, E0 + delta*.9, r'$\lambda_{crit}$', horizontalalignment = 'center')
ax1.text(tau+tau0/2, E0 - (E1-E0)*.25, r'$\tau_{low}$', horizontalalignment = 'center')
ax1.text(tau+tau0+tau1/2, E0 - (E1-E0)*.25, r'$\tau_{high}$', horizontalalignment = 'center')
ax1.text(tau+tau0+tau1/4, E0 + delta*.5, r'$\delta$')
ax1.text(tau-.08*tau, T0-.05*(E1-E0), r'$T_0$')

# figure axis
ax1.set_xlabel(r'Time, $t$')
ax1.set_ylabel(r'Environment, $E$', labelpad=20)
ax1.text(tau*.6 - 1.5, (E0-.4*(E1-E0) + 1.5*E1)/2, r'Tolerance, $T$', verticalalignment = 'center', rotation=90)
ax1.set_xlim((tau*.6,2.2*tau))
ax1.set_ylim((E0-.4*(E1-E0),1.5*E1))
# turn off ticks & labels
ax1.set_yticklabels([])
ax1.set_xticklabels([])
ax1.xaxis.set_ticks_position('none')
ax1.yaxis.set_ticks_position('none')


##-- Figure S1.b: Role of variance on growth rate
##- Parameters & settings
scalemin  = 0
scalemax  = 30
rmin      = 0
rmax      = .5
N         = 500

scale  = np.linspace(scalemin, scalemax, N+1)
Env0   = Eavg - scale*(Eavg-E0)*tau0/tau
Env1   = Eavg + scale*(E1 - Eavg)*tau1/tau
Sigma0 = np.sqrt(tau0/tau*(Env0-Eavg)**2 + tau1/tau*(Env1-Eavg)**2)
Delta0 = Env1 - T0

Lcrit0 = Delta0/tau0
Lcrit0[Env0>T0] = float('inf')

Lcrit1 = (np.log(Env1) - np.log(T0))/tau0
Lcrit1[Env0>T0] = float('inf')

Lcrit2 = (np.log(Env1 - T0/2) - np.log(T0 - T0/2))/tau0
Lcrit2[Env0>T0] = float('inf')

Lcrit3 = (np.log(Env1 + T0/2) - np.log(T0 + T0/2))/tau0
Lcrit3[Env0>T0] = float('inf')

##- Building up panel Figure S1.c
ax2.plot(Sigma0, Lcrit0, '-g', label=r'lin. T$_0$')
ax2.plot(Sigma0, Lcrit1, '-', color='darkblue', label=r'exp. $T_0$')
ax2.plot(Sigma0, Lcrit2, ':', color='darkblue', label=r'exp. $T_0$ - 50%')
ax2.plot(Sigma0, Lcrit3, '--', color='darkblue', label=r'exp. $T_0$ + 50%')
ax2.legend()

# add explaining text
ax2.text(np.max(Sigma0)/2 - 0.025*np.max(Sigma0)/2, rmin + 0.025*(rmax-rmin), 'No establishment', \
    horizontalalignment = 'right', verticalalignment = 'bottom')
ax2.text(0.05*np.max(Sigma0)/2, rmax - 0.05*(rmax-rmin), 'Establishment', \
    horizontalalignment = 'left', verticalalignment = 'top')

# No establishment at low variance (gray shaded)
ax2.fill_between([0,np.max(Sigma0[Env0>T0])], [0,0], [rmax,rmax], color='gray')

# figure axis
ax2.set_xlabel(r'Variance, $\sigma^2$(E)')
ax2.set_ylabel(r'Growth rate, $\lambda_{crit}$')
ax2.set_xlim((0, np.max(Sigma0)))
ax2.set_ylim((rmin, rmax))


##-- Figure 1.c: Role of period on growth rate
##- Parameters & settings
taumin  = 0
taumax  = 20
rmin    = 0
rmax    = .5

Tau0 = np.linspace(taumin, taumax, 501)

##- Building up panel Figure S1.c
# linear growth
ax3.plot(Tau0, delta/2/Tau0, ':g', lw=2, label=r'lin. $\delta$/2')
ax3.plot(Tau0, delta/Tau0, '--g', lw=2, label=r'lin. $\delta$')
ax3.plot(Tau0, 2*delta/Tau0, '-g', lw=2, label=r'lin. 2$\delta$')
# exponential growth
ax3.plot(Tau0, (np.log(E1) - np.log(T0))/Tau0, '-', lw=2, color='darkblue', label=r'exp. $\delta$')
ax3.plot(Tau0, (np.log(T0+.5*delta) - np.log(T0))/Tau0, ':', lw=2, color='darkblue', label=r'exp. $\delta$/2')
ax3.plot(Tau0, (np.log(T0+2*delta) - np.log(T0))/Tau0, '--', lw=2, color='darkblue', label=r'exp. 2$\delta$')
ax3.legend(loc='upper right', ncol=1, labelspacing=0.15)

# add explaining text
ax3.text(taumin + 0.025*(taumax-taumin), rmin + 0.025*(rmax-rmin), 'No establishment', \
    horizontalalignment = 'left', verticalalignment = 'bottom')
ax3.text(taumax - 0.5*(taumax-taumin), rmax - 0.035*(rmax-rmin), 'Establishment', \
    horizontalalignment = 'center', verticalalignment = 'top')
    
# figure axis
ax3.set_xlabel(r'Period, $\tau_{low}$')
ax3.set_ylabel(r'Growth rate, $\lambda_{crit}$')
ax3.set_xlim((taumin, taumax))
ax3.set_ylim((rmin, rmax))


# panel labels
ax1.text(0.6*tau-0.17*(2.2-.6)*tau, 1.5*E1, '(a)', fontsize=12, fontweight='bold')
ax2.text(-0.17*np.max(Sigma0), rmax, '(b)', fontsize=12, fontweight='bold')
ax3.text(-0.17*taumax, rmax, '(c)', fontsize=12, fontweight='bold')

# finish figure
plt.subplots_adjust(wspace = -0.5, hspace = .4)
plt.savefig('FigureS1.png', dpi = dpi, bbox_inches = 'tight')
plt.show()
plt.close()
print('---### Figure S1 created and saved ###---')


###=== FIGURE 3 ===###
#- set seed
np.random.seed(0)

# Figure panel properties
nrows   = 3
ncols   = 1
dpi     = 300
plt.rcParams['figure.figsize'] = [tmp[0] * .8, tmp[1] * 1.8]

##-- Build panel Figure 3
fig, (ax1, ax2, ax3) = plt.subplots(nrows, ncols)

##-- Figure 3.a: Dichotmous Markov Noise model
##- Parameters & settings
#- Environmental condition, E
E0      = 0.        # Emin
E1      = 1.        # Emax
tau0    = 10.       # waiting time in Emin
tau1    = 5.        # waiting time in Emax
k0      = 1./tau0   # transition rate from Emin to Emax
k1      = 1./tau1   # transition rate from Emin to Emax

tau     = tau0 + tau1   # total waiting time of the environmental condition
Eavg    = (tau0*E0 + tau1*E1)/tau

#- Tolerance of young organism, T
T0      = E0 + (Eavg-E0)*.5   # initial tolerance (for illustration it is in between Emin and Eavg)
delta   = E1 - T0

p0      = .9

lambmin = 0
lambmax = .5
nlamb   = 500
lambcrit= 0.1

#- Temporal settings
Tend    = 100
dt      = .01

# arrow setting
hl = tau/25     # head_length
hw = (E1-E0)/20  # head_width

##-- Data/results to show
# build Dichotomous Markov Noise sigal (Fig.3.a)
[E, t] = fDMN(E0, E1, k0, k1, Tend, dt)

# build Dichotomous Markov Noise sigal (Fig.3.b)
np.random.seed(100) #- reset seed
[Estats, tstats] = fDMN(E0, E1, k0, k1, Tend*100, dt)
w = fpotwin(Estats, T0, dt)

print ('\nStatistics of DMN displayed in Fig.3.b:')
print ('Mean lenght of time window =', np.mean(w))
print ('Minimum lenght of time window =', np.min(w))
print ('Maximum lenght of time window =', np.max(w))

taucrit0a = 1.*delta/lambcrit
taucrit1a = 2.*delta/lambcrit
taucrit2a = 4.*delta/lambcrit

##- Building up panel Figure 3.a
# set-up environmental signal
ax1.plot(t, E, color='black', lw=2)
ax1.plot([0, Tend],[Eavg, Eavg],'--', color='black', lw=1.5)

# add labels
ax1.text(-0.09*Tend, E0+.01*(E1-E0), r'$E_{low}$')
ax1.text(-0.09*Tend, Eavg+.01*(E1-E0), r'$E_{avg}$')
ax1.text(-0.09*Tend, E1+.01*(E1-E0), r'$E_{high}$')
ax1.text((23.2+39)/2, E0 - (E1-E0)*.25, r'$\tau_{low}$', horizontalalignment = 'center')
ax1.text((38.2+46.2)/2, E0 - (E1-E0)*.25, r'$\tau_{high}$', horizontalalignment = 'center')
ax1.annotate("", xy=(39, E0+0.1*(E1-E0)), xytext=(23.2, E0+0.1*(E1-E0)), \
    arrowprops=dict(arrowstyle="<->"))
ax1.annotate("", xy=(46.2, E0+0.1*(E1-E0)), xytext=(38.2, E0+0.1*(E1-E0)), \
    arrowprops=dict(arrowstyle="<->"))
ax1.text(49.7, (E1-E0)/2, r'$k_{low}$', horizontalalignment = 'center')
ax1.text(64, (E1-E0)/2, r'$k_{high}$', horizontalalignment = 'center')
ax1.annotate("", xy=(53.5, E1), xytext=(53.5, E0), \
    arrowprops=dict(arrowstyle="->"))
ax1.annotate("", xy=(60, E0), xytext=(60, E1), \
    arrowprops=dict(arrowstyle="->"))
    
# figure axis
ax1.set_xlabel(r'Time, $t$')
ax1.set_ylabel(r'Environment, $E$')
ax1.set_xlim((-0.1*Tend,0.8*Tend))
ax1.set_ylim((E0-.4*(E1-E0),E1+.4*(E1-E0)))
# turn off ticks & labels
ax1.set_yticklabels([])
ax1.yaxis.set_ticks_position('none')


##- Building up panel Figure 3.b
# the histogram of the data
bins = range(int(tau0*10))
n, bins, patches = ax2.hist(w, bins=bins, facecolor='g', alpha=0.75)

# instantiate a second axes that shares the same x-axis
ax2b = ax2.twinx()
ax2b.plot(bins, k0*np.exp(-k0*bins), '--k')


# figure axis
ax2.set_xlabel(r'Time window, $\tau_{low}$')
ax2.set_ylabel(r'Number of windows, $N(\tau_{low})$')
ax2.set_xlim((0,tau0*6.5))

ax2b.set_ylabel('Probability density')
ax2b.set_ylim((0,1.15*k0))

ax2b.annotate(r'$\delta$', xy=(taucrit0a, k0*np.exp(-k0*taucrit0a)), xytext=(taucrit0a, \
    k0*np.exp(-k0*taucrit0a)+k0/4), arrowprops=dict(arrowstyle="->", color='red'))
ax2b.plot([taucrit0a, taucrit0a], [0, k0*np.exp(-k0*taucrit0a)], '--r')
ax2b.annotate(r'2$\delta$', xy=(taucrit1a, k0*np.exp(-k0*taucrit1a)), xytext=(taucrit1a, \
    k0*np.exp(-k0*taucrit1a)+k0/4), arrowprops=dict(arrowstyle="->", color='red'))
ax2b.plot([taucrit1a, taucrit1a], [0, k0*np.exp(-k0*taucrit1a)], '--r')
ax2b.annotate(r'4$\delta$', xy=(taucrit2a, k0*np.exp(-k0*taucrit2a)), xytext=(taucrit2a, \
    k0*np.exp(-k0*taucrit2a)+k0/4), arrowprops=dict(arrowstyle="->", color='red'))
ax2b.plot([taucrit2a, taucrit2a], [0, k0*np.exp(-k0*taucrit2a)], '--r')


##- Estimate establishment probailities, Pest (Fig.3.c)
lamb = np.linspace(lambmin, lambmax, nlamb)
taucrit0b = 1.*delta/lamb
taucrit1b = 2.*delta/lamb
taucrit2b = 4.*delta/lamb

P_WoO0 = np.exp(-k0*taucrit0b)
P_WoO1 = np.exp(-k0*taucrit1b)
P_WoO2 = np.exp(-k0*taucrit2b)

xplt = np.array([lambmin, lambmax])
yplt = np.array([p0, p0])

##- Building up panel Figure 3.c
# set-up environmental signal
ax3.plot(lamb, P_WoO0, '-g', lw=2, label=r'$\delta$')
ax3.plot(lamb, P_WoO1, '--g', lw=2, label=r'2$\delta$')
ax3.plot(lamb, P_WoO2, ':g', lw=2, label=r'4$\delta$')

ax3.plot([lambcrit, lambcrit],[0,1], '--', color='red')
ax3.text(lambcrit-0.1*(lambmax-lambmin), 0.7, r'$\lambda_{crit}$')

# figure axis
ax3.set_xlabel(r'Growth rate, $\lambda$')
ax3.set_ylabel(r'Probability of WoO, $P_{WoO}$')
ax3.set_xlim((lambmin,lambmax))
ax3.set_ylim((0,1))
ax3.legend(loc='lower right', ncol=1, labelspacing=0.15)


# panel labels
ax1.text(-0.1*Tend-0.21*0.9*Tend, 1.5*E1, '(a)', fontsize=12, fontweight='bold')
ax2.text(-0.21*tau0*6.5, 89, '(b)', fontsize=12, fontweight='bold')
ax3.text(-0.21*(lambmax-lambmin), 1, '(c)', fontsize=12, fontweight='bold')

# finish figure
plt.subplots_adjust(wspace = -0.5, hspace = .4)
plt.savefig('Figure3.png', dpi = dpi, bbox_inches = 'tight')
plt.show()
plt.close()
print('---### Figure 3 created and saved ###---')


###=== FIGURE 4 ===###
# Figure panel properties
nrows   = 1
ncols   = 1
dpi     = 300
plt.rcParams['figure.figsize'] = [tmp[0] * .8, tmp[1] * .6]

##-- Build panel Figure 4
fig, ax1 = plt.subplots(nrows, ncols)

##-- Figure 4: Probability diagram
##- Parameters & settings (additional)
#- Environmental condition, E
E0      = 1         # E-
E1      = 2         # E+
tau0    = 10        # waiting time in E-
tau1    = 5         # waiting time in E+

k0      = 1./tau0   # transition rate from E0 to E1
k1      = 1./tau1   # transition rate from E1 to E0

tau     = tau0 + tau1   # total waiting time of the environmental condition
Eavg    = (tau0*E0 + tau1*E1)/tau

#- Tolerance of young organism, T
T0      = E0 + (Eavg-E0)*.5   # initial tolerance (for illustration it is in between E0 and Eavg)
delta   = E1 - T0
lcrit   = delta/tau0
p0      = .9

#- parameters phase-plain
nb_bins     = 50
nb_grpnts   = 500
nb_arrows   = 7

n_model     = 0     # choose model

levels1 = np.arange(-100, 100, 100)

#- initialize
# lambda along x-axis
xmin = -0.02#-0.05
xmax = 0.25
# T along x-axis
ymin = E0 - .25*(E1-E0)
ymax = E1 + .25*(E1-E0)

xplt = np.array([xmin, xmax])
yplt = np.array([Eavg, Eavg])
xplt2 = np.array([xmin, xmax])
yplt2 = np.array([E0, E0])
xplt3 = np.array([xmin, xmax])
yplt3 = np.array([E1, E1])

##- Phase-plain
#- contour
x = np.linspace(xmin, xmax, nb_grpnts)
y = np.linspace(ymin, ymax, nb_grpnts)

X, Y = np.meshgrid(x, y)    # create x- and y-grids

T = fTt(Y, X, tau0, n_model = 0) * (Y>E0)
T[T<=E1] = 0
DT = T - Y

levels2 = np.linspace(np.min(DT), np.max(DT), nb_bins)

#- probability
Delta = E1 - Y
Delta[(Delta<0)] = 0 # negative thresholds do not need to be overcome
Tau0  = Delta/X
P_WoO = np.exp(-k0*Tau0)
#P_WoO[(X<=0)] = 0
#P_WoO[(Y<E0)&(X>0)] = 0
P_WoO[(X<=0)] = np.nan
P_WoO[(Y<E0)&(X>0)] = np.nan
P_WoO[(Y>E1)&(X>0)] = 1

### stats
P_critWoO = np.exp(-k0*delta/lcrit)
print(r'P$_{WoO}$($T_0, \lambda_{crit}$) = ', P_critWoO)
f_critWoO = P_critWoO*tau0/tau
print(r'f$_{WoO}$',f_critWoO)
print(r'P$_{est}$ = ', 1 - (1 - p0)**f_critWoO)


#- quiver
x2 = np.linspace(xmin-0.025, xmax+0.025, nb_arrows)
y2 = np.linspace(ymin-0.025, ymax+0.4, nb_arrows)

X2, Y2 = np.meshgrid(x2, y2) # create x- and y-grids

T2 = fTt(Y2, X2, tau0, n_model = 0) * (Y2>E0)
T2[T2<=E1] = 0
DT2 = T2 - Y2

# Bifurcation line from Fig.2b
ax1.contour(X, Y, DT, levels1, colors='black', linewidths=1.5)

ax1.plot(xplt, yplt, linestyle='--', color='gray', lw=1.05)
ax1.plot(xplt2, yplt2, linestyle='--', color='gray', lw=1.05)
ax1.plot(xplt3, yplt3, linestyle='--', color='gray', lw=1.05)

cs = ax1.contourf(X, Y, P_WoO, np.linspace(0, 1, 101), cmap='YlGn', alpha = 0.95, vmin=0, vmax=1.1)
vt = np.linspace(0, 1, 5+1, endpoint=True)
cb = fig.colorbar(cs, ticks=vt, label=r'Probability for WoO, $P_{WoO}$')

# figure axis
ax1.set_xlabel(r'Growth rate, $\lambda$')
ax1.set_ylabel(r'Initial tolerance, $T_0$')
ax1.set_xlim((xmin,xmax))
ax1.set_ylim((ymin,ymax))
# turn off ticks & labels
ax1.set_xticks([0])
ax1.set_xticklabels(['0'])
ax1.set_yticks([E0, Eavg, E1])
ax1.set_yticklabels([r'$E_{min}$', r'$E_{avg}$', r'$E_{max}$'])

# add labels & arrow
ax1.text(xmin+(xmax-xmin)*.01, ymax + (ymax-ymin)*.03, \
    r'With variance, $\sigma^2(E)>0$')

    
# add explaining text
ax1.text(xmin+(xmax-xmin)*.29, ymin + (ymax-ymin)*0.02, \
    'No establishment')
ax1.text(xmin+(xmax-xmin)*.29, ymin + (ymax-ymin)*0.9, \
    'Establishment', color='white')


##- Finish & saving figure
plt.subplots_adjust(wspace = -0.5, hspace = .4)
plt.savefig('Figure4.png', dpi = dpi, bbox_inches = 'tight')
plt.show()
plt.close()
print('---### Figure 4 created and saved ###---')



###=== FIGURE 5 ===###
#- set seed
seed = 2

# Figure panel properties
nrows   = 2
ncols   = 2
dpi     = 300
plt.rcParams['figure.figsize'] = [tmp[0] * 1.5, tmp[1] * 1.07]

##-- Build panel Figure 5
fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(nrows, ncols)

##-- Figure 5.a: Example Hasselmann correlated noise signals
##- Parameters & settings
#- Environmental condition, E
mu      = 1     # average of stationary signal
tau0    = 1     # approximate period of the signal
tau1    = 10    # approximate period of the signal
sigma   = 0.5   # magnitude of Gaussian noise

#- parameters to rescale epsilon fro similar variances (sigma) of E with various tau
a = 0.74    # rate exponent
b = 0.493   # scaling exponent

#- Tolerance of young organism, T
T0 = mu/10   # initial tolerance, for illustration 10% of mu

#- Temporal settings
Tend    = 100
dt      = .01

#- estimate sigma of E, needed for rescaling
epsilon0 = sigma/(a*(tau0/dt)**b)
epsilon1 = sigma/(a*(tau1/dt)**b)
epsilon2 = sigma/(a*(tau0/10/dt)**b)

##-- Data/results to show
# build Dichotomous Markov Noise sigal (Fig.5.a)

#- Example environmental signal
np.random.seed(seed)
[E0, t] = fHRN(mu, tau0, epsilon0, Tend, dt)
np.random.seed(seed)
[E1, t] = fHRN(mu, tau1, epsilon1, Tend, dt)
np.random.seed(seed)
[E2, t] = fHRN(mu, tau0/10, epsilon2, Tend, dt)
np.random.seed(seed)
[E3, t] = fHRN(mu, tau0, epsilon0/2, Tend, dt)
np.random.seed(seed)
[E4, t] = fHRN(mu, tau0, epsilon0*2, Tend, dt)

print('\n### Statistics of example HRN signal ###')
print('\nmu E0 = %1.3f, mu E1 = %1.3f' % (np.mean(E0),np.mean(E1)))
print('median of E0 = %1.3f and E1 = %1.3f' % (np.quantile(E0, 0.5),np.quantile(E1, 0.5)))
print('sigma E0 = %1.3f, sigma E1 = %1.3f' % (np.std(E0),np.std(E1)))
print('\n#- Potential number of windows:')
W00 = fpotwin(E0, T0, dt)
W01 = fpotwin(E1, T0, dt)
print('\nmu N0(W) = %d and N1(W) = %d' % (len(W00),len(W01)))

path = '../Simulations/results/'
#- Load numerical results
fwoos1 = pd.read_csv(path + 'fwoos_tau_0.10_sigma_0.50.csv')
fwoos2 = pd.read_csv(path + 'fwoos_tau_1.00_sigma_0.50.csv')
fwoos3 = pd.read_csv(path + 'fwoos_tau_10.00_sigma_0.50.csv')
fwoos4 = pd.read_csv(path + 'fwoos_tau_1.00_sigma_0.25.csv')
fwoos5 = pd.read_csv(path + 'fwoos_tau_1.00_sigma_1.00.csv')

lambs   = np.arange(0, 2.6, 0.1)

#- Load E-stats
Estats1 = pd.read_csv(path + 'Estats_tau_0.10_sigma_0.50.csv')
Estats2 = pd.read_csv(path + 'Estats_tau_1.00_sigma_0.50.csv')
Estats3 = pd.read_csv(path + 'Estats_tau_10.00_sigma_0.50.csv')
Estats4 = pd.read_csv(path + 'Estats_tau_1.00_sigma_0.25.csv')
Estats5 = pd.read_csv(path + 'Estats_tau_1.00_sigma_1.00.csv')

print('\n### Statistics of replicate HRN signals ###')
print('tau = 0.1, sigma = 0.5 | N(E<T0) = %1.0f ± %1.0f' % (np.mean(Estats1.iloc[:, 5]),np.std(Estats1.iloc[:, 5])))
print('tau = 1, sigma = 0.5 | N(E<T0) = %1.1f ± %1.0f' % (np.mean(Estats2.iloc[:, 5]),np.std(Estats2.iloc[:, 5])))
print('tau = 10, sigma = 0.5 | N(E<T0) = %1.1f ± %1.0f' % (np.mean(Estats3.iloc[:, 5]),np.std(Estats3.iloc[:, 5])))
print('tau = 1, sigma = 0.25 | N(E<T0) = %1.1f ± %1.0f' % (np.mean(Estats4.iloc[:, 5]),np.std(Estats4.iloc[:, 5])))
print('tau = 1, sigma = 1 | N(E<T0) = %1.1f ± %1.0f' % (np.mean(Estats5.iloc[:, 5]),np.std(Estats5.iloc[:, 5])))


##- Building up panel Figure 5.a
# set-up environmental signal
ax1.plot(t, E2, color='blue', lw=2, label=r'$\langle\tau\rangle=0.1$')
ax1.plot(t, E0, color='black', lw=2, label=r'$\langle\tau\rangle=1$')
ax1.plot(t, E1, color='red', lw=2, label=r'$\langle\tau\rangle=10$')
ax1.plot([0, Tend],[mu, mu],'--', color='black', lw=1.5)
ax1.plot([0, Tend],[T0, T0],'--', color='green', lw=1.5)

# figure axis
ax1.set_xlabel(r'Time, $t$')
ax1.set_ylabel(r'Environment, $E$')
ax1.set_xlim((0,Tend))
ax1.set_ylim(-0.5*sigma, mu + 3*sigma)

##- Building up panel Figure 5.b
# set-up environmental signal
ax3.plot(lambs, np.mean(fwoos1.iloc[:, 1:], axis=0), '.-b', lw=2, label=r'$0.1\langle\tau\rangle, \sigma^2(E)$')
ax3.plot(lambs, np.mean(fwoos2.iloc[:, 1:], axis=0), '.-k', lw=2, label=r'$\langle\tau\rangle\, \sigma^2(E)$')
ax3.plot(lambs, np.mean(fwoos3.iloc[:, 1:], axis=0), '.-r', lw=2, label=r'$10\langle\tau\rangle, \sigma^2(E)$')

# error bars red, tau = 10
y = np.mean(fwoos3.iloc[:, 1:], axis=0)
yerr0 = np.quantile(fwoos3.iloc[:, 1:],0.9, axis=0) - y
yerr1 = y - np.quantile(fwoos3.iloc[:, 1:],0.1, axis=0)
ax3.errorbar(lambs, y, yerr=[yerr0, yerr1], ls='None', color='red')

# error bars black, tau = 1
y = np.mean(fwoos2.iloc[:, 1:], axis=0)
yerr0 = np.quantile(fwoos2.iloc[:, 1:],0.9, axis=0) - y
yerr1 = y - np.quantile(fwoos2.iloc[:, 1:],0.1, axis=0)
ax3.errorbar(lambs, y, yerr=[yerr0, yerr1], ls='None', color='black')

# error bars blue, tau = 0.1
y = np.mean(fwoos1.iloc[:, 1:], axis=0)
yerr0 = np.quantile(fwoos1.iloc[:, 1:],0.9, axis=0) - y
yerr1 = y - np.quantile(fwoos1.iloc[:, 1:],0.1, axis=0)
ax3.errorbar(lambs, y, yerr=[yerr0, yerr1], ls='None', color='blue')

# figure axis
ax3.set_xlabel(r'Growth rate, $\lambda$')
ax3.set_ylabel(r'Fraction attainable windows, $f_{WoO}$')
ax3.set_xlim((0,2.525))
ax3.set_ylim((0,1))
ax3.legend(loc='upper left', ncol=1, labelspacing=0.15)

##- Building up panel Figure 5.c
# set-up environmental signal
ax2.plot(t, E4, color='blue', lw=2, label=r'$2\sigma^2(E)$')
ax2.plot(t, E0, color='black', lw=2, label=r'$\sigma^2(E)$')
ax2.plot(t, E3, color='red', lw=2, label=r'$0.5\sigma^2(E)$')
ax2.plot([0, Tend],[mu, mu],'--', color='black', lw=1.5)
ax2.plot([0, Tend],[T0, T0],'--', color='green', lw=1.5)

# figure axis
ax2.set_xlabel(r'Time, $t$')
ax2.set_xlim((0,Tend))
ax2.set_ylim(-0.5*sigma, mu + 3*2*sigma)

##- Building up panel Figure 5.d
## error bars red, sigma = 0.25
y = np.mean(fwoos4.iloc[:, 1:], axis=0)
yerr0 = np.quantile(fwoos4.iloc[:, 1:],0.9, axis=0) - y
yerr1 = y - np.quantile(fwoos4.iloc[:, 1:],0.1, axis=0)
ax4.errorbar(lambs-1*1e-2, y, yerr=[yerr0, yerr1], ls='None', color='red')

# error bars black, sigma = 0.5
y = np.mean(fwoos2.iloc[:, 1:], axis=0)
yerr0 = np.quantile(fwoos2.iloc[:, 1:],0.9, axis=0) - y
yerr1 = y - np.quantile(fwoos2.iloc[:, 1:],0.1, axis=0)
ax4.errorbar(lambs, y, yerr=[yerr0, yerr1], ls='None', color='black')

# error bars blue, sigma = 1.0
y = np.mean(fwoos5.iloc[:, 1:], axis=0)
yerr0 = np.quantile(fwoos5.iloc[:, 1:],0.9, axis=0) - y
yerr1 = y - np.quantile(fwoos5.iloc[:, 1:],0.1, axis=0)
ax4.errorbar(lambs+1*1e-2, y, yerr=[yerr0, yerr1], ls='None', color='blue')

# numerical results
ax4.plot(lambs, np.mean(fwoos4.iloc[:, 1:], axis=0), '.-r', lw=2, label=r'$\langle\tau\rangle, 0.5\sigma^2(E)$')
ax4.plot(lambs, np.mean(fwoos2.iloc[:, 1:], axis=0), '.-k', lw=2, label=r'$\langle\tau\rangle, \sigma^2(E)$')
ax4.plot(lambs, np.mean(fwoos5.iloc[:, 1:], axis=0), '.-b', lw=2, label=r'$\langle\tau\rangle, 2\sigma^2(E)$')

# figure axis
ax4.set_xlabel(r'Growth rate, $\lambda$')
ax4.set_xlim((0,2.525))
ax4.set_ylim((0,1))
ax4.legend(loc='upper left', ncol=1, labelspacing=0.15)

ax1.text(Tend/2, 1.1*(mu + 3*sigma), r'effect of period, $\langle\tau\rangle$', fontsize=11, fontstyle='italic', horizontalalignment='center')
ax2.text(Tend/2, 1.1*(mu + 3*2*sigma), r'effect of variance, $\sigma^2(E)$', fontstyle='italic', fontsize=11, horizontalalignment='center')

# panel labels
ax1.text(-0.23*Tend, 2.5, '(a)', fontsize=12, fontweight='bold')
ax2.text(-0.18*Tend, 1.57*2.5, '(c)', fontsize=12, fontweight='bold')
ax3.text(-0.23*2.525, 1, '(b)', fontsize=12, fontweight='bold')
ax4.text(-0.18*2.525, 1, '(d)', fontsize=12, fontweight='bold')

# finish figure
plt.subplots_adjust(wspace = 0.2, hspace = .4)
plt.savefig('Figure5.png', dpi = dpi, bbox_inches = 'tight')
plt.show()
plt.close()
print('---### Figure 5 created and saved ###---')




###=== FIGURE 6 ===###
# Figure panel properties
nrows   = 1
ncols   = 1
dpi     = 300
plt.rcParams['figure.figsize'] = [tmp[0] * 1.2, tmp[1] * 1.07]

##-- Build panel Figure 6
fig, ax1 = plt.subplots(nrows, ncols)

##- Parameters & settings
#- Environmental condition, E
tau     = 10    # approximate period
k       = 1/tau # transition rate

# gradient in Eavg:
mu_min  = 0
mu_max  = 1.5

# gradient in variance, sigma^2(E):
sig_min = 0
sig_max = 1.5

#- Properties of the organism
T0      = 0.5   # Initial tolerance of young organism
lamb    = 0.05  # growth rate [/d]

#
N       = 500
Eavg    = np.linspace(mu_min, mu_max, N+1)
sigma   = np.linspace(sig_min, sig_max, N+1)

Mu, Sigma = np.meshgrid(Eavg, sigma)

E0      = Mu - Sigma
E1      = Mu + Sigma

Delta   = E1 - T0
Delta[Delta<0] = 0
Tau_crit= Delta/lamb
PWoO    = np.exp(-k*Tau_crit)
PWoO[T0<E0] = np.nan


##-- Figure 6: non-stationary
v = np.linspace(0, 1, 101, endpoint=True)
vt= np.linspace(0, 1, 11, endpoint=True)

# Fig 6: WoO
cs1 = ax1.contourf(Mu, Sigma, PWoO, v, cmap='YlGn', alpha=0.95)
cb1 = fig.colorbar(cs1, ax=ax1, ticks=vt)

ax1.plot([T0,T0], [sig_min,sig_max], '--g')
ax1.plot([T0,mu_max], [sig_min,1], '--r')
ax1.plot([T0,mu_min], [sig_min,0.5], '--w')

ax1.set_xlabel(r'Average, $E_{avg}$')
ax1.set_ylabel(r'Variance, $\sigma^2(E)$')
cb1.set_label(r'Probability of WoO, $P_{WoO}$', rotation=270, labelpad=15)

ax1.text(T0/2, 1.35, 'benign', fontweight='bold', horizontalalignment='center')
ax1.text(T0 + (mu_max - T0)/2, 1.35, 'harsh', fontweight='bold', horizontalalignment='center')
ax1.text(T0 + (mu_max - T0)/2, 0.1, 'no establishment', color='red', horizontalalignment='center')
ax1.text(T0, -0.065, r'T$_0$', horizontalalignment = 'center')

ax1.text(0.2, 0.1, 'regime i', horizontalalignment='center')
ax1.text(T0/2, 0.75, 'regime ii', horizontalalignment='center')
ax1.text(T0 + (mu_max - T0)/2, 0.75, 'regime iii', horizontalalignment='center')
ax1.text(T0 + (mu_max - T0)/2, 0.25, 'regime iv', horizontalalignment='center', color='red')

# finish figure
plt.subplots_adjust(wspace = 0.2, hspace = .4)
plt.savefig('Figure6.png', dpi = dpi, bbox_inches = 'tight')
plt.show()
plt.close()
print('---### Figure 6 created and saved ###---')




###=== FIGURE 7 ===###
#- set seed
np.random.seed(7)

# Figure panel properties
nrows   = 3
ncols   = 1
dpi     = 300

plt.rcParams['figure.figsize'] = [tmp[0] * .8, tmp[1] * 1.8]

##-- Build panel Figure 7
fig, (ax1, ax2, ax3) = plt.subplots(nrows, ncols)

##-- Figure 7.a: Conceptual model
##- Parameters & settings
#- Environmental condition, E
E0      = 0.        # Emin
E1      = 1         # Emax
tau0    = 10.       # waiting time in Emin
tau1    = 2.5       # waiting time in Emax
k0      = 1./tau0   # transition rate from Emin to Emax
k1      = 1./tau1   # transition rate from Emin to Emax

tau     = tau0 + tau1   # total waiting time of the environmental condition
Eavg    = (tau0*E0 + tau1*E1)/tau

#- Tolerance of young organism, T
T0      = E0 + (Eavg-E0)*1.15   # initial tolerance (for illustration it is in between Emin and Eavg)
delta   = E1 - T0

#- Temporal settings
Tend    = 300
dt      = .01

##-- Data/results to show
# build regular dichotomous sigal
[E, t] = fRDS(E0, E1, k0, k1, Tend, dt)

# arrow setting
hl = tau/25     # head_length
hw = (E1-E0)/20  # head_width

##- Building up panel Figure 7.a
x = np.array([tau+tau0-0.4*tau1, tau+tau0+0.4*tau1])
y = np.array([E1, E1])
ax1.fill_between(x, y, y2=E0, fc='#F5DEB3')

# set-up environmental signal
ax1.plot(t, E, color='gray', lw=2)
ax1.plot(t[int(tau/dt):int(2*tau/dt)], E[int(tau/dt):int(2*tau/dt)], '-k', lw=2)

# add establishing organisms
ax1.arrow(tau, T0, tau, delta/tau0*tau, head_width=hw, head_length=hl, lw=2, fc='g', ec='g')
ax1.annotate("", xy=(tau+tau0, E0-0.1*(E1-E0)), xytext=(tau, E0-0.1*(E1-E0)), \
    arrowprops=dict(arrowstyle="|-|", mutation_scale=2))
ax1.annotate("", xy=(tau+tau0-0.4*tau1, E1-delta*.5), xytext=(tau+tau0+0.4*tau1, E1-delta*.5), \
    arrowprops=dict(arrowstyle="<->"))
ax1.plot(tau, T0, 'og')

# add labels
ax1.text(tau+tau0/2, E0 - (E1-E0)*.25, r'$\tau_{low}$', horizontalalignment = 'center')

# figure axis
ax1.set_xlim((tau*.6,2.2*tau))
ax1.set_ylim((E0-.4*(E1-E0),1.5*E1))
# turn off ticks & labels
ax1.axis('off')


##- Building up panel Figure 7.a
x1 = np.array([tau-.7, tau+tau0+0.2*tau1])
y1 = np.array([T0+.15, T0+.15])
y2 = np.array([T0-.15, T0-.15])
ax2.fill_between(x1, y1, y2=y2, fc='#F5DEB3')
x1 = np.array([tau+tau0, 2*tau])
y1 = np.array([E1+.15, E1+.15])
y2 = np.array([E1-.15, E1-.15])
ax2.fill_between(x1, y1, y2=y2, fc='#F5DEB3')

# set-up environmental signal
ax2.plot(t, E, color='gray', lw=2)
ax2.plot(t[int(tau/dt):int(2*tau/dt)], E[int(tau/dt):int(2*tau/dt)], '-k', lw=2)

# add establishing organisms
ax2.arrow(tau, T0, tau, delta/tau0*tau, head_width=hw, head_length=hl, lw=2, fc='g', ec='g')
ax2.annotate("", xy=(tau+tau0+tau1/6, E1), xytext=(tau+tau0+tau1/6, T0), \
    arrowprops=dict(arrowstyle="|-|", mutation_scale=2))
ax2.plot(tau, T0, 'og')
ax2.annotate("", xy=(tau-.5, T0-.15), xytext=(tau-.5, T0+.15), \
    arrowprops=dict(arrowstyle="<->"))
ax2.annotate("", xy=(tau+tau0+tau1/2, E1-.15), xytext=(tau+tau0+tau1/2, E1+.15), \
    arrowprops=dict(arrowstyle="<->"))

# add labels
ax2.text(tau*.6 + 1.05*tau, E1+.05*(E1-E0), r'E$_{high}$')
ax2.text(tau+tau0+tau1/4, E1 - delta*.5, r'$\delta$')
ax2.text(tau-1.5, T0, r'T$_0$')

# figure axis
ax2.set_xlim((tau*.6,2.2*tau))
ax2.set_ylim((E0-.4*(E1-E0),1.5*E1))
# turn off ticks & labels
ax2.axis('off')


##- Building up panel Figure 7.c
x1 = np.array([tau, 2*tau+hl])
y1 = np.array([T0, T0+1.27*delta/tau0*(tau+hl)])
y2 = np.array([T0, T0+0.73*delta/tau0*(tau+hl)])
ax3.fill_between(x1, y1, y2=y2, fc='#F5DEB3')

# set-up environmental signal
ax3.plot(t, E, color='gray', lw=2)
ax3.plot(t[int(tau/dt):int(2*tau/dt)], E[int(tau/dt):int(2*tau/dt)], '-k', lw=2)

# add establishing organisms
ax3.arrow(tau, T0, tau, delta/tau0*tau, head_width=hw, head_length=hl, lw=2, fc='g', ec='g')
ax3.plot(tau, T0, 'og')
ax3.annotate("", xy=(tau+tau0*.75, T0+1.27*delta/tau0*(tau0*.75)+0.05), xytext=(tau+tau0*.75, 0.73*delta/tau0*(tau*.75)+0.05), arrowprops=dict(arrowstyle="<->"))

# add labels
ax3.text(tau+tau0/2, E0 + 1.3*delta, r'$\lambda_{crit}$', horizontalalignment = 'center')
ax3.text(tau+tau0/2, E0 + 1.05*delta, '\u2220', horizontalalignment = 'center')

# figure axis
ax3.set_xlim((tau*.6,2.2*tau))
ax3.set_ylim((E0-.4*(E1-E0),1.5*E1))
# turn off ticks & labels
ax3.axis('off')

# add titles to panels
ax1.set_title(r'Length of the benign window, $\tau_{low}$')
ax2.set_title(r'Height of the tolerance threshold, $\delta$')
ax3.set_title(r'Rate of tolerance growth, $\lambda_{crit}$')


##- Finish & saving figure
plt.subplots_adjust(wspace = -0.5, hspace = .4)
plt.savefig('Figure7.png', dpi = dpi, bbox_inches = 'tight')
plt.show()
plt.close()


print('\n\t###=== THE END ===###\n')
