/* Copyright (c) 2015 EPFL-BBP, All rights reserved.                             
                                                                                 
THIS SOFTWARE IS PROVIDED BY THE BLUE BRAIN PROJECT ``AS IS''                    
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO,            
THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR           
PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE BLUE BRAIN PROJECT                 
BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR           
CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF             
SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR                  
BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY,            
WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE             
OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN           
IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.                                    
                                                                                 
This work is licensed under a 
Creative Commons Attribution-NonCommercial-ShareAlike 4.0 International License.
To view a copy of this license, visit 
http://creativecommons.org/licenses/by-nc-sa/4.0/legalcode or send a letter to 
Creative Commons, 
171 Second Street, Suite 300, 
San Francisco, California, 94105, USA.                 
*/      

/*                                                                               
 * @file synapses.hoc                                                          
 * @brief Template that store the synapses                                
 * @author Werner Van Geit @ BBP                                                 
 * @date 2015                                                                    
*/        

begintemplate synapses_a28017c6c7                

public load_synapses, synapse_list, netcon_list, netstim_list, \
        pre_mtypes_bitmap, n_of_pre_mtypes, id_mtype_map, pre_mtypes, \
        active_pre_mtypes, update_synapses, pre_mtype_netconlists, \
        pre_mtype_freqs,  pre_mtype_netstimlists, reset_synapses, \
        weights, delays, synconf_list

objref synapse_list, netcon_list, netstim_list, pre_mtypes_bitmap, \
        id_mtype_map, pre_mtypes, active_pre_mtypes, pre_mtype_netconlists, \
        pre_cell_ids, pre_mtype_freqs, pre_mtype_netstimlists, rng_list, \
        pre_mtype_weightlists, pre_mtype_synapselists, \
        were_active_pre_mtypes, pre_mtypes_excinh, weights, delays, \
        synconf_list, stringtools, tmp_synapse, this

/** Constructor */                                                          
proc init() { local mtype_id localobj mtype_map_file, mtype_name_string
    strdef line, mtype_name

    // Number of m-types
    n_of_mtypes = 55

    synapse_list = new List()
    rng_list = new List()
    netcon_list = new List()
    netstim_list = new List()
    weights = new Vector()
    delays = new Vector()
    pre_mtype_netconlists = new List()
    pre_mtype_netstimlists = new List()
    pre_mtype_weightlists = new List()
    pre_mtype_synapselists = new List()
    pre_mtypes_bitmap = new Vector(n_of_mtypes, 0)
    pre_mtypes_excinh = new Vector(n_of_mtypes, 1)
    active_pre_mtypes = new Vector(n_of_mtypes, 1)
    were_active_pre_mtypes = new Vector(n_of_mtypes, 1)
    pre_mtype_freqs = new Vector(n_of_mtypes, 2.0)

    pre_mtypes = new Vector()
    pre_cell_ids = new Vector()

    for i=0, n_of_mtypes-1 {
        pre_mtype_netconlists.append(new List())
        pre_mtype_netstimlists.append(new List())
        pre_mtype_synapselists.append(new List())
        pre_mtype_weightlists.append(new Vector())
    }

    id_mtype_map = new List()

    // Read the mapping between m-type names and m-type IDs
    mtype_map_file = new File("synapses/mtype_map.tsv")
    mtype_map_file.ropen()
    for i=0, n_of_mtypes-1 {
        mtype_map_file.gets(line)
        sscanf(line, "%d\t%s", &mtype_id, mtype_name)
        mtype_name_string = new String(mtype_name)
        id_mtype_map.append(mtype_name_string)
    }
    mtype_map_file.close()

    stringtools = new StringFunctions()                                              
}

/** Update the active synapses in the GUI
                                                                                 
    Arguments:                                                                   
        synapse_plot: synapse plot object               
*/                                                          
proc update_synapses() { localobj pre_mtype_netcons, pre_mtype_netstims, \
                            pre_mtype_weights, pre_mtype_synapses, synapse_plot
   // synapse_plot = $o1

    for i=0, n_of_mtypes-1 {
        pre_mtype_netcons = pre_mtype_netconlists.o(i)
        pre_mtype_weights = pre_mtype_weightlists.o(i)
        pre_mtype_synapses = pre_mtype_synapselists.o(i)
		
		/*

        // Use different color for exc / inh synapses 
        if (pre_mtypes_excinh.x[i] == 0) {
            // Red, excitatory
            syn_color = 2
        } else {
            // Yellow, inhibitory
            syn_color = 5
        }
		
		*/

        // Loop over the synapses. 
        // enable and plot active synapses, 
        // disable and remove inactive ones
        // only do something if status of a synapse changed since last update
        for j=0, pre_mtype_netcons.count()-1 {
            if (active_pre_mtypes.x[i] == 1 && \
                    were_active_pre_mtypes.x[i] == 0) {
                // Set weight to correct value
                pre_mtype_netcons.o(j).weight = pre_mtype_weights.x[j]
                // Draw the synapses
            //    synapse_plot.point_mark(pre_mtype_synapses.o(j), \
            //                            syn_color, 4, 6)
            } else if (active_pre_mtypes.x[i] == 0 && \
                    were_active_pre_mtypes.x[i] == 1) {
                // Disable synapse by setting weight to 0
                pre_mtype_netcons.o(j).weight = 0.0
                // Remove synapse from plot
             //   synapse_plot.point_mark_remove(pre_mtype_synapses.o(j))
            }
        }

        // Update the list were_active_pre_mtypes, which is the list of 
        // mtypes that are active after this update
        if (active_pre_mtypes.x[i] == 1 && were_active_pre_mtypes.x[i] == 0) {
            were_active_pre_mtypes.x[i] = 1
        } else if (active_pre_mtypes.x[i] == 0 && \
                    were_active_pre_mtypes.x[i] == 1) {
            were_active_pre_mtypes.x[i] = 0
        }

        // Update the netstim frequencies
        pre_mtype_netstims = pre_mtype_netstimlists.o(i)
        for j=0, pre_mtype_netstims.count()-1 {
            pre_mtype_netstims.o(j).interval = 1000.0 / pre_mtype_freqs.x[i]
        }
    }

}

/** Disable all the synapses */
proc reset_synapses() {
    for i=0, n_of_mtypes-1 {
        were_active_pre_mtypes.x[i] = 0
    }
}



/** Load the synconf file
    Arguments:                                                                   
        nsynapses: number of synapses in cell
*/                                                          
proc load_synconf() { localobj synconf_file, synconf_gids
    strdef synconf_string
    nsynapses = $1

    synconf_list = new List()
    for isynapse=0, nsynapses-1 {
        synconf_list.append(new List())
    }

    synconf_file = new File("synapses/synconf.txt")                                                        
    {synconf_file.ropen()}
    
    synconf_gids = new Vector()                                                      
                                                                                                                                                                 
    while (synconf_file.gets(synconf_string) > 0) {                                  
        stringtools.left(synconf_string, stringtools.len(synconf_string)-1)          
        synconf_gids.scantil(synconf_file, -1e15)                                    
        for i=0, synconf_gids.size()-1 { 
            synconf_list.o(synconf_gids.x[i]).append(new String(synconf_string))
        }                                                 
    }                                                     
}

/** Load all the synapses
                                                                                 
    Arguments:                                                                   
        cell_ref: reference to the cell object               
*/                                                          


	    objref xx_p,yy_p,zz_p,length_p,range_p,xint_p,yint_p,zint_p,seg_cen,ref

proc load_synapses() { local isynapse, nsynapses, ncols, synapse_id, \
                            pre_cell_id, seg_x, synapse_type, \
                            dep, fac, use, tau_d, delay, weight, \
                            base_seed, gid, netstim_id \
                        localobj synapse_data, synapse_file, cell_ref, \
                            synapse, section, rng, netcon, netstim

    printf("Starting to add synapses\n")

    strdef sectionlist_name, synapse_type_name, synconf_string, head_string, \
        tail_string

    cell_ref = $o1
	
	rlx=$2
	rly=$3
	rlz=$4
	rlxang=$5
	
	ref=new Random(2)

    // Load the file that contains all the information about the synapses
    synapse_file = new File("synapses/synapses.tsv")
    {synapse_file.ropen()}

    // Data structure to store the data in synapses.tsv
    synapse_data = new Matrix()
    synapse_data.scanf(synapse_file)
    
    synapse_file.close()

    nsynapses = synapse_data.nrow
    ncols = synapse_data.ncol

    // There is only one cell in this simulation, so let's give it gid 1
    gid = 51423

    // Base seed for the rngs
    base_seed = 0

    // Load list of hoc commands that have to be execute on every synapse
    // to set certain parameters that are not specified in synapse.tsv
    load_synconf(nsynapses)

    for isynapse=0, nsynapses-1 {
        // Read the synapse parameters from the matrix
        synapse_id = synapse_data.x[isynapse][0]
        pre_cell_id = synapse_data.x[isynapse][1]
        pre_mtype = synapse_data.x[isynapse][2]
        sectionlist_id = synapse_data.x[isynapse][3]
        sectionlist_index = synapse_data.x[isynapse][4]
        seg_x = synapse_data.x[isynapse][5]
        synapse_type = synapse_data.x[isynapse][6]
        dep = synapse_data.x[isynapse][7] 
        fac = synapse_data.x[isynapse][8] 
        use = synapse_data.x[isynapse][9] 
        tau_d = synapse_data.x[isynapse][10] 
 //       delay = synapse_data.x[isynapse][11]
        delay = ref.uniform(0,2)  
//        weight = synapse_data.x[isynapse][12] 

        // Create sectionref to the section the synapse will be placed on
        if ( sectionlist_id == 0 ) {
            cell_ref.soma[sectionlist_index] section = new SectionRef()        
            sectionlist_name = "somatic"
			
		cell_ref.soma[sectionlist_index] nseg_p=nseg
		cell_ref.soma[sectionlist_index] nn = n3d()	
		
		xx_p = new Vector(nn) 
		yy_p = new Vector(nn)
		zz_p = new Vector(nn)
		length_p = new Vector(nn)
 		
		for ii = 0,nn-1 {					
		cell_ref.soma[sectionlist_index] xx_p.x[ii] = x3d(ii) 
		cell_ref.soma[sectionlist_index] yy_p.x[ii] = y3d(ii) 
		cell_ref.soma[sectionlist_index] zz_p.x[ii] = z3d(ii) 
		cell_ref.soma[sectionlist_index] length_p.x[ii] = arc3d(ii) 
		}

		
		length_p.div(length_p.x[nn-1])
        
		range_p = new Vector(nseg_p+2) 
		range_p.indgen(1/nseg_p) 
		range_p.sub(1/(2*nseg_p))
		range_p.x[0]=0
		range_p.x[nseg_p+1]=1

        
		xint_p = new Vector(nseg_p+2) 
		yint_p = new Vector(nseg_p+2)
		zint_p = new Vector(nseg_p+2)
		xint_p.interpolate(range_p, length_p, xx_p)
		yint_p.interpolate(range_p, length_p, yy_p)
		zint_p.interpolate(range_p, length_p, zz_p)									

		xint_p.remove(nseg_p+1)
		xint_p.remove(0)
		yint_p.remove(nseg_p+1)
		yint_p.remove(0)
		zint_p.remove(nseg_p+1)
		zint_p.remove(0)
		
		seg_cen=new Vector(nseg_p)
		
		for jj=1,nseg_p {
		      seg_cen.x[jj-1]=abs(((1/(nseg_p+1))*jj)-seg_x)
			  }
			  
		x_syn_real=xint_p.x[seg_cen.min_ind()]
		y_syn_real=yint_p.x[seg_cen.min_ind()]
		z_syn_real=zint_p.x[seg_cen.min_ind()]
		
			pnet_ad_x1=((rlx*cos(rlxang))-(rlz*sin(rlxang)))-(rlx)
            pnet_ad_y1=(rly)-(rly)
            pnet_ad_z1=(((rlx)*sin(rlxang))+((rlz)*cos(rlxang)))-(rlz)
			
			pxrot1=((rlx+x_syn_real)*cos(rlxang))-((rlz+z_syn_real)*sin(rlxang))
	        pyrot1=rly+y_syn_real
            pzrot1=((rlx+x_syn_real)*sin(rlxang))+((rlz+z_syn_real)*cos(rlxang))
	
	        px_fin=pxrot1-pnet_ad_x1
	        py_fin=pyrot1-pnet_ad_y1
	        pz_fin=pzrot1-pnet_ad_z1
		
		rdist1=sqrt(((px_fin-$6)^2)+((py_fin-$7)^2)+((pz_fin-$8)^2))
		rdist2=sqrt(((px_fin-$10)^2)+((py_fin-$11)^2)+((pz_fin-$12)^2))

		prob_act1=exp(-rdist1/$9)
		prob_act2=exp(-rdist2/$13)

		
		
		if ( ref.uniform(0,1) < prob_act1 || ref.uniform(0,1) < prob_act2 ) {
		        weight = synapse_data.x[isynapse][12] 
				} else {
		        weight = 0
				}
			
        } else if ( sectionlist_id == 1 ) {
            cell_ref.dend[sectionlist_index] section = new SectionRef()       
            sectionlist_name = "dend" 
			
		cell_ref.dend[sectionlist_index] nseg_p=nseg
		cell_ref.dend[sectionlist_index] nn = n3d()	
		
		xx_p = new Vector(nn) 
		yy_p = new Vector(nn)
		zz_p = new Vector(nn)
		length_p = new Vector(nn)
 		
		for ii = 0,nn-1 {					
		cell_ref.dend[sectionlist_index] xx_p.x[ii] = x3d(ii) 
		cell_ref.dend[sectionlist_index] yy_p.x[ii] = y3d(ii) 
		cell_ref.dend[sectionlist_index] zz_p.x[ii] = z3d(ii) 
		cell_ref.dend[sectionlist_index] length_p.x[ii] = arc3d(ii) 
		}

		
		length_p.div(length_p.x[nn-1])
        
		range_p = new Vector(nseg_p+2) 
		range_p.indgen(1/nseg_p) 
		range_p.sub(1/(2*nseg_p))
		range_p.x[0]=0
		range_p.x[nseg_p+1]=1

        
		xint_p = new Vector(nseg_p+2) 
		yint_p = new Vector(nseg_p+2)
		zint_p = new Vector(nseg_p+2)
		xint_p.interpolate(range_p, length_p, xx_p)
		yint_p.interpolate(range_p, length_p, yy_p)
		zint_p.interpolate(range_p, length_p, zz_p)									

		xint_p.remove(nseg_p+1)
		xint_p.remove(0)
		yint_p.remove(nseg_p+1)
		yint_p.remove(0)
		zint_p.remove(nseg_p+1)
		zint_p.remove(0)
		
		seg_cen=new Vector(nseg_p)
		
		for jj=1,nseg_p {
		      seg_cen.x[jj-1]=abs(((1/(nseg_p+1))*jj)-seg_x)
			  }
			  
        
	    x_syn_real=xint_p.x[seg_cen.min_ind()]
		y_syn_real=yint_p.x[seg_cen.min_ind()]
		z_syn_real=zint_p.x[seg_cen.min_ind()]
		
			pnet_ad_x1=((rlx*cos(rlxang))-(rlz*sin(rlxang)))-(rlx)
            pnet_ad_y1=(rly)-(rly)
            pnet_ad_z1=(((rlx)*sin(rlxang))+((rlz)*cos(rlxang)))-(rlz)
			
			pxrot1=((rlx+x_syn_real)*cos(rlxang))-((rlz+z_syn_real)*sin(rlxang))
	        pyrot1=rly+y_syn_real
            pzrot1=((rlx+x_syn_real)*sin(rlxang))+((rlz+z_syn_real)*cos(rlxang))
	
	        px_fin=pxrot1-pnet_ad_x1
	        py_fin=pyrot1-pnet_ad_y1
	        pz_fin=pzrot1-pnet_ad_z1
		
		rdist1=sqrt(((px_fin-$6)^2)+((py_fin-$7)^2)+((pz_fin-$8)^2))
		rdist2=sqrt(((px_fin-$10)^2)+((py_fin-$11)^2)+((pz_fin-$12)^2))

		prob_act1=exp(-rdist1/$9)
		prob_act2=exp(-rdist2/$13)

		
		
		if ( ref.uniform(0,1) < prob_act1 || ref.uniform(0,1) < prob_act2 ) {
		        weight = synapse_data.x[isynapse][12] 
				} else {
		        weight = 0
				}
			
        } else if ( sectionlist_id == 2 ) {
            cell_ref.apic[sectionlist_index] section = new SectionRef()        
            sectionlist_name = "apical"	
			
		cell_ref.apic[sectionlist_index] nseg_p=nseg
		cell_ref.apic[sectionlist_index] nn = n3d()	
		
		xx_p = new Vector(nn) 
		yy_p = new Vector(nn)
		zz_p = new Vector(nn)
		length_p = new Vector(nn)
 		
		for ii = 0,nn-1 {					
		cell_ref.apic[sectionlist_index] xx_p.x[ii] = x3d(ii) 
		cell_ref.apic[sectionlist_index] yy_p.x[ii] = y3d(ii) 
		cell_ref.apic[sectionlist_index] zz_p.x[ii] = z3d(ii) 
		cell_ref.apic[sectionlist_index] length_p.x[ii] = arc3d(ii) 
		}

		
		length_p.div(length_p.x[nn-1])
        
		range_p = new Vector(nseg_p+2) 
		range_p.indgen(1/nseg_p) 
		range_p.sub(1/(2*nseg_p))
		range_p.x[0]=0
		range_p.x[nseg_p+1]=1

        
		xint_p = new Vector(nseg_p+2) 
		yint_p = new Vector(nseg_p+2)
		zint_p = new Vector(nseg_p+2)
		xint_p.interpolate(range_p, length_p, xx_p)
		yint_p.interpolate(range_p, length_p, yy_p)
		zint_p.interpolate(range_p, length_p, zz_p)									

		xint_p.remove(nseg_p+1)
		xint_p.remove(0)
		yint_p.remove(nseg_p+1)
		yint_p.remove(0)
		zint_p.remove(nseg_p+1)
		zint_p.remove(0)
		
		seg_cen=new Vector(nseg_p)
		
		for jj=1,nseg_p {
		      seg_cen.x[jj-1]=abs(((1/(nseg_p+1))*jj)-seg_x)
			  }
			  
        
		x_syn_real=xint_p.x[seg_cen.min_ind()]
		y_syn_real=yint_p.x[seg_cen.min_ind()]
		z_syn_real=zint_p.x[seg_cen.min_ind()]
		
			pnet_ad_x1=((rlx*cos(rlxang))-(rlz*sin(rlxang)))-(rlx)
            pnet_ad_y1=(rly)-(rly)
            pnet_ad_z1=(((rlx)*sin(rlxang))+((rlz)*cos(rlxang)))-(rlz)
			
			pxrot1=((rlx+x_syn_real)*cos(rlxang))-((rlz+z_syn_real)*sin(rlxang))
	        pyrot1=rly+y_syn_real
            pzrot1=((rlx+x_syn_real)*sin(rlxang))+((rlz+z_syn_real)*cos(rlxang))
	
	        px_fin=pxrot1-pnet_ad_x1
	        py_fin=pyrot1-pnet_ad_y1
	        pz_fin=pzrot1-pnet_ad_z1
		
		rdist1=sqrt(((px_fin-$6)^2)+((py_fin-$7)^2)+((pz_fin-$8)^2))
		rdist2=sqrt(((px_fin-$10)^2)+((py_fin-$11)^2)+((pz_fin-$12)^2))

		prob_act1=exp(-rdist1/$9)
		prob_act2=exp(-rdist2/$13)

		
		
		if ( ref.uniform(0,1) < prob_act1 || ref.uniform(0,1) < prob_act2 ) {
		        weight = synapse_data.x[isynapse][12] 
				} else {
		        weight = 0    
				}
			
        } else if ( sectionlist_id == 3 ) {
            cell_ref.axon[sectionlist_index] section = new SectionRef()        
            sectionlist_name = "axonal" 
			
		cell_ref.axon[sectionlist_index] nseg_p=nseg
		cell_ref.axon[sectionlist_index] nn = n3d()	
		
		xx_p = new Vector(nn) 
		yy_p = new Vector(nn)
		zz_p = new Vector(nn)
		length_p = new Vector(nn)
 		
		for ii = 0,nn-1 {					
		cell_ref.axon[sectionlist_index] xx_p.x[ii] = x3d(ii) 
		cell_ref.axon[sectionlist_index] yy_p.x[ii] = y3d(ii) 
		cell_ref.axon[sectionlist_index] zz_p.x[ii] = z3d(ii) 
		cell_ref.axon[sectionlist_index] length_p.x[ii] = arc3d(ii) 
		}

		
		length_p.div(length_p.x[nn-1])
        
		range_p = new Vector(nseg_p+2) 
		range_p.indgen(1/nseg_p) 
		range_p.sub(1/(2*nseg_p))
		range_p.x[0]=0
		range_p.x[nseg_p+1]=1

        
		xint_p = new Vector(nseg_p+2) 
		yint_p = new Vector(nseg_p+2)
		zint_p = new Vector(nseg_p+2)
		xint_p.interpolate(range_p, length_p, xx_p)
		yint_p.interpolate(range_p, length_p, yy_p)
		zint_p.interpolate(range_p, length_p, zz_p)									

		xint_p.remove(nseg_p+1)
		xint_p.remove(0)
		yint_p.remove(nseg_p+1)
		yint_p.remove(0)
		zint_p.remove(nseg_p+1)
		zint_p.remove(0)
		
		seg_cen=new Vector(nseg_p)
		
		for jj=1,nseg_p {
		      seg_cen.x[jj-1]=abs(((1/(nseg_p+1))*jj)-seg_x)
			  }
			  
        
		x_syn_real=xint_p.x[seg_cen.min_ind()]
		y_syn_real=yint_p.x[seg_cen.min_ind()]
		z_syn_real=zint_p.x[seg_cen.min_ind()]
		
			pnet_ad_x1=((rlx*cos(rlxang))-(rlz*sin(rlxang)))-(rlx)
            pnet_ad_y1=(rly)-(rly)
            pnet_ad_z1=(((rlx)*sin(rlxang))+((rlz)*cos(rlxang)))-(rlz)
			
			pxrot1=((rlx+x_syn_real)*cos(rlxang))-((rlz+z_syn_real)*sin(rlxang))
	        pyrot1=rly+y_syn_real
            pzrot1=((rlx+x_syn_real)*sin(rlxang))+((rlz+z_syn_real)*cos(rlxang))
	
	        px_fin=pxrot1-pnet_ad_x1
	        py_fin=pyrot1-pnet_ad_y1
	        pz_fin=pzrot1-pnet_ad_z1
		
		rdist1=sqrt(((px_fin-$6)^2)+((py_fin-$7)^2)+((pz_fin-$8)^2))
		rdist2=sqrt(((px_fin-$10)^2)+((py_fin-$11)^2)+((pz_fin-$12)^2))

		prob_act1=exp(-rdist1/$9)
		prob_act2=exp(-rdist2/$13)

		
		
		if ( ref.uniform(0,1) < prob_act1 || ref.uniform(0,1) < prob_act2 ) {
		        weight = synapse_data.x[isynapse][12] 
				} else {
		        weight = 0
				}
			
			
        } else {                                                                
            printf("Sectionlist_id %d not support\n", sectionlist_id)           
            exit(1)                                                             
        }

        // If synapse_type < 100 the synapse is inhibitory, otherwise 
        // excitatory
        if ( synapse_type < 100 ) {
            synapse_type_name = "inhibitory"
            section.sec synapse = new ProbGABAAB_EMS(seg_x)
            synapse.tau_d_GABAA  = tau_d
            rng = new Random()                                                      
            rng.MCellRan4( isynapse*100000+100, gid+250+base_seed )                
            rng.lognormal(0.2, 0.1)                                                 
            synapse.tau_r_GABAA = rng.repick()
            pre_mtypes_excinh.x[pre_mtype] = 1                      
        } else {
            synapse_type_name = "excitatory"
            section.sec synapse = new ProbAMPANMDA_EMS(seg_x)
            synapse.tau_d_AMPA = tau_d
            pre_mtypes_excinh.x[pre_mtype] = 0                      
        }

        synapse.Use = abs( use )                                                  
        synapse.Dep = abs( dep )                                                  
        synapse.Fac = abs( fac )   

        // Execute all the extra synaptic configuration lines from synconf.txt
        tmp_synapse = synapse
        for isynconf=0,synconf_list.o(isynapse).count()-1 {
            synconf_string = synconf_list.o(isynapse).o(isynconf).s
            // Replacing all occurrences of %s with the temporary synapse name
            while( stringtools.substr( synconf_string, "%s" ) != -1 ) {
                stringtools.head(synconf_string, "%s", head_string)
                stringtools.tail(synconf_string, "%s", synconf_string)
                sprint(synconf_string, "%s%s%s", head_string, "tmp_synapse", synconf_string)
            }
            // Add {} around the string
            sprint(synconf_string, "{%s}", synconf_string)
            // Execute the statement 
            execute1(synconf_string, this)
        }

        // Create the random number generator for the synapse
        rng = new Random()                                                          
        rng.MCellRan4( isynapse*100000+100, gid+250+base_seed )                    
        rng.uniform(0,1)                                                            
        synapse.setRNG( rng )                                                       
        rng_list.append(rng)        
                                             
        synapse_list.append(synapse)

        // Check if there is already a netstim (spike generator) 
        // for the presynaptic cell
        // If it exists, use that one
        // Otherwise create a new one
        if (pre_cell_ids.contains(pre_cell_id)) {
            netstim_id = pre_cell_ids.indwhere("==", pre_cell_id)
            netstim = netstim_list.o(netstim_id)
        } else {
            section.sec netstim = new NetStim(seg_x)
            netstim.start = 201
            netstim.interval = 500
            netstim.number = 1e20
            netstim.noise = 0
            netstim_list.append(netstim)
            pre_mtype_netstimlists.o(pre_mtype).append(netstim)
            pre_cell_ids.append(pre_cell_id)
        }

        // Create a connection between the netstim and the synapse
        netcon = new NetCon(netstim, synapse)
        netcon.delay = delay
        netcon.weight = 0.0 
        netcon_list.append(netcon)  
        
        // Save all the objects we made to make them persistent        
        pre_mtype_netconlists.o(pre_mtype).append(netcon) 
        pre_mtype_weightlists.o(pre_mtype).append(weight) 
        weights.append(weight) 
        delays.append(delay) 
        pre_mtype_synapselists.o(pre_mtype).append(synapse) 
        pre_mtypes_bitmap.x[pre_mtype] = 1
        
		
		
        printf("Added %s synapse %d originating from cell %d of m-type %s on %s section %d(%f) and dep %f\n", \
           synapse_type_name, synapse_id, pre_cell_id, id_mtype_map.o(pre_mtype).s, sectionlist_name, \
            sectionlist_index, seg_x, dep)
 
    }    

    // Make a list of all the m-types presynaptic to this cell, used by the GUI
    for i=0, n_of_mtypes-1 {
        if (pre_mtypes_bitmap.x[i] == 1) {
            pre_mtypes.append(i)
        }
    } 

}

endtemplate synapses_a28017c6c7