// $Id: proc.hoc,v 1.22 1996/09/06 20:38:30 billl Exp $

//* crosstalk() calculates ndot and xdot
// ndot norm of vector against itself X.X
// xdot average cross product <X.Y>
// ddot = ndot-xdot
double szinp[1],ndot[1],xdot[1],ddot[1]
proc crosstalk () { local num,i,j,s1,s2,s3,ni,nf,szin,iflip
  num = $1
  szin=szinp[num] iflip=szin/2
  if (szin<100) {
    tot = fac(szin)/fac(iflip)^2 // the total number of vectors
  } else {
    tot = exp(logfac(szin) - 2*logfac(iflip))
  }
  s1 = s2 = s3 = 0
  nf = fac(iflip)
  for i=1,iflip { // let's not include X.X
    ni = iflip-i
    a = nf/fac(ni)/fac(i)
    a = a*a/tot
    s1 = s1 + a*(ni*BVBASE*BVBASE+ ni + 2*i*BVBASE)
  }
  ndot[num] = iflip*BVBASE*BVBASE + iflip // iflip*1*1 actually
  xdot[num] = s1
  ddot[num] = ndot[num]-xdot[num]
}

//* makemat(mat,inlist,outlist) recreates the outerproduct matrix from iovec's
proc makemat() { local p1,full,szinp
  if ($o2.count!=npatt || $o3.count!=npatt) { 
    printf("ERROR in makemat: list counts differ from npatt %s %s\n",$o2,$o3) return }
  szinp = $o2.object(0).size
  p1 = allocvecs(1)
  mso[p1].resize(szinp*szout)  // scratch matrix
  $o1.resize(szinp*szout)
  $o1.fill(0)
  // generate outer product matrix
  for ltr2(XO,YO,$o2,$o3) {
    mso[p1].outprod(YO,XO)
    $o1.add(mso[p1])
  }
  $o1.div(npatt)
  dealloc(p1)
}

//* makeinh() makes the inhibitory projection out of the output vectors
proc makeinh() { local ii
  ii = $3
  if ($o2.count!=npatt) { print "ERROR: wrong list count"  return }
  $o1.resize(szout)
  $o1.fill(0.)
  for ltr(XO,$o2) { $o1.add(XO) }
  $o1.mul(-xdot[ii]/npatt)
}

//* mkiovec(list,size)  put npatt vecs of size SIZE in LIST
proc mkiovec() { local i,flag,sz,flp
  sz=$2 flp=sz/2
  flag=0 // can't just use '||' because all get eval'ed and generate error
  if ($o1.count!=npatt) { flag=1
  } else if ($o1.object(0).size!=sz) {flag=1}
  if (flag) {
    $o1.remove_all
    for i=0,npatt-1 { 
      tmpvec = new Vector(sz,-1) // input vec
      tmpvec.randwd(flp)
      $o1.append(tmpvec)
    }
    tmpvec = nil
  } else for ltr(XO,$o1) {
    XO.randwd(flp)
  }
}

//* connmap(mat,pre,post) 
// map a connectivity matrix onto a list of pre's and posts'
proc connmap () { local ii,jj,isz,osz,num
  isz = $o2.count  osz = $o3.count
  if (numarg()==0) { print "connmap(mat,pre,post)" 
    return }
  if ($o1.size != isz*osz) {
    printf("ERROR: size mismatch in connmap: %d %d %d\n",$o1.size,isz,osz)
    return }
  ii=0
  for ltr(XO,$o3) {   // postsynaptic
    jj=0
    for ltr(YO,$o2) { //  presynaptic
      num = $o1.mget(ii,jj,isz)
      if (num != 0) {
        if (num<0) {
          printf("CONNMAP ERROR: gmax must not be <0.: %d %d %g.\n",ii,jj,num)
          return }
        XO.setlink(YO.link, YO.nsyn, YO.maxsyn)
        XO.gmax(-1,num)
      }
      jj=jj+1
    }
    ii=ii+1
  }
}

//* clearsyns() goes through postsyn list and reinitializes everything
proc clearsyns() { local i,j
  tmplist = new List("AMPA")
  for ltr(XO,tmplist) { XO.init_arrays(XO.maxsyn) }
  tmplist = new List("GABAA")
  for ltr(XO,tmplist) { XO.init_arrays(XO.maxsyn) }
  cleansyns()
}

//* delset ("AMPA") resets the delays for AMPA
proc delset () { local ii
  for sltr(XO,$s1) for ii=0,XO.nsyn-1 XO.delay(ii,0)
}