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)

# 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'])  

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"

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

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)