/************************************************************
'ca1' model code repository
Written by Marianne Bezaire, marianne.bezaire@gmail.com, www.mariannebezaire.com
In the lab of Ivan Soltesz, www.ivansolteszlab.org
Published and latest versions of this code are available online at:
ModelDB: 
Open Source Brain: http://www.opensourcebrain.org/projects/nc_ca1

Main code file: ../main.hoc
This file: Last updated on April 10, 2015


This file defines and calls a procedure for creating connections
between all of the cells of the network. Note that if it is not
desirable to have this procedure make connections from artificial
cells to other cells (real or artificial), perhaps because very 
specific connections need to be made that are defined within the 
stimulation protocol, then the user will either need to set the 
MakeArtConns flag to 0, or will need to ensure that the ConnData 
dataset (conndata_###.dat file) being used specifies 0 connections 
for the connections that they do not want managed by this algorithm. 

This algorithm actually uses a compiled mod file to select the
specific cells that will be connected, saving an order of magnitude
of time over the amount of time that would be required if the
decision were handled in hoc. However, the actual connection-making
process is still carried out by the hoc code below, using the
list of potential connections generated by the compiled mechanism.

For more information about how this works, see the following blog post:
http://mariannebezaire.com/2014/03/28/compile-connection-neuron/

Connectivity algorithm: repeatconn

This algorithm differs from the "fastconn" algorithm in that it
allows a presynaptic cell to make multiple connections to a single
postsynaptic cell if there are not enough potential unique presynaptic
cells available to satisfy the desired convergence level. In contrast, the
"fastconn" algorithm discourages multiple connections from the same
presynaptic cell, which for small networks may result in sparser connectivity
than was specified in the connection matrix dataset. Using this "repeatconn"
algorithm can sometimes result in a rigid, strongly interconnected network 
(for small networks with high convergence values).
************************************************************/

objref postcellgids
objref highIndices

strdef precellType_string, postcellType_string, memstr

proc makeConnections () {local precellType, postcellType, precellStart, precellEnd, postcellStart, postcellEnd, i, j, r

	install_repeatconn() 	// Install the compiled vector method used to choose which connections to make.
							//  By installing the vector method, it becomes available for use in hoc

	for postType = 0, numCellTypes-1 { // Iterate over potential POSTSYNAPTIC cell types
	
		// Load cell type string used to grab connection-specific properties later
		postcellType_string = cellType[postType].cellType_string
		numpost = cellType[postType].numThisHost 	// Track the number of cells of the postsynaptic
													//  cell type that are owned by this host
													
		postcellgids = new Vector(numpost)	// Create a vector to hold a list of gids of all postsynaptic
											//  cell type cells owned by this host
		highIndices = new Vector(numpost)	// Create a vector to hold a list of the last highIndex used
											//  by each cell of this postsynaptic type, owned by this host
		
		for r=0, numpost-1 {				// Fill the newly created vector with the gids
			postcellgids.x[r] = cellType[postType].CellList[r].gid*1.0
			highIndices.x[r] = 1+RandomSeedsConn
		}
		
		newhighindex = 1 + RandomSeedsConn	// Set the highIndex seed used by the random number generator within
											//  the compiled mechanism that decides the connections. By explicitly
											//  setting this seed to a known value that the user controls, (as well
											//  as setting the lowIndex seed, which happens within the compiled
											//  mechanism), this assures the simulation can be reproduced, with all
											//  the same random numbers being generated. Likewise, it becomes very
											//  straightforward to generate a statistically independent simulation.
		
		for preType = 0, numCellTypes-1 { // Iterate over potential PRESYNAPTIC cell types
			if ((MakeArtConns==1) || (cellType[preType].is_art()==0)) { // Only make incoming connections onto artificial cells if that's
																		 //  appropriate behavior for the given stimulation protocol, which
																		 //  will set the MakeArtConns flag.
			
				// Load cell type string used to grab connection-specific properties later
				precellType_string =  cellType[preType].cellType_string
				
				// Attempt to make all desired connections for this combination of presynaptic cell type and postsynaptic cell type.
				//  Return the total number of connections made into the cellType array for the presynaptic cell type.
				cellType[preType].numCons.x[postType] = connectCells(preType, postType, precellType_string, postcellType_string, numpost)
			}	
		}
		cellType[postType].LastHighIndex = 1+RandomSeedsConn	
		for r=0, numpost-1 {
			if (cellType[postType].LastHighIndex<highIndices.x[r]) {
				cellType[postType].LastHighIndex = highIndices.x[r]
			}
		}
		// Save the most recently used highIndex for each postsynaptic cell type
		//  (and for each post synaptic cell) since this info will need
		//  to be referenced by the random number generator compiled mechanism
		//  multiple times (for each cell, several calls are made to the random
		//  number generator, one per each presynaptic cell type), and we need
		//  to ensure that the random stream parts that we use each time do not
		//  overlap with each other.
	}															
	if ((PrintTerminal>0)) {print "Host ", pc.id, " connected cells."}
}

	objref params, conns2make
	objref randnum, synRand
	
	// This function does the meat of the work here. It gets called for each presynaptic cell type x postsynaptic cell type
	//  combination, and attempts to find specific presynaptic cells that can connect to each postsynaptic cell (each
	//  processor tries to find presynaptic inputs only for the postsynaptic cells that it owns). When looking for the
	//  specific presynaptic cells to connect to the postsynaptic cell, it takes into account the presynaptic cell type,
	//  the distance between the somata of each potential presynaptic cell and the postsynaptic cell, as well as the
	//  axonal distribution of the presynaptic cell type.  Arguments to this function are:
	//  $1 precellType: a numerical index
	//  $2 postcellType: a numerical index
	//  $s3 precellType_string: the string (word) for the presynaptic cell type, the 'friendly' name
	//  $s4 postcellType_string: the string (word) for the postsynaptic cell type, the 'friendly' name
	//  $5 numpost: the number of cells of type postcellType that are owned by this processor
	func connectCells () {local r, syn, distance, counter, precellType, postcellType, j, randSynNumber, pre_zpos, post_zpos , pre_xpos, post_xpos, pre_ypos, post_ypos, numSynTypes localobj cell
		// Indices of the cell types in the cellType array
		precellType = $1
		postcellType = $2
					
		// Properties of the connections between the presynaptic cell type and postsynaptic cell type
		synWeight = cellType[$1].wgtConns.x[$2]
		numSyns = cellType[$1].numSyns.x[$2]
		numConns = cellType[$1].numConns.x[$2]		

		numpost = $5	// Number of postsynaptic cells owned by this processor
		
		counter=0	// A variable that keeps track of the total number of connections made between 
					//  every presynaptic cell of type $1 and every postsynaptic cell of type $2

		if (PrintTerminal>1) {print "Host ", pc.id, " is connecting: ", $s3, "s to ", $s4, "s."}
		
		if (numConns != 0 && numpost !=0 && numSyns !=0) {	// If this processor owns any cells of the postsynaptic cell type
															//  AND there should be connections between this presynaptic cell
															//  type and this postsynaptic cell type, then proceed to make
															//  connections.
											
			params = new Vector(27)		// Create the vector of parameters that will be passed to the compiled
										//  mechanism to help it choose potential connections that satisfy the
										//  user's expectations.
										
			if (cellType[postcellType].numCells>=pc.nhost) {
				connlength = numConns*int(cellType[postcellType].numCells/pc.nhost+1.5)+cellType[postcellType].numCells
			} else {
				connlength = numConns+cellType[postcellType].numCells
			}
			
			conns2make = new Vector(1+connlength*6)

			params.x[0]=cellType[precellType].cellStartGid //gmin -- start gid of pre cell type
			params.x[1]=cellType[precellType].cellEndGid //gmax -- end gid of pre cell type
			params.x[2]= numConns // number of connections to make onto each postsynaptic cell
			params.x[3]= cellType[postcellType].numCells // number of cells of the postsynaptic type (total)
			params.x[4]= numpost  // number of cells of the postsynaptic type with gids on this computer

			params.x[5]= cellType[precellType].dist.x[2]*4 // distance
			//params.x[5]= sqrt(LongitudinalLength^2  + LayerVector.sum()^2 + TransverseLength^2) // + LayerVector.sum()^2
						// The distance of the dimension for which the fit equation was designed (should add all three...), in um
						//  this also confusing because, right now, the equation calculates the overall distance (through all dimensions).
						//  If per dimension, should compare to distance in that dimension only
			params.x[6]= distres // Steps - resolution of the fit , in number of steps to take			
			params.x[7]= cellType[precellType].dist.x[0] //a in the Gaussian fit equation
			params.x[8]= cellType[precellType].dist.x[1] //b in the Gaussian fit equation
			params.x[9]= cellType[precellType].dist.x[2] //c in the Gaussian fit equation
			params.x[10]= cellType[precellType].dentateXBins*1.0
			params.x[11]= cellType[precellType].dentateYBins*1.0
			params.x[12]= cellType[precellType].dentateZBins*1.0
			params.x[13]= cellType[precellType].dentateXBinSize
			params.x[14]= cellType[precellType].dentateYBinSize
			params.x[15]= cellType[precellType].dentateZBinSize
			addheight = 0
			if (cellType[precellType].layerflag>0) {addheight=LayerVector.sum(0,cellType[precellType].layerflag-1)}
			params.x[16]= addheight

			params.x[17]= cellType[postcellType].dentateXBins*1.0
			params.x[18]= cellType[postcellType].dentateYBins*1.0
			params.x[19]= cellType[postcellType].dentateZBins*1.0
			params.x[20]= cellType[postcellType].dentateXBinSize
			params.x[21]= cellType[postcellType].dentateYBinSize
			params.x[22]= cellType[postcellType].dentateZBinSize
			addheight = 0
			if (cellType[postcellType].layerflag>0) {addheight=LayerVector.sum(0,cellType[postcellType].layerflag-1)}
			params.x[23]= addheight
			params.x[24]= cellType[postcellType].cellStartGid
			params.x[25]= 1+RandomSeedsConn // newhighindex
			params.x[26]= connlength

			conns2make.repeatconn(params, postcellgids, highIndices)	// Use the newly installed vector method in the following manner:
																	//  1. Create a vector that is more than 3x as long as the number
																	//     of connections you expect to make from this presynaptic
																	//     cell type to the specific postsynaptic cells owned by this
																	//     host.
																	//  2. Create a second vector of length 27, which stores the
																	//     parameters necessary for the compiled mechanism to choose
																	//     the connections correctly.
																	//  3. Make a third vector that holds the gids of each potential
																	//     postsynaptic cell (of the postsynaptic cell type) owned by
																	//     by this processor.
																	//  4. Now, call the vector method on the vector created in step 1,
																	//     and pass in the vectors made in steps 2 and 3 as arguments.
																	//  5. The first vector, which had the vector method called on it,
																	//     will then be populated with info about the specific cells
																	//     to connect. The arrangement of the data within the vector
																	//     will be as follows, where connlength is the number of spaces
																	//     set aside for all connections between this pre- and post-
																	//     synaptic cell type combination, on this host:
																	//     vector.x[0]: the total number of potential connections found
																	//     vector.x[1]: numpost: the number of postsynaptic cells owned by this host
																	//     vector.x[2]: last highIndex used by postsynaptic cell #1
																	//     vector.x[3]: last highIndex used by postsynaptic cell #2
																	//     ...
																	//     vector.x[numpost+1]: last highIndex used by last postsynaptic cell
																	//     vector.x[numpost+2]: presynaptic cell gid for potential connection #1
																	//     vector.x[numpost+3]: presynaptic cell gid for potential connection #2															
																	//     ...
																	//     vector.x[numpost+2+vector.x[0]-1]: presynaptic cell gid for last potential connection														
																	//     ...
																	//     ...
																	//     vector.x[connlength+numpost+2]: postsynaptic cell gid for potential connection #1
																	//     vector.x[connlength+numpost+3]: postsynaptic cell gid for potential connection #2															
																	//     ...
																	//     vector.x[connlength+numpost+2+vector.x[0]-1]: postsynaptic cell gid for last potential connection															

			if (int(conns2make.x[0]/1)>connlength) {				// Make sure earlier data not overwritten in conns2make vector by later data
				print "Problem: vector only saved enough room for ", connlength, " conns but ", int(conns2make.x[0]/1) " were defined, overwriting some data! Simulation is compromised!"
				{pc.barrier()}
				{quit()}
			}

			sprint(memstr, "Defined %s to %s conns", $s3, $s4)	// Update memory usage after
			zzz = mallinfo(zzz, memstr)							//  connecting each cell type

			for r=2, numpost+1 {
				highIndices.x[r-2] = int(conns2make.x[r]/1)
			}
						
			for r=numpost+2, numpost+2+int(conns2make.x[0]/1)-1 {	// For each potential connection suggested by the compiled mechanism
				if (numpost>1) {
					postgid = conns2make.x[r+connlength]
				} else {
					postgid = postcellgids.x[0]
				}
				if (pc.gid_exists(postgid)) {		// Check that the postsynaptic cell exists on this processor
					cell = pc.gid2cell(postgid)		//  and create a reference to the cell object

													// Note that the compiled mechanism only suggests which cells
													//  to connect. It does not suggest which synapses on the
													//  postsynaptic cell should be used to make the desired number
													//  of synapses for the suggested connection. That happens below.
																		
					numSynTypes = cell.pre_list.o(precellType).count()	// Get total number of potential synapse types for this combination
																		//  of presynaptic cell type and postsynaptic cell type.
					if (numSynTypes > 0) {								// If at least one synapse type is defined, then proceed
					
						ransynlist.object(cell.randi).r.discunif(0,numSynTypes-1)	// Create a uniform random INTEGER distribution that will
																					//  allow the specification of a synapse type from all the
																					//  available synapse types available for that presynaptic
																					//  cell type (range = 0 to # synapse types-1)
					
						ranwgtlist.object(cell.randi).r.normal(synWeight, synWeight*synVar)	// Create a normal distribution that will
																							//  allow the specification of a synapse weight, where the
																							//  distribution has the specified mean and variance:
																							//  "r.normal(mean, variance)"
																					
						for s=1,numSyns {	// Randomly pick from the available synapses
											//  enough times to make the desired number 
											//  of synapses per connection.
											
							randSynNumber = ransynlist.object(cell.randi).repick	// Randomly pick a synapse type from the available synapse types
							if (AxConVel<=0) { // AxConVel is the axonal conduction velocity in um/ms
								conDelay =  3 	// Axonal conduction delay + synapse delay (time between when the presynaptic cell spikes
												//  and the synapse on the postsynaptic cell is activated)
							} else { // Don't use the z dimension when computing distances for now - we'll save that for another time
								xpos=xpos_algorithm(conns2make.x[r],cellType[precellType].numCells,cellType[precellType].cellStartGid,cellType[precellType].dentateXBins,cellType[precellType].dentateYBins*cellType[precellType].dentateZBins,cellType[precellType].dentateXBinSize)	// Algorithmically generate cell position
								ypos=ypos_algorithm(conns2make.x[r],cellType[precellType].numCells,cellType[precellType].cellStartGid,cellType[precellType].dentateYBins,cellType[precellType].dentateZBins,cellType[precellType].dentateYBinSize)	// Algorithmically generate cell position
								conDelay = int(10*sqrt((xpos - cell.x)*(xpos - cell.x) + (ypos - cell.y)*(ypos - cell.y))/AxConVel)/10 + 0.5 // use 0.5 for the synaptic cleft delay; this also ensures the delay will never be set to 0, which would error the program
								if (conDelay<dt || conDelay%dt>0 || conDelay/dt>=255) {
									print "BAD CONDELAY: host ", pc.id, " precell: ", conns2make.x[r], ", postcell: ", conns2make.x[r+connlength], ", conDelay = ", conDelay
									conDelay =  3 
								}
							}
																		
							// Connect the presynaptic cell to the chosen synapse on the postsynaptic cell.
							//  Note that the synapse weight parameter is also being used to pass the
							//  presynaptic cell gid to the synapse mechanism. Within the synapse mechanism,
							//  the presynaptic cell gid will be subtracted from the synapse weight, and the
							//  two parameters will be stored separately within the mechanism.
							if (RandomSynWgt==1) {
								nc_append(conns2make.x[r], postgid, precellType, randSynNumber, ranwgtlist.object(cell.randi).repick + (conns2make.x[r]+1)*1000, conDelay)
							} else {
								nc_append(conns2make.x[r], postgid, precellType, randSynNumber, synWeight + (conns2make.x[r]+1)*1000, conDelay)
							}
						}
						counter +=1	// After all of the synapses are formed for a single connection, increment the connection counter by 1
					}
				} else {	// Something is very wrong if the postsynaptic gid does not exist on this machine, print message to terminal
					print $s3, "s to ", $s4, "s:", " can't make gid ", postgid, " which is (r = ", r, " +  connlength = ", connlength, ") = ", r+connlength, " from x[r] = ", conns2make.x[r], " total conns = ", int(conns2make.x[0]/1)
				}
			}
		}
		return counter	// Return the total number of connections made between this presynaptic
						//  cell type and this postsynaptic cell type.
	}
makeConnections()

objref conns2make // This variable is likely to take up a large amount of memory, so redefine it to clear it at the end of use