import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
from scipy.optimize import minimize
from sklearn.linear_model import LinearRegression

H = np.array([85,75,70,85]*20)
T = np.array([65,70,85,85]*20)
#1
def eyring_gen(A,B,C,sigma):
  Y = np.exp(np.log(A) + B*H + C/(T+273)+sigma*np.random.normal(0,1,80))
  return Y

# 2

Y = eyring_gen(np.exp(-13.68),-.0422,8485,.0225)
plt.subplot(121)
plt.plot(H,Y,'ro')
plt.title('Humidity vs Failure Time')
plt.subplot(122)
plt.plot(T,Y,'bo')
plt.title('Temp vs Fail Time')
plt.show()

# 3

def logreg(Y,H,T):
  datmat = pd.DataFrame([H,1/(T+273)]).T
  X = datmat.as_matrix()
  mod = LinearRegression()
  mod.fit(X,np.log(Y))
  return(np.append(mod.intercept_,mod.coef_))

logreg(Y,H,T)

# 4 

n_run = 100
simvals = np.zeros((3,n_run))
for i in range(n_run):
  Y = eyring_gen(np.exp(-13.68),-.0422,8485,.0225)
  simvals[:,i] = logreg(Y,H,T)

simvals.mean(axis=1) - np.array([-13.68,-.0422,8485])

# 5

def target(coef):
  fitted = coef[0]*np.exp(coef[1]*H + coef[2]/(T+273))
  resids = Y - fitted
  return(np.sum(resids**2))

def gradient(coef):
  fitted = coef[0]*np.exp(coef[1]*H + coef[2]/(T+273))
  resids = Y - fitted
  ddA = -2*np.sum(resids*np.exp(coef[1]*H + coef[2]/(T+273)))
  ddB = -2*np.sum(resids*np.exp(coef[1]*H + coef[2]/(T+273))*coef[0]*H)
  ddC = -2*np.sum(resids*np.exp(coef[1]*H + coef[2]/(T+273))*coef[0]/(T+273))
  return(np.array([ddA,ddB,ddC]))

inits = np.array([np.exp(-13.68),-.0422,8485])

out = minimize(target,inits,jac=gradient,method='BFGS')

n_run = 100
simvals = np.zeros((3,n_run))
for i in range(n_run):
  Y = eyring_gen(np.exp(-13.68),-.0422,8485,.0225)
  simvals[:,i] = minimize(target,inits,jac=gradient,method='BFGS').x

simvals.mean(axis=1)
simvals.mean(axis=1) - np.array([np.exp(-13.68),-.0422,8485])

# 7

n_run = 100
pred_ll = np.zeros(n_run)
for i in range(n_run):
  Y = eyring_gen(np.exp(-13.68),-.0422,8485,.0225)
  coefs  = logreg(Y,H,T)
  pred_ll[i] = np.exp(coefs[0])*np.exp(coefs[1]*50+coefs[2]/(25+273))


pred_nl = np.zeros(n_run)
for i in range(n_run):
  Y = eyring_gen(np.exp(-13.68),-.0422,8485,.0225)
  coefs = minimize(target,inits,jac=gradient,method='BFGS').x
  pred_nl[i] = coefs[0]*np.exp(coefs[1]*50+coefs[2]/(25+273))

# 8 
pred_nl = np.zeros(n_run)
datmat = pd.DataFrame([Y,H,T]).T
X = datmat.as_matrix()
for i in range(n_run):
  newdat = X[np.random.choice(range(80),80),:]
  Y = newdat[:,0]
  H = newdat[:,1]
  T = newdat[:,2]
  coefs = minimize(target,inits,jac=gradient,method='BFGS').x
  pred_nl[i] = coefs[0]*np.exp(coefs[1]*50+coefs[2]/(25+273))




