//March 29th 2006.

// *** Set-up recording
objref igabaa, i_soma_br[nmitx], repere, cdIvec

i_soma_br = new Vector(nmitx)
repere = new Vector(10)
cdIvec = new Vector(5)

insert_iclamps_reg(tstop) // injection of the current
insert_iclamp_ipsc(tstop) // injection of the ipsc

mit[0][0].soma	i_soma_br[0] = new current_gauss(0.5)
i_soma_br[0].dur   = 700
i_soma_br[0].std0  = 0.9
i_soma_br[0].f0    = 4000
i_soma_br[0].tau_f = 1

xopen("fig4.ses")

input_ipsc[0][0].del = 4000
input_ipsc[0][0].tau = 5

a = startsw() // seed generator

// BOUCLE SUR L'INTENSITE DE STIMULATION //pour 3 ca marchait bien
input_reg[0][0].amp = 0.8
cdI = -0.5 	//default
nbnoise = 200

// button parameters cdI, amp ...
xpanel("Choose parameters")
xvalue("gGABAa only for abcd","cdI", 1,"change_ggaba()", 0, 1)
xvalue("Noise realizations number","nbnoise",1,"change_nbnoise()",0,1)
xvalue("Tau gGABAa","input_ipsc[0][0].tau", 1,"change_taugaba()", 0, 1)
xbutton("Run fig 4abcd", "run_fig4abcd()")
xbutton("Run fig 4ef", "run_fig4ef()")
xpanel(930,200)


xpanel("Hint")
    xlabel("To run the simulation, click on Run fig 4abcd or Run fig 4ef.")
    xlabel("Choose the number of noise realizations, 200 is a minimum")
    xlabel("200 may take few hours for fig4abcd, and a day for fig4ef")
xpanel()

proc change_ggaba() {
	cdI=cdI
}
proc change_nbnoise() {
	nbnoise = nbnoise
}
proc change_taugaba() {
	input_ipsc[0][0].tau = input_ipsc[0][0].tau
}

proc run_fig4abcd() {
	sprint(filename,"%s%f.dat","fig4_",cdI)
	outfile.wopen(filename)
	// Loop on the noise
	for w = 0, nbnoise-1 {
		//  1st simu without ipsc
		i_soma_br[0].seed(a+w)		//noise seed
		input_ipsc[0][0].amp = 0
		input_ipsc[0][0].del = 0
		run()
		mit[0][0].spiketimes.printf(outfile,"%10.3f")
		outfile.printf("\n")
		repere.x[0] = mit[0][0].spiketimes.x[9]
		print mit[0][0].spiketimes.x[9],mit[0][0].spiketimes.x[10],mit[0][0].spiketimes.x[10]-mit[0][0].spiketimes.x[9]

		//   with ipsc 
		for idel = 1, 6 {
			i_soma_br[0].seed(a+w)	//noise seed
			input_ipsc[0][0].amp = cdI
			input_ipsc[0][0].del = repere.x[0] + 5*idel
			run()
			mit[0][0].spiketimes.printf(outfile,"%10.3f")
			outfile.printf("\n")
			print mit[0][0].spiketimes.x[9],mit[0][0].spiketimes.x[10],mit[0][0].spiketimes.x[10]-mit[0][0].spiketimes.x[9]
		}
	print "w=", w
	}
	outfile.close()

	xpanel("Matlab analysis")
  	xlabel("In a matlab window, type figure4abcd(X)")
  	xlabel("X is the number of noise realizations you chose")
	xpanel()
}

proc run_fig4ef() {
	cdIvec.x[0] = -0.1
	cdIvec.x[1] = -0.2
	cdIvec.x[2] = -0.5
	cdIvec.x[3] = -1
	cdIvec.x[4] = -2
	for ncdI = 0, 4 {
		cdI = cdIvec.x[ncdI]
		sprint(filename,"%s%f.dat","fig4_",cdIvec.x[ncdI])
		outfile.wopen(filename)
		// Loop on the noise
		for w = 0, nbnoise-1 {
			//  1st simu without ipsc
			i_soma_br[0].seed(a+w)		//noise seed
			input_ipsc[0][0].amp = 0
			input_ipsc[0][0].del = 0
			run()
			mit[0][0].spiketimes.printf(outfile,"%10.3f")
			outfile.printf("\n")
			repere.x[0] = mit[0][0].spiketimes.x[9]
			print mit[0][0].spiketimes.x[9],mit[0][0].spiketimes.x[10],mit[0][0].spiketimes.x[10]-mit[0][0].spiketimes.x[9]

			//   with ipsc
			for idel = 1, 6 {
				i_soma_br[0].seed(a+w)	//noise seed
				input_ipsc[0][0].amp = cdI
				input_ipsc[0][0].del = repere.x[0] + 5*idel
				run()
				mit[0][0].spiketimes.printf(outfile,"%10.3f")
				outfile.printf("\n")
				print mit[0][0].spiketimes.x[9],mit[0][0].spiketimes.x[10],mit[0][0].spiketimes.x[10]-mit[0][0].spiketimes.x[9]
			}
			print "w=", w
		}
		outfile.close()
	}
	xpanel("Matlab analysis")
  	xlabel("In a matlab window, type figure4ef(X)")
  	xlabel("X is the number of noise realizations you chose")
	xpanel()
}