#  rophysteresis.py
#  
#  This script compute a Nonlinear PDE by calling XPPAUT, plot as many figures
#  are needed to be encoded into a movie, and play it --auxin.ode and auxin.set
#  must be provided.
#
#  Created by Victor Brena on 22/12/2011.
#

import os                         #For issuing commands to the OS.
import sys                        #For determining the Python version.

import matplotlib		  #Plot libraries
from matplotlib import rcParams   
from pylab import *
from numpy import *		  #Algebra library
import xppy			  #XPP Python interface library by Jakub Nowacki

rcParams['xtick.direction'] = 'out'
rcParams['ytick.direction'] = 'out'

# For plotting 3D graphs

#from mpl_toolkits.mplot3d import Axes3D
#from matplotlib.collections import PolyCollection
#from matplotlib.colors import colorConverter
#from matplotlib import *

# Print some useful info

print '\n Executing on', os.uname()
print '\n Python version', sys.version
print '\n matplotlib version', matplotlib.__version__

#Reaction parameter values

k1 = 0.01
b = 0.01
c = 0.1
r = 0.01
nu = 1.5

#Spatial parameter values

D1 = 0.1
D2 = 10.0
X = 50.0

#Varying parameter values

k2i = 0.0
k2f = 0.4
L1 = 40.0
L2 = 70.0

#-----

def ropxpprun(initialc,foutput):
	#load xppaut files
	ode_file = 'auxin.ode'
	set_file = 'auxin.set'
	
	#create temporal files
	xppy.createTmp(ode_file, set_file)
	
	if initialc == 0:
		#from null initial conditions data
		#asign parameter and initial condition data
		pars = [['par','b',b],
				['par','c',c],
				['par','k1',k1],
				['par','r',r],
				['par','du',D1],
				['par','dv',D2],
				['par','l1',L1],
				['par','l2',L2],
				['par','k2i',k2i],
				['par','k2f',k2f]]
			
		for i in range(201):
			pars.append(['init','U%i'%i,0])
			pars.append(['init','V%i'%i,0])
			
	elif initialc == 1: 
		#from xpp initial conditions data
		#asign parameter and initial condition data
		pars = [['par','b',b],
				['par','c',c],
				['par','k1',k1],
				['par','r',r],
				['par','du',D1],
				['par','dv',D2],
				['par','l1',L2],
				['par','l2',L1],
				['par','k2i',k2f],
				['par','k2f',k2i]]
				
		initialc = load('emergspikesup.npy')
		ics = initialc[:,1:803]
		icu = ics[:,0::2]
		icv = ics[:,1::2]
		icu0 = icu[40000,:]
		icv0 = icv[40000,:]
		
		for i in range(201):
			pars.append(['init','U%i'%i,icu0[i]])
			pars.append(['init','V%i'%i,icv0[i]])
			
	else:
		#from data already distributed (AUTO)
		#asign parameter and initial condition data
		pars = [['par','b',b],
				['par','c',c],
				['par','k1',k1],
				['par','r',r],
				['par','du',D1],
				['par','dv',D2],
				['par','l1',L1],
				['par','l2',L2],
				['par','k2i',k2i],
				['par','k2f',k2f]]
				
		icu0 = initialc[0::2]
		icv0 = initialc[1::2]
		
		for i in range(201):
			pars.append(['init','U%i'%i,icu0[i]])
			pars.append(['init','V%i'%i,icv0[i]])

	#change ode.set file
	xppy.changeSet(pars)
	
	#run xppaut
	print '\n Running xppaut...'
		
	o = xppy.run()

	print 'Done'

	#choose solutions of interest by columns
	so = o[:,1:803]
	uuo = so[:,0::2]
	vvo = so[:,1::2]
	
	solution = [uuo,vvo]
	
	print '\n Renaming data file...'

	hh = loadtxt('output.dat')
	newnamefile = str(foutput)
	save(newnamefile,hh)
	os.remove('output.dat')
	
	print 'Done'
	
	return solution

#-----
	
def plotcreation(solfile,way):	# 'way' is either increasing/0 or decreasing/1
	#spare solution vector components
	uuo = solfile[0]
	vvo = solfile[1]
	
	Lmax = max(L1,L2)
	
	#define time and space
	ddd = shape(uuo)
	time = ddd[0]
	posit = ddd[1]
	
	tt = linspace(0,1,time,endpoint=True)
	xx = linspace(0,1,posit,endpoint=True)
	
	#gradient function
	alpha = lambda z: exp(-nu*z)

	#build sweeping k2 and L parameters
	if way == 0:
		k2t = lambda t: k2i*(1.0 - t) + k2f*t
		Lt = lambda t: L1*(1.0 - t) + L2*t
	else:
		k2t = lambda t: k2f*(1.0 - t) + k2i*t
		Lt = lambda t: L2*(1.0 - t) + L1*t

	for yy in range(time):
		xt = yy*200
		if xt <= 40000:
			tvalue = str('Time = %03d'%xt)
	
			gra1 = subplot(221)
			gra1.plot(tt,Lt(tt),'g',lw=2)
			gra1.plot(xt/40000.0,Lt(xt/40000.0),'og', markersize=7)
			gra1.set_title(tvalue)
			ylim(0,Lmax)
			setp(gra1.get_xticklabels(), visible=False)
			ygra1 = gra1.set_ylabel(r'$L$',color='g')
			for tk in gra1.get_yticklabels():
				tk.set_color('g')
	
			gra2 = gra1.twinx()
			gra2.plot(tt,k2t(tt),'m',lw=2)
			gra2.plot(xt/40000.0,k2t(xt/40000.0),'om', markersize=7)
			grid(True)
			ygra2 = gra2.set_ylabel(r'$k_{20}$',color='m')
			for tk in gra2.get_yticklabels():
				tk.set_color('m')
				
			suf = subplot(222,autoscale_on = True)
			ima = suf.imshow(uuo.T,aspect='auto',origin='lower',extent=(0,Lmax,0,30),vmax=10)
			figi = gcf()
			figi.colorbar(ima)
			setp(suf.get_xticklabels(),visible=False)
			setp(suf.get_yticklabels(),visible=False)

			spi1 = subplot(212)
			spi1.plot(xx,uuo[xt,:],'b',xx,vvo[xt,:],'r',lw =2)
			grid(True)
			spi1.fill_between(xx,0,uuo[xt,:],color='b',alpha=0.2)
			spi1.fill_between(xx,0,vvo[xt,:],color='r',alpha=0.2)
			xspi0 = xlabel(r'$x$')
			yspi0 = ylabel('$U,V$')
			ylim(0, 15)
			xspi0.set_fontsize('large')
			yspi0.set_fontsize('large')
			for tk in spi1.get_yticklabels():
				tk.set_color('k')
				
			spi2 = spi1.twinx()
			spi2.plot(xx,alpha(xx),'--m')
			yspi2 = spi2.set_ylabel('Gradient',color='m')
			for tk in spi2.get_yticklabels():
				tk.set_color('m')
			
			if way == 0:
				fname = str('spikesframe%03d' % yy) + '.png'
				savefig(fname,dp1 = 300)
	
				clf()
	
				print 'Wrote file', fname
			
			else:
				zz = yy + 201
				fname = str('spikesframe%03d' % zz) + '.png'
				savefig(fname,dp1 = 300)
	
				clf()
	
				print 'Wrote file', fname
	print 'Done'

#-----

def encoding(filmname):
	command = ('mencoder',
		   '-mc',
		   '0',
		   '-noskip',
		   '-skiplimit',
		   '0',
		   'mf://*.png',
		   '-mf',
		   'type=png:w=1000:h=800:fps=25',
		   '-ovc',
		   'lavc',
		   '-lavcopts',
		   'vcodec=mpeg4:vhq:trell:mbd=2:vmax_b_frames=1:v4mv:vb_strategy=0:vlelim=0:vcelim=0:cmp=6:subcmp=6:precmp=6:predia=3:dia=3:vme=4:vqscale=1',
		   '-oac',
		   'copy',
		   '-o',
		    filmname)

	os.spawnvp(os.P_WAIT, 'mencoder', command)

	print 'Done'
	
#-----

def playfilm(filmname):
	command = ('mplayer',
		   '-loop',
		   '0',
		   filmname)
	   
	os.spawnvp(os.P_WAIT, 'mplayer', command)

#-----

#run and plot increasing way
print '\n Computing for increasing k_20 and L parameters...'

up = ropxpprun(0,'emergspikesup.npy')
plotcreation(up,0) 

#run and plot decreasing way
print '\n Computing for decreasing k_20 and L parameters...'

down = ropxpprun(1,'emergspikesdown.npy')
plotcreation(down,1)

#enconde all figures onto a movie
print '\n Encoding figures... \n'

encoding('emergspikes.avi')

#play the movie
print '\n Play an endless loop \n'

playfilm('emergspikes.avi')