# Gouldian2.py

# Import plotting routines
from pylab import *
import math as m
import numpy as np

# Generating stochasticity in epsilon
mean = 5  # average number of times for 1 success
lam = 1.0/mean
def randd(n):
	for i in range(n):
		a = np.random.rand()
		if a < 1 - m.exp(-1.0*lam): a = 1.0
		else: a = 0.0
	return a

# Simulation parameters


# Control parameter of the map:
delta = 0.1
c = 3.0
kslope = (c-1)/(pi/2)
eta = (sin(0.01 + 1.0/4) - sin(0.01)) / 200    # eta = Vmax/(100*2) = V(Smin, Wmax)/200
phi = 0.2

# Set up an array of iterates and set the initital condition
elist = [0.01+pi/4]
slist = [pi/2-0.01-pi/4]
klist = [c - slist[0]*kslope]

betalist = [0.01]


wlist = []
vlist = []
vplist = []
wplist = []
# The number of iterations to generate
N = 50000

	 

# The main loop that generates iterates and stores them
for n in xrange(0,N):
  
  epsilon = randd(1)
   
# Memory Box
  w = np.random.pareto(klist[n]) + 1  
  if w <= 2.0: w = 1/(2.0*pi)
  elif w > 2.0: w = 1.0/4
  wlist.append(w)
  
  if epsilon == 0.0: w_actualized = 0.0
  elif epsilon == 1.0: w_actualized = w
 
  v = sin(slist[n]+w_actualized)-sin(slist[n])
  if slist[n]+w_actualized >= pi/2: v = 1 - sin(slist[n])  
  vlist.append(v)
  
  vp = max(vlist)
  vplist.append(vp)

  if v >= vp-eta and v <= vp+eta: wp = w_actualized
  if max(vlist[0:n+1]) > 0.0 and v == 0.0: wp = wplist[-1]
  wplist.append(wp)
 
# Coupled 2-D Map 
  
  if w_actualized == wp:
  	e = elist[n]*exp(delta*(pi/2-elist[n])/(pi/2)) - w_actualized
	beta = 0.01
  elif w_actualized != wp:
  	e = elist[n]*exp(delta*(pi/2-elist[n])/(pi/2)) - w_actualized * betalist[n]
  	beta = betalist[n]*exp(phi)
  		
  if e <= 0.01: e = 0.01
  elif e >= pi/2-0.01: e = pi/2-0.01
  elist.append(e)  
  
  if beta <= 0.01: beta = 0.01
  elif beta >= 0.99: beta = 0.99
  betalist.append(beta)
  
  s = pi/2-e
  slist.append(s)
  
  k = c - s*kslope
  if k <= 1.0: k = 1.0
  elif k >= c: k = c
  klist.append(k)		
  


  
#  print 't:', n+1, 'eps', epsilon, 's:', slist[n+1], 'w:', w, 'w_act:', w_actualized, 'v:', v, 'v_choice:', vp, 'w_choice:', wp, 'beta:', betalist[n+1]

  


# Setup the plot
xlabel('Time step n') # set x-axis label
ylabel('s(n)') # set y-axis label
title('s_t at delta= ' + str(delta)) # set plot title
# Plot the time series: once with circles, once with lines
#plot(slist , 'b')
plot(slist, 'b')
plot(klist, 'r')
# Use command below to save figure
#savefig('LogisticMap', dpi=600)
#print 'u', u, '\nk', k
# Display the plot in a window
show()