Code by Scott B Miles.
import numpy as np
from pandas import *
from __future__ import division
from scipy.integrate import simps, trapz
import scipy.stats as stats
from matplotlib.font_manager import FontProperties
pandas.set_option('display.height', 500)
pandas.set_option('display.max_rows', 500)
pandas.set_option('display.max_columns', 500)
pandas.set_option('display.width', 500)
pandas.set_option('display.mpl_style', False)
rcParams['figure.figsize'] = 25, 25
t_run = 30
t_dt = linspace(1,t_run,t_run)
# Function/data for defining Q(t) -- "Quantity Resilience"
function_supply = "Binomial"
function_demand = "Binomial"
# Q_supply(t) defined by binomial function
if function_supply == "Binomial":
cdf = stats.binom.cdf
mu_supply = 90
sigma_supply = 0.0045
supply = Series(cdf(t_dt,mu_supply, sigma_supply), index=t_dt, name='Supply')
print "Binomial supply quantity curve"
print "Mu = " + str(mu_supply) + " Sigma = " + str(sigma_supply)
print
elif function_supply == "Empirical":
supply = read_excel('CanterburyHotels.xlsx', 'data_norm', parse_cols=[0,1], index_col=0, na_values=['..'])
#supply.index = t_dt
print "Empirical supply quantity curve: " + 'CanterburyHotels.xlsx'
else:
print "Supply quantity curve undefined"
print
# Q_demand(t) defined by binomial function
if function_demand == "Binomial":
cdf = stats.binom.cdf
mu_demand = 30
sigma_demand = 0.125
demand = Series(cdf(t_dt,mu_demand, sigma_demand), index=t_dt, name='Demand')
print "Binomial demand quantity curve"
print "Mu = " + str(mu_demand) + " Sigma = " + str(sigma_demand)
print
elif function_supply == "Empirical":
demand = read_excel('CanterburyHotels.xlsx', 'data_norm', parse_cols=[0,2], index_col=0, na_values=['..'])
#demand.index = t_dt
print "Empirical demand quantity curve: " + 'CanterburyHotels.xlsx'
else:
print "Demand quantity curve undefined"
recovery = concat([supply, demand], axis=1)
Binomial supply quantity curve Mu = 90 Sigma = 0.0045 Binomial demand quantity curve Mu = 30 Sigma = 0.125
Calculate adaptation resilience and speed resilience using the respective two formulas:
$A(t) =\frac{d^2Q(t)}{dt^2}$
$V(t) = \frac{dQ(t)}{dt}$
# Central finite difference function to estimate derivatives
def deriv(x, y):
"Finite difference derivative of the function f"
n = len(y)
d = zeros(n,'d') # assume double
# Use centered differences for the interior points, one-sided differences for the ends
for i in range(1,n-1):
d[i] = (y[i+1]-y[i])/(x[i+1]-x[i])
d[0] = (y[1]-y[0])/(x[1]-x[0])
d[n-1] = (y[n-1]-y[n-2])/(x[n-1]-x[n-2])
return d
# Calculate the first derivative of Q(t) to get V(t) -- "Speed Resilience"
recovery['Supply Speed'] = deriv(recovery.index, recovery['Supply'])
recovery['Demand Speed'] = deriv(recovery.index, recovery['Demand'])
# Calculate the second derivative of Q(t) to get A(t) -- "Adaptation Resilience"
recovery['Supply Adaptation'] = deriv(recovery.index, recovery['Supply Speed'])
recovery['Demand Adaptation'] = deriv(recovery.index, recovery['Demand Speed'])
Calculate Time and Sufficiency resilience using the following formulas:
$R(t) = \frac{1}{t} \int_{t_0}^t Q(t)\mathrm{d}t$
$S(t) = \frac{1}{2} [R_{s}(t) + R_{d}(t) - (R_{s}(t) + R_{d}(t)) \vert R_{supply}(t) - R_{d}(t) \vert]$
recovery['Demand Efficiency'] = recovery['Supply'] * 0.0
recovery['Supply Efficiency'] = recovery['Supply'] * 0.0
recovery['Demand Supply Ratio'] = recovery['Supply'] * 0.0
recovery['Sufficiency'] = recovery['Supply'] * 0.0
#recovery['Demand Supply Ratio'] = recovery['Supply'] / recovery['Demand']
for i in (recovery.index):
recovery['Demand Efficiency'][i] = simps(recovery['Demand'][0:i], dx=1, even='first') / i
recovery['Supply Efficiency'][i] = simps(recovery['Supply'][0:i], dx=1, even='first') / i
recovery['Sufficiency'] = 0.5*(recovery['Supply Efficiency'] + recovery['Demand Efficiency'] -
(recovery['Supply Efficiency'] + recovery['Demand Efficiency']) *
abs(recovery['Supply Efficiency'] - recovery['Demand Efficiency']))
recovery.fillna(0, inplace=True)
t_supply_adapt_max = int(recovery.index[argmax(recovery['Supply Adaptation'])])
t_demand_adapt_max = int(recovery.index[argmax(recovery['Demand Adaptation'])])
t_supply_speed_max = int(recovery.index[argmax(recovery['Supply Speed'])])
t_demand_speed_max = int(recovery.index[argmax(recovery['Demand Speed'])])
# Fraction defining "recovered" 1 = 100%
recovered = 0.99
if len(recovery['Supply'][recovery['Supply'] > recovered]) == 0:
t_supply_rec = t_run
else:
t_supply_rec = int(min(recovery['Supply'][recovery['Supply'] > recovered].index))
if len(recovery['Demand'][recovery['Demand'] > recovered]) == 0:
t_demand_rec = t_run
else:
t_demand_rec = int(min(recovery['Demand'][recovery['Demand'] > recovered].index))
supply_demand_diff = abs(recovery['Demand'] - recovery['Supply'])
if len(supply_demand_diff[supply_demand_diff > recovered]) == 0:
t_sufficiency = t_run
else:
t_sufficiency = int(min(supply_demand_diff[supply_demand_diff > recovered].index))
print "Maximum supply adaptation time = ", t_supply_adapt_max, "days"
print "Maximum demand adaptation time = ", t_demand_adapt_max, "days"
print "Maximum supply speed time = ", t_supply_speed_max, "days"
print "Maximum demand speed time = ", t_demand_speed_max, "days"
print "Supply quantity recovery time = ", t_supply_rec, "days"
print "Demand quantity recovery time = ", t_demand_rec, "days"
Maximum supply adaptation time = 2 days Maximum demand adaptation time = 2 days Maximum supply speed time = 2 days Maximum demand speed time = 3 days Supply quantity recovery time = 2 days Demand quantity recovery time = 8 days
r_supply_adapt = max(recovery['Supply Adaptation']) - min(recovery['Supply Adaptation'])
r_demand_adapt = max(recovery['Demand Adaptation']) - min(recovery['Demand Adaptation'])
r_supply_speed = max(recovery['Supply Speed'])
r_demand_speed = max(recovery['Demand Speed'])
r_supply_quant = recovery['Supply'][-1]
r_demand_quant = recovery['Demand'][-1]
r_supply = recovery['Supply Efficiency'][-1]
r_demand = recovery['Demand Efficiency'][-1]
r_sufficiency = recovery['Sufficiency'][-1]
print "Supply adaptation resilience = ", r_supply_adapt
print "Demand adaptation resilience = ", r_demand_adapt
print "Supply speed resilience = ", r_supply_speed
print "Demand speed resilience = ", r_demand_speed
print "Supply quantity resilience = ", r_supply_quant
print "Demand quantity resilience = ", r_demand_quant
print "Supply efficiency resilience = ", r_supply
print "Demand efficiency resilience = ", r_demand
print "Sufficiency resilience = ", r_sufficiency
Supply adaptation resilience = 0.101834969431 Demand adaptation resilience = 0.224120970033 Supply speed resilience = 0.054532927417 Demand speed resilience = 0.21551301576 Supply quantity resilience = 1.0 Demand quantity resilience = 1.0 Supply efficiency resilience = 0.965596430183 Demand efficiency resilience = 0.889789144955 Sufficiency resilience = 0.857366915817
subplot(511, axisbg='.97')
plot(recovery.index, recovery['Supply Adaptation'], label=r'$A_{s}(t)$', color='k')
plot(recovery.index, recovery['Demand Adaptation'], label=r'$A_{d}(t)$', color='k', linestyle='--')
tick_params(axis='both', which='major', labelsize=20)
ylabel('Adaptation', fontsize=25, family='serif')
xlim([1,t_run])
legend(loc='best', fontsize=20)
subplot(512, axisbg='.97')
plot(recovery.index, recovery['Supply Speed'], label=r'$V_{s}(t)$', color='k')
plot(recovery.index, recovery['Demand Speed'], label=r'$V_{d}(t)$', color='k', linestyle='--')
tick_params(axis='both', which='major', labelsize=20)
ylabel('Speed', fontsize=25, family='serif')
xlim([1,t_run])
legend(loc='best', fontsize=20)
subplot(513, axisbg='.97')
plot(recovery.index, recovery['Supply'], label=r'$Q_{s}(t)$', color='k')
plot(recovery.index, recovery['Demand'], label=r'$Q_{d}(t)$', color='k', linestyle='--')
tick_params(axis='both', which='major', labelsize=20)
ylabel('Quantity', fontsize=25, family='serif')
xlim([1,t_run])
legend(loc='best', fontsize=20)
subplot(514, axisbg='.97')
plot(recovery.index, recovery['Supply Efficiency'], label=r'$R_{s}(t)$', color='k')
plot(recovery.index, recovery['Demand Efficiency'], label=r'$R_{d}(t)$', color='k', linestyle='--')
tick_params(axis='both', which='major', labelsize=20)
ylabel('Efficiency', fontsize=25, family='serif')
xlim([1,t_run])
legend(loc='best', fontsize=20)
subplot(515, axisbg='.97')
plot(recovery.index, recovery['Sufficiency'], label=r'$S(t)$', color='k')
tick_params(axis='both', which='major', labelsize=20)
ylabel('Sufficiency', fontsize=25, family='serif')
xlabel('Time', fontsize=25)
xlim([1,t_run])
legend(loc='best', fontsize=20)
<matplotlib.legend.Legend at 0x825bd90>