import numpy as np
import scipy as sp
import matplotlib.pyplot as plt
from scipy.misc import imread
%matplotlib inline

# download the image to /tmp and open with scipy
!wget http://upload.wikimedia.org/wikipedia/commons/a/a9/Fuel_economy_vs_speed_1997.png -O /tmp/Fuel_economy_vs_speed_1997.png
pic = imread('/tmp/Fuel_economy_vs_speed_1997.png')

# the image we got
plt.imshow(pic,origin='upper');

# defining the model
beta = 2*(40)**3
alpha = 3*(40)**2*35
model = lambda v: alpha/( v**2 + 2*(40)**3/v )

# just checking that the model looks reasonable
vs = np.linspace(1,90,1000)
mpgs = model(vs)
plt.plot( vs, mpgs);

# looking at hte data, get the conversion factors from mph to pixels and mpg to pixels
vscale = ( 420-100) / 90.
mpgscale = ( 400 - 150 ) / 40.
# transform the model into pixel space
vpic = 100 + vscale*vs
mpgpic = 380 - (mpgscale*mpgs)

#Plot poth on top of one another
plt.figure(figsize=(22,16));
plt.imshow(pic, origin='upper');
plt.plot(vpic, mpgpic, linewidth=30, alpha=0.5);
plt.xlim((0,pic.shape[1]));
plt.ylim((pic.shape[0],0));
plt.axis('off');
#plt.savefig('modified.png')

gasprice = 3.753  # in dollars from wolfram alpha
costperhour = vs / mpgs  * gasprice 

plt.plot( vs, costperhour, lw=3 )
plt.grid()
plt.xlabel('speed (mph)')
plt.ylabel('cost ($/hour)');
#plt.savefig('costperhour.png');
plt.title("Operating cost per hour for average US vehicle");

# I get the change by taking a gradient
dcost5mph = np.gradient(costperhour, vs[1]-vs[0])*5

plt.plot(vs, dcost5mph, lw=3)
plt.grid()
plt.xlabel('speed (mph)')
plt.ylabel('dcost ($/hour/5 mph)');
#plt.savefig('dcostperhour.png');
plt.title("Change in operating cost for 5 mph difference in speed per hour");

costpermile = gasprice*10./mpgs
plt.plot(vs[100:], costpermile[100:], lw=3)
plt.grid()
plt.xlabel('speed (mph)')
plt.ylabel('$/10 miles');
plt.title("Cost per 10 miles");
#plt.savefig("costpermile.png");

dcostpermile = np.gradient(costpermile, vs[1]-vs[0])*5;
plt.plot(vs[100:], dcostpermile[100:], lw=3)
plt.grid()
plt.xlabel('speed (mph)')
plt.ylabel(r'delta $/10 miles/5 mph difference');
plt.title("Difference in cost per 10 miles for 5 mph change in speed");
#plt.savefig("dcostpermile.png");

block_size = 1./10  # in miles
speed = 30  # mph
#speed = 35
time_between_lights = block_size / speed * 60 * 60 # in seconds
off_time = 50  # seconds
on_time = 50  # seconds
p_on = on_time *1.0/(on_time + off_time)
n = 100

p_on

from random import random
def next_block(speed):
    time_between_lights = block_size / speed * 60 * 60
    if random() < p_on:
        return time_between_lights
    else:
        return time_between_lights + random() * off_time
    
    
def difftimes(nblocks, speed1, speed2):
    times = np.array([ sum( next_block(speed1) for i in xrange(nblocks) ) for i in xrange(10000) ])
    times2 = np.array([ sum( next_block(speed2) for i in xrange(nblocks) )  for i in xrange(10000) ])
    return times - times2

bins = np.linspace(-300,300,100)
plt.hist( difftimes(5, 45,50), bins, histtype='step', normed=True, alpha=0.6 ,lw=3, label='N=5');
plt.hist( difftimes(15, 45,50), bins, histtype='step', normed=True, alpha=0.6 ,lw=3, label='N=15');
plt.hist( difftimes(50, 45,50), bins, histtype='step', normed=True, alpha=0.6 ,lw=3, label='N=50');
plt.legend()
plt.xlim((-300,300));

# average speed
tau = 30 # seconds
vs = np.linspace(10,90,1000)
p = 0.65
d = 1./10

average_speed = 60* 60 * d / ( d/vs * 60 * 60 + tau/2 * (1-p) )
plt.plot(vs, average_speed, lw=3);
plt.grid()
plt.xlabel('target speed (mph)');
plt.ylabel('average speed (mph)');
#plt.savefig('averagecityspeed.png');

# Gains in time
plt.plot( vs, 5./vs , lw=3, label="5 mph gain");
plt.xlabel('target speed (mph)');
plt.ylabel('fractional time');

sigma = np.sqrt( tau**2 / 12. * ( 1 - p ) * ( 3 * p - 1) )
mu = tau/2. * ( 1 - p )
avgtime = d/vs * 60 * 60 + mu
plt.plot( vs, sigma/avgtime, ls='--', label='N=1' );
plt.plot( vs, sigma/avgtime/np.sqrt(10), ls='--', label='N=10' );
plt.plot( vs, sigma/avgtime/10.0, ls='--', label='N=100' );
plt.yscale('log');
plt.legend(loc=1);
ax = plt.gca()

plt.grid();
plt.tight_layout()
#plt.savefig('gainsversuslights.png');

plt.plot(vs, avgtime);