# TPL_CNN.py
# this program provides the infrastructure to build and manipulate arbitrarily large MxN arrays
# it is called on by programs like CNN_lattice_simulator.py, but can also be used on its own, although it is not interactive

from random import random
from numpy import cos, pi, array, mod, ones, zeros	#*	# must have cos, pi, array, mod
###################
## ** WARNING ** ## 	Do NOT use "from numpy import *" while running this program.  It will mess it up and give you a headache trying to figure out what is going on... trust me.
###################
print 'Warning:\nBe sure that you have not imported the entire numpy module before running this program.\n'
ThetaMax = pi



class SETJ:
	
	def __init__(self, param): #, neighbors):	# param = [a,b,gamma,theta_0] regardless of ind_var % a --> u == bias voltage at node
		self.typename = 'generic'
	
	def __repr__(self):
        	return 'SETJ %s cell' % self.typename	# with parameters [a,b,gamma,theta_0] = %s' % (self.typename, `self.param`)
	
	def __mul__(self, SomeTuple):	# SomeTuple should have 2 members like (M,N)
		"""
		eg. SETJinner*(2,3) == [[SETJinner, SETJinner, SETJinner],[SETJinner, SETJinner, SETJinner]]
		"""
		rows = int(SomeTuple[0])	# avoid using anything other than integers
		cols = int(SomeTuple[1])	# avoid using anything other than integers
		SETJ_list = [ ]
		for col in xrange(cols):
			SETJ_list.append(self)	
		SETJ_list = array(rows*[SETJ_list]) # eg. SETJinner*(2,3) == array([[SETJinner, SETJinner, SETJinner],[SETJinner, SETJinner, SETJinner]]) 	
		return SETJ_list

# below is a reminder of the indexing for the array

##	j -->
#   i	
#
#   |
#   |	
#   v


# SETJinnner.evolve will get this kind of call:
# self.k2p[i][j], self.k2q[i][j], self.k2u[i][j], self.k2x[i][j] = self.array[i][j].evolve(i,j,p,q,u,x,t)

class SETJisolated(SETJ):
	def __init__(self, param):
        	SETJ.__init__(self, param)
        	self.typename = 'isolated'
	def dxdt(self,i,j,p,q,u,x,t):
		return gamma * (u[i][j] - b*cos(t) - x[i][j]/(k*pi))
	def evolve(self,i,j,p,q,u,x,t):
		p_current = 0.	#p[i][j]	#self.dpdt(i,j,p,q,u,x,t)	
		q_current = 0.	#q[i][j]	#self.dqdt(i,j,p,q,u,x,t)	
		u_current = 0.	#u[i][j]	# keep u constant for now	
		x_current = self.dxdt(i,j,p,q,u,x,t)		
		return p_current, q_current, u_current, x_current
		
class SETJendL(SETJ):	# left end of a one-D lattice
	def __init__(self, param):
        	SETJ.__init__(self, param)
        	self.typename = '1D--left end'
	def dxdt(self,i,j,p,q,u,x,t):
		return gamma * (u[i][j] - b*cos(t) - (k+1)*x[i][j]/(k*pi) + x[i][j+1]/(k*pi) + q[i][j+1]/(k*pi))
	def evolve(self,i,j,p,q,u,x,t):
		p_current = 0.	#p[i][j]	#self.dpdt(i,j,p,q,u,x,t)	
		q_current = 0.	#q[i][j]	#self.dqdt(i,j,p,q,u,x,t)	
		u_current = 0.	#u[i][j]	# keep u constant for now	
		x_current = self.dxdt(i,j,p,q,u,x,t)		
		return p_current, q_current, u_current, x_current

class SETJendR(SETJ):	# right end of a one-D lattice	
	def __init__(self, param):
        	SETJ.__init__(self, param)
        	self.typename = '1D--right end'
	def dxdt(self,i,j,p,q,u,x,t):
		return gamma * (u[i][j] - b*cos(t) - (k+1)*x[i][j]/(k*pi) + x[i][j-1]/(k*pi) + q[i][j]/(k*pi))
	def dqdt(self,i,j,p,q,u,x,t):
		return (gamma / (k*epsilon*pi)) * (x[i][j-1] - x[i][j] - q[i][j])
	def evolve(self,i,j,p,q,u,x,t):	
		p_current = 0.	#self.dpdt(i,j,p,q,u,x,t)	#SETJ_P_integrator(i,j,p,q,u,x,t, self.dpdt)
		q_current = self.dqdt(i,j,p,q,u,x,t)	#SETJ_Q_integrator(i,j,p,q,u,x,t, self.dqdt)
		u_current = 0.	#u[i][j]	# keep u constant for now	#SETJ_U_integrator(i,j,p,q,u,x,t, self.dudt)
		x_current = self.dxdt(i,j,p,q,u,x,t)	#SETJ_X_integrator(i,j,p,q,u,x,t, self.dxdt)	
		return p_current, q_current, u_current, x_current
		
class SETJinnerLR(SETJ):	# inner cell of a horizontal 1-D lattice
	def __init__(self, param):
        	SETJ.__init__(self, param)
        	self.typename = 'horizontal 1D--inner'
	def dxdt(self,i,j,p,q,u,x,t):
		return gamma * (u[i][j] - b*cos(t) - (k+2)*x[i][j]/(k*pi) + x[i][j-1]/(k*pi) + x[i][j+1]/(k*pi) - q[i][j]/(k*pi) + q[i][j+1]/(k*pi))
	def dqdt(self,i,j,p,q,u,x,t):
		return (gamma / (k*epsilon*pi)) * (x[i][j-1] - x[i][j] - q[i][j])
	def evolve(self,i,j,p,q,u,x,t):	
		p_current = 0.	#self.dpdt(i,j,p,q,u,x,t)	#SETJ_P_integrator(i,j,p,q,u,x,t, self.dpdt)
		q_current = self.dqdt(i,j,p,q,u,x,t)	#SETJ_Q_integrator(i,j,p,q,u,x,t, self.dqdt)
		u_current = 0.	#u[i][j]	# keep u constant for now	#SETJ_U_integrator(i,j,p,q,u,x,t, self.dudt)
		x_current = self.dxdt(i,j,p,q,u,x,t)	#SETJ_X_integrator(i,j,p,q,u,x,t, self.dxdt)	
		return p_current, q_current, u_current, x_current	
		
class SETJendT(SETJ):	# top end of a 1D lattice
	def __init__(self, param):
        	SETJ.__init__(self, param)
        	self.typename = '1D--top'
	def dxdt(self,i,j,p,q,u,x,t):#param (=[a,b,gamma..., k?,t?]),t,k(in param?),u,x,p,q	# should t be an attribute? # could make k a fn. of i,j for more interesting network
		return gamma * (u[i][j] - b*cos(t) - (k+1)*x[i][j]/(k*pi) + x[i-1][j]/(k*pi) + x[i+1][j]/(k*pi) +  p[i+1][j]/(k*pi))
	def evolve(self,i,j,p,q,u,x,t):
		p_current = 0.	#self.dpdt(i,j,p,q,u,x,t)	#SETJ_P_integrator(i,j,p,q,u,x,t, self.dpdt)
		q_current = 0.	#self.dqdt(i,j,p,q,u,x,t)	#SETJ_Q_integrator(i,j,p,q,u,x,t, self.dqdt)
		u_current = 0.	#u[i][j]	# keep u constant for now	#SETJ_U_integrator(i,j,p,q,u,x,t, self.dudt)
		x_current = self.dxdt(i,j,p,q,u,x,t)	#SETJ_X_integrator(i,j,p,q,u,x,t, self.dxdt)	
		return p_current, q_current, u_current, x_current
		
class SETJendB(SETJ):	# bottom end of a 1D lattice
	def __init__(self, param):
        	SETJ.__init__(self, param)
        	self.typename = '1D--bottom'
	def dxdt(self,i,j,p,q,u,x,t):#param (=[a,b,gamma..., k?,t?]),t,k(in param?),u,x,p,q	# should t be an attribute? # could make k a fn. of i,j for more interesting network
		return gamma * (u[i][j] - b*cos(t) - (k+1)*x[i][j]/(k*pi) + x[i-1][j]/(k*pi) + x[i+1][j]/(k*pi) - p[i][j]/(k*pi))
	def dpdt(self,i,j,p,q,u,x,t):
		return (gamma / (k*epsilon*pi)) * (x[i-1][j] - x[i][j] - p[i][j])
	def evolve(self,i,j,p,q,u,x,t):
		p_current = self.dpdt(i,j,p,q,u,x,t)	#SETJ_P_integrator(i,j,p,q,u,x,t, self.dpdt)
		q_current = 0.	#self.dqdt(i,j,p,q,u,x,t)	#SETJ_Q_integrator(i,j,p,q,u,x,t, self.dqdt)
		u_current = 0.	#u[i][j]	# keep u constant for now	#SETJ_U_integrator(i,j,p,q,u,x,t, self.dudt)
		x_current = self.dxdt(i,j,p,q,u,x,t)	#SETJ_X_integrator(i,j,p,q,u,x,t, self.dxdt)	
		return p_current, q_current, u_current, x_current


class SETJinnerTB(SETJ):	# inner cell of a vertical 1-D lattice
	def __init__(self, param):
        	SETJ.__init__(self, param)
        	self.typename = 'vertical 1D--inner'
	def dxdt(self,i,j,p,q,u,x,t):#param (=[a,b,gamma..., k?,t?]),t,k(in param?),u,x,p,q	# should t be an attribute? # could make k a fn. of i,j for more interesting network
		return gamma * (u[i][j] - b*cos(t) - (k+2)*x[i][j]/(k*pi) + x[i-1][j]/(k*pi) + x[i+1][j]/(k*pi) - p[i][j]/(k*pi) +  p[i+1][j]/(k*pi))
	def dpdt(self,i,j,p,q,u,x,t):
		return (gamma / (k*epsilon*pi)) * (x[i-1][j] - x[i][j] - p[i][j])
	def evolve(self,i,j,p,q,u,x,t):
		p_current = self.dpdt(i,j,p,q,u,x,t)	#SETJ_P_integrator(i,j,p,q,u,x,t, self.dpdt)
		q_current = 0.	#self.dqdt(i,j,p,q,u,x,t)	#SETJ_Q_integrator(i,j,p,q,u,x,t, self.dqdt)
		u_current = 0.	#u[i][j]	# keep u constant for now	#SETJ_U_integrator(i,j,p,q,u,x,t, self.dudt)
		x_current = self.dxdt(i,j,p,q,u,x,t)	#SETJ_X_integrator(i,j,p,q,u,x,t, self.dxdt)	
		return p_current, q_current, u_current, x_current
#add elbows, etc. for complete set
class SETJinner(SETJ):		# inner cell of a 2-D lattice
	def __init__(self, param):
        	SETJ.__init__(self, param)
        	self.typename = 'inner'
	def dxdt(self,i,j,p,q,u,x,t):#param (=[a,b,gamma..., k?,t?]),t,k(in param?),u,x,p,q	# should t be an attribute? # could make k a fn. of i,j for more interesting network
		return gamma * (u[i][j] - b*cos(t) - (k+4)*x[i][j]/(k*pi) + x[i][j-1]/(k*pi) + x[i][j+1]/(k*pi) + x[i-1][j]/(k*pi) + x[i+1][j]/(k*pi) - q[i][j]/(k*pi) + q[i][j+1]/(k*pi) - p[i][j]/(k*pi) +  p[i+1][j]/(k*pi))
	def dpdt(self,i,j,p,q,u,x,t):
		return (gamma / (k*epsilon*pi)) * (x[i-1][j] - x[i][j] - p[i][j])
	def dqdt(self,i,j,p,q,u,x,t):
		return (gamma / (k*epsilon*pi)) * (x[i][j-1] - x[i][j] - q[i][j])	
	def evolve(self,i,j,p,q,u,x,t):
		p_current = self.dpdt(i,j,p,q,u,x,t)	#SETJ_P_integrator(i,j,p,q,u,x,t, self.dpdt)
		q_current = self.dqdt(i,j,p,q,u,x,t)	#SETJ_Q_integrator(i,j,p,q,u,x,t, self.dqdt)
		u_current = 0.	#u[i][j]	# keep u constant for now	#SETJ_U_integrator(i,j,p,q,u,x,t, self.dudt)
		x_current = self.dxdt(i,j,p,q,u,x,t)	#SETJ_X_integrator(i,j,p,q,u,x,t, self.dxdt)	
		return p_current, q_current, u_current, x_current
		# note that y_current is the tunneling phase output of the SETJ 
		#if binary == 1:		# binary should be global in top program 
		#	if mod(tau_i,2*pi) > pi:
		#		y_current = 0.999	# use 0.999 instead of 1.0 to make sure this is not cut to 0.0 later in the mod call for colormap
		#	# if...
		#	else:
		#		y_current = 0.0
		#else:
		#	y_current = mod(self.tau_i,2*pi)/(2*pi)		# --> y has a gray scale value from [0, 1) ;	
		

class SETJ_first_col(SETJ):
	
	def __init__(self, param):
        	SETJ.__init__(self, param)
        	self.typename = 'leftmost column'
	def dxdt(self,i,j,p,q,u,x,t):
		return gamma * (u[i][j] - b*cos(t) - (k+3)*x[i][j]/(k*pi) + x[i][j+1]/(k*pi) + x[i-1][j]/(k*pi) + x[i+1][j]/(k*pi) + q[i][j+1]/(k*pi) - p[i][j]/(k*pi) +  p[i+1][j]/(k*pi))
	def dpdt(self,i,j,p,q,u,x,t):
		return (gamma / (k*epsilon*pi)) * (x[i-1][j] - x[i][j] - p[i][j])	
	def evolve(self,i,j,p,q,u,x,t):
		p_current = self.dpdt(i,j,p,q,u,x,t)	
		q_current = 0.	#q[i][j]	#self.dqdt(i,j,p,q,u,x,t)	
		u_current = 0.	#u[i][j]	# keep u constant for now	
		x_current = self.dxdt(i,j,p,q,u,x,t)	
		return p_current, q_current, u_current, x_current
	# do we need a dummy function for dqdt? -- should be solved by evolve setting q_current = q ( = 0.0 forever...)
	

class SETJ_last_col(SETJ):
	
	def __init__(self, param):
        	SETJ.__init__(self, param)
        	self.typename = 'rightmost column'
	def dxdt(self,i,j,p,q,u,x,t):
		return gamma * (u[i][j] - b*cos(t) - (k+3)*x[i][j]/(k*pi) + x[i][j-1]/(k*pi) + x[i-1][j]/(k*pi) + x[i+1][j]/(k*pi) - q[i][j]/(k*pi) - p[i][j]/(k*pi) +  p[i+1][j]/(k*pi))
	def dpdt(self,i,j,p,q,u,x,t):
		return (gamma / (k*epsilon*pi)) * (x[i-1][j] - x[i][j] - p[i][j])
	def dqdt(self,i,j,p,q,u,x,t):
		return (gamma / (k*epsilon*pi)) * (x[i][j-1] - x[i][j] - q[i][j])
	def evolve(self,i,j,p,q,u,x,t):
		p_current = self.dpdt(i,j,p,q,u,x,t)	
		q_current = self.dqdt(i,j,p,q,u,x,t)	
		u_current = 0.	#u[i][j]	# keep u constant for now	
		x_current = self.dxdt(i,j,p,q,u,x,t)		
		return p_current, q_current, u_current, x_current	


class SETJ_first_row(SETJ):
	
	def __init__(self, param):
        	SETJ.__init__(self, param)
        	self.typename = 'top row'
	def dxdt(self,i,j,p,q,u,x,t):
		return gamma * (u[i][j] - b*cos(t) - (k+3)*x[i][j]/(k*pi) + x[i][j-1]/(k*pi) + x[i][j+1]/(k*pi) + x[i+1][j]/(k*pi) - q[i][j]/(k*pi) + q[i][j+1]/(k*pi) +  p[i+1][j]/(k*pi))	
	# dummy function for dpdt?
	def dqdt(self,i,j,p,q,u,x,t):
		return (gamma / (k*epsilon*pi)) * (x[i][j-1] - x[i][j] - q[i][j])
	def evolve(self,i,j,p,q,u,x,t):
		p_current = 0.	#p[i][j]	#self.dpdt(i,j,p,q,u,x,t)	
		q_current = self.dqdt(i,j,p,q,u,x,t)	
		u_current = 0.	#u[i][j]	# keep u constant for now	
		x_current = self.dxdt(i,j,p,q,u,x,t)		
		return p_current, q_current, u_current, x_current
		
		
class SETJ_last_row(SETJ):

	def __init__(self, param):
        	SETJ.__init__(self, param)
        	self.typename = 'bottom row'
	def dxdt(self,i,j,p,q,u,x,t):
		return gamma * (u[i][j] - b*cos(t) - (k+3)*x[i][j]/(k*pi) + x[i][j-1]/(k*pi) + x[i][j+1]/(k*pi) + x[i-1][j]/(k*pi) - q[i][j]/(k*pi) + q[i][j+1]/(k*pi) - p[i][j]/(k*pi))	
	def dpdt(self,i,j,p,q,u,x,t):
		return (gamma / (k*epsilon*pi)) * (x[i-1][j] - x[i][j] - p[i][j])
	def dqdt(self,i,j,p,q,u,x,t):
		return (gamma / (k*epsilon*pi)) * (x[i][j-1] - x[i][j] - q[i][j])
	def evolve(self,i,j,p,q,u,x,t):
		p_current = self.dpdt(i,j,p,q,u,x,t)	
		q_current = self.dqdt(i,j,p,q,u,x,t)	
		u_current = 0.	#u[i][j]	# keep u constant for now	
		x_current = self.dxdt(i,j,p,q,u,x,t)		
		return p_current, q_current, u_current, x_current


class SETJ_cornerTL(SETJ):

	def __init__(self, param):
        	SETJ.__init__(self, param)
        	self.typename = 'top left corner'
	def dxdt(self,i,j,p,q,u,x,t):
		return gamma * (u[i][j] - b*cos(t) - (k+2)*x[i][j]/(k*pi) + x[i][j+1]/(k*pi) + x[i+1][j]/(k*pi) + q[i][j+1]/(k*pi) + p[i+1][j]/(k*pi))
	# dummy for dpdt?
	# dummy for dqdt?
	def evolve(self,i,j,p,q,u,x,t):
		p_current = 0.	#p[i][j]	#self.dpdt(i,j,p,q,u,x,t)	
		q_current = 0.	#q[i][j]	#self.dqdt(i,j,p,q,u,x,t)	
		u_current = 0.	#u[i][j]	# keep u constant for now	
		x_current = self.dxdt(i,j,p,q,u,x,t)		
		return p_current, q_current, u_current, x_current
	
	
class SETJ_cornerTR(SETJ):
	
	def __init__(self, param):
        	SETJ.__init__(self, param)
        	self.typename = 'top right corner'
	def dxdt(self,i,j,p,q,u,x,t):
		return gamma * (u[i][j] - b*cos(t) - (k+2)*x[i][j]/(k*pi) + x[i][j-1]/(k*pi) + x[i+1][j]/(k*pi) - q[i][j]/(k*pi) + p[i+1][j]/(k*pi))
	# dummy for dpdt?
	def dqdt(self,i,j,p,q,u,x,t):
		return (gamma / (k*epsilon*pi)) * (x[i][j-1] - x[i][j] - q[i][j])
	def evolve(self,i,j,p,q,u,x,t):
		p_current = 0.	#p[i][j]	#self.dpdt(i,j,p,q,u,x,t)	
		q_current = self.dqdt(i,j,p,q,u,x,t)	
		u_current = 0.	#u[i][j]	# keep u constant for now	
		x_current = self.dxdt(i,j,p,q,u,x,t)		
		return p_current, q_current, u_current, x_current
		
		
class SETJ_cornerBL(SETJ):

	def __init__(self, param):
        	SETJ.__init__(self, param)
        	self.typename = 'bottom left corner'
	def dxdt(self,i,j,p,q,u,x,t):
		return gamma * (u[i][j] - b*cos(t) - (k+2)*x[i][j]/(k*pi) + x[i][j+1]/(k*pi) + x[i-1][j]/(k*pi) + q[i][j+1]/(k*pi) - p[i][j]/(k*pi))
	def dpdt(self,i,j,p,q,u,x,t):
		return (gamma / (k*epsilon*pi)) * (x[i-1][j] - x[i][j] - p[i][j])
	# dummy for dqdt?
	def evolve(self,i,j,p,q,u,x,t):
		p_current = self.dpdt(i,j,p,q,u,x,t)	
		q_current = 0.	#q[i][j]	#self.dqdt(i,j,p,q,u,x,t)	
		u_current = 0.	#u[i][j]	# keep u constant for now	
		x_current = self.dxdt(i,j,p,q,u,x,t)		
		return p_current, q_current, u_current, x_current
		
		
class SETJ_cornerBR(SETJ):

	def __init__(self, param):
        	SETJ.__init__(self, param)
        	self.typename = 'bottom right corner'
	def dxdt(self,i,j,p,q,u,x,t):
		return gamma * (u[i][j] - b*cos(t) - (k+2)*x[i][j]/(k*pi) + x[i][j-1]/(k*pi) + x[i-1][j]/(k*pi) - q[i][j]/(k*pi) - p[i][j]/(k*pi))
	def dpdt(self,i,j,p,q,u,x,t):
		return (gamma / (k*epsilon*pi)) * (x[i-1][j] - x[i][j] - p[i][j])
	def dqdt(self,i,j,p,q,u,x,t):
		return (gamma / (k*epsilon*pi)) * (x[i][j-1] - x[i][j] - q[i][j])	
	def evolve(self,i,j,p,q,u,x,t):
		p_current = self.dpdt(i,j,p,q,u,x,t)	
		q_current = self.dqdt(i,j,p,q,u,x,t)	
		u_current = 0.	#u[i][j]	# keep u constant for now	
		x_current = self.dxdt(i,j,p,q,u,x,t)		
		return p_current, q_current, u_current, x_current

#***********************************************#


class SETJarray:
	
	def __init__(self,M,N,param):		# this is the only call that will need to be made, although individual params can be varied afterward
		self.array = SETJinner(param)*(M,N)	# returns MxN array of SETJ cells (they are incorrectly all classified as 'inner' at this point)
		# change edge SETJ types from generic
		self.array[:,0] = SETJ_first_col(param)*(1,M)	# replace leftmost column with SETJ_first_col cells
		self.array[:,-1] = SETJ_last_col(param)*(1,M)	# replace rightmost column with SETJ_last_col cells
		self.array[0,:] = SETJ_first_row(param)*(1,N)	# replace top row with SETJ_first_row cells
		self.array[-1,:] = SETJ_last_row(param)*(1,N)	# replace bottom row with SETJ_last_row cells
		# change corner SETJ types
		self.array[0,0] = SETJ_cornerTL(param)		# replace top left corner with SETJ_cornerTL cell; note that this is an instance and not an array!
		self.array[0,-1] = SETJ_cornerTR(param)		# replace top right corner with SETJ_cornerTR cell; note that this is an instance and not an array!
		self.array[-1,0] = SETJ_cornerBL(param)		# replace bottom left corner with SETJ_cornerBL cell; note that this is an instance and not an array!
		self.array[-1,-1] = SETJ_cornerBR(param)	# replace bottom right corner with SETJ_cornerBR cell; note that this is an instance and not an array!
		# check for 0 or 1-D
		if M == 1:
			self.array = SETJinnerLR(param)*(M,N)
			self.array[:,0] = SETJendL(param)*(1,M)
			self.array[:,-1] = SETJendR(param)*(1,M)
		if N == 1:
			self.array = SETJinnerTB(param)*(M,N)
			self.array[0,:] = SETJendT(param)*(1,N)
			self.array[-1,:] = SETJendB(param)*(1,N)	
		if M*N < 2:
			self.array[:,:] = SETJisolated(param)*(1,1)
		#
		self.shape = (M,N)
		self.t = 0.0	# the array has its own reference time, so it is good to go for SR
		#
		self.binary = param.binary	#1
		self.randomIC = param.randomIC	#0
		self.b = param.b		#10.
		self.gamma = param.gamma	#0.1
		self.epsilon = param.epsilon	#0.1
		self.k = param.k		#10.
		self.dt = param.dt		#0.01
		global binary, randomIC, b, gamma, epsilon, k, dt	# these variables must be available globally in this module
		binary = param.binary		#1
		randomIC = param.randomIC	#0
		b = param.b			#10.
		gamma = param.gamma		#0.1
		epsilon = param.epsilon		#0.1
		k = param.k			#10.
		dt = param.dt			#0.01
		#
		self.u_current = ones((M,N),float)	# start with ones for ease of manipulation in CNN_lattice_simulator #float(param[0])*ones((M,N),float)	# initialize u as constant MxN array with value 'a'
		#self.u_last = self.u_current
		#consider putting u_distr() here
		self.x_current = self.u_current	#IC for x	# x_current takes on values from [-pi,pi]
		if randomIC == 1:
			for i in xrange(M):
				for j in xrange(N):
					self.x_current[i][j] = random()*(2*pi) - pi
		#self.x_last = self.x_current
		self.p_current = zeros((M,N),float) # however, only zeros((M-1,N),float) will be used, so this would be more appropriate if it were easier to implement in RK4
		#self.p_last = self.p_current
		self.q_current = zeros((M,N),float) # only zeros((M,N-1),float) used for array;	This is propaply fine though b/c these unused p's and q's will just stay at 0
		#self.q_last = self.q_current
		self.y_current = zeros((M,N),float)	# phase output currently takes on values from [0,1)
		self.tau_i = zeros((M,N),float)
		#
		#Runge-Kutta initialization, in alphabetical order
		self.k1p = zeros((M,N),float)	#self.p_current # replaces self.p_last
		self.k1q = zeros((M,N),float)	#self.q_current # replaces self.q_last
		self.k1u = zeros((M,N),float)	#self.u_current #
		self.k1x = zeros((M,N),float)	#self.x_current #
		#
		self.k2p = zeros((M,N),float)	#self.p_current 
		self.k2q = zeros((M,N),float)	#self.q_current 
		self.k2u = zeros((M,N),float)	#self.u_current 
		self.k2x = zeros((M,N),float)	#self.x_current # could be faster to create these all with zeros((_,_),float)
		#
		self.k3p = zeros((M,N),float)	#self.p_current 
		self.k3q = zeros((M,N),float)	#self.q_current 
		self.k3u = zeros((M,N),float)	#self.u_current 
		self.k3x = zeros((M,N),float)	#self.x_current
		#
		self.k4p = zeros((M,N),float)	#self.p_current 
		self.k4q = zeros((M,N),float)	#self.q_current 
		self.k4u = zeros((M,N),float)	#self.u_current 
		self.k4x = zeros((M,N),float)	#self.x_current

	def new_param(self,M,N,param):
		global binary, randomIC, b, gamma, epsilon, k, dt	# these variables must be available globally in this module
		binary = param.binary		#1
		randomIC = param.randomIC	#0
		b = param.b			#10.
		gamma = param.gamma		#0.1
		epsilon = param.epsilon		#0.1
		k = param.k			#10.
		dt = param.dt			#0.01
		#
		if randomIC == 1:
			for i in xrange(M):
				for j in xrange(N):
					self.x_current[i][j] = random()*(2*pi) - pi

		
	def evolvearray(self,M,N):	##***** t from self.t
		"""
		evolvearray(self,M,N)
		Evolves every variable of the M by N array (potentially a lot!) by a single timestep, using RK4
		"""
		# k1
		t = self.t
		for i in xrange(M):
			for j in xrange(N):
				self.k1p[i][j], self.k1q[i][j], self.k1u[i][j], self.k1x[i][j] = (dt * array(self.array[i][j].evolve(i,j,self.p_current,self.q_current,self.u_current,self.x_current,t))).tolist()
		# k2
		## change Runge-Kutta y & t outside of loop to save computation time
		t = t + .5*dt	# evolve t for k2 and k3;	# assuming dt is given in top module
		p = self.p_current + .5*self.k1p	# note that we are adding two arrays of the same size, element by element
		q = self.q_current + .5*self.k1q
		u = self.u_current + .5*self.k1u
		x = self.x_current + .5*self.k1x
		for i in xrange(M):
			for j in xrange(N):
				self.k2p[i][j], self.k2q[i][j], self.k2u[i][j], self.k2x[i][j] = (dt * array(self.array[i][j].evolve(i,j,p,q,u,x,t))).tolist()
		# k3		
		# t is still good for k3
		p = self.p_current + .5*self.k2p
		q = self.q_current + .5*self.k2q
		u = self.u_current + .5*self.k2u
		x = self.x_current + .5*self.k2x
		for i in xrange(M):
			for j in xrange(N):
				self.k3p[i][j], self.k3q[i][j], self.k3u[i][j], self.k3x[i][j] = (dt * array(self.array[i][j].evolve(i,j,p,q,u,x,t))).tolist()
		# k4
		t += .5*dt	# We are now using t = t0 + dt		
		p = self.p_current + self.k3p
		q = self.q_current + self.k3q
		u = self.u_current + self.k3u
		x = self.x_current + self.k3x
		for i in xrange(M):
			for j in xrange(N):
				self.k4p[i][j], self.k4q[i][j], self.k4u[i][j], self.k4x[i][j] = (dt * array(self.array[i][j].evolve(i,j,p,q,u,x,t))).tolist()
		# Now find y_(n+1) = y_n + (k1 + 2k2 + 2k3 + k4)/6.;	need to verify that this is correct: is another dt needed?
		# update p,q,u,x current values
		self.p_current = self.p_current + (self.k1p + 2.0*self.k2p + 2.0*self.k3p + self.k4p)/6.0
		self.q_current = self.q_current + (self.k1q + 2.0*self.k2q + 2.0*self.k3q + self.k4q)/6.0
		self.u_current = self.u_current + (self.k1u + 2.0*self.k2u + 2.0*self.k3u + self.k4u)/6.0
		self.x_current = self.x_current + (self.k1x + 2.0*self.k2x + 2.0*self.k3x + self.k4x)/6.0
		self.t = t	# note that t is really time normalized to the pump frequency
		# Check for tunneling events
		for i in xrange(M):
			for j in xrange(N):
				if abs(self.x_current[i][j]) > ThetaMax:
					if self.x_current[i][j] > ThetaMax:
						self.x_current[i][j] -= 2*pi
						self.tau_i[i][j] = self.t
						#self.array[i][j].tau_i = self.t
					if self.x_current[i][j] < -ThetaMax:	# Yang mentions this condition early on but does not use it for eq. 6 and on... is it necessary? appropriate? --perhaps exp(-c*theta/pi) always beats out theta*c(a-b*cos(t)) for neg. theta?
						self.x_current[i][j] += 2*pi	# probably this condition should never be needed
						self.tau_i[i][j] = self.t
						#self.array[i][j].tau_i = self.t
					#self.array[i][j].tau_i = self.t	# this holds the time of the last tunneling event					#.append(dt * step)# each entry holds the time of a tunneling event
				# Now is a good time to update phase output
				if binary == 1:		# binary should be global in top program 
					#if mod(self.array[i][j].tau_i,2*pi) > pi:
					if mod(self.tau_i[i][j],2*pi) > pi:
						self.y_current[i][j] = 0.999	# use 0.999 instead of 1.0 to make sure this is not cut to 0.0 later in the mod call for colormap
					else:
						self.y_current[i][j] = 0.0
				else:
					self.y_current[i][j] = mod(self.tau_i[i][j],2*pi)/(2*pi)		# --> y has a gray scale value from [0, 1) ;	
		
		# self.t = t	## Should we evolve self.t here or in main program??
		#***********
		
		# pass the SETJ the SETJ array info  ## the SETJ only uses the information from itself and nearest neighbors, so could possibly make program faster by only passing these?
				#p[i][j], q[i][j], u[i][j], x[i][j] = self.array[i][j].evolve(i,j,self.p,self.q,self.u,self.x)
				
		
		#firstline = [SETJ_cornerTL(param),] + (N-2)*SETJ_first_row(param) + [SETJ_cornerTR(param),]
		#self.array = array([SETJ_cornerTL(param),] + (N-2)*SETJ_first_row(param) + [SETJ_cornerTR(param),])
		#for line in xrange(2,M):	# index output is 2 - M-1, but we must not confuse the [0] index with the i0 = 1 index in actual programming matters
		#	nextline = [SETJ_first_row(param),] + (N-2)*SETJinner(param) + [SETJ_last_row(param),]
	
	# consider making def __repr__ return a plot that shows the arrangement of the array with different colors for different types