Synaptic scaling balances learning in a spiking model of neocortex (Rowan & Neymotin 2013)

 Download zip file   Auto-launch 
Help downloading and running models
Accession:147141
Learning in the brain requires complementary mechanisms: potentiation and activity-dependent homeostatic scaling. We introduce synaptic scaling to a biologically-realistic spiking model of neocortex which can learn changes in oscillatory rhythms using STDP, and show that scaling is necessary to balance both positive and negative changes in input from potentiation and atrophy. We discuss some of the issues that arise when considering synaptic scaling in such a model, and show that scaling regulates activity whilst allowing learning to remain unaltered.
Reference:
1 . Rowan MS,Neymotin SA (2013) Synaptic Scaling Balances Learning in a Spiking Model of Neocortex Adaptive and Natural Computing Algorithms, Tomassini M, Antonioni A, Daolio F, Buesser P, ed. pp.20
Model Information (Click on a link to find other models with that property)
Model Type: Realistic Network;
Brain Region(s)/Organism: Neocortex;
Cell Type(s): Neocortex V1 L6 pyramidal corticothalamic GLU cell; Neocortex V1 L2/6 pyramidal intratelencephalic GLU cell; Neocortex V1 interneuron basket PV GABA cell; Neocortex fast spiking (FS) interneuron; Neocortex spiny stellate cell; Neocortex spiking regular (RS) neuron; Neocortex spiking low threshold (LTS) neuron; Abstract integrate-and-fire adaptive exponential (AdEx) neuron;
Channel(s):
Gap Junctions:
Receptor(s): GabaA; AMPA; NMDA;
Gene(s):
Transmitter(s): Gaba; Glutamate;
Simulation Environment: NEURON; Python;
Model Concept(s): Synaptic Plasticity; Long-term Synaptic Plasticity; Learning; STDP; Homeostasis;
Implementer(s): Lytton, William [bill.lytton at downstate.edu]; Neymotin, Sam [Samuel.Neymotin at nki.rfmh.org]; Rowan, Mark [m.s.rowan at cs.bham.ac.uk];
Search NeuronDB for information about:  Neocortex V1 L6 pyramidal corticothalamic GLU cell; Neocortex V1 L2/6 pyramidal intratelencephalic GLU cell; Neocortex V1 interneuron basket PV GABA cell; GabaA; AMPA; NMDA; Gaba; Glutamate;
/
stdpscalingpaper
mod
infot.mod
intf6.mod
intfsw.mod *
misc.mod *
myfft.mod *
nstim.mod *
place.mod *
sampen.mod *
staley.mod *
stats.mod *
tsa.mod *
updown.mod *
vecst.mod *
bpf.h *
misc.h *
mkmod *
parameters.multi *
                            
: $Id: tsa.mod,v 1.18 2011/10/25 21:45:39 billl Exp $ 


NEURON {
 SUFFIX tsa
 GLOBAL INSTALLED
 GLOBAL verbose
}

PARAMETER {
  INSTALLED=0
  verbose=0
}

VERBATIM
#include "misc.h"
static double *x1x, *y1y, *z1z;
//void spctrm(double data[], double p[], int m, int k);
extern double dfftpow(double* x,int n,double* ppow,int powlen,int* fftlen);

#define WINDOW(j,a,b) (1.0-fabs((((j)-1)-(a))*(b)))       /* Bartlett */

double *nrvector(long nl, long nh)
/* allocate a double vector with subscript range v[nl..nh] */
{
	double *v;

	v=(double *)malloc((size_t) ((nh-nl+1+1)*sizeof(double)));
	if (!v) { nrerror("allocation failure in vector()"); return 0x0;}
	return v-nl+1;
}

void nrfree_vector(double *v, long nl, long nh)
/* free a double vector allocated with vector() */
{
  free((char*) (v+nl-1));
}

#define SWAP(a,b) tempr=(a);(a)=(b);(b)=tempr

static double nrsqrarg;
#define SQR(a) ((nrsqrarg=(a)) == 0.0 ? 0.0 : nrsqrarg*nrsqrarg)


void nrfour1(double mdata[], unsigned long nn, int isign)
{
	unsigned long n,mmax,m,j,istep,i;
	double wtemp,wr,wpr,wpi,wi,theta;
	double tempr,tempi;

	n=nn << 1;
	j=1;
	for (i=1;i<n;i+=2) {
		if (j > i) {
			SWAP(mdata[j],mdata[i]);
			SWAP(mdata[j+1],mdata[i+1]);
		}
		m=n >> 1;
		while (m >= 2 && j > m) {
			j -= m;
			m >>= 1;
		}
		j += m;
	}
	mmax=2;
	while (n > mmax) {
		istep=mmax << 1;
		theta=isign*(6.28318530717959/mmax);
		wtemp=sin(0.5*theta);
		wpr = -2.0*wtemp*wtemp;
		wpi=sin(theta);
		wr=1.0;
		wi=0.0;
		for (m=1;m<mmax;m+=2) {
			for (i=m;i<=n;i+=istep) {
				j=i+mmax;
				tempr=wr*mdata[j]-wi*mdata[j+1];
				tempi=wr*mdata[j+1]+wi*mdata[j];
				mdata[j]=mdata[i]-tempr;
				mdata[j+1]=mdata[i+1]-tempi;
				mdata[i] += tempr;
				mdata[i+1] += tempi;
			}
			wr=(wtemp=wr)*wpr-wi*wpi+wr;
			wi=wi*wpr+wtemp*wpi+wi;
		}
		mmax=istep;
	}
}

void nrspctrm(double mdata[], double p[], int m, int k)
{
  // assume overlap = 1
	void nrfour1(double mdata[], unsigned long nn, int isign);
	int mm,m44,m43,m4,kk,joffn,joff,j2,j,c=0;
	double w,facp,facm,*w1,*w2,sumw=0.0,den=0.0;

	mm=m+m;
	m43=(m4=mm+mm)+3;
	m44=m43+1;
	w1=nrvector(1,m4);
	w2=nrvector(1,m);
	facm=m;
	facp=1.0/m;
	for (j=1;j<=mm;j++) sumw += SQR(WINDOW(j,facm,facp));
	for (j=1;j<=m;j++) p[j]=0.0;
	for (j=1;j<=m;j++) w2[j] = mdata[c++];
	for (kk=1;kk<=k;kk++) {
		for (joff = -1;joff<=0;joff++) {
			for (j=1;j<=m;j++) w1[joff+j+j]=w2[j];
			for (j=1;j<=m;j++) w2[j] = mdata[c++];
			joffn=joff+mm;
			for (j=1;j<=m;j++) w1[joffn+j+j]=w2[j];
		}
		for (j=1;j<=mm;j++) {
			j2=j+j;
			w=WINDOW(j,facm,facp);
			w1[j2] *= w;
			w1[j2-1] *= w;
		}
		nrfour1(w1,mm,1);
		p[1] += (SQR(w1[1])+SQR(w1[2]));
		for (j=2;j<=m;j++) {
			j2=j+j;
			p[j] += (SQR(w1[j2])+SQR(w1[j2-1])
				+SQR(w1[m44-j2])+SQR(w1[m43-j2]));
		}
		den += sumw;
	}
	den *= m4;
	for (j=1;j<=m;j++) p[j] /= den;
	nrfree_vector(w2,1,m);
	nrfree_vector(w1,1,m4);
}

double mymax(double x, double y){
  return x > y ? x : y;
}

/* ********************************************************************* */

void mysvd(int m, int n, double** u, double w[], double** v, int* ierr)
/*
 *   This subroutine is a translation of the Algol procedure svd,
 *   Num. Math. 14, 403-420(1970) by Golub and Reinsch.
 *   Handbook for Auto. Comp., Vol II-Linear Algebra, 134-151(1971).
 *
 *   This subroutine determines the singular value decomposition
 *        t
 *   A=usv  of a real m by n rectangular matrix, where m is greater
 *   then or equal to n.  Householder bidiagonalization and a variant
 *   of the QR algorithm are used.
 *  
 *
 *   On input.
 *
 *      m is the number of rows of A (and u).
 *
 *      n is the number of columns of A (and u) and the order of v.
 *
 *      u contains the rectangular input matrix A to be decomposed.
 *
 *   On output.
 *
 *      w contains the n (non-negative) singular values of a (the
 *        diagonal elements of s).  they are unordered.  if an
 *        error exit is made, the singular values should be correct
 *        for indices ierr+1,ierr+2,...,n.
 *
 *      u contains the matrix u (orthogonal column vectors) of the
 *        decomposition.
 *        if an error exit is made, the columns of u corresponding
 *        to indices of correct singular values should be correct.
 *
 *      v contains the matrix v (orthogonal) of the decomposition.
 *        if an error exit is made, the columns of v corresponding
 *        to indices of correct singular values should be correct.
 *
 *      ierr is set to
 *        zero       for normal return,
 *        k          if the k-th singular value has not been
 *                   determined after 30 iterations,
 *        -1         if memory allocation fails.
 *
 *   Questions and comments should be directed to B. S. Garbow,
 *   Applied Mathematics division, Argonne National Laboratory
 *
 *   Modified to eliminate machep
 *
 *   Translated to C by Michiel de Hoon, Human Genome Center,
 *   University of Tokyo, for inclusion in the C Clustering Library.
 *   This routine is less general than the original svd routine, as
 *   it focuses on the singular value decomposition as needed for
 *   clustering. In particular,
 *     - We require m >= n
 *     - We calculate both u and v in all cases
 *     - We pass the input array A via u; this array is subsequently
 *       overwritten.
 *     - We allocate for the array rv1, used as a working space,
 *       internally in this routine, instead of passing it as an
 *       argument. If the allocation fails, svd sets *ierr to -1
 *       and returns.
 *   2003.06.05
 */
{ int i, j, k, i1, k1, l1, its;
  double c,f,h,s,x,y,z;
  int l = 0;
  double g = 0.0;
  double scale = 0.0;
  double anorm = 0.0;
  double* rv1 = (double*)malloc(n*sizeof(double));
  if (!rv1)
  { *ierr = -1;
    return;
  }
  *ierr = 0;
  /* Householder reduction to bidiagonal form */
  for (i = 0; i < n; i++)
  { l = i + 1;
    rv1[i] = scale * g;
    g = 0.0;
    s = 0.0;
    scale = 0.0;
    for (k = i; k < m; k++) scale += fabs(u[k][i]);
    if (scale != 0.0)
    { for (k = i; k < m; k++)
      { u[k][i] /= scale;
        s += u[k][i]*u[k][i];
      }
      f = u[i][i];
      g = (f >= 0) ? -sqrt(s) : sqrt(s);
      h = f * g - s;
      u[i][i] = f - g;
      if (i < n-1)
      { for (j = l; j < n; j++)
        { s = 0.0;
          for (k = i; k < m; k++) s += u[k][i] * u[k][j];
          f = s / h;
          for (k = i; k < m; k++) u[k][j] += f * u[k][i];
        }
      }
      for (k = i; k < m; k++) u[k][i] *= scale;
    }
    w[i] = scale * g;
    g = 0.0;
    s = 0.0;
    scale = 0.0;
    if (i<n-1)
    { for (k = l; k < n; k++) scale += fabs(u[i][k]);
      if (scale != 0.0)
      { for (k = l; k < n; k++)
        { u[i][k] /= scale;
          s += u[i][k] * u[i][k];
        }
        f = u[i][l];
        g = (f >= 0) ? -sqrt(s) : sqrt(s);
        h = f * g - s;
        u[i][l] = f - g;
        for (k = l; k < n; k++) rv1[k] = u[i][k] / h;
        for (j = l; j < m; j++)
        { s = 0.0;
          for (k = l; k < n; k++) s += u[j][k] * u[i][k];
          for (k = l; k < n; k++) u[j][k] += s * rv1[k];
        }
        for (k = l; k < n; k++)  u[i][k] *= scale;
      }
    }
    anorm = mymax(anorm,fabs(w[i])+fabs(rv1[i]));
  }
  /* accumulation of right-hand transformations */
  for (i = n-1; i>=0; i--)
  { if (i < n-1)
    { if (g != 0.0)
      { for (j = l; j < n; j++) v[j][i] = (u[i][j] / u[i][l]) / g;
        /* double division avoids possible underflow */
        for (j = l; j < n; j++)
        { s = 0.0;
          for (k = l; k < n; k++) s += u[i][k] * v[k][j];
          for (k = l; k < n; k++) v[k][j] += s * v[k][i];
        }
      }
    }
    for (j = l; j < n; j++)
    { v[i][j] = 0.0;
      v[j][i] = 0.0;
    }
    v[i][i] = 1.0;
    g = rv1[i];
    l = i;
  }
  /* accumulation of left-hand transformations */
  for (i = n-1; i >= 0; i--)
  { l = i + 1;
    g = w[i];
    if (i!=n-1)
      for (j = l; j < n; j++) u[i][j] = 0.0;
    if (g!=0.0)
    { if (i!=n-1)
      { for (j = l; j < n; j++)
        { s = 0.0;
          for (k = l; k < m; k++) s += u[k][i] * u[k][j];
          /* double division avoids possible underflow */
          f = (s / u[i][i]) / g;
          for (k = i; k < m; k++) u[k][j] += f * u[k][i];
        }
      }
      for (j = i; j < m; j++) u[j][i] /= g;
    }
    else
      for (j = i; j < m; j++) u[j][i] = 0.0;
    u[i][i] += 1.0;
  }
  /* diagonalization of the bidiagonal form */
  for (k = n-1; k >= 0; k--)
  { k1 = k-1;
    its = 0;
    while(1)
    /* test for splitting */
    { for (l = k; l >= 0; l--)
      { l1 = l-1;
        if (fabs(rv1[l]) + anorm == anorm) break;
        /* rv1[0] is always zero, so there is no exit
         * through the bottom of the loop */
        if (fabs(w[l1]) + anorm == anorm)
        /* cancellation of rv1[l] if l greater than 0 */
        { c = 0.0;
          s = 1.0;
          for (i = l; i <= k; i++)
          { f = s * rv1[i];
            rv1[i] *= c;
            if (fabs(f) + anorm == anorm) break;
            g = w[i];
            h = sqrt(f*f+g*g);
            w[i] = h;
            c = g / h;
            s = -f / h;
            for (j = 0; j < m; j++)
            { y = u[j][l1];
              z = u[j][i];
              u[j][l1] = y * c + z * s;
              u[j][i] = -y * s + z * c;
            }
          }
          break;
        }
      }
      /* test for convergence */
      z = w[k];
      if (l==k) /* convergence */
      { if (z < 0.0)
        /*  w[k] is made non-negative */
        { w[k] = -z;
          for (j = 0; j < n; j++) v[j][k] = -v[j][k];
        }
        break;
      }
      else if (its==30)
      { *ierr = k;
        break;
      }
      else
      /* shift from bottom 2 by 2 minor */
      { its++;
        x = w[l];
        y = w[k1];
        g = rv1[k1];
        h = rv1[k];
        f = ((y - z) * (y + z) + (g - h) * (g + h)) / (2.0 * h * y);
        g = sqrt(f*f+1.0);
        f = ((x - z) * (x + z) + h * (y / (f + (f >= 0 ? g : -g)) - h)) / x;
        /* next qr transformation */
        c = 1.0;
        s = 1.0;
        for (i1 = l; i1 <= k1; i1++)
        { i = i1 + 1;
          g = rv1[i];
          y = w[i];
          h = s * g;
          g = c * g;
          z = sqrt(f*f+h*h);
          rv1[i1] = z;
          c = f / z;
          s = h / z;
          f = x * c + g * s;
          g = -x * s + g * c;
          h = y * s;
          y = y * c;
          for (j = 0; j < n; j++)
          { x = v[j][i1];
            z = v[j][i];
            v[j][i1] = x * c + z * s;
            v[j][i] = -x * s + z * c;
          }
          z = sqrt(f*f+h*h);
          w[i1] = z;
          /* rotation can be arbitrary if z is zero */
          if (z!=0.0)
          { c = f / z;
            s = h / z;
          }
          f = c * g + s * y;
          x = -s * g + c * y;
          for (j = 0; j < m; j++)
          { y = u[j][i1];
            z = u[j][i];
            u[j][i1] = y * c + z * s;
            u[j][i] = -y * s + z * c;
          }
        }
        rv1[l] = 0.0;
        rv1[k] = f;
        w[k] = x;
      }
    }
  }
  free(rv1);
  return;
}

int mycompare( const void* a, const void* b ) {
  double* arg1 = (double*) a;
  double* arg2 = (double*) b;
  if( *arg1 < *arg2 ) return -1;
  else if( *arg1 == *arg2 ) return 0;
  else return 1;
}  

double* myspectrum(double* v1,int dc,int* anslen){
  // n data pts will be divided into k partitions of size m
  // the spectrum will have m values from 0 to m/2 cycles/dt.
  //  n = (2*k+1)*m

  // data set
  //  Vect* v1 = vector_arg(1);

  //  int dc = v1->capacity();
  int mr;
  //  if (ifarg(2)) mr = (int)(*getarg(2)); else mr = dc/8;
  mr = dc / 8;

  // make sure the length of partitions is integer power of 2
  int m = 1;
  while(m<mr) m*=2;
  *anslen = m;

  int k = (int) (ceil(((double)(dc)/m-1.)/2.));
  int n = (2*k+1)*m;

  int i;

  double *mdata = (double *)calloc(n,(unsigned)sizeof(double));
  for (i=0;i<dc;++i) mdata[i] = v1[i]; //v1->elem(i);

  //if (ans->capacity() < m) ans->resize(m);
  double* ans = (double *)calloc(m,(unsigned)sizeof(double));
  //spctrm(&mdata[0], &ans->elem(0)-1, m, k);
  nrspctrm(&mdata[0], &ans[0], m, k);

  free((char *)mdata);

  return ans;
}

double myspct(void* v){
  double* x;
  int n = vector_instance_px(v,&x) , anslen = 0 , i = 0;
  double* ans = myspectrum(x,n,&anslen);
  for(i=0;i<anslen;i++) printf("%g ",ans[i]);
  printf("\n");
  double* p;
  int outlen = vector_arg_px(1,&p);
  for(i=0;i<outlen && i<anslen;i++) p[i]=ans[i];
  free(ans);
  return 1.0;
}

int iinrange(double* x,int sz,double dmin,double dmax){
  int i = 0, cnt = 0;
  for(i=0;i<sz;i++) if(x[i]>=dmin && x[i]<=dmax) cnt++;
  return cnt;
}

double inrange(void* v){
  double* x , dmin, dmax;
  int sz = vector_instance_px(v,&x) , i , cnt = 0;
  return (double) iinrange(x,sz,*getarg(1),*getarg(2));
}

extern double dfftpow(double* x,int n,double* ppow,int powlen,int* fftlen);

ENDVERBATIM

: get averagecorrelations between channels, but first bandpass filter them
: tscoravgband(output vector,vec-list of eeg time series,window size,increment,vec:60/120Hz power,vec:saturation,
: vec:delta power,vec:theta,vec:beta,vec:gamma,vec:ripple,sampling rate,vec:low-cutoffs,vec:high-cutoffs, [optional
: output list of vectors for filtered time-series]
:  --
: low and high cutoffs are different bands that will be kept, so can do theta+gamma, i.e. lo=4,31 , hi=12,100
: or can just do theta : lo=4 , hi=12
FUNCTION tscoravgband () {
  VERBATIM

  int i;

  double dRet = 0.0;
  ListVec* pTS = 0;

  double* pcvin = 0;
  double** pfilter = 0x0;

  double** pTSN = 0; // normalized time series

  double* vTCor = 0; //correlation matrix abs sum avg for a time point
  int corsz = vector_arg_px(1,&vTCor);

  double* pTMP = 0x0;

  ListVec* pTSOut = 0;

  pTS = AllocListVec(*hoc_objgetarg(2)); //list of eeg time series
  
  if(!pTS){
    printf("tscoravg ERRA: problem initializing 1st 2 args!\n");
    goto CLEANUP;
  }
  int iChans = pTS->isz; // # of channels
  if(iChans < 2) {
    printf("tscoravg ERRB: must have at least 2 EEG channels!\n");
    goto CLEANUP;
  }
  int iWinSz = (int) *getarg(3); //window size for analysis
  int iSamples = pTS->plen[0]; //make sure all time series have same length
  int iINC = ifarg(4)?(int)*getarg(4):1; 
  double *phz = 0, *psat = 0, *pfft = 0; int sz = 0 , iNoiseTest = 0;
  if(ifarg(5) && ifarg(6)){ //do saturation and 60/120 Hz Frequency test
    if((sz=vector_arg_px(5,&phz))!=corsz){
      printf("tscoravg ERRG: invalid size for phz should be %d, is %d!\n",corsz,sz);
      goto CLEANUP;
    } else if((sz=vector_arg_px(6,&psat))!=corsz){
      printf("tscoravg ERRH: invalid size for psat should be %d, is %d!\n",corsz,sz);
      goto CLEANUP;
    } else {	      
      iNoiseTest = 1;
      pfft = (double*)calloc(iWinSz+1,sizeof(double));
    }
  }
  int iOutSz = (iSamples - iWinSz + 1) / iINC;//output List size

  double *pdelta=0,*ptheta=0,*pbeta=0,*pgamma=0,*pripple=0;
  if(ifarg(7) && vector_arg_px(7,&pdelta)<iOutSz){ printf("bad pdelta size\n"); goto CLEANUP; } //delta 0-3
  if(ifarg(8) && vector_arg_px(8,&ptheta)<iOutSz){ printf("bad ptheta size\n"); goto CLEANUP; } //theta 4-12
  if(ifarg(9)&& vector_arg_px(9,&pbeta)<iOutSz){ printf("bad pbeta size\n"); goto CLEANUP; } //beta 13-29
  if(ifarg(10)&& vector_arg_px(10,&pgamma)<iOutSz){ printf("bad pgamma size\n"); goto CLEANUP; } //gamma 30-100
  if(ifarg(11)&& vector_arg_px(11,&pripple)<iOutSz){ printf("bad pripple size\n"); goto CLEANUP; } //ripple 101-300

  double sampr = *getarg(12);

  double *plohz=0x0,*phihz=0x0;
  int iFilters = 1;
  if((iFilters=vector_arg_px(13,&plohz))!=(i=vector_arg_px(14,&phihz))) {
    printf("invalid arg 13,14 lohz/hihz filter sizes: %d %d\n",iFilters,i); goto CLEANUP;
  } else if(iFilters<1) { printf("no filters specified!\n"); goto CLEANUP; }
  int N = 1, M = 1025;
  while(N < iWinSz) N*=2;
  if(N>M) M=N+1;
  pfilter = getdouble2D(iFilters,N);
  //normalize low/high cutoff frequencies by sampling rate (should be between 0 and 0.5)
  for(i=0;i<iFilters;i++){ plohz[i]/=sampr; phihz[i]/=sampr; }
  for(i=0;i<iFilters;i++){ wsfirBP(pfilter[i], M, 1, plohz[i], phihz[i], N); wrap(pfilter[i],N,M); }
  pcvin = (double*) calloc(N,sizeof(double));

  if(verbose) printf("iWinSz=%d iINC=%d iSamples=%d\n",iWinSz,iINC,iSamples);
  for(i=1;i<iChans;i++) {
    if(pTS->plen[i] != iSamples){
      printf("tscoravg ERRC: time series of unequal size %d %d!\n",iSamples,pTS->plen[i]);
      goto CLEANUP;
    }
  }

  if(corsz < iOutSz){
    printf("tscoravg ERRD: need output vector of at least %d as arg 1, have only %d!\n",iOutSz,corsz);
    goto CLEANUP;
  }
  if(verbose) printf("iOutSz=%d\n",iOutSz);

  pTSN = getdouble2D(iChans,2*N); if(verbose) printf("got pTSN\n");
  if(!pTSN){
    printf("tscor ERRD: out of memory!\n");
    goto CLEANUP;
  }

  if(ifarg(15) && (pTSOut = AllocListVec(*hoc_objgetarg(15)))) { //optional - copy filtered output
    if(pTSOut->isz!=iChans) { printf("incorrect output vec-list size: %d %d\n",pTSOut->isz,iChans); goto CLEANUP; }
    for(i=0;i<iChans;i++) if(pTSOut->plen[i] < iSamples) { printf("incorrect output vec-size %d %d\n",pTSOut->plen[i],iSamples);
      goto CLEANUP; }}

  pTMP = (double*) calloc(2*N,sizeof(double));

  double* psums = (double*)calloc(iChans,sizeof(double)); //stores sum of values in channel's time window
  double* psums2 = (double*)calloc(iChans,sizeof(double));//stores sum of squared values in channel's time window

  double dVal = 0.0;
  int iStartIDX = 0 , ierr = 0 , iOutIDX = 0 , iFilt = 0;
  for(iOutIDX=0;iOutIDX<corsz;iOutIDX++) vTCor[iOutIDX] = 0.;
  if( phz ) for(iOutIDX=0;iOutIDX<corsz;iOutIDX++) phz[iOutIDX]=0.;
  if( psat ) for(iOutIDX=0;iOutIDX<corsz;iOutIDX++) psat[iOutIDX]=0.;
  iOutIDX=0;
  double CP = iChans*(iChans-1.0)/2.0;
  for(iStartIDX=0;iOutIDX<iOutSz && iStartIDX+iWinSz<iSamples;iStartIDX+=iINC,iOutIDX++){
    if(iOutIDX%1000==0) printf("%d\n",iOutIDX);
    int IDX = iStartIDX , iChan = 0 , iEndIDX = iStartIDX + iWinSz;
    int iCurSz = iEndIDX - iStartIDX;
    if(phz) phz[iOutIDX]=0.; if(psat)psat[iOutIDX]=0.;
    if(pdelta) pdelta[iOutIDX]=0.; if(ptheta) ptheta[iOutIDX]=0.;
    if(pbeta) pbeta[iOutIDX]=0.; if(pgamma) pgamma[iOutIDX]=0.; if(pripple) pripple[iOutIDX]=0.;

    double dsatavg = 0., dhzavg = 0.;
    for(iChan=0;iChan<iChans;iChan++){       
      i=0;      
      int iSat = 0; 
      double* pc = &pTS->pv[iChan][iStartIDX]; //saturation test
      for(IDX=iStartIDX;IDX<iEndIDX;IDX++,pc++)if((*pc>=900. && *pc<=1010.) || (*pc>=-1010. && *pc<=-900.)) iSat++;
      dsatavg += (double)iSat/ (double)iWinSz;
      int fftlen = 0; memset(pfft,0,sizeof(double)*(iWinSz+1));
      if(dfftpow(&pTS->pv[iChan][iStartIDX],iWinSz,pfft,iWinSz+1,&fftlen)){       //60,120 Hz Frequency test
	int f = 0; pfft[0]=0.0; double fsum = 0.0; double* ppf = pfft;
	for(f=0;f<fftlen && f<=sampr/2;f++) fsum += *ppf++; //normalize power btwn 0-sampr/2Hz to 1.0
	for(f=55;f<=65 && f<fftlen;f++) dhzavg += pfft[f] / fsum; //60Hz contribution
	for(f=115;f<=125 && f<fftlen;f++) dhzavg += pfft[f] / fsum; //120Hz contribution
	if(pdelta) {
	  fsum=0.0;
	  for(f=0;f<fftlen;f++) fsum+=pfft[f];
	  for(f=0;f<=3;f++) pdelta[iOutIDX] += pfft[f]/fsum;
	  if(ptheta) for(f=4;f<=12;f++) ptheta[iOutIDX] += pfft[f]/fsum;
	  if(pbeta) for(f=13;f<=29;f++) pbeta[iOutIDX] += pfft[f]/fsum;
	  if(pgamma) for(f=30;f<=100;f++) pgamma[iOutIDX] += pfft[f]/fsum;
	  if(pripple) for(f=101;f<=300 && f<fftlen;f++) pripple[iOutIDX] += pfft[f]/fsum;
	}
      }	
    }
    phz[iOutIDX] = dhzavg / (double)iChans;
    psat[iOutIDX] = dsatavg / (double)iChans;
    if(pdelta) pdelta[iOutIDX] /= (double)iChans;
    if(ptheta) ptheta[iOutIDX] /= (double)iChans;
    if(pbeta) pbeta[iOutIDX] /= (double)iChans;
    if(pgamma) pgamma[iOutIDX] /= (double)iChans;
    if(pripple) pripple[iOutIDX] /= (double)iChans;

    //bandpass filter and normalize
    for(iChan=0;iChan<iChans;iChan++){
      memset(pTMP,0,sizeof(double)*2*N);
      double* pC = &pTS->pv[iChan][iStartIDX], *pIN = pcvin;
      if(0) for(IDX=iStartIDX;IDX<iEndIDX;IDX++) *pIN++ = *pC++; //copy input
      for(iFilt=0;iFilt<iFilters;iFilt++) {        
        if(1) {
          memset(pTSN[iChan],0,sizeof(double)*N);//zero it out first...
          pIN = pcvin; pC = &pTS->pv[iChan][iStartIDX];
          for(IDX=iStartIDX;IDX<iEndIDX;IDX++) *pIN++ = *pC++; //copy input
        }
        convlv(pcvin-1,N,pfilter[iFilt]-1,M,1,pTSN[iChan]-1);//do filtering
        for(IDX=0;IDX<N;IDX++) pTMP[IDX] += pTSN[iChan][IDX]; 
      }
      for(IDX=0;IDX<N;IDX++) pTSN[iChan][IDX] = pTMP[IDX]; //copy filtered output to pTSN
      if(pTSOut)for(IDX=0;IDX<N && IDX+iStartIDX<pTSOut->plen[iChan];IDX++) pTSOut->pv[iChan][iStartIDX+IDX]=pTSN[iChan][IDX];

      psums[iChan]=psums2[iChan]=0.;
      pC = pTSN[iChan]; //pointer to start of channel
      for(i=0;i<N;i++,pC++){ 
	psums[iChan] += *pC;//for average
	psums2[iChan] += (*pC * *pC); //for std-dev
      }
      double davg = psums[iChan] / N;
      double dstdev = psums2[iChan]/N - davg*davg; 
      if(dstdev > 0. ) dstdev = sqrt(dstdev); else dstdev=1.; //make sure no problems with nan
      if(dstdev <= 0.) dstdev = 1.0;
      pC = pTSN[iChan]; //pointer to start of channel
      for(i=0;i<N;i++,pC++) *pC = (*pC - davg)/dstdev; //normalization
    }
    //get average of pairwise correlations across channels
    int iC1,iC2;
    double dpsum = 0.0;
    for(iC1=0;iC1<iChans;iC1++){
      for(iC2=0;iC2<iC1;iC2++){
        double r = 0.0, *p1 = pTSN[iC1], *p2 = pTSN[iC2];
	for(i=0;i<N;i++) r += (*p1++ * *p2++);
        r /= (double) N;
        dpsum += r;
      }
    }
    dpsum /= CP;
    if(dpsum > 1.0 || dpsum < -1.0) printf("out of bounds = %g!\n",dpsum);
    vTCor[iOutIDX] = dpsum;
  }
  dRet = 1.0;
  printf("\n");

CLEANUP:
  if(pTS) FreeListVec(&pTS);                       if(pTSN) freedouble2D(&pTSN,iChans);
  free(psums);                                     free(psums2); 
  if(pfft) free(pfft);
  if(pfilter) freedouble2D(&pfilter,iFilters);     if(pcvin) free(pcvin);
  if(pTMP) free(pTMP);                             if(pTSOut) FreeListVec(&pTSOut);
  return dRet;
  ENDVERBATIM
}


: normalize values in list of vectors
FUNCTION tsnorm () {
  VERBATIM
  double dRet = 0.0;
  ListVec* pEIG = AllocListVec(*hoc_objgetarg(1));
  if(!pEIG || pEIG->isz < 1 || pEIG->plen[0]<1){
    printf("tsnorm ERRA: problem initializing 1st arg!\n");
    goto TSNCLEANUP;
  }
  double* dmean = 0, *dstdev = 0;
  int cols = vector_arg_px(2,&dmean);
  if(cols!=vector_arg_px(3,&dstdev) || cols!=pEIG->plen[0]){
    printf("tsnorm ERRB: problem initializing 2nd,3rd arg!\n");
    goto TSNCLEANUP;
  }
  int i,j;
  double val;
  for(i=0;i<pEIG->isz;i++)for(j=0;j<cols;j++)pEIG->pv[i][j]=(pEIG->pv[i][j]-dmean[j])/dstdev[j];
  dRet = 1.0;
TSNCLEANUP:
  if(pEIG) FreeListVec(&pEIG);
  return dRet;
  ENDVERBATIM
}

VERBATIM
int ddiff(double* pin,double* pout,int sz){
  int i;
  for(i=0;i<sz-1;i++) pout[i]=pin[i+1]-pin[i];
  return 1;
}
ENDVERBATIM

: get averages from tscorfull into Vectors
: tsgetmeans( 1 cor vec list, 2 phzlist, 3 psatlist, 4 pderivlist, 5 avg cor vec A, 6 avg cor vec B, 7 avg cor vec C, 8 phzt, 9 psatt, 10 pderivt , 11 numchannels , 12 phza factor, 13 vfid, 14 fid, 15 vprctbadch)
: arg 5 will have average correlation without getting rid of any channels on any epochs
: arg 6 will have average correlation with getting rid of saturated and flat channels
: arg 7 will have avg. cor. with getting rid of saturated,flat channels and channels with 60/120Hz power > phzt
: pfctr not used currently
: 13 vector of file ids so can selectively clean up particular files
: 14 file id to clean up -- iff == -1, do all
: 15 % of 'bad' channels for each time period -- channels excluded for "avgC" / # of channels
FUNCTION tsgetmeans () {
VERBATIM
  int iC1,iC2,iChans,idx,jdx;
  double dRet = 0.0;
  ListVec *plvcor = 0, *plvhz = 0, *plvsat = 0, *plvderiv = 0;
  double* pavgA=0,*pavgB=0,*pavgC=0,*pprctbadch=0,*pfid=0,*pbadc=0;
  int vsz = 0, fid=-1; 
  ////////////////////////////////////////////////////////////////////////////////////////////////
  //
  //                                     set up input args.
  //
  if(!(plvcor = AllocListVec(*hoc_objgetarg(1)))){
    printf("tsgetmeans ERRA: couldn't get cor vec list!\n");
    goto MCLEANUP;
  }
  if(!(plvhz = AllocListVec(*hoc_objgetarg(2)))){
    printf("tsgetmeans ERRB: couldn't get phz vec list!\n");
    goto MCLEANUP;
  }
  if(plvhz->isz < plvcor->isz){
    printf("tsgetmeans ERRC: phz vec list size < cor vec list size : %d %d\n",plvhz->isz,plvcor->isz);
    goto MCLEANUP;
  }
  if(!(plvsat = AllocListVec(*hoc_objgetarg(3)))){
    printf("tsgetmeans ERRD: couldn't get psat vec list!\n");
    goto MCLEANUP;
  }
  if(plvsat->isz < plvcor->isz){
    printf("tsgetmeans ERRE: psat vec list size < cor vec list size : %d %d\n",plvsat->isz,plvcor->isz);
    goto MCLEANUP;
  }  
  if(!(plvderiv = AllocListVec(*hoc_objgetarg(4)))){
    printf("tsgetmeans ERRF: couldn't get pderiv vec list!\n");
    goto MCLEANUP;
  }
  if(plvderiv->isz < plvcor->isz){
    printf("tsgetmeans ERRG: pderiv vec list size < cor vec list size : %d %d\n",plvderiv->isz,plvcor->isz);
    goto MCLEANUP;
  }  
  if(vector_arg_px(5,&pavgA) < plvcor->isz || 
     vector_arg_px(6,&pavgB) < plvcor->isz || 
     vector_arg_px(7,&pavgC) < plvcor->isz ||
     vector_arg_px(13,&pfid) < plvcor->isz ||
     vector_arg_px(15,&pprctbadch)< plvcor->isz){
    printf("tsgetmeans ERRH: output vecs must have size >= %d!\n",plvcor->isz);
    goto MCLEANUP;
  }
  double phzt = *getarg(8), psatt = *getarg(9) , pderivt = *getarg(10);
  iChans = (int)*getarg(11);
  double pfctr = *getarg(12);
  fid = (int)*getarg(14);
  pbadc = (double*) malloc(iChans*sizeof(double));
  //
  //
  ////////////////////////////////////////////////////////////////////////////////////////////////

  /////////
  // do the calculations  
  double dSumA,dSumB,dSumC;
  int cntA,cntB,cntC,cntBadCh;
  for(idx=0;idx<plvcor->isz;idx++){
    if(fid>=0 && pfid[idx]!=fid) continue;
    cntA=cntB=cntC=jdx=0;
    dSumA=dSumB=dSumC=0.;
    for(iC1=0;iC1<iChans;iC1++){ //mark bad channels
      if(plvsat->pv[idx][iC1] >= psatt ||
         plvderiv->pv[idx][iC1] >= pderivt)
        pbadc[iC1]=1;
      else
        pbadc[iC1]=0;
    }
    for(iC1=0;iC1<iChans;iC1++){
      if(pbadc[iC1])continue;
      for(iC2=0;iC2<iC1;iC2++,jdx++){
        dSumA += plvcor->pv[idx][jdx];
        cntA++;
        if(pbadc[iC2]) continue;
        dSumB += plvcor->pv[idx][jdx];
        cntB++;
      }
    }
    pavgA[idx] = dSumA / (double) cntA;
    if(cntB) pavgB[idx] = dSumB / (double) cntB; else pavgB[idx] = -666.;
    jdx=0;
    for(iC1=0;iC1<iChans;iC1++) if(plvhz->pv[idx][iC1] >= phzt) pbadc[iC1]=1; //mark bad channels
    for(iC1=0;iC1<iChans;iC1++){
      if(pbadc[iC1]) continue;
      for(iC2=0;iC2<iC1;iC2++,jdx++){
        if(pbadc[iC2]) continue;
        dSumC += plvcor->pv[idx][jdx];
        cntC++;
      }
    }
    if(cntC) pavgC[idx] = dSumC / (double) cntC; else pavgC[idx] = -666.;
    cntBadCh=0;
    for(iC1=0;iC1<iChans;iC1++) if(pbadc[iC1]) cntBadCh++;
    pprctbadch[idx]= (double)cntBadCh/(double)iChans;
  }

  dRet = 1.0;
MCLEANUP:
 if(pbadc) free(pbadc);
 return dRet;
ENDVERBATIM
}

: get full time series correlation matrix for each time window into list of vectors
: tscorfull( 1 cor vec list, 2 time series list, 3 winsz, 4 incsz, 5 phzlist, 6 psatlist, 7 pderivlist)
FUNCTION tscorfull () {
  VERBATIM

  int i;

  double dRet = 0.0;
  ListVec* pTS = 0, *plvcor = 0, *plvhz = 0, *plvsat = 0, *plvderiv = 0;

  double** pTSN = 0; // normalized time series
  double* pdiff = 0; // for straight line test
  double* psums = 0, *psums2 = 0, *pfft = 0;

  int corsz=0; //# of vectors in output correlation vector list
  plvcor = AllocListVec(*hoc_objgetarg(1)); //each vec is correlation mat in vector form for a time point
  if(!plvcor || (corsz=plvcor->isz)<1){
    printf("tscorfull ERRA: problem initializing arg1 plvcor!\n");
    goto FCLEANUP;
  }

  pTS = AllocListVec(*hoc_objgetarg(2)); //list of eeg time series  
  if(!pTS){
    printf("tscorfull ERRA: problem initializing arg2 pTS!\n");
    goto FCLEANUP;
  }

  int iChans = pTS->isz; // # of channels
  if(iChans < 2) {
    printf("tscorfull ERRB: must have at least 2 timeseries!\n");
    goto FCLEANUP;
  }

  int iMatSz = iChans*(iChans-1)/2;

  for(i=0;i<corsz;i++){   //make sure all output vectors have right size
    if(plvcor->plen[i]<iMatSz){
      printf("tscorfull ERRC: each cor vec must have sz >= %d, [%d] has only %d\n",iMatSz,i,plvcor->plen[i]);
      goto FCLEANUP;
    }
  }

  int iWinSz = (int) *getarg(3); //window size for analysis
  int iINC = ifarg(4)?(int)*getarg(4):1;  //increment size for analysis

  plvhz = AllocListVec(*hoc_objgetarg(5));
  plvsat = AllocListVec(*hoc_objgetarg(6));
  plvderiv = AllocListVec(*hoc_objgetarg(7));

  pfft = (double*)calloc(iWinSz+1,sizeof(double));
  pdiff = (double*)calloc(iWinSz,sizeof(double));

  int iSamples = pTS->plen[0]; //make sure all time series have same length
  int iOutSz = (iSamples - iWinSz + 1) / iINC;//output List size

  if(verbose) printf("iWinSz=%d iINC=%d iSamples=%d\n",iWinSz,iINC,iSamples);
  for(i=1;i<iChans;i++) {
    if(pTS->plen[i] != iSamples){
      printf("tscorfull ERRD: time series of unequal size %d %d!\n",iSamples,pTS->plen[i]);
      goto FCLEANUP;
    }
  }

  if(corsz < iOutSz){
    printf("tscorfull ERRE: need output list vector of at least sz %d as arg 1, have only %d!\n",iOutSz,corsz);
    goto FCLEANUP;
  }   if(verbose) printf("iOutSz=%d\n",iOutSz);
  if(plvhz->isz<iOutSz || plvsat->isz<iOutSz || plvderiv->isz<iOutSz){
    printf("tscorfull ERRF: need output plvhz,plvsat,plvderiv size of at least %d!\n",iOutSz);
    goto FCLEANUP;
  }
  for(i=0;i<iOutSz;i++){
    if(plvhz->plen[i]<iChans || plvsat->plen[i]<iChans || plvderiv->plen[i]<iChans){
      printf("tscorfull ERRG: each vec in plvhz,plvsat,plvderiv must be sz >= %d!\n",iChans);
      goto FCLEANUP;
    }
  }

  pTSN = getdouble2D(iChans,iWinSz); if(verbose) printf("got pTSN\n"); //normalized time series
  if(!pTSN){
    printf("tscor ERRD: out of memory!\n");
    goto FCLEANUP;
  }

  psums = (double*)calloc(iChans,sizeof(double));   //for sum(X)
  psums2 = (double*)calloc(iChans,sizeof(double));  //for sum(X^2)
  double* pc = 0;
  double dVal = 0.0;
  int iStartIDX = 0 , ierr = 0 , iOutIDX = 0;
  FillListVec(plvcor,0.0); FillListVec(plvhz,0.); FillListVec(plvsat,0.); FillListVec(plvderiv,0.);//zero results
  for(iStartIDX=0;iOutIDX<iOutSz && iStartIDX+iWinSz<iSamples;iStartIDX+=iINC,iOutIDX++){
    if(iOutIDX%1000==0) printf("%d\n",iOutIDX);
    int IDX = iStartIDX , iChan = 0 , iEndIDX = iStartIDX + iWinSz;
    int iCurSz = iEndIDX - iStartIDX;
    for(iChan=0;iChan<iChans;iChan++){ //get sum(X) , sum(X^2) for ALL channels
      psums[iChan]=psums2[iChan]=0.;
      pc = &pTS->pv[iChan][iStartIDX];
      for(IDX=iStartIDX;IDX<iEndIDX;IDX++,pc++){
        psums[iChan] += *pc; 
        psums2[iChan] += (*pc * *pc);
      }
    }
    for(iChan=0;iChan<iChans;iChan++){       //do noise checks and normalize for ALL channels
      i=0;
      double davg = psums[iChan] / iCurSz;
      double dstdev = psums2[iChan]/iCurSz - davg*davg; 
      if(dstdev > 0. ) dstdev = sqrt(dstdev); else dstdev=1.; //make sure no problems with nan
      if(dstdev <= 0.) dstdev = 1.0;
      //saturation test
      int iSat = 0;
      pc = &pTS->pv[iChan][iStartIDX];
      for(IDX=iStartIDX;IDX<iEndIDX;IDX++,i++,pc++){	
        if((*pc>=900. && *pc<=1010.) || (*pc>=-1010. && *pc<=-900.)) iSat++;
        pTSN[iChan][i] = (*pc - davg) / dstdev;	//normalization
      }
      plvsat->pv[iOutIDX][iChan] = (double)iSat/ (double)iWinSz;
      //60,120 Hz Frequency test
      int fftlen = 0; memset(pfft,0,sizeof(double)*(iWinSz+1));
      if(dfftpow(&pTS->pv[iChan][iStartIDX],iWinSz,pfft,iWinSz+1,&fftlen)){
        int f = 0; pfft[0]=0.0; double fsum = 0.0;
        for(f=0;f<fftlen && f<=1000;f++) fsum += pfft[f]; //normalize power btwn 0-1000Hz to 1.0
        pc = &plvhz->pv[iOutIDX][iChan];
        for(f=55;f<=65 && f<fftlen;f++)   *pc += pfft[f]; //60Hz contribution
        for(f=115;f<=125 && f<fftlen;f++) *pc += pfft[f]; //120Hz contribution
        *pc /= fsum;
      }
      //discrete derivative test to measure how flat of a line signal looks like
      ddiff(&pTS->pv[iChan][iStartIDX],pdiff,iWinSz);
      plvderiv->pv[iOutIDX][iChan] = (double)iinrange(pdiff, iWinSz-1, -1.0, 1.0) / (double)(iWinSz-1.0);
    }

    //form correlation vector
    int iC1,iC2; pc = &plvcor->pv[iOutIDX][0];
    for(iC1=0;iC1<iChans;iC1++){
      for(iC2=iC1+1;iC2<iChans;iC2++){
	for(i=0;i<iWinSz;i++) 
          *pc += pTSN[iC1][i]*pTSN[iC2][i];
        *pc /= (double) iWinSz;
        pc++;
      }
    }
  }
  dRet = 1.0;
  printf("\n");

FCLEANUP:
  if(pTS) FreeListVec(&pTS);
  if(pTSN) freedouble2D(&pTSN,iChans);
  if(psums) free(psums);
  if(psums2) free(psums2);
  if(pfft) free(pfft);
  if(pdiff) free(pdiff);
  return dRet;
  ENDVERBATIM
}

: doesnt use eigenvalues 
FUNCTION tscoravg () {
  VERBATIM

  int i;

  double dRet = 0.0;
  ListVec* pTS = 0;

  double** pTSN = 0; // normalized time series
  double** pCorrel = 0; //correlation matrix

  double* vTCor = 0; //correlation matrix abs sum avg for a time point
  int corsz = vector_arg_px(1,&vTCor);

  pTS = AllocListVec(*hoc_objgetarg(2)); //list of eeg time series
  
  if(!pTS){
    printf("tscoravg ERRA: problem initializing 1st 2 args!\n");
    goto CLEANUP;
  }

  int iChans = pTS->isz; // # of channels
  if(iChans < 2) {
    printf("tscoravg ERRB: must have at least 2 EEG channels!\n");
    goto CLEANUP;
  }
  int iWinSz = (int) *getarg(3); //window size for analysis


  int iSamples = pTS->plen[0]; //make sure all time series have same length

  int iINC = ifarg(4)?(int)*getarg(4):1; 

  double *phz = 0, *psat = 0, *pfft = 0; int sz = 0 , iNoiseTest = 0;
  if(ifarg(5) && ifarg(6)){ //do saturation and 60/120 Hz Frequency test
    if((sz=vector_arg_px(5,&phz))!=corsz){
      printf("tscoravg ERRG: invalid size for phz should be %d, is %d!\n",corsz,sz);
      goto CLEANUP;
    } else if((sz=vector_arg_px(6,&psat))!=corsz){
      printf("tscoravg ERRH: invalid size for psat should be %d, is %d!\n",corsz,sz);
      goto CLEANUP;
    } else {	      
      iNoiseTest = 1;
      pfft = (double*)calloc(iWinSz+1,sizeof(double));
    }
  }

  int idoabs = ifarg(7)?(int)*getarg(7):0; //do abs or not

  int iOutSz = (iSamples - iWinSz + 1) / iINC;//output List size

  double *pdelta=0,*ptheta=0,*pbeta=0,*pgamma=0,*pripple=0;
  if(ifarg(8) && vector_arg_px(8,&pdelta)<iOutSz){ printf("bad pdelta size\n"); goto CLEANUP; } //delta 0-3
  if(ifarg(9) && vector_arg_px(9,&ptheta)<iOutSz){ printf("bad ptheta size\n"); goto CLEANUP; } //theta 4-12
  if(ifarg(10)&& vector_arg_px(10,&pbeta)<iOutSz){ printf("bad pbeta size\n"); goto CLEANUP; } //beta 13-29
  if(ifarg(11)&& vector_arg_px(11,&pgamma)<iOutSz){ printf("bad pgamma size\n"); goto CLEANUP; } //gamma 30-100
  if(ifarg(12)&& vector_arg_px(12,&pripple)<iOutSz){ printf("bad pripple size\n"); goto CLEANUP; } //ripple 101-300

  if(verbose) printf("iWinSz=%d iINC=%d iSamples=%d\n",iWinSz,iINC,iSamples);
  for(i=1;i<iChans;i++) {
    if(pTS->plen[i] != iSamples){
      printf("tscoravg ERRC: time series of unequal size %d %d!\n",iSamples,pTS->plen[i]);
      goto CLEANUP;
    }
  }

  if(corsz < iOutSz){
    printf("tscoravg ERRD: need output vector of at least %d as arg 1, have only %d!\n",iOutSz,corsz);
    goto CLEANUP;
  }
  if(verbose) printf("iOutSz=%d\n",iOutSz);

  pCorrel = getdouble2D(iChans,iChans); if(verbose) printf("got pCorrel\n");
  pTSN = getdouble2D(iChans,iWinSz); if(verbose) printf("got pTSN\n");
  if(!pCorrel || !pTSN){
    printf("tscor ERRD: out of memory!\n");
    goto CLEANUP;
  }

  double* psums = (double*)calloc(iChans,sizeof(double));
  double* psums2 = (double*)calloc(iChans,sizeof(double));

  double dVal = 0.0;
  int iStartIDX = 0 , ierr = 0 , iOutIDX = 0;
  for(iOutIDX=0;iOutIDX<corsz;iOutIDX++) vTCor[iOutIDX] = 0.;
  if( phz ) for(iOutIDX=0;iOutIDX<corsz;iOutIDX++) phz[iOutIDX]=0.;
  if( psat ) for(iOutIDX=0;iOutIDX<corsz;iOutIDX++) psat[iOutIDX]=0.;
  iOutIDX=0;
  for(iStartIDX=0;iOutIDX<iOutSz && iStartIDX+iWinSz<iSamples;iStartIDX+=iINC,iOutIDX++){
    if(iOutIDX%1000==0) printf("%d\n",iOutIDX);
    int IDX = iStartIDX , iChan = 0 , iEndIDX = iStartIDX + iWinSz;
    int iCurSz = iEndIDX - iStartIDX;
    if(phz) phz[iOutIDX]=0.; if(psat)psat[iOutIDX]=0.;
    if(pdelta) pdelta[iOutIDX]=0.; if(ptheta) ptheta[iOutIDX]=0.;
    if(pbeta) pbeta[iOutIDX]=0.; if(pgamma) pgamma[iOutIDX]=0.; if(pripple) pripple[iOutIDX]=0.;
    if(iStartIDX==0 || iINC > 1){
      for(iChan=0;iChan<iChans;iChan++){ //normalize data
	psums[iChan]=psums2[iChan]=0.;
	for(IDX=iStartIDX;IDX<iEndIDX;IDX++){
	  dVal = pTS->pv[iChan][IDX];
	  psums[iChan] += dVal;
	  psums2[iChan] += dVal*dVal;
	}
      }
    } else {
      for(iChan=0;iChan<iChans;iChan++){
	dVal = pTS->pv[iChan][iStartIDX-1];
	psums[iChan] -= dVal;
	psums2[iChan] -= dVal*dVal;
	dVal = pTS->pv[iChan][iEndIDX-1];
	psums[iChan] += dVal;
	psums2[iChan] += dVal*dVal;
      }
    }
    double dsatavg = 0., dhzavg = 0.;
    for(iChan=0;iChan<iChans;iChan++){       
      i=0;
      double davg = psums[iChan] / iCurSz;
      double dstdev = psums2[iChan]/iCurSz - davg*davg; 
      if(dstdev > 0. ) dstdev = sqrt(dstdev); else dstdev=1.; //make sure no problems with nan
      if(dstdev <= 0.) dstdev = 1.0;
      if(iNoiseTest){	
	//saturation test
	int iSat = 0;
	for(IDX=iStartIDX;IDX<iEndIDX;IDX++,i++){	
	  dVal = pTS->pv[iChan][IDX];
	  if((dVal>=900. && dVal<=1010.) || (dVal>=-1010. && dVal<=-900.)) iSat++;
	  pTSN[iChan][i] = (dVal - davg) / dstdev;	
	}
	dsatavg += (double)iSat/ (double)iWinSz;
	//60,120 Hz Frequency test
	int fftlen = 0; memset(pfft,0,sizeof(double)*(iWinSz+1));
	if(dfftpow(&pTS->pv[iChan][iStartIDX],iWinSz,pfft,iWinSz+1,&fftlen)){
	  int f = 0; pfft[0]=0.0; double fsum = 0.0;
	  for(f=0;f<fftlen && f<=1000;f++) fsum += pfft[f]; //normalize power btwn 0-1000Hz to 1.0
	  for(f=55;f<=65 && f<fftlen;f++) dhzavg += pfft[f] / fsum; //60Hz contribution
	  for(f=115;f<=125 && f<fftlen;f++) dhzavg += pfft[f] / fsum; //120Hz contribution
	  if(pdelta) {
	    for(f=0;f<=3;f++) pdelta[iOutIDX] += pfft[f]/fsum;
	    if(ptheta) for(f=4;f<=12;f++) ptheta[iOutIDX] += pfft[f]/fsum;
	    if(pbeta) for(f=13;f<=29;f++) pbeta[iOutIDX] += pfft[f]/fsum;
	    if(pgamma) for(f=30;f<=100;f++) pgamma[iOutIDX] += pfft[f]/fsum;
	    if(pripple) for(f=101;f<=300 && f<fftlen;f++) pripple[iOutIDX] += pfft[f]/fsum;
	  }
	}	
      } else {
	for(IDX=iStartIDX;IDX<iEndIDX;IDX++,i++){	
	  dVal = pTS->pv[iChan][IDX];
	  pTSN[iChan][i] = (dVal - davg) / dstdev;	
	}
      }
    }
    if(iNoiseTest){
      phz[iOutIDX] = dhzavg / (double)iChans;
      psat[iOutIDX] = dsatavg / (double)iChans;
      if(pdelta) pdelta[iOutIDX] /= (double)iChans;
      if(ptheta) ptheta[iOutIDX] /= (double)iChans;
      if(pbeta) pbeta[iOutIDX] /= (double)iChans;
      if(pgamma) pgamma[iOutIDX] /= (double)iChans;
      if(pripple) pripple[iOutIDX] /= (double)iChans;
    }

    //form correlation matrix
    int iC1,iC2;
    for(iC1=0;iC1<iChans;iC1++){
      for(iC2=0;iC2<=iC1;iC2++){
	if(iC1==iC2){
	  pCorrel[iC1][iC2]=1.0; //diagonals == 1.0
	  continue;
	} else pCorrel[iC1][iC2]=0.; //init to 0.
	for(i=0;i<iWinSz;i++){
	  pCorrel[iC1][iC2] += pTSN[iC1][i]*pTSN[iC2][i];
	}
	pCorrel[iC1][iC2] /= (double)iWinSz;
	pCorrel[iC2][iC1] = pCorrel[iC1][iC2];
      }
    }
    double dpsum = 0.0; int icc = 0; //get avg of abs val of correlations
    if(idoabs){
      for(iC1=0;iC1<iChans;iC1++){
	for(iC2=0;iC2<iC1;iC2++){
	  dpsum += fabs(pCorrel[iC1][iC2]);
	  icc++;
	}
      }
    } else {
      for(iC1=0;iC1<iChans;iC1++){
	for(iC2=0;iC2<iC1;iC2++){
	  dpsum += pCorrel[iC1][iC2];
	  icc++;
	}
      }
    }
    dpsum /= (double) icc;
    if(dpsum > 1. || dpsum < -1.){
      printf("out of bounds = %g!\n",dpsum);
    }
    vTCor[iOutIDX] = dpsum;
  }
  dRet = 1.0;
  printf("\n");

CLEANUP:
  if(pTS) FreeListVec(&pTS);
  if(pCorrel) freedouble2D(&pCorrel,iChans);
  if(pTSN) freedouble2D(&pTSN,iChans);
  free(psums);
  free(psums2);
  if(pfft) free(pfft);
  return dRet;
  ENDVERBATIM
}


FUNCTION tscor () {
  VERBATIM

  int i;

  double dRet = 0.0;

  ListVec* pEIG = 0, *pTS = 0;

  double** pTSN = 0; // normalized time series
  double** pCorrel = 0; //correlation matrix
  double** pV = 0;
  double* pW = 0; // eigs

  pEIG = AllocListVec(*hoc_objgetarg(1)); //list of eigenvalue vectors, one for each time
  pTS = AllocListVec(*hoc_objgetarg(2)); //list of eeg time series
  
  if(!pEIG || !pTS){
    printf("tscor ERRA: problem initializing 1st 2 args!\n");
    goto CLEANUP;
  }

  int iChans = pTS->isz; // # of channels
  if(iChans < 2) {
    printf("tscor ERRB: must have at least 2 EEG channels!\n");
    goto CLEANUP;
  }
  int iWinSz = (int) *getarg(3); //window size for analysis
  int iINC = (int) *getarg(4); //increment for analysis

  int iSamples = pTS->plen[0]; //make sure all time series have same length

  if(verbose) printf("iWinSz=%d iINC=%d iSamples=%d\n",iWinSz,iINC,iSamples);
  for(i=1;i<iChans;i++) {
    if(pTS->plen[i] != iSamples){
      printf("tscor ERRC: time series of unequal size %d %d!\n",iSamples,pTS->plen[i]);
      goto CLEANUP;
    }
  }
  int iOutSz = (iSamples - iWinSz + 1) / iINC;//output List size
  if(pEIG->isz < iOutSz){
    printf("tscor ERRD: need List of size at least %d as arg 1, have only %d!\n",iOutSz,pEIG->isz);
    goto CLEANUP;
  }
  if(verbose) printf("iOutSz=%d\n",iOutSz);

  pCorrel = getdouble2D(iChans,iChans); if(verbose) printf("got pCorrel\n");
  pTSN = getdouble2D(iChans,iWinSz); if(verbose) printf("got pTSN\n");
  pV = getdouble2D(iChans,iChans);   if(verbose) printf("got pV\n");
  pW = (double*) malloc(sizeof(double)*iChans); if(verbose) printf("got pW\n");
  if(!pCorrel || !pTSN || !pV || !pW){
    printf("tscor ERRD: out of memory!\n");
    goto CLEANUP;
  }

  double dVal = 0.0;
  int iStartIDX = 0 , ierr = 0 , iTIDX = 0;
  for(iStartIDX=0;iStartIDX+iWinSz<iSamples;iStartIDX+=iINC,iTIDX++){
    if(iStartIDX%100==0) printf("%d\n",iStartIDX);
    int IDX = iStartIDX , iChan = 0 , iEndIDX = iStartIDX + iWinSz;
    int iCurSz = iEndIDX - iStartIDX;
    for(iChan=0;iChan<iChans;iChan++){ //normalize data
      double dSum = 0.0, dSum2 = 0.0;
      for(IDX=iStartIDX;IDX<iEndIDX;IDX++){
        dVal = pTS->pv[iChan][IDX];
	dSum += dVal;
	dSum2 += dVal*dVal;
      }
      dSum /= iCurSz; //mean
      dSum2 /= iCurSz; 
      dSum2 -= dSum*dSum;
      dSum2 = sqrt(dSum2); //standard deviation
      i=0;
      for(IDX=iStartIDX;IDX<iEndIDX;IDX++,i++){
	pTSN[iChan][i] = (pTS->pv[iChan][IDX] - dSum) / dSum2;
      }
    }
    //form correlation matrix
    int iC1,iC2;
    for(iC1=0;iC1<iChans;iC1++){
      for(iC2=0;iC2<=iC1;iC2++){
	if(iC1==iC2){
	  pCorrel[iC1][iC2]=1.0; //diagonals == 1.0
	  continue;
	} else pCorrel[iC1][iC2]=0.;
	for(i=0;i<iWinSz;i++){
	  pCorrel[iC1][iC2] += pTSN[iC1][i]*pTSN[iC2][i];
	}
	pCorrel[iC1][iC2] /= (double)iWinSz;
	pCorrel[iC2][iC1] = pCorrel[iC1][iC2];
      }
    }
    //perform SVD to get eigenvalues into pW
    mysvd(iChans,iChans,pCorrel,pW,pV,&ierr); 
    qsort(pW,iChans,sizeof(double),mycompare); //sort eigenvalues in ascending order
//    for(i=0;i<iChans;i++) pEIG->pv[iTIDX][i] = pW[i]; //store sorted eigenvalues
    memcpy(pEIG->pv[iTIDX],pW,sizeof(double)*iChans); //store sorted eigenvalues
  }
  dRet = 1.0;
  printf("\n");

CLEANUP:
  if(pEIG) FreeListVec(&pEIG);
  if(pTS) FreeListVec(&pTS);
  if(pCorrel) freedouble2D(&pCorrel,iChans);
  if(pTSN) freedouble2D(&pTSN,iChans);
  if(pW) free(pW);
  if(pV) freedouble2D(&pV,iChans);
  return dRet;
  ENDVERBATIM
}

VERBATIM

double getmean(double* p,int n) {
  if(n<1) return 0.0;
  int i = 0;
  double sum = 0.0, *pp=p;
  for(i=0;i<n;i++,pp++) sum += *pp;
  return sum / (double) n;
}

double getstd(double* p,int n,double X) {
  if(n<1) return 0.;
  int i;
  double X2=0., *pp=p;
  for(i=0;i<n;i++,pp++) X2 += (*pp * *pp);
  X2 = X2/(double)n - X*X;
  if(X2>0.) return sqrt(X2);
  return 0.;
}

double dnormv(double* p,int n) {
  if(n<1) return 0.0;
  double X = getmean(p,n);
  double s = getstd(p,n,X);
  double* pp = p;
  int i;
  if(s>0.)
    for(i=0;i<n;i++,pp++) *pp = (*pp - X)/s; 
  else
    for(i=0;i<n;i++,pp++) *pp -= X;
  return 1.0;
}

static double normv(void* v){ 
  double* x;
  int sz = vector_instance_px(v,&x);
  if(sz<1){
    printf("normv ERRA: empty input size!\n");
    return 0.;
  }
  return dnormv(x,sz);
}

double dcopynzidx(double* pin,double* pidx,double* pout,int sz) {
  int szout,i;
  szout=0;
  for(i=0;i<sz;i++) if(pidx[i]) pout[szout++]=pin[i];
  return (double) szout;
}

static double copynzidx (void* v) {
  double *x,*y,*z,ret; int sz;
  sz = vector_instance_px(v,&x);
  if(sz!=vector_arg_px(1,&y)) y=vector_newsize(vector_arg(1),sz);
  if(sz!=vector_arg_px(2,&z)) z=vector_newsize(vector_arg(2),sz);
  ret=dcopynzidx(x,y,z,sz);
  vector_resize(vector_arg(2),(int)ret);
  return ret;
}
ENDVERBATIM

: assumes v1,v2 are normalized
FUNCTION tsmul () {
  VERBATIM
  double* x, *y, *px, *py, sum = 0.;
  int xsz,ysz,i;
  xsz = vector_arg_px(1,&x);
  ysz = vector_arg_px(2,&y);
  px=x; 
  py=y;
  if(xsz!=ysz || xsz<1 || ysz<1){
    printf("tsmul ERRB: input vecs are invalid sizes: %d %d!\n",xsz,ysz);
    return -2.;
  }
  for(i=0;i<xsz;i++,px++,py++) sum += *px * *py;
  return sum / (double) xsz;
  ENDVERBATIM
}

PROCEDURE install () {
  if (INSTALLED==1) {
    printf("already installed tsa.mod")
  } else {
    INSTALLED=1
    VERBATIM
      //    install_vector_method("myspct", myspct);
    install_vector_method("inrange",inrange);
    install_vector_method("normv",normv);
    install_vector_method("copynzidx",copynzidx);
    ENDVERBATIM
  }
}

Loading data, please wait...