/************************************************************
'ca1' model code repository
Written by Ivan Raikov, ivan.g.raikov@gmail.com and Marianne
Bezaire, marianne.bezaire@gmail.com, www.mariannebezaire.com
In the lab of Ivan Soltesz, www.ivansolteszlab.org
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 9, 2015

This file computes an approximation to the LFP generated by
the pyramidal cells in the network. The approximate LFP is 
calculated as the difference in voltage between apical distal 
and basal dendrites. It also contains functions for computing
an average LFP across all pyramidal cells in the network or
across pyramidal cells within a given distance of a particular
point in the network (ie, where an electrode may be positioned).
The time resolution of the LFP calculation may be lower than
that of the simulation by setting lfp_dt in this file.
************************************************************/

func dipole_lfp() { local avglfp, vlfp, ncells, celltype, pci, gid, i, typei localobj cell
	// Calculate the average LFP of all pyramidal cells in the network
	
	ncells = 0			// Initialize the cell count
 	vlfp = 0			// Initialize the LFP variable
   for celltype=0, numCellTypes-1 {											// Iterate over all cell types
		if (strcmp(cellType[celltype].cellType_string, "pyramidalcell")==0) {	// Choosing pyramidal cells	
		   ncells = cellType[celltype].cellEndGid - cellType[celltype].cellStartGid	// Get number of pyramidal cells
		   for pcitr(&i, &typei, &gid, cellType[celltype].cellStartGid, cellType[celltype].cellEndGid) {
			   cell = pc.gid2cell(gid)
			   vlfp = vlfp + cell.apical[cell.NumApical-1].v(0.5) - cell.basal[0].v(0.5)	// For each pyramidal cell
																							//  on this processor, sum
																							//  the LFP contribution
		   }
	   }      
   }
   
   pc.barrier()			// Ensure all processors have completed their calculation before proceeding.
   
   if (ncells > 0) {
       avglfp = pc.allreduce(vlfp, 1)/ncells									// Sum all the LFP contributions across
																				//  all processors using allreduce with
																				//  argument 1, so that the sum of the
																				//  LFP contribution from all pyramidal
																				//  cells in the network is given, and
																				//  then divide by the total number of
																				//  cells for the average.
   }
   
   return avglfp
}

func dipole_pos_lfp() { local avgplfp, plfp, ncellsp, celltype, pci, gid, i, typei, MaxEDist localobj cell
 	// Calculate the average LFP of select pyramidal cells in the network,
	//  only including pyramidal cells whose somata are within MaxEDist
	//  microns of the (x,y,z) location passed into this function as
	//  arguments: $1 = x, $2 = y, $3 = z in microns
	
	MaxEDist = 500 // microns
	sumcell = 0
	ncellsp = 0			// Initialize pyramidal cell counter
	plfp = 0				// Initialize the LFP variable
   
   for celltype=0, numCellTypes-1 {											// Iterate over all cell types
		if (strcmp(cellType[celltype].cellType_string, "pyramidalcell")==0) {	// Choosing pyramidal cells	

		   for pcitr(&i, &typei, &gid, cellType[celltype].cellStartGid, cellType[celltype].cellEndGid) {
				cell = pc.gid2cell(gid)
			   
				// Use one of the following if conditions to constrain position of cells
				//  that contribute to LFP signal:
			   
				//Absolute position based
				//if ((cell.x<1000) && (cell.y<700) && (cell.y>200)) {
				// OR
				// Relative to a point
				if (sqrt((cell.x-$1)*(cell.x-$1) + (cell.y-$2)*(cell.y-$2))<MaxEDist) { 	// Ignoring z location in this example, use
																							//  larger value for MaxEDist if z included
				ncellsp += 1	// Update the number of pyramidal cells
								//  included in the calculation
				plfp = plfp + cell.apical[cell.NumApical-1].v(0.5) - cell.basal[0].v(0.5)	// For each pyramidal cell
																							//  on this processor that
																							//  is within the desired
																							//  location, sum the LFP
																							//  contribution
			   }
		   }
	   }      
   }
   
   
	pc.barrier()			// Ensure all processors have completed their calculation before proceeding.
   
	sumcell=pc.allreduce(ncellsp, 1)
   if (sumcell > 0) {
       avgplfp = pc.allreduce(plfp, 1)/sumcell									// Sum all the LFP contributions across
																				//  all processors using allreduce with
																				//  argument 1, so that the sum of the
																				//  LFP contribution from all pyramidal
																				//  cells in the network is given, and
																				//  then divide by the number of included
																				//  cells for the average.
   }
   
   return avgplfp
}



func dipole_lfp_traceidx() { local vlfp, gid, startgid, endgid, i, j, typei localobj traceidxlist, lfptracelist, lfptrace, cell, vec
	// Compute the lfp for the cells whose voltage is recorded,
	//  based on the list of gids of cells to record that is
	//  passed in as argument 1 ($o1), and record the lfp for
	//  each of those cells into a vector appended to a
	//  results list passed in as argument 2 ($o2).
    
    traceidxlist = $o1		// List if gids to record intracellularly
							//  (the subset of these that are pyramidal
							//  cells will have their lfp recorded here)
							
    lfptracelist = $o2		// List of results vectors into which the
							//  LFP will be recorded for each cell of
							//  interest
    
    for celltype=0, numCellTypes-1 {
		if (strcmp(cellType[celltype].cellType_string, "pyramidalcell")==0) {
            startgid = cellType[celltype].cellStartGid			// Get the gid range for
            endgid = cellType[celltype].cellEndGid				//  pyramidal cells
            break
        }
    }
    
    j = 0
	
    for i=0, traceidxlist.size()-1 {
        
        gid = traceidxlist.x(i)
        
        if (pc.gid_exists(gid) && (gid >= startgid) && (gid <= endgid)) {	// Only include pyramidal cells
																			//  that are on the list of cells
																			//  we are recording from intra-
																			//  cellularly
            
            lfptrace = lfptracelist.o(j)	// Each list object is itself a list (dedicated to a particular
											//  cell, where the gid is given by using the same index into 
											//  into the traceidxlist
            
            cell = pc.gid2cell(gid)
            vlfp = cell.apical[cell.NumApical-1].v(0.5) - cell.basal[0].v(0.5)
            vec = new Vector()
            vec.append(t, vlfp)				// For this time step, create a vector with entries of time
											//  and LFP for this cell
            lfptrace.append(vec.c)			// Append the vector for this time step to the list for this cell
			
            j = j + 1						// Increment the list index
        }
    }
    
    return 0
}

objref lfplist, lfptracelist
lfplist = new List()
lfptracelist = new List()

objref mytraceidxlist, lfptracelist
proc sample_lfp() { local i, avglfp, startgid, endgid, gid, celltype localobj vec, lst
    // Execute this code regularly throughout the simulation
	//  to add datapoints to all the LFP recordings (average
	//  and per-cell).
	
	// Either compute average LFP across all pyramidal cells
	//  or across a subset of cells defined by their position:
	// All pyramidal cells
    avglfp = dipole_lfp()
	// OR
	// Subset of pyramidal cells defined by position
	//avglfp = dipole_pos_lfp(mypoint.x[0],mypoint.x[1]) // x, y, z position of electrode in microns
    
    if (pc.id == 0) {
        vec = new Vector()
		vec.append(t, avglfp)			// For this time step, create a vector with entries of time
										//  and average LFP
		lfplist.append(vec.c)			// Append the vector for this time step to the list
    }
    
    
	// Create a list to hold the LFP data for each recorded pyramidal cell
    if ((t == 0) && name_declared("traceidxlist")) {
        
        execute("~mytraceidxlist = traceidxlist")
        lfptracelist = new List()
        
        for celltype=0, numCellTypes-1 {
	    if (strcmp(cellType[celltype].cellType_string, "pyramidalcell")==0) {
                startgid = cellType[celltype].cellStartGid
                endgid = cellType[celltype].cellEndGid
                break
            }
        }
        
        for i=0, mytraceidxlist.size()-1 {
            gid = mytraceidxlist.x(i)
            if (pc.gid_exists(gid) && (gid >= startgid) && (gid <= endgid)) {
                lst = new List()
                lfptracelist.append(lst)
            }
        }
    }
    
	// For each recorded pyramidal cell, calculate the LFP t this time step
	//  and add it to the cell's LFP recording list.
    if (name_declared("traceidxlist")) {
        execute("~mytraceidxlist = traceidxlist")	// Because traceidxlist is only available
													//  at the top level and not within this
													//  function, execute this code to make
													//  a copy of the variable available within
													//  this function
													
        dipole_lfp_traceidx(mytraceidxlist,lfptracelist)
    }
    cvode.event(t + lfp_dt, "sample_lfp()")		// Add another event to the event queue, to 
												//  execute sample_lfp again, lfp_dt ms from now
}