from __future__ import print_function
import numpy as np
import sys
import os
import datetime
import matplotlib.pyplot as plt
from matplotlib.backends.backend_pdf import PdfPages
from matplotlib.ticker import MultipleLocator
from matplotlib.ticker import ScalarFormatter
from matplotlib.ticker import NullFormatter
import matplotlib.image as image
import matplotlib.patches as patches
import matplotlib.lines as lines
import matplotlib.ticker as ticker
from scipy.optimize import curve_fit

from Aux import *

def compute_results (N_times,dt,AvgPvalues):
	ProbPhasesDuration = np.zeros((3,N_times))
	normBPphase = 0.0
	normABPphase = 0.0
	avgBPphase = 0.0
	avgABPphase = 0.0
	ProbPhasesDuration[0] = AvgPvalues[0]
	ProbPhasesDuration[1][0] = AvgPvalues[1][0]/dt
	ProbPhasesDuration[2][0] = AvgPvalues[2][0]/dt
	avgBPphase+=ProbPhasesDuration[0][0]*ProbPhasesDuration[1][0]*dt
	avgABPphase+=ProbPhasesDuration[0][0]*ProbPhasesDuration[2][0]*dt
	normBPphase+=ProbPhasesDuration[1][0]*dt
	normABPphase+=ProbPhasesDuration[2][0]*dt
	for t in range(1,N_times):
		ProbPhasesDuration[1][t] = ProbPhasesDuration[1][t-1]/AvgPvalues[1][t-1]  *(1.-AvgPvalues[1][t-1])*AvgPvalues[1][t]
		ProbPhasesDuration[2][t] = ProbPhasesDuration[2][t-1]/AvgPvalues[2][t-1]  *(1.-AvgPvalues[2][t-1])*AvgPvalues[2][t]
		avgBPphase+=ProbPhasesDuration[0][t]*ProbPhasesDuration[1][t]*dt
		avgABPphase+=ProbPhasesDuration[0][t]*ProbPhasesDuration[2][t]*dt
		normBPphase+=ProbPhasesDuration[1][t]*dt
		normABPphase+=ProbPhasesDuration[2][t]*dt
	
	soglia = 0.01
	Comululative = np.zeros((3,N_times))
	Comululative[0] = AvgPvalues[0]
	Comululative[1][0] = 0.0
	Comululative[2][0] = 0.0
	sBP=0
	sABP=0
	logBP = 1
	logABP = 1
	for t in range(1,N_times):
		Comululative[1][t] = Comululative[1][t-1] + ProbPhasesDuration[1][t-1]*dt
		Comululative[2][t] = Comululative[2][t-1] + ProbPhasesDuration[2][t-1]*dt
		if (1.0-Comululative[1][t] < soglia and logBP==1): 
			sBP=t
			logBP=0
		if (1.0-Comululative[2][t] < soglia and logABP==1): 
			sABP=t
			logABP=0
			
	print (normBPphase,normABPphase)
	return [AvgPvalues,ProbPhasesDuration,Comululative,sBP,sABP,avgBPphase,avgABPphase]

def make_figure (Results_ini,Results_Pe1,Results_Pe2,Results_Pe5,Results_Pe10,Results_Pe20,Results_Pe50,Results_Pe100,Results_Pe200,Results_Pe500,Results_Pe1000,Results_Pe2000):
	# extract results
	Pes = np.asarray([1,2,5,10,20,50,100,200,500,1000,2000])
	avgBPphaseduration_ini = np.asarray([Results_ini[5] for Pe in Pes])
	avgABPphaseduration_ini = np.asarray([Results_ini[6] for Pe in Pes])
	avgBPphaseduration = np.asarray([Results_Pe1[5],Results_Pe2[5],Results_Pe5[5],Results_Pe10[5],Results_Pe20[5],Results_Pe50[5],Results_Pe100[5],Results_Pe200[5],Results_Pe500[5],Results_Pe1000[5],Results_Pe2000[5]])
	avgABPphaseduration = np.asarray([Results_Pe1[6],Results_Pe2[6],Results_Pe5[6],Results_Pe10[6],Results_Pe20[6],Results_Pe50[6],Results_Pe100[6],Results_Pe200[6],Results_Pe500[6],Results_Pe1000[6],Results_Pe2000[6]])
	
	
	print(Pes)
	print(avgBPphaseduration)
	print(avgABPphaseduration)
	
	# setto alcune variabili comuni
	axisticslabelfontsize=9
	axisticslabelfontsizeinset=7
	axislabelfontsize=11 
	axislabelfontsizeinset=9
	
	xfig,yfig=7.2,4.6
	factor = xfig/yfig
	
	xgap = 0.082
	ygap = 0.075
	dy_bottom_panel = 0.32
	dy_upper_panel = 0.42
	dx_3panels = (0.98-3*xgap)/3
	dx_2panels = (0.98-2*xgap)/2
	
	
	#~ xd1,yd1,dxd1,dyd1 = xgap, ygap, dx_2panels*0.6, dy_bottom_panel
	#~ xd2,yd2,dxd2,dyd2 = xd1+dxd1, ygap, dx_2panels*0.4, dy_bottom_panel
	
	#~ xe1,ye1,dxe1,dye1 = xd1+dx_2panels+xgap, ygap, dx_2panels*0.6, dy_bottom_panel/2
	#~ xe2,ye2,dxe2,dye2 = xe1, ye1+dye1, dxe1, dye1
	#~ xe3,ye3,dxe3,dye3 = xe1+dxe1, ye1, dx_2panels*0.4, dye1
	#~ xe4,ye4,dxe4,dye4 = xe3, ye3+dye3, dxe3, dye3
	
	xe1,ye1,dxe1,dye1 = xgap, ygap, dx_2panels*0.6, dy_bottom_panel/2
	xe2,ye2,dxe2,dye2 = xe1, ye1+dye1, dxe1, dye1
	xe3,ye3,dxe3,dye3 = xe1+dxe1, ye1, dx_2panels*0.4, dye1
	xe4,ye4,dxe4,dye4 = xe3, ye3+dye3, dxe3, dye3
	
	xe5,ye5,dxe5,dye5 = xgap, ygap, dx_2panels, dy_bottom_panel/2
	xe6,ye6,dxe6,dye6 = xe5, ye5+dye5, dxe5, dye5
	
	xd1,yd1,dxd1,dyd1 = xe1+dx_2panels+xgap, ygap, dx_2panels*0.6, dy_bottom_panel
	xd2,yd2,dxd2,dyd2 = xd1+dxd1, ygap, dx_2panels*0.4, dy_bottom_panel
	
	xd3,yd3,dxd3,dyd3 = xe1+dx_2panels+xgap, ygap, dx_2panels, dy_bottom_panel
	
	
	
	xa1,ya1,dxa1,dya1 = xgap, yd1+dyd1+ygap+0.04, dx_3panels, dy_upper_panel/2
	xa2,ya2,dxa2,dya2 = xa1, ya1+dya1, dxa1, dya1
	
	xb1,yb1,dxb1,dyb1 = xa1+dxa1+xgap, ya1, dxa1, dya1
	xb2,yb2,dxb2,dyb2 = xb1, yb1+dyb1, dxb1, dyb1
	
	xc1,yc1,dxc1,dyc1 = xb1+dxb1+xgap, ya1, dxa1, dya1
	xc2,yc2,dxc2,dyc2 = xc1, yc1+dyc1, dxc1, dyc1
	
	xalogscale = 1
	xalim_left,xalim_right = 0,0.04
	xatics_main,xatics_minor = 0.01,0.005
	xalimlog_left,xalimlog_right = 0.00007,1.0
	
	yalogscale = 1
	yblogscale = 1
	
	xdlim_left,xdlim_middle,xdlim_right = -1,50,1000
	ydlim_low,ydlim_up = 0.0, 1.0
	xd1tics_main,xd1tics_minor = 10,2
	xd2tics_main,xd2tics_minor = 200,50
	
	xelim_left,xelim_right = 0.8,3000
	yelim_low,yelim_up = 0.00006, 0.6
	
	
	
	
	
	log_epanel_Pe = 1
	log_dpanel_Pe = 1
	
	
	
	with PdfPages('fig2.pdf') as pdf:
		fig = plt.figure(figsize=(xfig,yfig))
		plt.rc('text', usetex=True)
		plt.rc('text.latex', preamble = ','.join('''\usepackage{txfonts} \usepackage{lmodern}'''.split()))
		
		
		
		######################## panel A1 ########################################
		panel = fig.add_axes([xa1,ya1,dxa1,dya1])
		panel.tick_params(axis='both',which='both',direction='in',bottom=True,top=True,left=True,right=True)
		for tick in panel.xaxis.get_major_ticks(): tick.label.set_fontsize(axisticslabelfontsize)
		for tick in panel.yaxis.get_major_ticks(): tick.label.set_fontsize(axisticslabelfontsize)
		panel.set_xlabel(r'$\omega / \tau$',fontsize=axislabelfontsize)
		panel.set_ylabel(r'$p(\phi \!=\! 0,\omega)$',fontsize=axislabelfontsize-1)
		if (xalogscale == 0):
			panel.set_xlim(xalim_left,xalim_right)
			panel.xaxis.set_major_locator(MultipleLocator(xatics_main))
			panel.xaxis.set_minor_locator(MultipleLocator(xatics_minor))
		else:
			panel.set_xscale('log')
			panel.set_xlim(xalimlog_left,xalimlog_right)
			panel.set_xticks([0.0001,0.001,0.01,0.1,1.0])
			x_minor = ticker.LogLocator(base=10.0,subs=(0.2,0.4,0.6,0.8),numticks=12)
			panel.xaxis.set_minor_locator(x_minor)
			panel.xaxis.set_minor_formatter(ticker.NullFormatter())
		
		if (yalogscale == 0):
			panel.set_ylim(-0.01, 0.49)
			panel.yaxis.set_major_locator(MultipleLocator(0.1))
			panel.yaxis.set_minor_locator(MultipleLocator(0.05))
		else:
			panel.set_yscale('log')
			panel.set_ylim(0.00005, 2.0)
			panel.set_yticks([0.0001,0.001,0.01,0.1,1.0])
			y_minor = ticker.LogLocator(base=10.0,subs=(0.2,0.4,0.6,0.8),numticks=12)
			panel.yaxis.set_minor_locator(y_minor)
			panel.yaxis.set_minor_formatter(ticker.NullFormatter())
		
		panel.text(0.77,0.68,r'BP',fontsize=axislabelfontsize+1,transform=panel.transAxes)
		
		# plot avg line
		#~ panel.plot(Results_Pe1[0][0],Results_Pe1[0][1],'b-',markersize=2,linewidth=1.5,label=r'Pe = $1$')
		panel.plot(Results_Pe2[0][0],Results_Pe2[0][1],'b-',markersize=2,linewidth=1.5,label=r'Pe = $2$')
		panel.plot(Results_Pe20[0][0],Results_Pe20[0][1],'r-',markersize=2,linewidth=1.5,label=r'Pe = $20$')
		panel.plot(Results_Pe100[0][0],Results_Pe100[0][1],'g-',markersize=2,linewidth=1.5,label=r'Pe = $100$')
		panel.plot(Results_ini[0][0],Results_ini[0][1],'k--',markersize=2,linewidth=1.5,label=r'initial policy')
		#~ panel.legend(loc='upper right', bbox_to_anchor=(0.99, 0.99),ncol=1,fontsize=8,handlelength=1.5,labelspacing=0.2)
		
		######################## panel A2 ########################################
		panel = fig.add_axes([xa2,ya2,dxa2,dya2])
		panel.tick_params(axis='both',which='both',direction='in',bottom=True,top=True,left=True,right=True)
		for tick in panel.xaxis.get_major_ticks(): tick.label.set_fontsize(axisticslabelfontsize)
		for tick in panel.yaxis.get_major_ticks(): tick.label.set_fontsize(axisticslabelfontsize)
		panel.set_ylabel(r'$p(\phi \!=\! 1,\omega)$',fontsize=axislabelfontsize-1)
		if (xalogscale == 0):
			panel.set_xticklabels([])
			panel.set_xlim(xalim_left,xalim_right)
			panel.xaxis.set_major_locator(MultipleLocator(xatics_main))
			panel.xaxis.set_minor_locator(MultipleLocator(xatics_minor))
		else:
			panel.set_xscale('log')
			panel.set_xticklabels([])
			panel.set_xlim(xalimlog_left,xalimlog_right)
			panel.set_xticks([0.0001,0.001,0.01,0.1,1.0])
			x_minor = ticker.LogLocator(base=10.0,subs=(0.2,0.4,0.6,0.8),numticks=12)
			panel.xaxis.set_minor_locator(x_minor)
			panel.xaxis.set_minor_formatter(ticker.NullFormatter())
			
		if (yalogscale == 0):
			panel.set_ylim(-0.01, 0.49)
			panel.yaxis.set_major_locator(MultipleLocator(0.1))
			panel.yaxis.set_minor_locator(MultipleLocator(0.05))
		else:
			panel.set_yscale('log')
			panel.set_ylim(0.00005, 2.0)
			panel.set_yticks([0.0001,0.001,0.01,0.1,1.0])
			y_minor = ticker.LogLocator(base=10.0,subs=(0.2,0.4,0.6,0.8),numticks=12)
			panel.yaxis.set_minor_locator(y_minor)
			panel.yaxis.set_minor_formatter(ticker.NullFormatter())
		
		panel.text(0.76,0.55,r'ABP',fontsize=axislabelfontsize+1,transform=panel.transAxes)
		panel.text(0.85,0.8,r'(a)',fontsize=axislabelfontsize,transform=panel.transAxes)
		# plot avg line
		#~ panel.plot(Results_Pe1[0][0],Results_Pe1[0][2],'b-',markersize=2,linewidth=1.5,label=r'Pe = $1$')
		panel.plot(Results_Pe2[0][0],Results_Pe2[0][2],'b-',markersize=2,linewidth=1.5,label=r'Pe = $2$')
		panel.plot(Results_Pe20[0][0],Results_Pe20[0][2],'r-',markersize=2,linewidth=1.5,label=r'Pe = $20$')
		panel.plot(Results_Pe100[0][0],Results_Pe100[0][2],'g-',markersize=2,linewidth=1.5,label=r'Pe = $100$')
		#~ panel.plot(Results_Pe100[0][0][:Results_Pe100[4]],Results_Pe100[0][2][:Results_Pe100[4]],'g-',markersize=2,linewidth=1.5,label=r'Pe = $100$')
		#~ panel.plot(Results_Pe100[0][0][Results_Pe100[4]-1:],Results_Pe100[0][2][Results_Pe100[4]-1:],'-',color='lightgreen', markersize=2,linewidth=1.0)
		panel.plot(Results_ini[0][0],Results_ini[0][2],'k--',markersize=2,linewidth=1.5,label=r'initial policy')
		
		#~ panel.legend(loc='upper right', bbox_to_anchor=(0.99, 0.99),ncol=1,fontsize=8,handlelength=1.5,labelspacing=0.2)
		
		
		######################## panel B1 ########################################
		panel = fig.add_axes([xb1,yb1,dxb1,dyb1])
		panel.tick_params(axis='both',which='both',direction='in',bottom=True,top=True,left=True,right=True)
		for tick in panel.xaxis.get_major_ticks(): tick.label.set_fontsize(axisticslabelfontsize)
		for tick in panel.yaxis.get_major_ticks(): tick.label.set_fontsize(axisticslabelfontsize)
		panel.set_xlabel(r'$\tilde{\omega} / \tau$',fontsize=axislabelfontsize)
		panel.set_ylabel(r'$P(\omega \!=\! \tilde{\omega}|\phi \!=\! 0)$',fontsize=axislabelfontsize-1)
		if (xalogscale == 0):
			panel.set_xlim(xalim_left,xalim_right)
			panel.xaxis.set_major_locator(MultipleLocator(xatics_main))
			panel.xaxis.set_minor_locator(MultipleLocator(xatics_minor))
		else:
			panel.set_xscale('log')
			panel.set_xlim(xalimlog_left,xalimlog_right)
			panel.set_xticks([0.0001,0.001,0.01,0.1,1.0])
			x_minor = ticker.LogLocator(base=10.0,subs=(0.2,0.4,0.6,0.8),numticks=12)
			panel.xaxis.set_minor_locator(x_minor)
			panel.xaxis.set_minor_formatter(ticker.NullFormatter())
			
		if (yblogscale == 0):
			panel.set_ylim(-100, 2400)
			panel.yaxis.set_major_locator(MultipleLocator(500))
			panel.yaxis.set_minor_locator(MultipleLocator(100))
		else:
			panel.set_yscale('log')
			panel.set_ylim(0.2, 8000.0)
			panel.set_yticks([1.0,10,100,1000])
			y_minor = ticker.LogLocator(base=10.0,subs=(0.2,0.4,0.6,0.8),numticks=12)
			panel.yaxis.set_minor_locator(y_minor)
			panel.yaxis.set_minor_formatter(ticker.NullFormatter())
			
		panel.text(0.77,0.68,r'BP',fontsize=axislabelfontsize+1,transform=panel.transAxes)
		
		# plot avg line
		#~ panel.plot(Results_Pe1[1][0],Results_Pe1[1][1],'b-',markersize=2,linewidth=1.5,label=r'Pe = $1$')
		panel.plot(Results_Pe2[1][0],Results_Pe2[1][1],'b-',markersize=2,linewidth=1.5,label=r'Pe = $2$')
		panel.plot(Results_Pe20[1][0],Results_Pe20[1][1],'r-',markersize=2,linewidth=1.5,label=r'Pe = $20$')
		panel.plot(Results_Pe100[1][0],Results_Pe100[1][1],'g-',markersize=2,linewidth=1.5,label=r'Pe = $100$')
		#~ panel.plot(Results_Pe100[1][0][:Results_Pe100[3]],Results_Pe100[1][1][:Results_Pe100[3]],'g-',markersize=2,linewidth=1.5,label=r'Pe = $100$')
		#~ panel.plot(Results_Pe100[1][0][Results_Pe100[3]-1:],Results_Pe100[1][1][Results_Pe100[3]-1:],'-',color='lightgreen', markersize=2,linewidth=1.0)
		panel.plot(Results_ini[1][0],Results_ini[1][1],'k--',markersize=2,linewidth=1.5,label=r'initial policy')
		
		#~ panel.legend(loc='upper right', bbox_to_anchor=(0.99, 0.99),ncol=1,fontsize=8,handlelength=1.5,labelspacing=0.2)
		
		
		######################## panel B2 ########################################
		panel = fig.add_axes([xb2,yb2,dxb2,dyb2])
		panel.tick_params(axis='both',which='both',direction='in',bottom=True,top=True,left=True,right=True)
		for tick in panel.xaxis.get_major_ticks(): tick.label.set_fontsize(axisticslabelfontsize)
		for tick in panel.yaxis.get_major_ticks(): tick.label.set_fontsize(axisticslabelfontsize)
		panel.set_ylabel(r'$P(\omega \!=\! \tilde{\omega}|\phi \!=\! 1)$',fontsize=axislabelfontsize-1)
		if (xalogscale == 0):
			panel.set_xticklabels([])
			panel.set_xlim(xalim_left,xalim_right)
			panel.xaxis.set_major_locator(MultipleLocator(xatics_main))
			panel.xaxis.set_minor_locator(MultipleLocator(xatics_minor))
		else:
			panel.set_xscale('log')
			panel.set_xticklabels([])
			panel.set_xlim(xalimlog_left,xalimlog_right)
			panel.set_xticks([0.0001,0.001,0.01,0.1,1.0])
			x_minor = ticker.LogLocator(base=10.0,subs=(0.2,0.4,0.6,0.8),numticks=12)
			panel.xaxis.set_minor_locator(x_minor)
			panel.xaxis.set_minor_formatter(ticker.NullFormatter())
			
		if (yblogscale == 0):
			panel.set_ylim(-100, 1700)
			panel.yaxis.set_major_locator(MultipleLocator(500))
			panel.yaxis.set_minor_locator(MultipleLocator(100))
		else:
			panel.set_yscale('log')
			panel.set_ylim(0.2, 8000.0)
			panel.set_yticks([1.0,10,100,1000])
			y_minor = ticker.LogLocator(base=10.0,subs=(0.2,0.4,0.6,0.8),numticks=12)
			panel.yaxis.set_minor_locator(y_minor)
			panel.yaxis.set_minor_formatter(ticker.NullFormatter())
		panel.text(0.85,0.8,r'(b)',fontsize=axislabelfontsize,transform=panel.transAxes)
		panel.text(0.76,0.55,r'ABP',fontsize=axislabelfontsize+1,transform=panel.transAxes)
		# plot avg line
		#~ panel.plot(Results_Pe1[1][0],Results_Pe1[1][2],'b-',markersize=2,linewidth=1.5,label=r'Pe = $1$')
		panel.plot(Results_Pe2[1][0],Results_Pe2[1][2],'b-',markersize=2,linewidth=1.5,label=r'Pe = $2$')
		panel.plot(Results_Pe20[1][0],Results_Pe20[1][2],'r-',markersize=2,linewidth=1.5,label=r'Pe = $20$')
		panel.plot(Results_Pe100[1][0],Results_Pe100[1][2],'g-',markersize=2,linewidth=1.5,label=r'Pe = $100$')
		#~ panel.plot(Results_Pe100[1][0][:Results_Pe100[4]],Results_Pe100[1][2][:Results_Pe100[4]],'g-',markersize=2,linewidth=1.5,label=r'Pe = $100$')
		#~ panel.plot(Results_Pe100[1][0][Results_Pe100[4]-1:],Results_Pe100[1][2][Results_Pe100[4]-1:],'-',color='lightgreen', markersize=2,linewidth=1.0)
		panel.plot(Results_ini[1][0],Results_ini[1][2],'k--',markersize=2,linewidth=1.5,label=r'initial policy')
		
		panel.legend(loc='center', bbox_to_anchor=(0.45, 1.2),ncol=4,fontsize=8,handlelength=1.5,labelspacing=0.2)
		
		
		######################## panel C1 ########################################
		panel = fig.add_axes([xc1,yc1,dxc1,dyc1])
		panel.tick_params(axis='both',which='both',direction='in',bottom=True,top=True,left=True,right=True)
		for tick in panel.xaxis.get_major_ticks(): tick.label.set_fontsize(axisticslabelfontsize)
		for tick in panel.yaxis.get_major_ticks(): tick.label.set_fontsize(axisticslabelfontsize)
		panel.set_xlabel(r'$\tilde{\omega} / \tau$',fontsize=axislabelfontsize)
		panel.set_ylabel(r'$P(\omega \!<\! \tilde{\omega}|\phi\!=\!0)$',fontsize=axislabelfontsize-1)
		if (xalogscale == 0):
			panel.set_xlim(xalim_left,xalim_right)
			panel.xaxis.set_major_locator(MultipleLocator(xatics_main))
			panel.xaxis.set_minor_locator(MultipleLocator(xatics_minor))
		else:
			panel.set_xscale('log')
			panel.set_xlim(xalimlog_left,xalimlog_right)
			panel.set_xticks([0.0001,0.001,0.01,0.1,1.0])
			x_minor = ticker.LogLocator(base=10.0,subs=(0.2,0.4,0.6,0.8),numticks=12)
			panel.xaxis.set_minor_locator(x_minor)
			panel.xaxis.set_minor_formatter(ticker.NullFormatter())
			
		panel.set_ylim(-0.07, 1.07)
		panel.yaxis.set_major_locator(MultipleLocator(0.2))
		panel.yaxis.set_minor_locator(MultipleLocator(0.1))
		
		panel.text(0.77,0.18,r'BP',fontsize=axislabelfontsize+1,transform=panel.transAxes)
		
		# plot avg line
		#~ panel.plot(Results_Pe1[2][0],Results_Pe1[2][1],'b-',markersize=2,linewidth=1.5,label=r'Pe = $1$')
		panel.plot(Results_Pe2[2][0],Results_Pe2[2][1],'b-',markersize=2,linewidth=1.5,label=r'Pe = $2$')
		panel.plot(Results_Pe20[2][0],Results_Pe20[2][1],'r-',markersize=2,linewidth=1.5,label=r'Pe = $20$')
		panel.plot(Results_Pe100[2][0],Results_Pe100[2][1],'g-',markersize=2,linewidth=1.5,label=r'Pe = $100$')
		panel.plot(Results_ini[2][0],Results_ini[2][1],'k--',markersize=2,linewidth=1.5,label=r'initial policy')
		
		#~ panel.legend(loc='upper right', bbox_to_anchor=(0.99, 0.99),ncol=1,fontsize=8,handlelength=1.5,labelspacing=0.2)
		
		
		######################## panel C2 ########################################
		panel = fig.add_axes([xc2,yc2,dxc2,dyc2])
		panel.tick_params(axis='both',which='both',direction='in',bottom=True,top=True,left=True,right=True)
		for tick in panel.xaxis.get_major_ticks(): tick.label.set_fontsize(axisticslabelfontsize)
		for tick in panel.yaxis.get_major_ticks(): tick.label.set_fontsize(axisticslabelfontsize)
		panel.set_ylabel(r'$P(\omega \!<\! \tilde{\omega}|\phi\!=\!1)$',fontsize=axislabelfontsize-1)
		if (xalogscale == 0):
			panel.set_xticklabels([])
			panel.set_xlim(xalim_left,xalim_right)
			panel.xaxis.set_major_locator(MultipleLocator(xatics_main))
			panel.xaxis.set_minor_locator(MultipleLocator(xatics_minor))
		else:
			panel.set_xscale('log')
			panel.set_xticklabels([])
			panel.set_xlim(xalimlog_left,xalimlog_right)
			panel.set_xticks([0.0001,0.001,0.01,0.1,1.0])
			x_minor = ticker.LogLocator(base=10.0,subs=(0.2,0.4,0.6,0.8),numticks=12)
			panel.xaxis.set_minor_locator(x_minor)
			panel.xaxis.set_minor_formatter(ticker.NullFormatter())
			
		panel.set_ylim(-0.07, 1.07)
		panel.yaxis.set_major_locator(MultipleLocator(0.2))
		panel.yaxis.set_minor_locator(MultipleLocator(0.1))
		panel.text(0.07,0.8,r'(c)',fontsize=axislabelfontsize,transform=panel.transAxes)
		panel.text(0.76,0.18,r'ABP',fontsize=axislabelfontsize+1,transform=panel.transAxes)
		# plot avg line
		#~ panel.plot(Results_Pe1[2][0],Results_Pe1[2][2],'b-',markersize=2,linewidth=1.5,label=r'Pe = $1$')
		panel.plot(Results_Pe2[2][0],Results_Pe2[2][2],'b-',markersize=2,linewidth=1.5,label=r'Pe = $2$')
		panel.plot(Results_Pe20[2][0],Results_Pe20[2][2],'r-',markersize=2,linewidth=1.5,label=r'Pe = $20$')
		panel.plot(Results_Pe100[2][0],Results_Pe100[2][2],'g-',markersize=2,linewidth=1.5,label=r'Pe = $100$')
		panel.plot(Results_ini[2][0],Results_ini[2][2],'k--',markersize=2,linewidth=1.5,label=r'initial policy')
		
		#~ panel.legend(loc='upper right', bbox_to_anchor=(0.99, 0.99),ncol=1,fontsize=8,handlelength=1.5,labelspacing=0.2)
		
		
		if (log_dpanel_Pe == 0):
			######################## panel E1 ########################################
			panel = fig.add_axes([xd1,yd1,dxd1,dyd1])
			panel.tick_params(axis='both',which='both',direction='in',bottom=True,top=True,left=True,right=False)
			for tick in panel.xaxis.get_major_ticks(): tick.label.set_fontsize(axisticslabelfontsize)
			for tick in panel.yaxis.get_major_ticks(): tick.label.set_fontsize(axisticslabelfontsize)
			panel.set_xlabel(r'episodes',fontsize=axislabelfontsize)
			panel.xaxis.set_label_coords(.82, -.1)
			panel.set_ylabel(r'fraction of time in BP phase',fontsize=axislabelfontsize-2)
			panel.set_xlim(xdlim_left,xdlim_middle)
			panel.set_ylim(ydlim_low,ydlim_up)
			panel.xaxis.set_major_locator(MultipleLocator(xd1tics_main))
			panel.xaxis.set_minor_locator(MultipleLocator(xd1tics_minor))
			panel.yaxis.set_major_locator(MultipleLocator(0.5))
			panel.yaxis.set_minor_locator(MultipleLocator(0.1))
			
			######################## panel E2 ########################################
			panel = fig.add_axes([xd2,yd2,dxd2,dyd2])
			panel.tick_params(axis='both',which='both',direction='in',bottom=True,top=True,left=False,right=True)
			for tick in panel.xaxis.get_major_ticks(): tick.label.set_fontsize(axisticslabelfontsize)
			for tick in panel.yaxis.get_major_ticks(): tick.label.set_fontsize(axisticslabelfontsize)
			panel.set_xlim(xdlim_middle,xdlim_right)
			panel.set_ylim(ydlim_low,ydlim_up)
			panel.xaxis.set_major_locator(MultipleLocator(xd2tics_main))
			panel.xaxis.set_minor_locator(MultipleLocator(xd2tics_minor))
			panel.yaxis.set_major_locator(MultipleLocator(0.5))
			panel.yaxis.set_minor_locator(MultipleLocator(0.1))
			panel.set_yticklabels([])
			panel.text(0.75,0.88,r'(e)',fontsize=axislabelfontsize,transform=panel.transAxes)
		else:
			######################## panel E3 ########################################
			panel = fig.add_axes([xd3,yd3,dxd3,dyd3])
			panel.tick_params(axis='both',which='both',direction='in',bottom=True,top=True,left=True,right=True)
			for tick in panel.xaxis.get_major_ticks(): tick.label.set_fontsize(axisticslabelfontsize)
			for tick in panel.yaxis.get_major_ticks(): tick.label.set_fontsize(axisticslabelfontsize)
			panel.set_xlabel(r'Pe',fontsize=axislabelfontsize)
			panel.set_xscale('log')
			panel.set_xlim(xelim_left,xelim_right)
			
			panel.set_ylabel(r'fraction of time in BP phase',fontsize=axislabelfontsize-2)
			panel.set_ylim(ydlim_low,ydlim_up)
			panel.yaxis.set_major_locator(MultipleLocator(0.5))
			panel.yaxis.set_minor_locator(MultipleLocator(0.1))
			
			panel.plot(Pes,avgBPphaseduration_ini/(avgBPphaseduration_ini+avgABPphaseduration_ini),'k--',markersize=2,linewidth=1.5,label=r'initial policy')
			panel.plot(Pes,avgBPphaseduration/(avgBPphaseduration+avgABPphaseduration),'o-',color='brown',markersize=3,linewidth=1.8)
			
			panel.text(0.9,0.88,r'(e)',fontsize=axislabelfontsize,transform=panel.transAxes)
		
		
		
		
		
		if (log_epanel_Pe == 0):
			######################## panel D1 ########################################
			panel = fig.add_axes([xe1,ye1,dxe1,dye1])
			panel.tick_params(axis='both',which='both',direction='in',bottom=True,top=True,left=True,right=False)
			for tick in panel.xaxis.get_major_ticks(): tick.label.set_fontsize(axisticslabelfontsize)
			for tick in panel.yaxis.get_major_ticks(): tick.label.set_fontsize(axisticslabelfontsize)
			panel.set_xlabel(r'episodes',fontsize=axislabelfontsize)
			panel.xaxis.set_label_coords(.82, -.2)
			panel.set_xlim(xdlim_left,xdlim_middle)
			panel.xaxis.set_major_locator(MultipleLocator(xd1tics_main))
			panel.xaxis.set_minor_locator(MultipleLocator(xd1tics_minor))
			
			panel.set_ylabel(r'average duration of phase',fontsize=axislabelfontsize-2)
			panel.yaxis.set_label_coords(-0.15, 0.9)
			panel.set_ylim(yelim_low,yelim_up)
			panel.yaxis.set_major_locator(MultipleLocator(0.5))
			panel.yaxis.set_minor_locator(MultipleLocator(0.1))
			
			######################## panel D2 ########################################
			panel = fig.add_axes([xe2,ye2,dxe2,dye2])
			panel.tick_params(axis='both',which='both',direction='in',bottom=True,top=True,left=True,right=False)
			for tick in panel.xaxis.get_major_ticks(): tick.label.set_fontsize(axisticslabelfontsize)
			for tick in panel.yaxis.get_major_ticks(): tick.label.set_fontsize(axisticslabelfontsize)
			panel.set_xlim(xdlim_left,xdlim_middle)
			panel.xaxis.set_major_locator(MultipleLocator(xd1tics_main))
			panel.xaxis.set_minor_locator(MultipleLocator(xd1tics_minor))
			panel.set_xticklabels([])
			
			panel.set_ylim(yelim_low,yelim_up)
			panel.yaxis.set_major_locator(MultipleLocator(0.5))
			panel.yaxis.set_minor_locator(MultipleLocator(0.1))
			
			
			######################## panel D3 ########################################
			panel = fig.add_axes([xe3,ye3,dxe3,dye3])
			panel.tick_params(axis='both',which='both',direction='in',bottom=True,top=True,left=False,right=True)
			for tick in panel.xaxis.get_major_ticks(): tick.label.set_fontsize(axisticslabelfontsize)
			for tick in panel.yaxis.get_major_ticks(): tick.label.set_fontsize(axisticslabelfontsize)
			panel.set_xlim(xdlim_middle,xdlim_right)
			panel.xaxis.set_major_locator(MultipleLocator(xd2tics_main))
			panel.xaxis.set_minor_locator(MultipleLocator(xd2tics_minor))
			
			panel.set_ylim(yelim_low,yelim_up)
			panel.yaxis.set_major_locator(MultipleLocator(0.5))
			panel.yaxis.set_minor_locator(MultipleLocator(0.1))
			panel.set_yticklabels([])
			
			
			
			######################## panel D4 ########################################
			panel = fig.add_axes([xe4,ye4,dxe4,dye4])
			panel.tick_params(axis='both',which='both',direction='in',bottom=True,top=True,left=False,right=True)
			for tick in panel.xaxis.get_major_ticks(): tick.label.set_fontsize(axisticslabelfontsize)
			for tick in panel.yaxis.get_major_ticks(): tick.label.set_fontsize(axisticslabelfontsize)
			panel.set_xlim(xdlim_middle,xdlim_right)
			panel.xaxis.set_major_locator(MultipleLocator(xd2tics_main))
			panel.xaxis.set_minor_locator(MultipleLocator(xd2tics_minor))
			panel.set_xticklabels([])
			
			panel.set_ylim(yelim_low,yelim_up)
			panel.yaxis.set_major_locator(MultipleLocator(0.5))
			panel.yaxis.set_minor_locator(MultipleLocator(0.1))
			panel.set_yticklabels([])
			
			panel.text(0.75,0.76,r'(d)',fontsize=axislabelfontsize,transform=panel.transAxes)
			
			
		else:
			######################## panel D5 ########################################
			panel = fig.add_axes([xe5,ye5,dxe5,dye5])
			panel.tick_params(axis='both',which='both',direction='in',bottom=True,top=True,left=True,right=True)
			for tick in panel.xaxis.get_major_ticks(): tick.label.set_fontsize(axisticslabelfontsize)
			for tick in panel.yaxis.get_major_ticks(): tick.label.set_fontsize(axisticslabelfontsize)
			panel.set_xlabel(r'Pe',fontsize=axislabelfontsize)
			panel.set_xscale('log')
			panel.set_xlim(xelim_left,xelim_right)
			
			panel.set_ylabel(r'average duration of phase',fontsize=axislabelfontsize-2)
			panel.yaxis.set_label_coords(-0.12, 0.96)
			panel.set_yscale('log')
			panel.set_ylim(yelim_low,yelim_up)
			panel.set_yticks([0.0001,0.001,0.01,0.1])
			y_minor = ticker.LogLocator(base=10.0,subs=(0.2,0.4,0.6,0.8),numticks=12)
			panel.yaxis.set_minor_locator(y_minor)
			panel.yaxis.set_minor_formatter(ticker.NullFormatter())
			
			panel.plot(Pes,avgBPphaseduration_ini,'k--',markersize=2,linewidth=1.5,label=r'initial policy')
			panel.plot(Pes,avgBPphaseduration,'o-',color='brown',markersize=3,linewidth=1.8)
			
			panel.text(0.82,0.27,r'BP',fontsize=axislabelfontsize+1,transform=panel.transAxes)
			
			######################## panel D6 ########################################
			panel = fig.add_axes([xe6,ye6,dxe6,dye6])
			panel.tick_params(axis='both',which='both',direction='in',bottom=True,top=True,left=True,right=True)
			for tick in panel.xaxis.get_major_ticks(): tick.label.set_fontsize(axisticslabelfontsize)
			for tick in panel.yaxis.get_major_ticks(): tick.label.set_fontsize(axisticslabelfontsize)
			panel.set_xscale('log')
			panel.set_xlim(xelim_left,xelim_right)
			panel.set_xticklabels([])
			
			panel.set_yscale('log')
			panel.set_ylim(yelim_low,yelim_up)
			panel.set_yticks([0.0001,0.001,0.01,0.1])
			y_minor = ticker.LogLocator(base=10.0,subs=(0.2,0.4,0.6,0.8),numticks=12)
			panel.yaxis.set_minor_locator(y_minor)
			panel.yaxis.set_minor_formatter(ticker.NullFormatter())
			
			panel.plot(Pes,avgABPphaseduration_ini,'k--',markersize=2,linewidth=1.5,label=r'initial policy')
			panel.plot(Pes,avgABPphaseduration,'o-',color='brown',markersize=3,linewidth=1.8)
			
			t = panel.text(0.9,0.76,r'(d)',fontsize=axislabelfontsize,transform=panel.transAxes)
			t.set_bbox(dict(facecolor='white', alpha=1.0, edgecolor='none',pad=1.5))
			panel.text(0.8,0.38,r'ABP',fontsize=axislabelfontsize+1,transform=panel.transAxes)
		
		
		
		pdf.savefig(fig)
	return

def main():
	N_runs = 5000				# number of completely independent runs
	N_episodes = 1000		# number of episodes. Each episode lasts for a time_single_episode
	N_times = 10000			# number of different phase durations
	dt = 0.0001				# integration time step
	time_single_episode = 20.0
	
	
	# initial policy
	inip0 = 0.01
	inip1 = 0.001
	AvgPvalues = np.zeros((3,N_times))
	for t in range(0,N_times):
		AvgPvalues[0][t]=(1+t)*dt
		AvgPvalues[1][t]=inip0
		AvgPvalues[2][t]=inip1
	Results_ini = compute_results(N_times,dt,AvgPvalues)
		
	# Pe100
	foldername = '../RESULTS/Pe100/'
	filename = foldername + 'AvgHvalues_bunch'
	AvgPvalues = compute_AvgPvalues_from_AvgH_bunches(filename,dt,N_times,1,N_episodes)
	#~ filename = foldername + 'AvgPvalues_bunch'
	#~ AvgPvalues = read_AvgPvalues_bunches(filename,dt,N_times,1,N_episodes)
	Results_Pe100 = compute_results(N_times,dt,AvgPvalues)
		
	# read final probabilities of switching Pe=1
	foldername = '../RESULTS/Pe1_Kaur/'
	#~ filename = foldername + 'AvgPvaluesEpisode.dat'
	#~ AvgPvaluesEpisode = read_AvgPvaluesEpisode(filename,dt,N_times)
	#~ AvgPvalues = AvgPvaluesEpisode[N_episodes-1]
	filename = foldername + 'AvgHvalues_bunch'
	#~ AvgPvalues = compute_AvgPvalues_from_AvgH_bunches(filename,dt,N_times,2,N_episodes)
	AvgPvalues = compute_AvgPvalues_from_AvgH_bunches(filename,dt,N_times,1,N_episodes)
	#~ filename = foldername + 'AvgPvalues_bunch'
	#~ AvgPvalues = read_AvgPvalues_bunches(filename,dt,N_times,2,N_episodes)
	Results_Pe1 = compute_results(N_times,dt,AvgPvalues)
		
	# read final probabilities of switching Pe=20
	foldername = '../RESULTS/Pe20/'
	filename = foldername + 'AvgHvalues_bunch'
	AvgPvalues = compute_AvgPvalues_from_AvgH_bunches(filename,dt,N_times,1,N_episodes)
	#~ filename = foldername + 'AvgPvalues_bunch'
	#~ AvgPvalues = read_AvgPvalues_bunches(filename,dt,N_times,1,N_episodes)
	Results_Pe20 = compute_results(N_times,dt,AvgPvalues)
	
	
	# read final probabilities of switching Pe=2
	foldername = '../RESULTS/Pe2/'
	filename = foldername + 'AvgHvalues_bunch'
	AvgPvalues = compute_AvgPvalues_from_AvgH_bunches(filename,dt,N_times,2,N_episodes)
	#~ filename = foldername + 'AvgPvalues_bunch'
	#~ AvgPvalues = read_AvgPvalues_bunches(filename,dt,N_times,2,N_episodes)
	Results_Pe2 = compute_results(N_times,dt,AvgPvalues)
	
	# read final probabilities of switching Pe=5
	foldername = '../RESULTS/Pe5/'
	filename = foldername + 'AvgHvalues_bunch'
	AvgPvalues = compute_AvgPvalues_from_AvgH_bunches(filename,dt,N_times,2,N_episodes)
	#~ filename = foldername + 'AvgPvalues_bunch'
	#~ AvgPvalues = read_AvgPvalues_bunches(filename,dt,N_times,2,N_episodes)
	Results_Pe5 = compute_results(N_times,dt,AvgPvalues)
	
	# read final probabilities of switching Pe=10
	foldername = '../RESULTS/Pe10/'
	filename = foldername + 'AvgHvalues_bunch'
	AvgPvalues = compute_AvgPvalues_from_AvgH_bunches(filename,dt,N_times,2,N_episodes)
	#~ filename = foldername + 'AvgPvalues_bunch'
	#~ AvgPvalues = read_AvgPvalues_bunches(filename,dt,N_times,2,N_episodes)
	Results_Pe10 = compute_results(N_times,dt,AvgPvalues)
	
	# read final probabilities of switching Pe=1000
	foldername = '../RESULTS/Pe1000/'
	filename = foldername + 'AvgHvalues_bunch'
	AvgPvalues = compute_AvgPvalues_from_AvgH_bunches(filename,dt,N_times,1,N_episodes)
	#~ filename = foldername + 'AvgPvalues_bunch'
	#~ AvgPvalues = read_AvgPvalues_bunches(filename,dt,N_times,1,N_episodes)
	Results_Pe1000 = compute_results(N_times,dt,AvgPvalues)
	
	# read final probabilities of switching Pe=500
	foldername = '../RESULTS/Pe500/'
	filename = foldername + 'AvgHvalues_bunch'
	AvgPvalues = compute_AvgPvalues_from_AvgH_bunches(filename,dt,N_times,1,N_episodes)
	#~ filename = foldername + 'AvgPvalues_bunch'
	#~ AvgPvalues = read_AvgPvalues_bunches(filename,dt,N_times,1,N_episodes)
	Results_Pe500 = compute_results(N_times,dt,AvgPvalues)
	
	# read final probabilities of switching Pe=50
	foldername = '../RESULTS/Pe50/'
	filename = foldername + 'AvgHvalues_bunch'
	AvgPvalues = compute_AvgPvalues_from_AvgH_bunches(filename,dt,N_times,1,N_episodes)
	#~ filename = foldername + 'AvgPvalues_bunch'
	#~ AvgPvalues = read_AvgPvalues_bunches(filename,dt,N_times,1,N_episodes)
	Results_Pe50 = compute_results(N_times,dt,AvgPvalues)
	
	# read final probabilities of switching Pe=200
	foldername = '../RESULTS/Pe200/'
	filename = foldername + 'AvgHvalues_bunch'
	AvgPvalues = compute_AvgPvalues_from_AvgH_bunches(filename,dt,N_times,1,N_episodes)
	#~ filename = foldername + 'AvgPvalues_bunch'
	#~ AvgPvalues = read_AvgPvalues_bunches(filename,dt,N_times,1,N_episodes)
	Results_Pe200 = compute_results(N_times,dt,AvgPvalues)
	
	# read final probabilities of switching Pe=2000
	foldername = '../RESULTS/Pe2000/'
	filename = foldername + 'AvgHvalues_bunch'
	AvgPvalues = compute_AvgPvalues_from_AvgH_bunches(filename,dt,N_times,5,N_episodes)
	#~ filename = foldername + 'AvgPvalues_bunch'
	#~ AvgPvalues = read_AvgPvalues_bunches(filename,dt,N_times,1,N_episodes)
	Results_Pe2000 = compute_results(N_times,dt,AvgPvalues)
	
	
	# plots
	make_figure(Results_ini,Results_Pe1,Results_Pe2,Results_Pe5,Results_Pe10,Results_Pe20,Results_Pe50,Results_Pe100,Results_Pe200,Results_Pe500,Results_Pe1000,Results_Pe2000)

main()