proc connect_cells() {local i, j  localobj cvec
	for j=0, ngranule-1 {
		// g2m_connections require previously specified distribution
		// parameters
		cvec = g2m_connections(j)
		for i=0, cvec.size-1 {
			mgconnect(cvec.x[i], j, 0)
		}
	}
}

// connect mitral $1 to granule $2 based on mitral position and granule
// position. Only connect if secden long enough to reach granule
proc mgconnect() {local mp, gp, delta, x, gpdarc  localobj mgrs
	mp = mitral_x.x[$1]
	gp = granule_x.x[$2]
            delta = gp - mp
	if (ring) { 
               if (delta > Lsec) { 
                       delta -= len 
               }else if (delta < -Lsec) { 
                       delta += len 
               } 
            } 
	if (abs(delta) > Lsec) { return }
	x = delta/Lsec
	gpdarc = g2m_priden_position($1, $2, x)
	if (delta >= 0) { // use secden[0]
		mgrs = mk_mgrs($1, "secden[0]", x, $2, "priden2[0]", gpdarc, $3)
	}else{ // use secden[1]
		mgrs = mk_mgrs($1, "secden[1]", -x, $2, "priden2[0]", gpdarc, $3)
	}
	// store position and relation info in synapse to aid in
	// selection and movie making
	if (object_id(mgrs)) {
		if (object_id(mgrs.fi)) {
			mgrs.fi.x = gp
		}
		if (object_id(mgrs.ampanmda)) {
			mgrs.ampanmda.x = gp
		}
	}
}

// given a granule cell index, return a vector of mitral cell indices which
// specify the existence of a mgrs.
// The primary distribution parameter is the fractional connectivity, g2m_mean
// Other aspects of the distribution are coded in the implementation.
// In particular, in the distribution below, a specified number of mitral
// indices are chosen randomly from the set of possible mitral indices. The
// number chosen is specified from a uniform random distribution
// ranging from (g2m_mean +- g2m_var)*nmitral but of course bounded by
// 1 and the possible number (nmitral-1 if the secondary dendrites span
// the domain).
// Note that with this algorithm, the distribution of number of granules
// per mitral cell is unknown. Possibly gaussian with some variance.
// It is possible some mitrals may have a lot of granule connections and
// some mitrals may have only a few. I would think that since there are
// many more granules than mitrals, the number of granules per mitral would
// be narrowly distributed about the mean of g2m_mean*ngranule
// Note that g2m_connections may be called several times with the same
// granule index and it will return the same "random" mitral index vector.
// This allows computation of load balance prior to construction of the
// network as well as simulation independence with respect to cpu location
// of cells.  For a statistically independent connectivity matrix, choose
// a different value for g2m_ranseed.

/* see param.hoc
default_var("g2m_mean", 0.2)
default_var("g2m_var" , 0.1)
default_var("g2m_ranseed", 1)
*/

obfunc g2m_connections() {local i, j, xm, xg, saveseed, n, nmin, nmax \
 localobj r, set, c
	saveseed = mcell_ran4_init(g2m_ranseed)
	nmin = (g2m_mean - g2m_var)*nmitral
	nmax = (g2m_mean + g2m_var)*nmitral
	if (nmin < 1) nmin = 1
	if (nmax > nmitral) nmax = nmitral
	r = new Random()
	r.MCellRan4($1*1000 + 1)
	if (nmin == nmax) {
		n = nmin
	}else{
		n = r.discunif(nmin, nmax)
	}

	// what is the set of mitrals with a secden domain that overlaps
	// the granule? Usually it is all of them but check to make sure
	set = new Vector(nmitral) set.resize(0)
	xg = granule_x.x[$1]
	for i=0, nmitral-1 {
		xm = mitral_x.x[i]
		if (abs(xm - xg) < Lsec) {
			set.append(i)
		}
                        if (ring) {
                            // mitral on far right looping to beginning
                            if ((len - xm + xg) < Lsec) {
                                set.append(i)
                            }
                            // mitral at beginning looking left
                            if ((len - xg + xm) < Lsec) {
                                set.append(i)
                            }
                        }
	}

	// now pick n integers from at set without replacement
	// there are several ways which are better or worse depending on
	// n relative to nmitral
	c = new Vector(n)
	for i=0, n-1 {
		j = r.discunif(i, set.size-1)
		c.x[i] = set.x[j]
		set.x[j] = set.x[i]
		set.x[i] = c.x[i]
	}
	mcell_ran4_init(saveseed)
	return c
}

// return a sparse matrix specifying the mitral,granule connections
// space expensive but need for load balance of mitral pieces
objref sparse_connection_matrix_
proc sparse_connections() {local g, i  localobj c, sm
	sm = new Matrix(nmitral, ngranule, 2)
	for g = 0, ngranule-1 {
		c = g2m_connections(g)
		for i=0, c.size-1 {
			sm.x[c.x[i]][g] = 1
		}
	}
	sparse_connection_matrix_ = sm
}

// return arc position on priden2[0] of granule cell for m=$1 and g=$2
// and secondary dendrite arc position $3
// following implementation is all at 0.5 or random 0-1 independent of distance
// see param.hoc
//default_var("g2m_priden_seed_offset", 0) // 0 means all at 0.5
func g2m_priden_position() {local saveseed, x  localobj r
	if (g2m_priden_seed_offset == 0) {
		return 0.5
	}
	saveseed = mcell_ran4_init(g2m_priden_seed_offset + g2m_ranseed)
	r = new Random()
	r.MCellRan4(1 + $1 + $2*1000)
	x = r.uniform(0,1)
	mcell_ran4_init(saveseed)
	return x
}