numTrial = 1
StoreEveryTrial = 1           // Determines how often StoreW is called

//PLAST=0	// NOT USED
 
ScaleExEx =0     //1 //0.5 : 0=no plasticity, 1=plasticity: SCALING
ScaleExInh=0
ScaleInhEx=0
STDP = 0         //1
STDPINH = 0

GRAPHICS = 1
COMPUTE_SLOPE = 1
t_slope_start = 20
t_slope_window = 0.6

alfa_ScaleExEx = 0.5
alfa_ScaleExInh = 0.5

WRITE_VOLTAGE_FILE = 0	// Write voltage.dat file (can get big really easy...)
WRITE_W_FILE = 0	// Write W.dat file (can get big really easy...)
WRITE_RASTER_FILE = 0

xopen("NETWORK.oc")




FILE = 1 // Write output files


xopen("PANEL_PARAMETERS.oc")	// Loads panel for PARAMETERS

objectvar fSynSpace_Out
fSynSpace_Out = new File()
fSynSpace_Out.wopen("SynSpace_Out.dat")

P1_start = 0.02		// Input->Ex 0.02
P1_step_size = 0.00075	//0.00055
P1_steps = 35		//48

P2_start = 0		// 0.1 Inh->Ex
P2_step_size = 0.0175	//0.0127
P2_steps = 35		//48

nr_iterations = 50	// Number of times the whole sampling of synaptic space will run

xopen("ParamExtract.oc")	// functions to get some stats on the voltage traces

proc MULTI() {	local p1, p2, ii, j, input_nr, run_t, input_block, input_idx

			// Sets EXCITATORY synaptic strengths ONTO Inh neurons  <--  This synapse is FIXED W throughout
			for (ii = 0; ii < nInh; ii = ii + 1) {		
				j=0
				sExInh[ii][j].gmaxAMPA = (totEx_Inh/nInh)*1.7
			}
			for (ii = 0; ii < nInh; ii = ii + 1) {		
				for (j = 1; j < nEx_Inh; j = j + 1) {
					// As this is, all the inputs WITHIN 'A' Inh will be the same, but then one Inh will receive all the weak inputs, the other all the strong, etc...
					sExInh[ii][j].gmaxAMPA = ((totEx_Inh/nInh)/3.8) * (1 + (ii/nInh)*1.4) 
				}
			}
			
	
	for (p1 = 1; p1 <= P1_steps; p1 = p1 + 1) {		// Runs for parameter 1

		for (p2 = 1; p2 <= P2_steps; p2 = p2 + 1) {	// Runs for parameter 2

			print "P1: ",p1,"P2: ",p2

			for (ii = ExINPUT; ii < nEx; ii = ii + 1) {		// Sets EXCITATORY synaptic strengths ONTO true Ex neurons (not input)
				for (j = 0; j < nEx_Ex; j = j + 1) {
					sExEx[ii][j].gmaxAMPA = (P1_start + P1_step_size * (p1 - 1)) / ExINPUT
				}
			}
			for (ii = ExINPUT; ii < nEx; ii = ii + 1) {
				for (j = 0; j < nInh_Ex; j = j + 1) {
					sInhEx[ii][j].gmaxGABA = (P2_start + P2_step_size * (p2 - 1)) / nInh	// 18 Oct 2007, replaced ExINPUT here
				}
			}

			// Selects which input cells will fire
			for (ipt_idx = 0; ipt_idx < ExINPUT; ipt_idx = ipt_idx + 1) {	// RESET: Inhibits all the inputs from firing
				Ex[ipt_idx].soma.end_INPUT = 0
			}
			for (input_block = 0; input_block < ExINPUT_end; input_block = input_block + ExINPUT_block_size) {	// Runs code for different inputs firing
				for (input_nr = 0; input_nr < ExINPUT_block_size; input_nr = input_nr + 1) {		// Gets INPUTS firing in block fashion
					input_idx = input_block + input_nr
					Ex[input_block + input_nr].soma.end_INPUT = 20
				}
							
				for (run_t = 1; run_t <= nr_iterations; run_t = run_t + 1) {	// Repeats the same conditions "nr_iterations" times

					InitCond()
					for (trial=1;trial<=numTrial;trial=trial+1) {				// This are the trials, if learning is implemented


						SCALE_ExIAF = alfa_ScaleExEx * ScaleExEx         //0.5
						SCALE_InhIAF = alfa_ScaleExInh * ScaleExInh

						//print ">>>>>>>>>>>>>>>Trial=",trial
						if (trial%StoreEveryTrial==0 || trial==1) {
							if (WRITE_W_FILE) StoreW()
						}

						run()

						//get_EPSP_amp()
						max_Ex_EPSP_amp = 0
						max_Inh_EPSP_amp = 0
						total_Inh_firing = get_how_many_Inh_fired()

						slope = 0
						if (COMPUTE_SLOPE) {
							slope = get_EPSP_slope(Ex[0].soma.start_INPUT, t_slope_start, t_slope_window, dt)
						}
				
						// Write synaptic space OUTPUT FILE
						// fOutParam.printf("%4d\t%8.4f\t%4d\n",p1,sExEx[1][0].gmaxAMPA,Ex[1].soma.Ca_ExIAF) // JUST BACKUP LINE
						
						total_Inh_firing = 0
						for (Inh_counter = 0; Inh_counter < nInh; Inh_counter = Inh_counter + 1) {
							if (Inh[Inh_counter].soma.Ca_InhIAF > 1  ) {
								my_single_Inh_firing = 1
							} else {
								my_single_Inh_firing = Inh[Inh_counter].soma.Ca_InhIAF
							}
							total_Inh_firing = total_Inh_firing + my_single_Inh_firing
						}
						
						//fSynSpace_Out.printf("%d\t%f\t%d\t%f\t%d\t%f\t%d\t%f\t%d\t%d\t%f\n", p1, sExEx[ExINPUT][0].gmaxAMPA, Ex[ExINPUT].soma.Ca_ExIAF, max_Ex_EPSP_amp, p2, sInhEx[ExINPUT][0].gmaxGABA, Inh[0].soma.Ca_InhIAF, max_Inh_EPSP_amp, total_Inh_firing, input_idx+1, slope)
						fSynSpace_Out.printf("%d\t%d\t%d\t%d\t%d\t%.2e\n", p1, Ex[ExINPUT].soma.Ca_ExIAF, p2, total_Inh_firing, input_idx+1, slope)
					}	// End TRIALS loop
				}	// End nr_iterations
			}	// End variable inputs
		}	// End P1
	}	// End times_it_will_run
    //CLOSEFILES()

	fSynSpace_Out.close()

	print "Code finished running successfully!"
}



// Writes parameters used to scan parameter space

objectvar fSynSpace_Param
fSynSpace_Param = new File()
fSynSpace_Param.wopen("SynSpace_Param.dat")


fSynSpace_Param.printf("P1_start = %f\n", P1_start)
fSynSpace_Param.printf("P1_step_size = %f\n", P1_step_size)
fSynSpace_Param.printf("P1_steps = %d\n", P1_steps)
fSynSpace_Param.printf("P2_start = %f\n", P2_start)
fSynSpace_Param.printf("P2_step_size = %f\n", P2_step_size)
fSynSpace_Param.printf("P2_steps = %d\n", P2_steps)
fSynSpace_Param.printf("nr_repetitions = %d\n", nr_iterations)
//fSynSpace_Param.printf("alfa_ScaleExEx = %f\n", alfa_ScaleExEx)
//fSynSpace_Param.printf("alfa_ScaleExInh = %f\n", alfa_ScaleExInh)
fSynSpace_Param.printf("AmpaMaxExEx = %f\n", AmpaMaxExEx)
fSynSpace_Param.printf("AmpaMaxExInh = %f\n", AmpaMaxExInh)
fSynSpace_Param.printf("GabaMax = %f\n", GabaMax)
fSynSpace_Param.printf("AMPANMDARATIO_EPlasSyn = %f\n", AMPANMDARATIO_EPlasSyn)
fSynSpace_Param.printf("totEx_Ex = %f\n", totEx_Ex)
fSynSpace_Param.printf("totEx_Inh = %f\n", totEx_Inh)
fSynSpace_Param.printf("totInh_Ex = %f\n", totInh_Ex)
fSynSpace_Param.printf("ExNoise = %f\n", ExNoise)
fSynSpace_Param.printf("InhNoise = %f\n", InhNoise)
fSynSpace_Param.printf("Thr_IAF_noise = %f\n", Thr_IAF_noise)
fSynSpace_Param.printf("nEx_Ex = %d\n", nEx_Ex)
fSynSpace_Param.printf("nEx_Inh = %d\n", nEx_Inh)
fSynSpace_Param.printf("nInh_Ex = %d\n", nInh_Ex)
fSynSpace_Param.printf("nEx = %d\n", nEx)
fSynSpace_Param.printf("nInh = %d\n", nInh)
fSynSpace_Param.printf("ExINPUT = %d\n", ExINPUT)
fSynSpace_Param.printf("ExINPUT_block_size = %d\n", ExINPUT_block_size)	
fSynSpace_Param.printf("ExINPUT_end = %d\n", ExINPUT_end)
fSynSpace_Param.printf("ExEx_Delay = %f\n", Ex[ExINPUT].soma.Delay_EPlasSom)
fSynSpace_Param.printf("ExInh_Delay = %f\n", Ex[ExINPUT-1].soma.Delay_EtoIPlasSom)
fSynSpace_Param.printf("InhEx_Delay = %f\n", Inh[nInh-1].soma.Delay_IPlasSom)
//fSynSpace_Param.printf("ElapsedTime = %f\n", ElapsedTime)


fSynSpace_Param.printf("INPUT_diam = %d\n", Ex[0].soma.diam)
fSynSpace_Param.printf("INPUT_L = %d\n", Ex[0].soma.L)
fSynSpace_Param.printf("INPUT_gpas = %f\n", Ex[0].soma.gPAS_ExIAF)

fSynSpace_Param.printf("Ex_diam = %d\n", Ex[ExINPUT].soma.diam)
fSynSpace_Param.printf("Ex_L = %d\n", Ex[ExINPUT].soma.L)
fSynSpace_Param.printf("Ex_soma_e_pas = %f\n", Ex[ExINPUT].soma.ePAS_ExIAF)
fSynSpace_Param.printf("Ex_soma_g_pas = %f\n", Ex[ExINPUT].soma.gPAS_ExIAF)
fSynSpace_Param.printf("Thr_ExIAF = %f\n", Ex[ExINPUT].soma.Thr_ExIAF)
fSynSpace_Param.printf("ePAS_ExIAF = %f\n", Ex[ExINPUT].soma.ePAS_ExIAF)
fSynSpace_Param.printf("gPAS_ExIAF = %f\n", Ex[ExINPUT].soma.gPAS_ExIAF)
fSynSpace_Param.printf("eOFF_ExIAF = %f\n", Ex[ExINPUT].soma.eOFF_ExIAF)


fSynSpace_Param.printf("Inh_diam = %d\n", Inh[0].soma.diam)
fSynSpace_Param.printf("Inh_L = %d\n", Inh[0].soma.L)
fSynSpace_Param.printf("Inh_gpas = %f\n", Inh[0].soma.gPAS_InhIAF)

fSynSpace_Param.printf("WRITE_VOLTAGE_FILE = %d\n", WRITE_VOLTAGE_FILE)
//fSynSpace_Param.printf("WRITE_CONDUCTANCE_FILE = %d\n", WRITE_CONDUCTANCE_FILE)
fSynSpace_Param.printf("GRAPHICS = %d\n", GRAPHICS)
fSynSpace_Param.printf("WRITE_RASTER_FILE = %d\n", WRITE_RASTER_FILE)
//fSynSpace_Param.printf("PLOT_CONDUCTANCES = %d\n", PLOT_CONDUCTANCES)
//fSynSpace_Param.printf("WRITE_IW_FILE = %d\n", WRITE_IW_FILE)

fSynSpace_Param.printf("COMPUTE_SLOPE = %d\n", COMPUTE_SLOPE)
fSynSpace_Param.printf("t_slope_start = %.2f\n", t_slope_start)
fSynSpace_Param.printf("t_slope_window = %.2f\n", t_slope_window)


fSynSpace_Param.printf("dt = %.3f\n", dt)
fSynSpace_Param.printf("tstop = %d\n", tstop)


print "Wrote Synaptic Space Parameters file successsfuly!\n"
fSynSpace_Param.close()





// Calculates how many Inh fired
func get_how_many_Inh_fired() { local i, global_Ca_Inh
	
	global_Ca_Inh = 0
	for(i = 0; i < nInh; i = i + 1) {
		if (Inh[i].soma.Ca_InhIAF > 1) {
			binary_counter = 1
		} else {
			binary_counter = Inh[i].soma.Ca_InhIAF
		}
		global_Ca_Inh = global_Ca_Inh + binary_counter
	}
	return global_Ca_Inh
}




proc InitCond() {	local i, j

	// ExEx
	for(i=ExINPUT; i<nEx; i=i+1) {                // ExINPUT CELLS
		for(j=0;j<nEx_Ex;j=j+1) {
                   sExEx[i][j].scale = ScaleExEx
                   sExEx[i][j].stdp = STDP
		}                               //end nInh_Ex
	}
	
	// ExInh
	for(i=0; i<nInh; i=i+1) {            // ExINPUT CELLS
		for(j=0;j<nEx_Inh;j=j+1) {
			sExInh[i][j].scale =  ScaleExInh
            sExInh[i][j].stdp =   STDP
           }                               //end nInh_Ex
	}

	// InhEx
	for(i=ExINPUT; i<nEx; i=i+1) {
		for(j=0;j<nInh_Ex;j=j+1) {
			sInhEx[i][j].scale = ScaleInhEx
			sInhEx[i][j].stdp = STDPINH
		}
	}
}                                       //InitCond()