From 7539856cf3cc82896ae5d614b3c993049c483fcb Mon Sep 17 00:00:00 2001 From: Marianne Bezaire <6456664+mbezaire@users.noreply.github.com> Date: Sat, 30 Jul 2022 19:31:12 -0400 Subject: [PATCH 1/9] add gitignore --- .gitignore | 3 +++ 1 file changed, 3 insertions(+) create mode 100644 .gitignore diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..686d33e --- /dev/null +++ b/.gitignore @@ -0,0 +1,3 @@ +*.c +*.o +*.dll From 0c3706240cf1a3049d87a0c790674150e7845550 Mon Sep 17 00:00:00 2001 From: Marianne Bezaire <6456664+mbezaire@users.noreply.github.com> Date: Sat, 30 Jul 2022 19:31:29 -0400 Subject: [PATCH 2/9] rm sys/time.h includes --- misc.h | 2 +- misc.mod | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/misc.h b/misc.h index 08704be..350efa0 100644 --- a/misc.h +++ b/misc.h @@ -4,7 +4,7 @@ #include #include /* contains LONG_MAX */ #include -#include +//#include #include #include diff --git a/misc.mod b/misc.mod index 9344209..730aa9e 100644 --- a/misc.mod +++ b/misc.mod @@ -46,7 +46,7 @@ VERBATIM #include /* F_OK */ #include /* errno */ #include -#include /* MUST REMEMBER THIS */ +//#include /* MUST REMEMBER THIS */ #include #include #include From 2d00e0d022e048284d825268cd967e66403ac950 Mon Sep 17 00:00:00 2001 From: Marianne Bezaire <6456664+mbezaire@users.noreply.github.com> Date: Fri, 5 Aug 2022 14:47:46 -0400 Subject: [PATCH 3/9] added place to add gabapentin expression --- geom.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/geom.py b/geom.py index ec5ae23..d48fc8e 100644 --- a/geom.py +++ b/geom.py @@ -136,6 +136,8 @@ def set_conductances(self): self.soma.insert('Iholmw') self.soma.insert('Caolmw') self.soma.insert('ICaolmw') + for seg in self.soma: + seg.gca_ICaolmw = 1 # adjust the conductance when gabapentin present self.soma.insert('KCaolmw') def set_synapses(self): From a7fe95668da292e52b75be6125791b690815cf3f Mon Sep 17 00:00:00 2001 From: Marianne Bezaire <6456664+mbezaire@users.noreply.github.com> Date: Fri, 5 Aug 2022 14:55:42 -0400 Subject: [PATCH 4/9] compiles on windows but didn't address integer overflow and unsure about rand changes, etc. --- misc.mod => misc.mignoreod | 4 +- stats.mod | 1932 +++------------------------------- vecst.mod | 2026 +++++++++--------------------------- 3 files changed, 611 insertions(+), 3351 deletions(-) rename misc.mod => misc.mignoreod (99%) diff --git a/misc.mod b/misc.mignoreod similarity index 99% rename from misc.mod rename to misc.mignoreod index 730aa9e..63bbfbe 100644 --- a/misc.mod +++ b/misc.mignoreod @@ -42,12 +42,12 @@ THREADSAFE } VERBATIM -#include "misc.h" +//#include "misc.h" #include /* F_OK */ #include /* errno */ #include //#include /* MUST REMEMBER THIS */ -#include +//#include #include #include ENDVERBATIM diff --git a/stats.mod b/stats.mod index a84a031..31b0d63 100644 --- a/stats.mod +++ b/stats.mod @@ -1,4 +1,4 @@ -: $Id: stats.mod,v 1.215 2010/08/18 19:20:23 samn Exp $ +: $Id: stats.mod,v 1.8 2006/06/30 19:03:36 billl Exp $ :* COMMENT COMMENT @@ -11,56 +11,48 @@ fac not vec related - returns factorial logfac not vec related - returns log factorial vseed set some C level randomizer seeds slope(num) does a linear regression to find the slope, assuming num=timestep of vector -vslope(v2) does a linear regression to find the slope, assuming num=timestep of vector -stats(num,[out]) does a linear regression, assuming num=timestep of vector -vstats(v2,[out]) does a linear regression, using v2 as the x-coords -setrnd(v,flag) does set rand using 1:rand, 2:drand48 -v.hash(list) // make a hash out values in vecs in list -v.unnan([nan_value,][inf_value]) // remove nan's and inf's from a vector +vslope(num) does a linear regression to find the slope, assuming num=timestep of vector +stats(num) does a linear regression, assuming num=timestep of vector +vstats(v2) does a linear regression, using v2 as the x-coords ENDCOMMENT NEURON { -THREADSAFE - SUFFIX stats - GLOBAL INSTALLED,seed,kmeasure,verbose,self_ok_combi,hretval,flag,transpose,newline + SUFFIX nothing + GLOBAL STATS_INSTALLED } PARAMETER { - : BVBASE = 0. : defined in vecst.mod - INSTALLED=0 - kmeasure=0 - verbose=0 - self_ok_combi=0 - hretval=0 - transpose=0 - newline=90 - flag=0 : flag can be used by any of the routines for different things -} - -ASSIGNED { seed } - -VERBATIM -#include "misc.h" -#define MIN_MERGESORT_LIST_SIZE 32 - -union dblint { - int i[2]; - double d; -}; - -unsigned int valseed; -static double *x1x, *y1y, *z1z; -static void vprpr(); - -static int compare_ul(const void* l1, const void* l2) { - int retval; - unsigned long d; - d = (*((const unsigned long*) l1)) - (*((const unsigned long*) l2)); - if(d==0) return 1; - if(d < 0) return -1; - return 0; -} - + : BVBASE = -1. : defined in vecst.mod + STATS_INSTALLED=0 +} + +ASSIGNED { } + +VERBATIM +#include +#include +//#include +extern double BVBASE; +extern double* hoc_pgetarg(); +extern double hoc_call_func(Symbol*, int narg); +extern FILE* hoc_obj_file_arg(int narg); +extern Object** hoc_objgetarg(); +extern void vector_resize(); +extern int vector_instance_px(); +extern void* vector_arg(); +extern double* vector_vec(); +extern double hoc_epsilon; +extern void set_seed(); +extern int ivoc_list_count(Object*); +extern Object* ivoc_list_item(Object*, int); +static void hxe() { hoc_execerror("",0); } + +typedef struct BVEC { + int size; + int bufsize; + short *x; + Object* o; +} bvec; ENDVERBATIM :* v1.slope(num) does a linear regression to find the slope, assuming num=timestep of vector @@ -87,44 +79,14 @@ static double slope(void* vv) { sigy += y[i]; sigx2 += x[i]*x[i]; } - free(x); return (n*sigxy - sigx*sigy)/(n*sigx2 - sigx*sigx); } ENDVERBATIM - -:* v1.moment(v2) stores moments: -VERBATIM -static double moment (void* vv) { - int i, j, n, fl; - double *mdata, *y; - double ave,adev,sdev,svar,skew,curt,s,p; - n = vector_instance_px(vv, &mdata); - fl=0; - if (n<=1) {printf("n must be at least 2 in stats:moment()"); hxe();} - if(ifarg(1)) { - if (hoc_is_object_arg(1)) { - y=vector_newsize(vector_arg(1), 6); fl=1; - } else { - printf("vec.moment(ovec) stores in ovec: ave,adev,sdev,svar,skew,kurt\n"); - return 0; - } - } - for (j=0,s=0;j*maxsqerr) *maxsqerr = val; - *meansqerr += val; - } - *meansqerr=*meansqerr/(double)n; - return *meansqerr; -} -ENDVERBATIM :* v1.stats(num) does a linear regression, assuming num=timestep of vector + VERBATIM static double stats(void* vv) { int i, n; - double *x, *y, *out; + double *x, *y; double timestep, sigxy, sigx, sigy, sigx2, sigy2; - double r, m, b, dmeansqerr,dmaxsqerr; + double r, m, b; /* how to get the instance data */ n = vector_instance_px(vv, &y); @@ -203,571 +139,24 @@ static double stats(void* vv) { m = (n*sigxy - sigx*sigy)/(n*sigx2 - sigx*sigx); b = (sigy*sigx2 - sigx*sigxy)/(n*sigx2 - sigx*sigx); r = (n*sigxy - sigx*sigy)/(sqrt(n*sigx2-sigx*sigx) * sqrt(n*sigy2-sigy*sigy)); - getsqerr(x,y,m,b,n,&dmeansqerr,&dmaxsqerr); //mean,max squared error - if(ifarg(2)){ //save results to output - out=vector_newsize(vector_arg(2),5); - out[0]=m; out[1]=b; out[2]=r; out[3]=dmeansqerr; out[4]=dmaxsqerr; - } else { - printf("Examined %d data points\n", n); - printf("slope = %f\n", m); - printf("intercept = %f\n", b); - printf("R = %f\n", r); - printf("R-squared = %f\n", r*r); - printf("MeanSQErr = %f\n",dmeansqerr); - printf("MaxSQErr = %f\n",dmaxsqerr); - } - free(x); - return 1; -} - -typedef struct pcorst_ { - int pidse[2]; - double* X; - double* Y; - double sigx; - double sigy; - double sigx2; - double sigy2; - double sigxy; -} pcorst; -void* PCorrelTHFunc(void *arg) { - pcorst* p; - int i; - double *X,*Y; - p=(pcorst*)arg; -// X=&p->X[p->pidse[0]]; Y=&p->Y[p->pidse[1]]; - X = p->X; Y = p->Y; - p->sigx=p->sigy=p->sigxy=p->sigx2=p->sigy2=0.0; - for(i=p->pidse[0]; ipidse[1]; i++) { -// p->sigxy += *X * *Y; -// p->sigx += *X; -// p->sigy += *Y; -// p->sigx2 += *X * *X; -// p->sigy2 += *Y * *Y; -// X++; Y++; - p->sigxy += X[i] * Y[i]; - p->sigx += X[i]; - p->sigy += Y[i]; - p->sigx2 += X[i] * X[i]; - p->sigy2 += Y[i] * Y[i]; - } - return NULL; -} - -/* v1.pcorrels2(v2) does a Pearson correlation*/ -#if defined(t) -static double pcorrelsmt(double *x, double* y, int n,int nth) { - int i,nperth,idx,rc; - double sigxy, sigx, sigy, sigx2, sigy2, ret; - pcorst** pp; - pthread_t* pth; - pthread_attr_t attr; - ret=sigxy=sigx=sigy=sigx2=sigy2=0.0; // initialize these - nperth = n / nth; - //allocate thread args - pp = (pcorst**)malloc(sizeof(pcorst*)*nth); - idx=0; - for(i=0;iX = x; - pp[i]->Y = y; - pp[i]->pidse[0] = idx; - pp[i]->pidse[1] = idx + nperth; - idx += nperth; - } - i--; if(pp[i]->pidse[1] < n || - pp[i]->pidse[1] > n) pp[i]->pidse[1] = n; //make sure all values used - //allocate thread IDs - pth=(pthread_t*)malloc(sizeof(pthread_t)*nth); - pthread_attr_init(&attr); - pthread_attr_setdetachstate(&attr, PTHREAD_CREATE_JOINABLE); - //start threads - for(i=0;isigx; - sigy += pp[i]->sigy; - sigxy += pp[i]->sigxy; - sigx2 += pp[i]->sigx2; - sigy2 += pp[i]->sigy2; - } - sigxy -= (sigx * sigy) / n; - sigx2 -= (sigx * sigx) / n; - sigy2 -= (sigy * sigy) / n; - if(sigx2 <= 0) goto PCMTFREE; - if(sigy2 <= 0) goto PCMTFREE; - ret = sigxy / sqrt(sigx2*sigy2); -PCMTFREE: - //free memory - for(i=0;i 1.0 ){ - printf("ppval ERRA: r=%g must be : -1.0 <= r <= 1.0\n",r); - return -1.0; - } - if( n < 3 ){ - printf("ppval ERRB: n too small, can't calc probability on samples with < 3 values!\n"); - return -1.0; - } - df = n-2; // degres of freedom - // Use a small floating point value to prevent divide-by-zero nonsense - // fixme: TINY is probably not the right value and this is probably not - // the way to be robust. The scheme used in spearmanr is probably better. - TINY = 1.0e-20; - ts = r*sqrt(df/((1.0-r+TINY)*(1.0+r+TINY))); - mpval = betai(_threadargscomma_ 0.5*df,0.5,df/(df+ts*ts)); - return mpval; - ENDVERBATIM -} - -VERBATIM - - -static const double* sortdata = NULL; /* used in the quicksort algorithm */ - -/* Helper function for sort. Previously, this was a nested function under - * sort, which is not allowed under ANSI C. - */ -static -int compare(const void* a, const void* b) -{ const int i1 = *(const int*)a; - const int i2 = *(const int*)b; - const double term1 = sortdata[i1]; - const double term2 = sortdata[i2]; - if (term1 < term2) return -1; - if (term1 > term2) return +1; - return 0; -} - -void csort (int n, const double mdata[], int index[]) -/* Sets up an index table given the data, such that mdata[index[]] is in - * increasing order. Sorting is done on the indices; the array mdata - * is unchanged. - */ -{ int i; - sortdata = mdata; - for (i = 0; i < n; i++) index[i] = i; - qsort(index, n, sizeof(int), compare); -} - -static double* getrank (int n, double mdata[]) -/* Calculates the ranks of the elements in the array mdata. Two elements with - * the same value get the same rank, equal to the average of the ranks had the - * elements different values. The ranks are returned as a newly allocated - * array that should be freed by the calling routine. If getrank fails due to - * a memory allocation error, it returns NULL. - */ -{ int i; - double* rank; - int* index; - rank = (double*)malloc(n*sizeof(double)); - if (!rank) return NULL; - index = (int*)malloc(n*sizeof(int)); - if (!index) - { free(rank); - return NULL; - } - /* Call csort to get an index table */ - csort (n, mdata, index); - /* Build a rank table */ - for (i = 0; i < n; i++) rank[index[i]] = i; - /* Fix for equal ranks */ - i = 0; - while (i < n) - { int m; - double value = mdata[index[i]]; - int j = i + 1; - while (j < n && mdata[index[j]] == value) j++; - m = j - i; /* number of equal ranks found */ - value = rank[index[i]] + (m-1)/2.; - for (j = i; j < i + m; j++) rank[index[j]] = value; - i += m; - } - free (index); - return rank; -} - -/* -The spearman routine calculates the Spearman rank correlation between two vectors. -n (input) int The number of elements in a data vector -data1 (input) double array -- the first vector -data2 (input) double array -- the second vector -*/ -static double spearman(int n, double* data1, double* data2) -{ int i; - int m = 0; - double* rank1; - double* rank2; - double result = 0.; - double denom1 = 0.; - double denom2 = 0.; - double avgrank; - double* tdata1; - double* tdata2; - tdata1 = (double*)malloc(n*sizeof(double)); - if(!tdata1) return 0.0; /* Memory allocation error */ - tdata2 = (double*)malloc(n*sizeof(double)); - if(!tdata2) /* Memory allocation error */ - { free(tdata1); - return 0.0; - } - for (i = 0; i < n; i++) - { tdata1[m] = data1[i]; - tdata2[m] = data2[i]; - m++; - } - if (m==0) return 0; - rank1 = getrank(m, tdata1); - free(tdata1); - if(!rank1) return 0.0; /* Memory allocation error */ - rank2 = getrank(m, tdata2); - free(tdata2); - if(!rank2) /* Memory allocation error */ - { free(rank1); - return 0.0; - } - avgrank = 0.5*(m-1); /* Average rank */ - for (i = 0; i < m; i++) - { const double value1 = rank1[i]; - const double value2 = rank2[i]; - result += value1 * value2; - denom1 += value1 * value1; - denom2 += value2 * value2; - } - /* Note: denom1 and denom2 cannot be calculated directly from the number - * of elements. If two elements have the same rank, the squared sum of - * their ranks will change. - */ - free(rank1); - free(rank2); - result /= m; - denom1 /= m; - denom2 /= m; - result -= avgrank * avgrank; - denom1 -= avgrank * avgrank; - denom2 -= avgrank * avgrank; - if (denom1 <= 0) return 0; /* include '<' to deal with roundoff errors */ - if (denom2 <= 0) return 0; /* include '<' to deal with roundoff errors */ - result = result / sqrt(denom1*denom2); - return result; -} - -static double scorrel(void* vv) { - int i, n; - double *x, *y; - n = vector_instance_px(vv, &x); - if ((i=vector_arg_px(1, &y)) != n ) {printf("scorrelERRA: %d %d\n",n,i); hxe();} - return spearman(n,x,y); -} - -//* Kendall's correlation routines -//** Erfcc() Returns the complementary error function erfc(x) with fractional error -// everywhere less than 1.2 x 10^-7. -// from place.mod, is that from numerical recipes?? -double Erfcc (double x) { - double mt,z,ans; - z=fabs(x); - mt=1.0/(1.0+0.5*z); - ans=mt*exp(-z*z-1.26551223+mt*(1.00002368+mt*(0.37409196+mt*(0.09678418+\ - mt*(-0.18628806+mt*(0.27886807+mt*(-1.13520398+mt*(1.48851587+\ - mt*(-0.82215223+mt*0.17087277))))))))); - return x >= 0.0 ? ans : 2.0-ans; -} - -//** Rktau() R version of kendall tau, doesnt have huge memory footprint -//function returns kendall's tau -double Rktau (double* x, double* y, int n){ - int i,j; double c,vx,vy,sx,sy,var,z,tau; - c = vx = vy = 0.0; - for(i = 0; i < n; i++) { - for(j = 0; j < i; j++) { - sx = (x[i] - x[j]); - sx = ((sx > 0) ? 1 : ((sx == 0)? 0 : -1)); - sy = (y[i] - y[j]); - sy = ((sy > 0) ? 1 : ((sy == 0)? 0 : -1)); - vx += sx * sx; - vy += sy * sy; - c += sx * sy; - } - } - if(vx>0 && vy>0) { - tau = c / sqrt(vx*vy); - return tau; - } - return 0.; -} - -//** vec1.kcorrel(vec2,[fast version -- useful for large arrays, vec of size 1 holding p-value]) -// kendall's tau correlation -static double kcorrel (void* vv) { - int i, n; - double *x, *y, *prob, *i1d, *i2d, *ps, var, z, tau; - n = vector_instance_px(vv, &x); - if ((i=vector_arg_px(1, &y)) != n ) {printf("kcorrel ERRA: %d %d\n",n,i); hxe();} - if(ifarg(2) && *getarg(2)) { - i1d=dcrset(n*3); i2d=&i1d[n]; ps=&i2d[n]; tau=kcorfast(x,y,i1d,i2d,n,ps); - } else { - tau = Rktau(x,y,n); - } - if(!(ifarg(3) && vector_arg_px(3,&prob))) prob = 0x0; //does user want to store p-value? - if(prob) { //get p-value - var = (4.0 * n + 10.0) / (9.0 * n * (n - 1.0)); - z = tau / sqrt(var); - *prob = Erfcc(fabs(z)/1.4142136); //when prob small, chance of tau having its value by chance is small - } - return tau; -} - -//** mycompare() comparison function for qsort -- sorts in ascending order -static int mycompare (const void* a, const void* b) -{ double d1,d2; - d1 = *(double*)a; - d2 = *(double*)b; - if (d1 < d2) return -1; - if (d1 > d2) return +1; - return 0; -} - -//** mergesort_array() -// recursive mergesort -- sorts a by splitting into sublists, sorting, and recombining -void mergesort_array (double a[], int size, double temp[],unsigned long* swapcount) { - int i1, i2, i, tempi, j, vv; - double *right,*left; - if(size<=1) return; //base case -- 1 element is sorted by definition - if (0 && size = 0; j--) { - if (a[j] <= vv) break; - a[j + 1] = a[j]; - } - a[j + 1] = vv; - } - return; - } - mergesort_array(a, size/2, temp,swapcount); //sort left half - mergesort_array(a + size/2, size - size/2, temp,swapcount); //sort right half - //merge halves together - i=tempi=0; i1 = size/2; i2 = size - size/2; - left = a; right = &a[size/2]; - while(i1>0 && i2>0) { - if(*right < *left) { - *swapcount += i1; - temp[i] = *right++; - i2--; - } else { - temp[i] = *left++; - i1--; - } - i++; - } - if(i2>0) { - while(i2-->=0 && i=0 && i2) printf("nPair=%lu\n",nPair); - if(verbose>3){printf("i1d after qsort2: "); for(i=0;i 0) { - qsort(&i2d[i-tieCount-1],tieCount+1,sizeof(double),mycompare); - m1 += tieCount * (tieCount + 1) / 2; - s += getMs(&i2d[i-tieCount-1],tieCount+1); - tieCount = 0; - } - } - if(verbose>2) printf("tieCount=%lu\n",tieCount); - if(tieCount > 0) { - qsort(&i2d[n-tieCount-1],tieCount+1,sizeof(double),mycompare); - m1 += tieCount * (tieCount + 1) / 2; - s += getMs(&i2d[n-tieCount-1],tieCount+1); - } - if(verbose>2) printf("tieCount=%lu\n",tieCount); - swapCount = 0; - - mergesort_array(i2d,n,ps,&swapCount); //sort input2 & count # of swaps to get into sorted order - if(verbose>3) { printf("i2d after mergesort: "); for(i=0;i2) printf("swapCount=%lu\n",swapCount); - - m2 = getMs(i2d,n); if(verbose>2) printf("s=%lu m1=%lu m2=%lu\n",s,m1,m2); - s -= (m1 + m2) + 2 * swapCount; - denom1=nPair-m1; denom2=nPair-m2; if(verbose>2) printf("s=%lu d1=%g d2=%g\n",s,denom1,denom2); - if(denom1>0. && denom2>0.) return s / sqrt(denom1*denom2); else return 0.; -} -//root mean square of vector's elements -static double rms (void* vv) { - int i,n; - double *x,sum; - if(!(n=vector_instance_px(vv, &x))) {printf("rms ERRA: 0 sized vector!\n"); hxe();} - sum=0.0; - for(i=0;i0.) return sqrt(sum); else return 0.0; -} - -//cumulative sum of vector's elements -static double cumsum (void* vv) { - int i,n; - double *x,*y; - if(!(n=vector_instance_px(vv, &x))) {printf("cumsum ERRA: 0 sized vector!\n"); hxe();} - if(vector_arg_px(1, &y) != n) {printf("cumsum ERRB: output vec size needs size of %d\n",n); hxe();} - memcpy(y,x,sizeof(double)*n); - for(i=1;i=nx) { hoc_execerror("# of flips exceeds (or ==) vector size", 0); } for (i=0; i < nx; i++) { x[i] = BVBASE; } for (i=0,jj=0; i < flip; i++) { /* flip these bits */ - ii = (int) ((nx+1)*drand48()); + ii = (int) ((nx+1)*((double)rand()/RAND_MAX)); // replaced drand48() with (double)rand()/RAND_MAX if (x[ii]==BVBASE) { x[ii] = 1.; if (flag) { y[jj] = ii; jj++; } @@ -836,682 +218,52 @@ static double randwd(void* vv) { } ENDVERBATIM -:* v1.smash(veclist,base) -: smash squeezes a set of numbers into a single double by considering them as digits -: in base base -- x[i]+=vvo[i][j]*wt; where wt is base^i -: note that handles transpose -- ie can smash on (transpose==1) or across each vec in a veclist -VERBATIM -static double smash (void* vv) { - int i, j, nx, nv[VRRY], num; - Object* ob; - double *x, *vvo[VRRY], wt, wtj; - nx = vector_instance_px(vv, &x); - ob=*hoc_objgetarg(1); - if (ifarg(2)) wtj=*getarg(2); else wtj=10.; - num = ivoc_list_count(ob); - if (num>VRRY) {printf("vecst:smash ERRA: can only handle %d vecs: %d\n",VRRY,num); hxe();} - if (transpose) if (nx!=num) { printf("vecst:smash ERRB %d %d %d\n",i,nx,nv[i]);hxe(); } - for (i=0;iVRRY) {printf("stats:dpro ERR: can only handle %d vecs: %d\n",VRRY,num); hxe();} - for (i=0;i= randvar and mul by step -: v1.setrnd(4.5,min,max[,seed]) // [min,max) -: v1.setrnd(5[,n,seed]) -- integers [0,100) or [0,n] -: v1.setrnd(5[,min,max,seed]) -- integers [min,max] -- if seed=0 it's not reset -: v1.setrnd(5,ind[,seed]) -- random values from ind -: v1.setrnd(6) -- unique integers as follows: -: v1.setrnd(6,min,max[,seed]) -- unique in [min,max] -: v1.setrnd(6,min,max,exclude_vec[,seed]) -- unique in [min,max] excluding values in exclude_vec -VERBATIM -static double setrnd (void* vv) { - int flag, i,j,k,n,cnt; unsigned int nx, nx1, nex, lt, rt, mid; - double *x, y, *ex, *ex2, min, max, dfl, tmp, step, num; - unsigned long value; - value=1; - nx = vector_instance_px(vv, &x); - flag = (int)(dfl=*getarg(1)); - if (flag==1) { - for (i=0; i < nx; i++) x[i] = (double)rand()/RAND_MAX; - } else if (flag==2) { - for (i=0; i < nx; i++) x[i] = drand48(); - } else if (flag==3) { // scop_random()'s cheap and dirty rand - unsigned long a = 2147437301, c = 453816981, m = ~0; - value = (unsigned long) seed; - for (i=0; i < nx; i++) { - value = a * value + c; - x[i] = (fabs((double) value / (double) m)); - } - seed=(double)value; - } else if (flag==4) { // mcell_ran4() doubles - ex=0x0; i=2; - if (ifarg(i)) { - if (hoc_is_object_arg(i)) { - nex=vector_arg_px(i++,&ex); // vector to look in - step=ifarg(i)?*getarg(i):1.0; - max=1.0; i++; - } else { - if (dfl==4.5 || ifarg(4)) { // flag 4.5 resolves ambiguity of arg3 max or seed - min=*getarg(i++); - max=*getarg(i++)-min; - dfl=4.5; - } else { - max=*getarg(i++); - } - } - } else max=1.0; // default - if (ifarg(i)) { y=*getarg(i++); if (y) valseed=(unsigned int)y; } // look for seed - if (max==0) { for (i=0;iex[mid]) { - if (num1) { // sort the exclude vector - scrset(nex); - x1x = (double *)realloc(x1x,sizeof(double)*nx*4); - for (i=0;imax || ex[j]0 && ex[j]<=ex[j-1])) { - printf("%g in exclusion list -- out of range [%g,%g] or a repeat\n",ex[j],min,max);hxe();} - if (max-min+1-nex==nx) { - for (i=min,k=0;i<=max;i++) { - y=(double)i; - for (j=0;j= to an excluded value must shift by # its >= to - for (i=0;i=ex[j]) k++; // will move it up by k - x[i]=z1z[i]+k; - } - } else for (i=0;i1;) { - mcell_ran4(&valseed, y, 1, n); - n--; - k=(int)y[0]; // random int(n) // 0 <= k < n. - temp = x[n]; - x[n] = x[k]; - x[k] = temp; - } -} - -//shuffle array of unsigned ints -void uishuffle(unsigned int* x,int nx) { - int n,k; unsigned int temp; double y[1]; - for (n=nx;n>1;) { - mcell_ran4(&valseed, y, 1, n); - n--; - k=(int)y[0]; // random int(n) // 0 <= k < n. - temp = x[n]; - x[n] = x[k]; - x[k] = temp; - } -} - -//shuffle array of unsigned ints -void ishuffle(int* x,int nx) { - int n,k,temp; double y[1]; - for (n=nx;n>1;) { - mcell_ran4(&valseed, y, 1, n); - n--; - k=(int)y[0]; // random int(n) // 0 <= k < n. - temp = x[n]; - x[n] = x[k]; - x[k] = temp; - } -} - -unsigned long choose (int n, int k) { - int i,delta; - unsigned long ret; - //assert ((n >= 0) && (k >= 0)); - if (n < k) return 0; - if (n==k) return 1; - if (k < n - k) { - delta = n - k; - } else { - delta = k; - k = n - k; - } - ret = delta + 1; - for (i = 2; i <= k; ++i) - ret = (ret * (delta + i)) / i; - return ret; -} - -// computes combinadic given combination vector -unsigned long syncci (int nn,int kk,int* ccvv) { - unsigned long c = 0; - while ((kk > 0) && (*ccvv >= kk)) { - c += choose (*ccvv++, kk--); - } - return c; -} - - -// computes combinadic given combination vector -unsigned long syncc (int nn,int kk,double* ccvv) { - unsigned long c = 0; - while ((kk > 0) && (*ccvv >= kk)) { - c += choose (*ccvv++, kk--); - } - return c; -} - -// computes combination vector given combinadic -void synccv (int nn,int kk,int cc,double* ccvv) { - unsigned long n_k; - while (--nn >= 0) { - n_k = choose (nn, kk); - if (cc >= n_k) { - cc -= n_k; - *ccvv++ = nn; - --kk; - } - } -} - -ENDVERBATIM - -: returns the # of combinations for choosing k from a set of n -: combs(n,k) -FUNCTION combs () { - VERBATIM - unsigned int n,k; - n=(unsigned int)*getarg(1); - k=(unsigned int)*getarg(2); - return choose(n,k); - ENDVERBATIM -} - -:* vec.comb(nn,kk,cc) -: returns combination indexed by cc, from selecting kk elements from set of nn elements -: element ids are from 0..nn-1 -VERBATIM -static double comb (void* vv) { - double* x; - int nn,kk,cc,sz,i; - sz=vector_instance_px(vv,&x); - nn=(int)*getarg(1); - kk=(int)*getarg(2); - cc=(int)*getarg(3); - if(sz= %d , is %d!\n",kk,sz); - return 0.0; - } - memset(x,0,sizeof(double)*kk); - synccv(nn,kk,cc,x); - vector_resize((IvocVect*)vv,kk); - return 1.0; -} -ENDVERBATIM - -:* vec.combid(nn,kk) -: returns combination index corresponding to data in vector -: indices from selecting kk elements from set of nn elements -: element ids are from 0..nn-1 -VERBATIM -static double combid (void* vv) { - double* x; - int nn,kk,sz,i; - sz=vector_instance_px(vv,&x); - nn=(int)*getarg(1); - kk=(int)*getarg(2); - if(sz= %d , is %d!\n",kk,sz); - return 0.0; - } - return syncc(nn,kk,x); -} -ENDVERBATIM - -VERBATIM -int findlong (unsigned long* p,unsigned long val,int istart,int iend) { - int i; - for(i=istart;i<=iend;i++) if(p[i]==val) return 1; - return 0; -} - -ENDVERBATIM - -:* rsampsig(vIN0,vIN1,prc,"vhfmeasure","compfunc",vhsout,[onesided,nocmbchk,inorder]) -: vIN0 - elements of group 0 -: vIN1 - elements of group 1 -: prc - fraction of subsets to measure with "vhfmeasure" , or total # of combinations to test iff prc > 1.0 -: "vhfmeasure" - hoc function that takes 1 vector as arg and returns a double, must set hretval_stats -: to the return value so can be read from here -: "compfunc" - hoc function that takes 2 doubles as args & returns 1 iff a boolean condition is satisfied, must set -: hretval_stats to return value so can be read from here -: vhsout - vector used in hocmeasure , must have size >= vIN0.size + vIN1.size -: onesided - optional, use one sided p-value or two-sided (default = 1) , onsided==0 means use two-sided -: nocmbchk - skip combination ID check, useful for groups with large size ( > 30000 ) , to avoid overflow, def=1 == skip combination ID check -: returns -1.0 on error -VERBATIM -static double rsampsig(void* vv){ - int n0,n1,na,nn,kk,cc,i,j,*pm,szthis,onesided,nocmbchk,bti,*pids; - unsigned long nruncombs,nallcombs,*pcombids; - double *g0,*g1,*ga,prc,*g0t,*g1t,dmobs,dm0,dm1,*phso,nmatch,*pthis,dret; - IvocVect* vhso; //vector * for changing size at end - Symbol* pHocVecFunc,*pHocCompFunc; //hoc function pointers - dret=-1.0; - g0t=g1t=NULL; pm=pids=NULL; pcombids=NULL;//init arrays to null - szthis=vector_instance_px(vv, &pthis); //size of calling vector - n0=vector_arg_px(1,&g0);//group 0 size - n1=vector_arg_px(2,&g1);//group 1 size - na=n0+n1;//total # of elements - prc=*getarg(3);//% of combinations to try, or total # of combinations to try iff > 1.0 - if(!(pHocVecFunc=hoc_lookup(gargstr(4)))){ - printf("rsampsig ERRA: couldn't find hoc vec func %s\n",gargstr(4)); - goto CLEANRSAMPSIG; - } - if(!(pHocCompFunc=hoc_lookup(gargstr(5)))){ - printf("rsampsig ERRB: couldn't find hoc comp func %s\n",gargstr(5)); - goto CLEANRSAMPSIG; - } - if(vector_arg_px(6,&phso)= %d!\n",na); - goto CLEANRSAMPSIG; - } - if(prc<=0.0) { - printf("rsampsig ERRD: invalid value for arg 3, must be > 0.0!\n"); - goto CLEANRSAMPSIG; - } - vhso=vector_arg(6);//get vector arg for resize - ga=(double*)malloc(sizeof(double)*na);//all elements - memcpy(ga,g0,sizeof(double)*n0);//copy elements of group 0 to ga - memcpy(ga+n0,g1,sizeof(double)*n1);//append elements of group 1 to ga - //form nruncombs combinations, select elements into groups, and compare measure on groups - g0t=(double*)malloc(sizeof(double)*n0);//temp storage for group 0 - g1t=(double*)malloc(sizeof(double)*n1);//temp storage for group 1 - if(n01) printf("choose(%d,%d)=%ld\n",na,kk,choose(na,kk)); - nallcombs=choose(na,kk);//all possible combinations selecting kk elements from a set of size na - nruncombs=prc>1.0?prc:prc*nallcombs; nmatch=0.0;//# of combinations to run - if(szthis1) printf("na=%d , kk=%d, n0=%d, n1=%d\n",na,kk,n0,n1); - if(verbose>1) printf("nruncombs = %ld\n",nruncombs); - if(verbose>2) pm=(int*)malloc(sizeof(int)*na); - for(i=0;i0 && findlong(pcombids,pcombids[i],0,i-1)); //make sure pcombids[i] wasnt used yet - if(verbose) if(i%100==0) { printf("."); fflush(stdout); } - for(j=0;j2){ for(j=0;j observed test statistic, to get p-value - for(i=0;i dmobs) nmatch++; - dret=nmatch/(double)nruncombs;//this output value doesnt mean anything yet -CLEANRSAMPSIG: //free memory and return - if(ga) free(ga); if(g0t) free(g0t); if(g1t) free(g1t); if(pm) free(pm); - if(pcombids) free(pcombids); if(pids) free(pids); - return dret; -} -ENDVERBATIM - - -VERBATIM -// rantran(index_vec1,value_vec1[,index_vec2,value_vec2,...]) -// index_vec can be replaced with a step which gives regular pointers or base values -// in presence of a step value can be eg 10 for 0..9 [or 10.3 for 0,3,6,9 -- not implemented] -// rantran() translates a random index vec through a series of mappings from a value onto -// one or more possible values: eg rand col indices to rand minicols within that col to -// rand cell within the minicol to rand locations on the cell dend -static double rantran (void* vv) { - int i,j,ix,ixe,ixvn,nvn,rvn,na,xj; - double *ixv, *nv, *x, y[1], ixn,step,indx; - rvn=vector_instance_px(vv, &x); - for (na=1;ifarg(na);na++) {} na--; // count args - for (i=1;i1.) { - x1x = (double *)realloc(x1x,sizeof(double)*rvn); - mcell_ran4(&valseed, x1x, rvn, indx); - } - } - for (j=0;j-1) { // step in by xi[type] and then augment by random val - x[j]=step + x[j]*indx + ((indx>1.)?(floor(x1x[j])):0.0); - } else { - xj=(int)x[j]; // prior level index - ix=(int)ixv[xj]; ixn=ixv[xj+1]-ix; // possible next level indices; pick 1 in [ix,ixe] - if (ixn==1.) { - x[j]=nv[ix]; - } else { - mcell_ran4(&valseed, y, 1, ixn); // pick 1 rand value [0,ix) - x[j]=nv[(int)y[0]+ix]; +static double hamming(void* vv) { + int i, nx, ny, nz, prflag; + double* x, *y, *z,sum; + sum = 0.; + nx = vector_instance_px(vv, &x); + ny = vector_arg_px(1, &y); + if (ifarg(2)) { /* write a diff vector to z */ + prflag = 1; nz = vector_arg_px(2, &z); + } else { prflag = 0; } + if (nx!=ny || (prflag && nx!=nz)) { + hoc_execerror("Vectors must be same size", 0); + } + for (i=0; i < nx; ++i) { + if (x[i] != y[i]) { sum++; + if (prflag) { z[i] = 1.; } + } else if (prflag) { z[i] = 0.; } } - } - } - } - return (double)rvn; // number of substitutions performed -} -ENDVERBATIM - -:* vec.shuffle() Fisher-Yates shuffle (from wikipedia) -: pointlessly extended to augment the vector n-fold (used for cell projections in columns) -: makes more sense to just have a parallel random vector of intra-columnar indices -VERBATIM -static double shuffle (void* vv) { - int i,j,k,n,nx,augfac; double *x, y[1], temp, augstep; - nx=vector_instance_px(vv, &x); - if (ifarg(1)) { - augfac=(int)*getarg(1); - if (ifarg(2)) augstep=*getarg(2); else augstep=1.0/augfac; - x=vector_newsize((IvocVect*)vv,nx*augfac); - for (i=1;inx) {printf("stats vpr ERRA OOB: %d %d\n",min,max); hxe();} if (ifarg(1)) { - if (hoc_is_double_arg(1)) { - flag=(int) *getarg(1); - if (flag==2) { // binary which is default - for (i=min; iBVBASE)?1:0); - } else if (flag==0) { - for (i=min; i=10?"|":"",(int)x[i],x[i]>=10?"|":""); - } else if (flag==-1) { - for (i=min; iBVBASE) { fprintf(f,"%d",1); - } else { fprintf(f,"%d",0); } - } - fprintf(f,"\n"); + f = hoc_obj_file_arg(1); + for (i=0; iBVBASE) { fprintf(f,"%d",1); + } else { fprintf(f,"%d",0); } } + fprintf(f,"\n"); } else { - for (i=min; iBVBASE)?1:0); - printf("\n"); - } - return 1.; -} -ENDVERBATIM - -:* v1.vpr2(v2,[BASE]) prints out 2 vectors using = where same -- optional arg allows base 10,16,64 -: use VERBOSE to adjust printing verbose=1 for no spaces; verbose=2 no spaces and abbrev 00->_ -VERBATIM -static double vpr2 (void* vv) { - int i, flag, n, nx, ny, colc, base, min,max, fl2; - double *x, *y, cnt, ign; char c; - nx = vector_instance_px(vv, &x); - cnt=min=0; max=nx; - ny = vector_arg_px(1, &y); - if (nx!=ny) {printf("vpr2 diff sizes %d %d\n",nx,ny); hxe();} - base=(int)*getarg(2); - flag=(ifarg(3)?(int)*getarg(3):2); - for (i=min,fl2=0,colc=0; i0 && x[i]==BVBASE && y[i]==BVBASE) { - if (!fl2) {printf(" _ "); colc+=3;} - fl2=1; - } else { - fl2=0; colc++; - vprpr(x[i],base); - if (colc>(int)newline){printf("\n "); colc=0;} - } - } - printf("\n"); - for (i=min,fl2=0,colc=0; i0 && x[i]==BVBASE && y[i]==BVBASE) { - if (!fl2) {printf(" _ "); colc+=3;} - fl2=1; - } else { fl2=0; - vprpr(y[i],base); - colc++; - if (colc>(int)newline){printf("\n ",colc); colc=0;} - } - } - printf("\n"); - if (flag==1) { - for (i=min,n=0,fl2=0,colc=0; i(int)newline){printf("\n ",colc); colc=0;} + for (i=0; iBVBASE) { printf("%d",1); + } else { printf("%d",0); } } printf("\n"); } - return 0; -} - -static void vprpr (double x, int base) { - int xx; - xx=(int)x; - if (base==0) { printf("%5.2f",x); - } else if (xx>=base && base!=0) { printf("+"); - } else if (base==64) { // base 64 - if (xx<16) { printf("%x",xx); // 0-f 0-15 - } else if (xx<36) { printf("%c",xx+87); // g-z 16-35 - } else if (xx<62) { printf("%c",xx+29); // A-Z 36-61 - } else if (xx<63) { printf("@"); // @ 62 - } else if (xx<64) { printf("="); // = 63 - } else printf("ERROR"); - } else if (base==10) { printf("%d",xx); - } else if (base==16) { printf("%x",xx); - } + return 1.; } ENDVERBATIM -:* v1.bin(targ,invl[min,max]) place counts for each interval -:* v1.bin(list,invl[min,max]) using NQS("count","index") for easy sorting -: like .hist() but doesn't throw away values max -: note that optional max denotes the start of an interval not end +:* v1.bin(targ,invl) place counts for each interval VERBATIM -static double bin (void* vv) { - int i, j, nx, maxsz, lfl; - double* x, *y, *ix, invl, min, max, maxf, jj; - Object* ob; - IvocVect* voi[2]; +static double bin(void* vv) { + int i, j, nx, ny, next, maxsz; + double* x, *y, invl, loc; - min=0; max=1e9; maxf=-1e9; nx = vector_instance_px(vv, &x); - ob = *hoc_objgetarg(1); - if (strncmp(hoc_object_name(ob),"Vector",6)==0) { lfl=0; // lfl is list flag - if ((maxsz=openvec(1, &y))==0) hxe(); - voi[0]=vector_arg(1); - } else { // list of 2 - lfl=1; - maxsz = list_vector_px3(ob, 0, &y, &voi[0]); - if (maxsz!=(i=list_vector_px3(ob,1,&ix,&voi[1]))){printf("binERRA %d %d\n",maxsz,i); hxe();} - } - invl = *getarg(2); - if (ifarg(4)) {min=*getarg(3); max=*getarg(4); - } else if (ifarg(3)) max=*getarg(3); + ny = vector_arg_px(1, &y); + invl = (int)*getarg(2); + + vv=vector_arg(1); + maxsz=vector_buffer_size(vv); + vector_resize(vv, maxsz); + if (x[nx-1]/invl>(double)(maxsz-1)) { + printf("Need size %d in target vector (%d)\n",(int)(x[nx-1]/invl+1),maxsz); + hoc_execerror("",0); } for (j=0; j=max) { y[(j=(int)(jj=(max-min)/invl))]++; - } else if (x[i]<=min) { y[(j=0)]++; jj=0.; - } else { - if (x[i]>maxf) maxf=x[i]; // max found - j=(int)(jj=(x[i]-min)/invl); - if (j>maxsz-1) {printf("bin ERRB OOB: %d>%d\n",j,maxsz-1); hxe();} + for (i=0,j=0,loc=invl; iloc) { loc+=invl; j++; } y[j]++; } - if (lfl) ix[j]=jj+min; - } - maxsz=(max==1e9)?(int)(maxf/invl+1):(int)((max-min)/invl+1); - vector_resize(voi[0], maxsz); - if (lfl) vector_resize(voi[1], maxsz); - return (double)maxsz; -} -ENDVERBATIM - -:* ind.ihist(vec,list,tmin,tmax,binsz,[,ioff]) does binning for individual indices -: into list -- ioff gives an offset such that index N maps to list.o(0) -: normally used for a raster pair of cell indices (ind) and spike time (tvec) -: eg nq=new NQS(#cells) -: ind.ihist(tvec,nq.vl,tmin,tmax,100) -VERBATIM -static double ihist (void* vv) { - unsigned int i, j, k, n, nx, c; - double *x, *tv, ioff, min, max, nbin, binsz; - ListVec* pL; Object* obl; - nx = vector_instance_px(vv, &x); // vector of indices - i = vector_arg_px(1, &tv); // vector of times - if (i!=nx) {printf("vecst:ihist()ERR0: diff size %d %d\n",nx,i); hxe();} - if (!flag && !ismono1(tv,nx,1)){ - printf("vecst:ihist()ERR0A: set flag_stats for non-monotonic time vec\n",nx,i); hxe();} - pL = AllocListVec(obl=*hoc_objgetarg(2)); - min=*getarg(3); max=*getarg(4); binsz=*getarg(5); - if (binsz<=0) {printf("stats:ihist()ERR0B: binsz must be >0 (%g)\n",binsz); hxe();} - nbin=floor((max-min)/binsz); - max=min+binsz*nbin; // new max - if (verbose) printf("%g-%g in %g bins of %g\n",min,max,nbin,binsz); - if (ifarg(6)) ioff=*getarg(6); else ioff=0.; - if (pL->isz<2) {printf("stats:ihist()ERRA: %d\n",pL->isz); FreeListVec(&pL); hxe();} - c=pL->isz; // number of columns (list length) - // pL->pv[i][j] -- i goes through the list and j goes through each vector - for (i=0;ipv[i]=list_vector_resize(obl, i, (int)nbin); - for (j=0;j<(int)nbin;j++) pL->pv[i][j]=0.; - } - i=n=0; - if (!flag) for (;i=min) break; // tvec sorted so can move forward to begin - for (;i=max || tv[i]=max) break; - k=(int)(x[i]-ioff); // which cell index - if (k>=c || k<0) continue; // ignore outside of index range - j=(int)((tv[i]-min)/binsz); // which time bin - // if(j>=nbin||j<0){printf("INTERR %d %d %d %g %g %g %g\n",k,c,j,nbin,tv[i],min,binsz);hxe();} - pL->pv[k][j]++; - n++; - } - flag=0; - FreeListVec(&pL); - return (double)n; -} -ENDVERBATIM - -:* vec.irate(vhist,binsz) - sets contents of vec to instantaneous rate using inverse -: of interspike intervals -VERBATIM -static double irate (void* vv) { - unsigned int i, j, n, nx; - double *prate,*phist,binsz,t1,t2; - nx = vector_arg_px(1, &phist); - vector_resize((IvocVect*)vv,nx); - vector_instance_px(vv, &prate); - binsz = *getarg(2); - for(i=0;i0) { - t1=t2=i; - prate[i]=phist[i]*1e3/binsz; - i++; - break; - } - } - for(;i0) { - t2=i; - prate[i]=phist[i]*1e3/binsz; - break; - } - } - if(verbose>1) if(t1==t2) printf("t1==t2!\n"); - for(i=t1;i0) { - prate[i] = phist[i]*1e3/binsz; - t1=t2; t2=i; - if(verbose>2) printf("t1 %g t2 %g\n",t1,t2); - } else { - if(verbose>1) if(t1==t2) printf("t1==t2!\n"); - prate[i] = 1e3 / (binsz*(t2-t1)); - } - } - return 1.0; -} -ENDVERBATIM - -VERBATIM -//cholesky decomposition from numerical recipies chapter 2.9 -//Given a positive-dfinite symmetric matrix a[1..n][1..n], this routine constructs its Cholesky -//decomposition, A = L * LT . On input, only the upper triangle of a need be given; it is not -//modified. The Cholesky factor L is returned in the lower triangle of a, except for its diagonal -//elements which are returned in p[1..n] -int choldc(double **a, int n, double p[]) -{ - int i,j,k; - double sum; - for (i=0;i=0;k--) sum -= a[i][k]*a[j][k]; - if (i == j) { - if (sum <= 0.0) return 0; - p[i]=sqrt(sum); - } else a[j][i]=sum/p[i]; - } } - return 1; + vector_resize(vv, j); + return (double)j; } ENDVERBATIM -: mktscor(pearson r-value,list of time-series) -: each time-series should be normally distributed -: upon return, the values will be changed and each pair of time-series will be correlated with pcorrel = ~r -: same method as here http://comisef.wikidot.com/tutorial:correlation -FUNCTION mktscor () { -VERBATIM - double dret,*ptmp,**cF,**Y,**X,**cFT,r; - ListVec *pTS; - int i,j,k,nrows,ncols,nrcf; - dret=0.0; - pTS=0x0; ptmp=0x0; cF=cFT=0x0; X=Y=0x0; - r = *getarg(1); - pTS = AllocListVec(*hoc_objgetarg(2)); - nrows=pTS->plen[0]; ncols=pTS->isz; - X = getdouble2D(nrows,ncols); - for(i=0;ipv[i][j]; - nrcf=ncols; - cF = getdouble2D(nrcf,nrcf); - for(i=0;ii) cF[i][j]=0.0; else if(j==i) cF[i][j]=ptmp[i]; - cFT = getdouble2D(nrcf,nrcf); - for(i=0;ipv[j][i]=Y[i][j]; -MKTSCORFREE: - dret=1.0; - if(pTS) FreeListVec(&pTS); - if(ptmp) free(ptmp); - if(cF) freedouble2D(&cF,nrcf); - if(cFT) freedouble2D(&cFT,nrcf); - if(Y) freedouble2D(&Y,nrows); - if(X) freedouble2D(&X,nrows); - return dret; -ENDVERBATIM -} - :* PROCEDURE install_stats() -PROCEDURE install () { - if (INSTALLED==1) { - printf("$Id: stats.mod,v 1.215 2010/08/18 19:20:23 samn Exp $") - } else { - INSTALLED=1 +PROCEDURE install_stats () { + STATS_INSTALLED=1 VERBATIM - x1x=y1y=z1z=0x0; - valseed=1; install_vector_method("slope", slope); - install_vector_method("moment", moment); install_vector_method("vslope", vslope); install_vector_method("stats", stats); - install_vector_method("pcorrel", pcorrel); - install_vector_method("scorrel",scorrel); - install_vector_method("kcorrel",kcorrel); - install_vector_method("rms",rms); install_vector_method("vstats", vstats); install_vector_method("randwd", randwd); install_vector_method("hamming", hamming); install_vector_method("flipbits", flipbits); install_vector_method("flipbalbits", flipbalbits); install_vector_method("vpr", vpr); - install_vector_method("vpr2", vpr2); install_vector_method("bin", bin); - install_vector_method("ihist", ihist); - install_vector_method("setrnd", setrnd); - install_vector_method("rantran", rantran); - install_vector_method("distance", distance); - install_vector_method("ndprd", ndprd); - install_vector_method("smash", smash); - install_vector_method("smash1", smash1); - install_vector_method("dpro", dpro); - install_vector_method("unnan", unnan); - install_vector_method("combi", combi); - install_vector_method("shuffle", shuffle); - install_vector_method("comb", comb); - install_vector_method("combid", combid); - install_vector_method("rsampsig",rsampsig); - install_vector_method("irate",irate); - install_vector_method("cumsum",cumsum); -ENDVERBATIM - } -} - -PROCEDURE prhash (x) { -VERBATIM { - long unsigned int xx; - xx=*(long unsigned int*)(&_lx); - printf("%16lx\n",xx); -} ENDVERBATIM } @@ -1972,55 +430,34 @@ VERBATIM { ENDVERBATIM } -FUNCTION getseed () { - VERBATIM - seed=(double)valseed; - return seed; - ENDVERBATIM -} - : unable to get the drand here to recognize the same fseed used in rand FUNCTION vseed () { VERBATIM #ifdef WIN32 + double seed; if (ifarg(1)) seed=*getarg(1); else { printf("TIME ACCESS NOT PRESENT IN WINDOWS\n"); hxe(); } - srand48((unsigned)seed); + srand((unsigned)seed); // changed from srand48 to srand - mjb set_seed(seed); return seed; #else struct timeval tp; struct timezone tzp; + double seed; if (ifarg(1)) seed=*getarg(1); else { gettimeofday(&tp,&tzp); seed=tp.tv_usec; } - srand48((unsigned)seed); + srand((unsigned)seed); // srand48 set_seed(seed); srandom(seed); - valseed=(unsigned int)seed; return seed; #endif ENDVERBATIM } -: mc4seed() set valseed for mccell_ran4() as series of factors so as not to lose low order -: digits before casting double to unsigned int, then call mcell_ran4_init with new valseed -FUNCTION mc4seed () { - VERBATIM - int i; - valseed=(unsigned int)(*getarg(1)); - for (i=2;ifarg(i);i++) { - valseed*=(unsigned int)(*getarg(i)); - } - mcell_ran4_init(valseed); // do initialization - return valseed; - ENDVERBATIM -} - - : from Numerical Recipes in C FUNCTION gammln (xx) { VERBATIM { @@ -2043,15 +480,14 @@ FUNCTION gammln (xx) { FUNCTION betai(a,b,x) { VERBATIM { double bt; - if (_lx < 0.0 || _lx > 1.0) {printf("Bad x in routine BETAI\n"); hxe();} if (_lx == 0.0 || _lx == 1.0) bt=0.0; else - bt=exp(gammln(_threadargscomma_ _la+_lb)-gammln(_threadargscomma_ _la)-gammln(_threadargscomma_ _lb)+_la*log(_lx)+_lb*log(1.0-_lx)); + bt=exp(gammln(_la+_lb)-gammln(_la)-gammln(_lb)+_la*log(_lx)+_lb*log(1.0-_lx)); if (_lx < (_la+1.0)/(_la+_lb+2.0)) - return bt*betacf(_threadargscomma_ _la,_lb,_lx)/_la; + return bt*betacf(_la,_lb,_lx)/_la; else - return 1.0-bt*betacf(_threadargscomma_ _lb,_la,1.0-_lx)/_lb; + return 1.0-bt*betacf(_lb,_la,1.0-_lx)/_lb; } ENDVERBATIM } @@ -2093,91 +529,3 @@ VERBATIM { } ENDVERBATIM } - -FUNCTION symval() { -VERBATIM { - Symbol *sym; - sym = hoc_get_symbol(* hoc_pgargstr(1)); - // should do type check eg sym->type == VAR - return *(hoc_objectdata[sym->u.oboff]._pval); // is this ._pval safe?? - } -ENDVERBATIM -} - -:* tval(r,N) , computes t-statistic , r == pearson correlation coefficient, N == size of sample -FUNCTION tstat() { - VERBATIM - double r = fabs(*getarg(1)); - double N = *getarg(2); - if(N < 2) { printf("tstat ERRA: N must be > 2!\n"); return -1; } - return r * sqrt(N-2.)/sqrt(1.0-(r*r)); - ENDVERBATIM -} - -FUNCTION tdistrib() { - VERBATIM - double x = *getarg(1); - double dof = *getarg(2); - double res = (gammln(_threadargscomma_ (dof+1.0) / 2.0 ) / gammln(_threadargscomma_ dof / 2.0 ) ); - double pi = 3.14159265358979323846; - res *= (1.0 / sqrt( dof * pi ) ); - res *= pow((1 + x*x/dof),-1.0*((dof+1.0)/2.0)); - return res; - ENDVERBATIM -} - -: based on code from http://www.psychstat.missouristate.edu/introbook/rdist.htm -: return estimate of critical correlation coefficient (r) that >= to it would be considered significant -: at a given sample size of $1 (N) . probability of null hypothesis p < 0.05 -: iff $2 == 1, probability of null hypothesis p < 0.01 ( r will then be higher ) -: null hypothesis == correlation between variables , r , == 0 -: this function uses a lookup table so is limited in estimates -- max sample size , N, == 1000 , but -: if your sample is bigger and your r value is also bigger, the correlation is significant because -: in general, as the sample size , N , increases, the r value above which is considered significant decreases -FUNCTION rcrit () { - VERBATIM - double rtbl[68][3]={//for each row: index0==N, index1==p0.05 value of r, index2==p0.01 value of r - 1 , 0.997 , 0.999 , 2 , 0.950 , 0.990 , 3 , 0.878 , 0.959 , 4 , 0.811 , 0.917 , 5 , 0.754 , 0.874 , - 6 , 0.707 , 0.834 , 7 , 0.666 , 0.798 , 8 , 0.632 , 0.765 , 9 , 0.602 , 0.735 , 10 , 0.576 , 0.708 , - 11 , 0.553 , 0.684 , 12 , 0.532 , 0.661 , 13 , 0.514 , 0.641 , 15 , 0.482 , 0.606 , 16 , 0.468 , 0.590 , - 17 , 0.456 , 0.575 , 18 , 0.444 , 0.561 , 19 , 0.433 , 0.549 , 20 , 0.423 , 0.537 , 21 , 0.413 , 0.526 , - 22 , 0.404 , 0.515 , 23 , 0.396 , 0.505 , 24 , 0.388 , 0.496 , 25 , 0.331 , 0.487 , 26 , 0.374 , 0.478 , - 27 , 0.367 , 0.470 , 28 , 0.361 , 0.463 , 29 , 0.355 , 0.456 , 30 , 0.349 , 0.449 , 31 , 0.344 , 0.442 , - 32 , 0.339 , 0.436 , 33 , 0.334 , 0.430 , 34 , 0.329 , 0.424 , 35 , 0.325 , 0.418 , 36 , 0.320 , 0.413 , - 37 , 0.316 , 0.408 , 38 , 0.312 , 0.403 , 39 , 0.308 , 0.398 , 40 , 0.304 , 0.393 , 41 , 0.301 , 0.389 , - 42 , 0.297 , 0.384 , 43 , 0.294 , 0.380 , 44 , 0.291 , 0.376 , 45 , 0.288 , 0.372 , 46 , 0.284 , 0.368 , - 47 , 0.281 , 0.364 , 48 , 0.279 , 0.361 , 58 , 0.254 , 0.330 , 63 , 0.244 , 0.317 , 68 , 0.235 , 0.306 , - 73 , 0.227 , 0.296 , 78 , 0.220 , 0.286 , 83 , 0.213 , 0.278 , 88 , 0.207 , 0.270 , 93 , 0.202 , 0.263 , - 98 , 0.195 , 0.256 , 123 , 0.170 , 0.230 , 148 , 0.159 , 0.210 , 173 , 0.148 , 0.194 , - 198 , 0.138 , 0.181 , 298 , 0.113 , 0.148 , 398 , 0.098 , 0.128 , 498 , 0.088 , 0.115 , - 598 , 0.080 , 0.105 , 698 , 0.074 , 0.097 , 798 , 0.070 , 0.091 , 898 , 0.065 , 0.086 , - 998 , 0.062 , 0.081 - }; - double N; - int get99, i , tablen, df; - N = *getarg(1); - get99 = ifarg(2) ? (int)*getarg(2) : 0; - tablen = 68; - - if(N < 3){ - printf("rcrit ERRA: N must be >= 3!\n"); - return -1.0; - } - - if( N > 998+2 ) { printf("rcrit WARNA: Using N=1000 as estimate.\n"); N = 998+2; } - - df = (int)N - 2; - - for(i=0;i df) { - if(get99) - return rtbl[i-1][2] + ((rtbl[i][2] - rtbl[i-1][2])*((df - rtbl[i-1][0] )/(rtbl[i][0] - rtbl[i-1][0]))); - else - return rtbl[i-1][1] + ((rtbl[i][1] - rtbl[i-1][1])*((df - rtbl[i-1][0] )/(rtbl[i][0] - rtbl[i-1][0]))); - } - } - return -1.0; - ENDVERBATIM -} diff --git a/vecst.mod b/vecst.mod index 7ee49ec..5bc7839 100644 --- a/vecst.mod +++ b/vecst.mod @@ -1,4 +1,4 @@ -: $Id: vecst.mod,v 1.494 2010/12/03 17:26:01 billl Exp $ +: $Id: vecst.mod,v 1.326 2006/06/24 23:38:36 billl Exp $ :* COMMENT COMMENT @@ -7,22 +7,19 @@ triplet return location of a triplet in a vector onoff turns state vec on or off depending on values in other vecs bpeval used by backprop: vo=outp*(1-outp)*del w like where but sets chosen nums // modified from .wh in 1.23 -whi find indices where vec equal to given values dest.xing(src,tvec,thresh) determines where vector crosses in a positive direction +src.updown(thresh,destlist) determines crossings in both directions+peaks dest.snap(src,tvec,dt) interpolate src with tvec to prior dt step, saves only highest value xzero count 0 crossings by putting in a 1 or -1 negwrap wrap negative values around to pos (with flag just set to 0) indset(ind,val) sets spots indexed by ind to val ismono([arg]) checks if a vector is monotonically increaseing (-1 decreasing) count(num) count the number of nums in vec -muladd(mul,add) mul*x+add -binfind(num) find index for num in a sorted vector scr.fewind(ind,veclist) // uses ind as index into other vecs ind.findx(vecAlist,vecBlist) // uses ind as index into vecA's to load in vecB's ind.sindx(vecAlist,vecBlist) // replace ind elements in vecA's with vecB's values ind.sindv(vecAlist,valvec) // replace ind elements in vecA's with vals scr.nind(ind,vec1,vec2[,vec3,vec4]) // uses ind to elim elems in other vecs -yv.circ(xv,x0,y0,rad) // put coords of circle into xv and yv ind.keyind(key,vec1,vec2[,vec3,vec4]) // pick out bzw. values from vectors ind.slct(key,args,veclist) // pick out bzw. values from vectors ind.slor(key,args,veclist) // do OR rather than AND function of slct @@ -53,25 +50,17 @@ rnd round off to nearest integer ENDCOMMENT NEURON { -THREADSAFE SUFFIX nothing : BVBASE is bit vector base number (typically 0 or -1) - GLOBAL BVBASE, RES, VECST_INSTALLED, DEBUG_VECST, VERBOSE_VECST, INTERP_VECST, LOOSE + GLOBAL BVBASE, RES, VECST_INSTALLED, SHM_UPDOWN, DEBUG_VECST } PARAMETER { - BVBASE = 0. + BVBASE = -1. VECST_INSTALLED=0 DEBUG_VECST=0 - VERBOSE_VECST=1 - INTERP_VECST=1 : interpolation - LOOSE=1e-6 - : code words - ERR= -1.3480e121 - GET= -1.3479e121 - SET= -1.3478e121 - OK= -1.3477e121 - NOP= -1.3476e121 + : misc + ERR=-1.3479e121 : 0 args ALL=-1.3479e120 @@ -87,236 +76,145 @@ PARAMETER { EQU=-1.3470e120 EQV=-1.3469e120 : value equal to same row value in parallel vector EQW=-1.3468e120 : value found in other vector - EQX=-1.3467e120 : value found in sorted vector - EQY=-1.34665e120 : value found in sorted vector -- within LOOSE - NEQ=-1.3466e120 - SEQ=-1.3465e120 - RXP=-1.3464e120 - : 2 args - IBE=-1.3463e120 - EBI=-1.3462e120 - IBI=-1.3461e120 - EBE=-1.3460e120 + NEQ=-1.3467e120 + SEQ=-1.3466e120 + RXP=-1.3465e120 + : 2 args + IBE=-1.3464e120 + EBI=-1.3463e120 + IBI=-1.3462e120 + EBE=-1.3461e120 + + SHM_UPDOWN=4 : used in updown() for measuring sharpness } ASSIGNED { RES } VERBATIM -#include "misc.h" -ENDVERBATIM - -VERBATIM -// Maintain parallel int vector to avoid slowness of repeated casts -int cmpdfn (double a, double b) {return ((a)<=(b))?(((a) == (b))?0:-1):1;} -static unsigned int bufsz=0; -unsigned int scrsz=0; -unsigned int *scr=0x0; -unsigned int dcrsz=0; -double *dcr=0x0; - -int *iscr=0x0; -unsigned int iscrsz=0; - -int *iscrset (int nx) { - if (nx>iscrsz) { - iscrsz=nx+10000; - if (iscrsz>0) { iscr=(int *)realloc((void*)iscr,(size_t)iscrsz*sizeof(int)); - } else { iscr=(int *)ecalloc(iscrsz, sizeof(int)); } - } - return iscr; -} - -unsigned int *scrset (int nx) { - if (nx>scrsz) { - scrsz=nx+10000; - if (scrsz>0) { scr=(unsigned int *)realloc((void*)scr,(size_t)scrsz*sizeof(int)); - } else { scr=(unsigned int *)ecalloc(scrsz, sizeof(int)); } - } - return scr; -} +#ifndef NRN_VERSION_GTEQ_8_2_0 +#include +#include +#include +extern double* hoc_pgetarg(); +extern double hoc_call_func(Symbol*, int narg); +extern FILE* hoc_obj_file_arg(int narg); +extern Object** hoc_objgetarg(); +extern void vector_resize(); +extern int vector_instance_px(); +extern void* vector_arg(); +extern double* vector_vec(); +extern double hoc_epsilon; +extern double chkarg(); +extern void set_seed(); +extern int ivoc_list_count(Object*); +extern Object* ivoc_list_item(Object*, int); +extern int hoc_is_double_arg(int narg); +extern char* hoc_object_name(Object*); +char ** hoc_pgargstr(); +#endif +static int ismono1 (double *x, int n, int flag); +int list_vector_px(Object *ob, int i, double** px); +int list_vector_px2 (Object *ob, int i, double** px, void** vv); +int list_vector_resize (Object *ob, int i, int sz); +static double sc[6]; +static void hxe() { hoc_execerror("",0); } + +typedef struct BVEC { + int size; + int bufsize; + short *x; + Object* o; +} bvec; + +#define VRRY 50 + +#define BYTEHEADER int _II__; char *_IN__; char _OUT__[16]; int BYTESWAP_FLAG=0; +#define BYTESWAP(_X__,_TYPE__) \ + if (BYTESWAP_FLAG == 1) { \ + _IN__ = (char *) &(_X__); \ + for (_II__=0;_II__dcrsz) { - dcrsz=nx+10000; - if (dcrsz>0) { dcr=(double*) realloc((void*)dcr,(size_t)dcrsz*sizeof(double)); - } else { dcr=(double*)ecalloc(dcrsz, sizeof(double)); } - } - return dcr; -} +#define UNCODE(_X_,_J_,_Y_) {(_Y_)=floor((_X_)/sc[(_J_)])/sc[4]; \ + (_Y_)=floor(sc[4]*((_Y_)-floor(_Y_))+0.5);} ENDVERBATIM :* v1.ident() gives addresses and sizes VERBATIM -static double ident (void* vv) { +static double ident(void* vv) { int nx,bsz; double* x; nx = vector_instance_px(vv, &x); - bsz=vector_buffer_size((IvocVect*)vv); + bsz=vector_buffer_size(vv); printf("Obj*%x Dbl*%x Size: %d Bufsize: %d\n",vv,x,nx,bsz); - return (double)nx; + return 0.0; } ENDVERBATIM :* v1.indset(ind,x[,y]) sets indexed values to x and other values to optional y -: v1.indset(ind,vvec[,y]) sets indexed values to vvec values -: v1.indset(ind,"INC",1) increments values -: v1.indset(ind,"INC",-1) decrements values -: v1.indset(v2, "EQU",val,x[,y]) checks if v2 is EQU to val VERBATIM -static double indset (void* vv) { - int i, nx, ny, nz, flag, equ; char *op; - double *x, *y, *z, val, val2, inc; - nx = vector_instance_px(vv, &x); - ny = vector_arg_px(1, &y); - val2=flag=equ=inc=0; - if (hoc_is_object_arg(2)) { - flag=1; - nz = vector_arg_px(2, &z); - if (ny!=nz) z=vector_newsize((IvocVect*)vector_arg(2),ny); - } else if (hoc_is_double_arg(2)) { - val=*getarg(2); - } else if (hoc_is_str_arg(2)) { - op = gargstr(2); - if (strcmp(op,"EQU")==0) equ=1; else if (strcmp(op,"INC")==0) inc=1; else { - printf("indset %s not recog\n",op); hxe(); } - } - if (equ) { - val2 = *getarg(3); val=*getarg(4); - for (i=0; i nx) {printf("vecst:inset ERRB Index exceeds vector size %g %d\n",y[i],nx); hxe();} - if (inc!=0) x[(int)y[i]]+=inc; else if (flag) x[(int)y[i]]=z[i]; else x[(int)y[i]]=val; - } - } - return (double)i; -} -ENDVERBATIM - -:* v1.mkind(ind) creates a simple index for starts of a sorted integer vector v1 -VERBATIM -static double mkind(void* vv) { - int i, j, nx, ny, flag; - double *x, *y, flox, last, min, max; - nx = vector_instance_px(vv, &x); - ny = vector_arg_px(1, &y); - flag=ifarg(2)?(int)*getarg(2):0; - min=x[0]; max=x[nx-1]; - ny=max-min+4+(flag?0:min); - if (ny>1e6) printf("vecst:mkind() WARNING: index of size %d being built\n",ny); - y=vector_newsize((IvocVect*)vector_arg(1),ny); - y[0]=0.; y[ny-3]=nx; y[ny-2]=min; y[ny-1]=max; // last 2 record min,max - if (min==max) {printf("vecst:mkind() ERRA: min==max %g %g\n",min,max); hxe();} - if (!flag) for (j=1;j<=min;j++) y[j]=0.; else j=1; - for (i=1,last=floor(min); i=last+1;j++,last++) y[j]=(double)i; - last=flox; - } - return min; -} -ENDVERBATIM - -:* yv.circ(xv,x0,y0,rad) sets cartesian x,y coords for circle with center x0,y0; radius rad -VERBATIM -static double circ (void* vv) { - int i, nx, ny, flag, lnew; - double *x, *y, x0, y0, x1, y1, rad, theta; - lnew=0; - nx = vector_instance_px(vv, &x); - ny = vector_arg_px(1, &y); - if (ny!=nx) { hoc_execerror("v.circ: Vector sizes don't match.", 0); } - x0=*getarg(2); y0=*getarg(3); rad=*getarg(4); - if (ifarg(6)) { // gives center and a point on the circle instead of radius - x1=rad; y1=*getarg(5); - rad=sqrt((x0-x1)*(x0-x1) + (y0-y1)*(y0-y1)); - lnew=*getarg(6); // resize the vectors - } else if (ifarg(5)) lnew=*getarg(5); // resize the vectors - if (lnew) { nx=lnew; x=vector_newsize((IvocVect*)vv,nx); y=vector_newsize((IvocVect*)vector_arg(1),ny=nx); } - for (i=0,theta=0; i nx) { hoc_execerror("v.indset: Index exceeds vector size", 0); } + if (flag) x[(int)y[i]]=z[i]; else x[(int)y[i]]=val; } - if (nz>nx) {printf("v.roton: Can't rotate %d vals on vec of size %d\n",nz,nx); hxe();} - for (i=nz,j=0; iVRRY) hoc_execerror("ERR: fewind can only handle VRRY vectors", 0); - if (flag) ix=(char*)ecalloc(nx,sizeof(char)); - if (!flag && nxscrsz) { + if (ni>scrsz) { if (scrsz>0) { free(scr); scr=(unsigned int *)NULL; } - scrsz=nx+10000; + scrsz=ni+10000; scr=(unsigned int *)ecalloc(scrsz, sizeof(int)); } - if (flag) { - for (i=0;i=nx || j<0) { printf("fewind ERR1A %d %d\n",j,nx); free(ix);hxe();} - scr[k]=j; - ix[j]=1; // indicate that has already been used - k++; - } - } - ni=nx; // index up to max - } else for (i=0;i=nx || scr[i]<0) { printf("fewind ERR1 %d %d\n",scr[i],nx); hoc_execerror("Index vector out-of-bounds", 0); } } - if (flag) for (i=0;iscrsz) { @@ -366,55 +263,6 @@ static double findx (void* vv) { } ENDVERBATIM -:* ind.lma(vecAlist,vecBlist,beg,end[,mul,add]) -: copy each vector in Alist to Blist indices from beg to end; with sequential mul/add -: ind picks out specific vectors if not doing all of them -VERBATIM -static double lma (void* vv) { - IvocVect* vv_tmp = (IvocVect*)vv; - int i, j, k, ia, ib, ni, nj, nx, av[VRRY], bv[VRRY], num, numb, beg, end, *xx; - Object *ob1, *ob2; - double *ind, *avo[VRRY], *bvo[VRRY], mul,mmul,add,madd; - ni = vector_instance_px(vv, &ind); - xx=iscrset(ni); - ob1 = *hoc_objgetarg(1); ob2 = *hoc_objgetarg(2); - beg = (int)*getarg(3); end= (int)*getarg(4); - mul=ifarg(5)?*getarg(5):1; add=ifarg(6)?*getarg(6):0; - num = ivoc_list_count(ob1); - numb = ivoc_list_count(ob2); - if ((ni==0 && numb!=num) || (ni>0 && ni!=num)) { - printf("lma ERRA: wrong# of outvecs: %d (inlist:%d,inds:%d)\n",numb,num,ni); hxe();} - for (i=0,j=0;inum){printf("lma ERRA1 OOB: %d %d\n",j,num); hxe();} - } - nj=j; - if (num>VRRY){printf("lma ****ERRB****: can only handle %d vectors\n",VRRY); hxe();} - for (i=0;i=end || beg<0 || end>nx) {printf("lma ERRC1 OOB %d - %d (%d)\n",beg,end,nx); hxe();} - for (i=0;i0) { - for (ia=0,ib=0,mmul=1.,madd=0.;iaVRRY) hoc_execerror("ERR: vecst::slct can only handle VRRY vectors", 0); - for (i=0,j=0;i=EQV && key[j]<=EQY) { // EQV-EQY take extra vec arg of any size - i++; - nv[i] = list_vector_px(lob, i, &vvo[i]);// pick up extra vector of any size - } else if (ni!=nv[i]) { + if (ni!=nv[i] && (key[j]!=EQW || k!=1)) { printf("vecst::slct ERR %d %d %d %d %d\n",i,j,k,ni,nv[i]); hoc_execerror("index and searched vectors must all be same size: ", 0); } + if (key[j]==EQW || key[j]==EQV) if (k==0){j--;k++;} else k=0; // EQW,EQV take 2 vector args } for (j=0;j=EQV && key[i]<=EQY) n++; // special cases take 2 vec args + for (i=0,n=0;ivvo[i+1][p]) lt=p+1; else if (valvvo[i+1][p]+LOOSE) lt=p+1; else if (val=arg[m+1])) { fl=0; break; } else continue; // IBE="[)" include-bracket-exclude @@ -581,12 +409,9 @@ static double slct (void* vv) { fl=0; break; } else continue; // "()" : exclude-bracket-exclude } else {printf("vecst::slct ERR4 %g\n",key[n]); hoc_execerror("Unknown key",0);} } - if (fl) { - ind[k++]=j; // all equal - if (flag==1) break; - } + if (fl) ind[k++]=j; // all equal } - vector_resize((IvocVect*)vv, k); + vector_resize(vv, k); return (double)k; } ENDVERBATIM @@ -594,8 +419,8 @@ ENDVERBATIM :* ind.slor(key,args,vec1,vec2[,vec3,vec4,...]) : picks out indices of numbers in key from multiple vectors VERBATIM -static double slor (void* vv) { - int i, j, k, m, n, p, ni, nk, na, nv[VRRY], num, fl, field[VRRY], lt, rt; +static double slor(void* vv) { + int i, j, k, m, n, p, ni, nk, na, nv[VRRY], num, fl, field[VRRY]; Object* lob; double *ind, *key, *arg, *vvo[VRRY], val; @@ -606,15 +431,13 @@ static double slor (void* vv) { num = ivoc_list_count(lob); if (num>VRRY) hoc_execerror("ERR: vecst::slor can only handle VRRY vectors", 0); - for (i=0,j=0;i=EQV && key[j]<=EQY) { // EQV-EQY take extra vec arg of any size - i++; - nv[i] = list_vector_px(lob, i, &vvo[i]);// pick up extra vector of any size - } else if (ni!=nv[i]) { - printf("vecst::slor ERR %d %d %d %d %d\n",i,j,k,ni,nv[i]); + if (ni!=nv[i] && (key[j]!=EQW || k!=1)) { + printf("vecst::slct ERR %d %d %d %d %d\n",i,j,k,ni,nv[i]); hoc_execerror("index and searched vectors must all be same size: ", 0); } + if (key[j]==EQW || key[j]==EQV) if (k==0){j--;k++;} else k=0; // k counts 2 vecs } for (j=0;j=EQV && key[i]<=EQY) n++; // special case takes 2 vec args + for (i=0,n=0;ivvo[i+1][p]) lt=p+1; else if (val=arg[m+1])) { continue; } else {fl=1; break;} // IBE="[)" include-bracket-exclude @@ -674,78 +487,53 @@ static double slor (void* vv) { } if (fl) ind[k++]=j; // all equal } - vector_resize((IvocVect*)vv, k); + vector_resize(vv, k); return (double)k; } ENDVERBATIM -:* src.whi(valvec,indvec[,&i0,&i1,...]) returns index of first one for each val +:* v.iwr() VERBATIM -static double whi (void* vv) { - int i, j, nx, na, nb, cnt; double *x, *val, *ind; - nx = vector_instance_px(vv, &x); - i = vector_arg_px(1, &val); - na = vector_arg_px(2, &ind); - if (i!=na) {printf("vecst:whi() takes 2 eq length vecs:%d %d\n",i,na); hxe();} - scrset(na); - for (i=0;i2) { - if (nb-2!=na) {printf("vecst:whi() wrong #args:%d %d\n",na,nb-2); hxe();} - for (i=3,j=0;jscrsz) { + if (scrsz>0) { free(scr); scr=(unsigned int *)NULL; } + scrsz=nx+10000; + scr=(unsigned int *)ecalloc(scrsz, sizeof(int)); + } for (i=0;iscrsz) { if (scrsz>0) { free(scr); scr=(unsigned int *)NULL; } scrsz=n+10000; scr=(unsigned int *)ecalloc(scrsz, sizeof(int)); } if (n!=nx) { - nx=vector_buffer_size((IvocVect*)vv); + nx=vector_buffer_size(vv); if (n<=nx) { - vector_resize((IvocVect*)vv, n); nx=n; + vector_resize(vv, n); nx=n; } else { printf("%d > %d :: ",n,nx); hoc_execerror("Vector max capacity too small for ird ", 0); } } - r=fread(scr,sizeof(int),n,f); + fread(scr,sizeof(int),n,f); for (i=0;imaxsz) { printf("%d > %d :: ",n,maxsz); hoc_execerror("Vector max capacity too small for fread2 ", 0); } else { - vector_resize((IvocVect*)vv, n); + vector_resize(vv, n); } if (type==6 || type==16) { // unsigned ints unsigned int *xs; @@ -778,7 +566,7 @@ static double fread2 (void* vv) { scr=(unsigned int *)ecalloc(scrsz, sizeof(int)); } xs=(unsigned int*)scr; - r=fread(xs,sizeof(int),n,fp); + fread(xs,sizeof(int),n,fp); if (type==16) BYTESWAP_FLAG=1; for (i=0;i=nx) x=vector_newsize((IvocVect*)vv,nx*=4); - x[k]=*getarg(i); - } else { - ny=vector_arg_px(i, &y); - if (k+ny>=nx) x=vector_newsize((IvocVect*)vv,nx=2*(nx+ny)); - for (j=0;j=ny) y=vector_newsize(vc,ny*=4); - y[j++]=(double)i; - } else { - *(hoc_pgetarg(2)) = (double)i; - return 1.0; - } - } else return 1; - } - if (j>0) vector_resize(vc,j); // only fall through here when doing indvwhere - return (double)j; + return 0.0; } ENDVERBATIM @@ -857,60 +594,19 @@ static double insct (void* vv) { int i, j, k, nx, nv1, nv2, maxsz; double *x, *v1, *v2; nx = vector_instance_px(vv, &x); - if (maxsz==0) maxsz=1000; - maxsz=vector_buffer_size((IvocVect*)vv); - x=vector_newsize((IvocVect*)vv, maxsz); + maxsz=vector_buffer_size(vv); + vector_resize(vv, maxsz); nv1 = vector_arg_px(1, &v1); nv2 = vector_arg_px(2, &v2); for (i=0,k=0;imaxsz) { + printf("\tinsct WARNING: ran out of room: %d<%d\n",maxsz,k); + } else { vector_resize(vv, k); } return (double)k; } ENDVERBATIM -:* ind.linsct(vlist,this_many) -: return indices that show up this_many times or more in the vectors in vlist -VERBATIM -static double linsct (void* vv) { - int i,j,k,nx,maxsz,min,cnt,lt,rt,p,iv,jv,jj,fl; ListVec* pL; - double *x, val; - nx = vector_instance_px(vv, &x); - maxsz=vector_buffer_size((IvocVect*)vv); - if (maxsz==0) maxsz=1000; - x=vector_newsize((IvocVect*)vv, maxsz); - pL = AllocListVec(*hoc_objgetarg(1)); - min= (int)*getarg(2); - for (iv=0;ivisz;iv++) { - if (!ismono1(pL->pv[iv],pL->plen[iv],2)){printf("linsct() ERRA not monotonic %d\n",iv); hxe();} - if (pL->pv[iv][0]<0){printf("linsct() ERRB neg value in %d is %g\n",iv,pL->pv[iv][0]);hxe();} - } - for (iv=0,k=0;iv<=pL->isz-min+1;iv++) for (j=0;jplen[iv];j++) { - val=pL->pv[iv][j]; - for (i=0,fl=0;iisz;jv++) { - fl=lt=0; rt=pL->plen[jv]-1; - while (lt <= rt) { - p = (lt+rt)/2; - if (val>pL->pv[jv][p]) lt=p+1; else if (valpv[jv][p]) rt=p-1; else { - fl=1; break;} - } - if (fl) cnt++; - if (cnt>=min) { - if (k==maxsz) x=vector_newsize((IvocVect*)vv, maxsz*=2); - x[k++]=val; - break; - } - } - } - vector_resize((IvocVect*)vv, k); - return (double)k; -} -ENDVERBATIM - :* vdest.vfill(vsrc) : fill vdest with multiple instances of vsrc until reach size VERBATIM @@ -920,7 +616,7 @@ static double vfill (void* vv) { nx = vector_instance_px(vv, &x); nv1 = vector_arg_px(1, &v1); for (i=0;imaxsz) { + printf("\tcull WARNING: ran out of room: %d<%d\n",maxsz,k); + } else { vector_resize(vv, k); } return (double)k; } ENDVERBATIM @@ -952,7 +650,7 @@ VERBATIM static double redundout (void* vv) { int i, j, nx, nv1, maxsz, ncntr, indflag, cntflag; double *x, *v1, *cntr, val; - IvocVect *vc; + void *vc; if (ifarg(2)) indflag=(int)*getarg(2); else indflag=0; if (ifarg(3)) { cntflag=1; ncntr = vector_arg_px(3, &cntr); vc=vector_arg(3); @@ -960,7 +658,7 @@ static double redundout (void* vv) { for (i=0;i=maxsz) { printf("\tredundout WARNING: ran out of room: %d=ncntr) { printf("\tredundout WARNING: cntr ran out of room: %dVRRY) hoc_execerror("mredundout ****ERRA****: can only handle VRRY vectors", 0); for (i=0;imaxsz){printf("mredundout****ERRE**** vec overflow %d>%d\n",k,maxsz);hxe();} for (i=0;i0 && k0.) hi=mid-1; else if (tmp<0) lo=mid+1; else break; - } - dist=fabs(x[mid]-targ); minind=mid; - for (i=-1;i<=1;i+=2) { kk=mid+i; // check the flanking values - if (kk>0 && kk=nv[1]) { printf("nearall WARN: oor %d %d\n",k,nv[1]); return -1.; } - vvo[0][k]=dist; vvo[2][k]=targ; vvo[3][k]=x[minind]; - k++; - } - } - if (k>scrsz) { - if (scrsz>0) { free(scr); scr=(unsigned int *)NULL; } - scrsz=k+10000; - scr=(unsigned int *)ecalloc(scrsz, sizeof(int)); - } - for (kk=0,i=0;iy[i]+epsilon) return 0; - return 1; -} -ENDVERBATIM - -:* v1.samp(VEC,rate) does something like a resample -VERBATIM -//linear interpolation -static double samp (void* vv) { - int i, nx, cnt, iOrigSz, maxsz, iNewSz, isrcidx,isrcidx1; - double *x, *y, dNewSz, scale, dsrcidx, frac; - nx = vector_instance_px(vv, &x); // dest - iOrigSz = vector_arg_px(1, &y); // source - dNewSz = *getarg(2); // new size - maxsz=vector_buffer_size((IvocVect*)vv); - iNewSz = (int)dNewSz; - if (iNewSz>maxsz) {printf("VECST samp ERRA: dest vec too small: %d %d\n",iNewSz,maxsz); hxe();} - vector_resize((IvocVect*)vv,iNewSz); - - scale = (double) iOrigSz / (double) iNewSz; - for(i=0;i")) { - if (f==1) {for (i=0; ith) { x[(int)ind[i]]=val;c++;}} - } else {for (i=0; ith) { x[i]=val;c++;}}} + if (f==1) {for (i=0; ith) { x[(int)ind[i]] = val;}} + } else {for (i=0; ith) { x[i] = val;}}} } else if (!strcmp(op,"<")) { - if (f==1) {for (i=0; i=")) { - if (f==1) {for (i=0; i=th) { x[(int)ind[i]]=val;c++;}} - } else {for (i=0; i=th) { x[i]=val;c++;}}} + if (f==1) {for (i=0; i=th) { x[(int)ind[i]] = val;}} + } else {for (i=0; i=th) { x[i] = val;}}} } else if (!strcmp(op,"<=")) { - if (f==1) {for (i=0; ith) { // ? passing thresh @@ -1637,36 +1213,164 @@ static double xing (void* vv) { hoc_execerror("Dest vec too small in xing ", 0); } if (i>0) { // don't record if first value is above thresh - if (INTERP_VECST) { - if (tvf) { - dest[j++] = tvec[i-1] + (tvec[i]-tvec[i-1])*(th-src[i-1])/(src[i]-src[i-1]); - } else { - dest[j++] = (i-1) + (th-src[i-1])/(src[i]-src[i-1]); - } + if (tvf) { + dest[j++] = tvec[i-1] + (tvec[i]-tvec[i-1])*(th-src[i-1])/(src[i]-src[i-1]); } else { - if (tvf) dest[j++]=tvec[i]; else dest[j++]=i; + dest[j++] = (i-1) + (th-src[i-1])/(src[i]-src[i-1]); } } f=1; } } else { // below thresh if (f==1) { f=0; } // just passed going down - if (d2f) { - if (INTERP_VECST) { - if (tvf) { - dest2[j++] = tvec[i-1] + (tvec[i]-tvec[i-1])*(th-src[i-1])/(src[i]-src[i-1]); - } else { - dest2[j++] = (i-1) + (th-src[i-1])/(src[i]-src[i-1]); + } + } + vector_resize(vv, j); + return (double)i; +} +ENDVERBATIM + +:* src.updown(thresh,dlist) -- returns indices, change into time by .mul(tstep) +: dest.updown(src) -- default thresh=0; returns indices +: a .where that looks for threshold crossings and then doesn't pick another till +: comes back down again; places values from tvec in dest; interpolates +VERBATIM +#define UDSL 50 +#define UDNQ 10 +static double updown (void* vv) { + int i, k, m, n, nsrc, jj[UDSL], f[UDSL], lc, dsz[UDSL], nqsz[UDNQ], thsz, lc2, done, dbn; + double *src, *tvec, *th, *dest[UDSL], *nq[UDNQ], *tmp, *dbx, lt, thdist; + Object *ob, *ob2; + void *vvd[UDSL], *vvth, *vnq[UDNQ]; + nsrc = vector_instance_px(vv, &src); + thsz = vector_arg_px(1, &th); + ob = *hoc_objgetarg(2); + ob2 = *hoc_objgetarg(3); + tmp = (double *)ecalloc(nsrc, sizeof(double)); + lc = ivoc_list_count(ob); + lc2 = ivoc_list_count(ob2); + if (lc>UDSL) {printf("updown ERRF mismatch: max slice list:%d %d\n",UDSL,lc); hxe();} + if (lc2>UDNQ) {printf("updown ERRG mismatch: max nq sz:%d %d\n",UDNQ,lc2); hxe();} + // nq=new NQS("LOC","PEAK","WIDTH","BASE","HEIGHT","START","SLICES","SHARP","INDEX","FILE") + if (lc2!=10) {printf("updown ERRB mismatch: nq sz should be 10:%d\n",lc2); hxe();} + if (nsrcth[k];i++) {} // make sure start below thresh + for (; ith[k]) { + if (f[k]==0) { // ? passing thresh + if (jj[k]>=dsz[k]){printf("(%d,%d,%d) :: ",k,jj[k],dsz[k]); + hoc_execerror("Dest vec too small in updown ", 0); } + dest[k][jj[k]++] = (i-1) + (th[k]-src[i-1])/(src[i]-src[i-1]); + f[k]=1; + tmp[k]=-1e9; dest[k][jj[k]]=-1.; + } + if (f[k]==1 && src[i]>tmp[k]) { // use tmp[] temporarily + tmp[k]=src[i]; // pick out max + dest[k][jj[k]] = (double)i; // location of this peak + } + } else { // below thresh + if (f[k]==1) { // just passed going down + jj[k]++; + dest[k][jj[k]++] = (i-1) + (src[i-1]-th[k])/(src[i-1]-src[i]); + f[k]=0; + } + } + } + } + // truncate vectors to multiples of 3: + for (k=0;k=nqsz[0]) { printf("updown OOR: %d %d\n",k,nqsz[0]); hxe(); } + nq[0][k]=(double)i; // approx location of the peak of the spike + nq[2][k]=tmp[i+1]; // location of right side + nq[5][k]=tmp[i-1]; // start of spike + nq[6][k]=tmp[i]; // neg of # of slices + k++; + } + if (DEBUG_VECST && ifarg(4)) { dbn = vector_arg_px(4, &dbx); // DEBUG + if (dbn0 && nq[5][i]=0 && !done;m--) { // go from bottom to top + for (n=1;nnq[0][i-1]) { + // nq[6][i]=nq[5][i]; // overwrite #of slices with end of overlap + nq[5][i]=dest[m][n-1]; nq[2][i]=dest[m][n+1]; done=1; } - } else { - if (tvf) dest2[j++]=tvec[i]; else dest2[j++]=i; } } } + if ((i+1)nq[0][i+1]) { // look at RT side + if (DEBUG_VECST) printf("RT problem %d %g %g>%g\n",i,nq[0][i],nq[2][i],nq[0][i+1]); + for (m=lc-1,done=0;m>=0 && !done;m--) { // go from bottom to top + for (n=1;nmaxsz) { - printf("%d > %d\n",size,maxsz); + printf("%g > %g\n",size,maxsz); hoc_execerror("v.snap: insufficient room in dest", 0); } - vector_resize((IvocVect*)vv, size); + vector_resize(vv, size); if (nsrc!=ntvec) hoc_execerror("v.snap: src and tvec not same size", 0); - for (tt=0,i=0;itt) dest[i]=src[j-1]; else { for (;jval) val=src[j]; if (val==-1e9) printf("vecst:snap() internal ERROR\n"); @@ -1701,64 +1405,30 @@ static double snap (void* vv) { } ENDVERBATIM -:* v1.xzero(vec[,thresh]) finds indices of zero [or thresh] crossings in vec +:* v1.xzero() looks for zero crossings VERBATIM -static double xzero (void* vv) { - int i, n, nv, up, cnt; - double *x, *vc, th; +static double xzero(void* vv) { + int i, n, nv, up; + double *x, *vc, th, cnt=0.; n = vector_instance_px(vv, &x); nv = vector_arg_px(1, &vc); if (ifarg(2)) { th = *getarg(2); } else { th=0.0;} if (vc[0]th) { - up=1; - if (cnt>=nv) x=vector_newsize((IvocVect*)vv,(n+=100)); - x[cnt++]=(double)i; - } - } - x=vector_newsize((IvocVect*)vv,cnt); - return (double)cnt; -} -ENDVERBATIM - -:* v1.peak(vec[,vamp]) puts indices of zero [or thresh] crossings of vec into v1 -: optional vamp gives amplitues -VERBATIM -static double peak (void* vv) { - int i, n, nc, ny, up, cnt; - double *x, *y, *vc, last; IvocVect* vy; - n = vector_instance_px(vv, &x); // for indices - nc = vector_arg_px(1, &vc); // source vectors - if (n==0) x=vector_newsize((IvocVect*)vv,n=100); - if (ifarg(2)) y=vector_newsize(vy=vector_arg(2),ny=n); else ny=0; // same size as v1 - if (vc[1]-vc[0]<0) up=0; else up=1; // F or T - last=vc[1]; - for (i=2,cnt=0; i=n) { x=vector_newsize((IvocVect*)vv,(n+=100)); - if (ny) y=vector_newsize(vy,n);} - x[cnt]=(double)(i-1); - if (ny) y[cnt]=last; - cnt++; - } - } else if (vc[i]-last>0) up=1; - last=vc[i]; + if (vc[i]th) { up=x[i]=1; cnt++; } } - x=vector_newsize((IvocVect*)vv,cnt); - if (ny) y=vector_newsize(vy,cnt); - return (double)cnt; + return cnt; } ENDVERBATIM :* v1.negwrap([FLAG]) wrap neg values to pos, FLAG==0 set them to 0, FLAG!=0 wrap : above FLAG VERBATIM -static double negwrap (void* vv) { +static double negwrap(void* vv) { int i, n; double *x, cnt, sig; n = vector_instance_px(vv, &x); @@ -1785,7 +1455,7 @@ ENDVERBATIM :* v1.sw(FROM,TO) switchs all FROMs to TO VERBATIM -static double sw (void* vv) { +static double sw(void* vv) { int i, n; double *x, fr, to; n = vector_instance_px(vv, &x); @@ -1800,7 +1470,7 @@ ENDVERBATIM :* v.b2v(bytevec) copies from vector to bytevec VERBATIM -static double b2v (void* vv) { +static double b2v(void* vv) { int i, n, num; double *x; bvec* to; Object *ob; n = vector_instance_px(vv, &x); @@ -1814,7 +1484,7 @@ ENDVERBATIM :* v.v2d(&x) copies from vector to double area -- a seg error waiting to happen VERBATIM -static double v2d (void* vv) { +static double v2d(void* vv) { int i, n, num; double *x, *to; n = vector_instance_px(vv, &x); @@ -1826,177 +1496,9 @@ static double v2d (void* vv) { } ENDVERBATIM -:* v.v2p(&x0[,&x1,&x2...]) copies from vector or vectors to doubles -VERBATIM -static double v2p (void* vv) { - int i, j, n, cnt; - double *x; - n = vector_instance_px(vv, &x); - for (i=0,j=1,cnt=0;ifarg(j) && i=cnt) {printf("vecst:l2p() ERRA: %d %d %d\n",i,ix,cnt); hxe();} - x[i]=y[ix]; - } - for (i=0,j=3,cnt=0;ifarg(j) && ipL->isz) vector_resize((IvocVect*)vector_arg(3),pL->isz); // don't make bigger if only want a few - for (i=0,j=0,cnt=0;iisz && jpL->plen[i]) {printf("vecst:fetch()ERRB: %d %d %x\n",i,ix,pL->pv[i]); - FreeListVec(&pL); hxe();} - y[j]=pL->pv[i][ix]; - cnt++; - } - ret=y[j-1]; // final value - } else { - for (i=0,j=3,cnt=0;iisz && ifarg(j);i++,j++) { - if (hoc_is_double_arg(j)) continue; // skip this one - if (ix>pL->plen[i]) {printf("vecst:fetch()ERRB1: %d %d %x\n",i,ix,pL->pv[i]); - FreeListVec(&pL); hxe();} - *hoc_pgetarg(j)=ret=pL->pv[i][ix]; - cnt++; - } - } - FreeListVec(&pL); - return ret; -} -ENDVERBATIM - - -:* v.covar(list,vec) generates covariance matrix in vec form -: generally data are in the columns; here data are in the vectors -- ie as if in the rows -VERBATIM -static double covar (void* vv) { - int ix, i, j, j2, n, m; - double *x, *y, *mean; ListVec* pL; - n = vector_instance_px(vv, &x); // number of data vectors - if (n==0) { pL = AllocListVec(*hoc_objgetarg(1)); // get all of them - } else pL = AllocILV(*hoc_objgetarg(1),n,x); - if (pL->isz<2) {printf("vecst:covar()ERRA: %d\n",pL->isz); FreeListVec(&pL); hxe();} - n=pL->isz; // number of data points - m=pL->plen[0]; // dimensionality of data - for (i=1;iisz;i++) if (m!=pL->plen[i]) { - printf("vecst:covar()ERRB: sz mismatch %d %d@%d\n",m,pL->plen[i],i);FreeListVec(&pL);hxe();} - y=vector_newsize((IvocVect*)vector_arg(2),m*m); - // pL->pv[i][j] -- i goes through the list and j goes through each vector - mean=(double*)malloc(sizeof(double)*m); - for (j=0;jpv[i][j]; - mean[j]/=(double)n; - } - for (i=0;ipv[i][j] -= mean[j]; // center the vectors - for (j=0;jpv[i][j]*pL->pv[i][j2]; - y[j*m+j2]/=(n-1); - y[j2*m+j]=y[j*m+j2]; - } - for (i=0;ipv[i][j] += mean[j]; // restore vectors - free(mean); - FreeListVec(&pL); - return m; -} -ENDVERBATIM - -:* ind.vlxpose(src_list,dest_list) does 'transpose' -VERBATIM -static double vlxpose (void* vv) { - int i, j, k, n, c, c2, sz, err; - double *x; ListVec *pL, *pL2; Object* obl; - err=0; - n = vector_instance_px(vv, &x); // vector of indices - if (n==0) { pL = AllocListVec(*hoc_objgetarg(1)); // get all of them - } else pL = AllocILV(*hoc_objgetarg(1),n,x); - pL2 = AllocListVec(obl=*hoc_objgetarg(2)); - c=pL->isz; c2=pL2->isz; // number of columns (list length) - // pL->pv[i][j] -- i goes through the list and j goes through each vector - for (j=0;jpbuflen[j]); - n=pL->plen[0]; // length of vector - if (n!=c2) err=1; - for (j=1;jplen[j] || err) { - printf("vecst:vlxpose()ERRA: %d %d %d\n",n,pL->plen[j],c2); - FreeListVec(&pL); FreeListVec(&pL2); hxe(); } - for (j=0,k=0;j=pL2->pbuflen[i]) { sz=pL2->pbuflen[i]; // need to grow vector - sz=(sz<10)?100:(sz*2); - pL2->pv[i]=list_vector_resize(obl, i, pL2->pbuflen[i]=sz); - } - pL2->pv[i][k]=pL->pv[j][i]; - } - for (j=0;jisz<2) {printf("vecst:ixsort()ERRA: %d\n",pL->isz); FreeListVec(&pL); hxe();} - c=pL->isz; // number of columns (list length) - // pL->pv[i][j] -- i goes through the list and j goes through each vector - for (j=0;jpbuflen[j]); - for (j=0;j=c) {printf("vecst:ixsort()ERRB: OOB %d %d\n",i,c); FreeListVec(&pL); hxe();} - if (pL->plen[i]>=pL->pbuflen[i]) { - if (pL->pbuflen[i]) pL->pbuflen[i]*=2; else pL->pbuflen[i]=100; - pL->pv[i]=list_vector_resize(obl, i, pL->pbuflen[i]); - } - pL->pv[i][pL->plen[i]++]=tv[j]; - } - for (j=0;jplen[j]); - FreeListVec(&pL); - return (double)n; -} -ENDVERBATIM - :* v.d2v(&x) copies from double area to vector -- a seg error waiting to happen VERBATIM -static double d2v (void* vv) { +static double d2v(void* vv) { int i, n, num; double *x, *fr; n = vector_instance_px(vv, &x); @@ -2010,13 +1512,13 @@ ENDVERBATIM :* v.lcat(LIST) VERBATIM -static double lcat (void* vv) { +static double lcat(void* vv) { int i, j, k, n, lc, cap, maxsz; Object *ob1; double *x, *fr; - IvocVect *vw; + void *vw; n = vector_instance_px(vv, &x); - vector_resize((IvocVect*)vv,maxsz=vector_buffer_size((IvocVect*)vv)); // open it up fully + vector_resize(vv,maxsz=vector_buffer_size(vv)); // open it up fully ob1 = *hoc_objgetarg(1); lc = ivoc_list_count(ob1); for (i=0,j=0;i5) hoc_execerror("uncode ****ERRA****: can only handle 5 vectors", 0); - for (i=0;iu.this_pointer); - *px = vector_vec((IvocVect*)obv->u.this_pointer); + sz = vector_capacity(obv->u.this_pointer); + *px = vector_vec(obv->u.this_pointer); return sz; } //* list_vector_px2(LIST,ITEM#,DOUBLE PTR ADDRESS,VEC POINTER ADDRESS) // returns the vector pointer as well as the double pointer -int list_vector_px2 (Object *ob, int i, double** px, IvocVect** vv) { +int list_vector_px2 (Object *ob, int i, double** px, void** vv) { Object* obv; int sz; obv = ivoc_list_item(ob, i); - if (! ISVEC(obv)) return -1; - sz = vector_capacity((IvocVect*)obv->u.this_pointer); - *px = vector_vec((IvocVect*)obv->u.this_pointer); - *vv = (IvocVect*) obv->u.this_pointer; + sz = vector_capacity(obv->u.this_pointer); + *px = vector_vec(obv->u.this_pointer); + *vv = (void*) obv->u.this_pointer; return sz; } //* list_vector_px3(LIST,ITEM#,DOUBLE PTR ADDRESS,VEC POINTER ADDRESS) // same as px2 but returns max vec size instead of current vecsize // side effect -- increase vector size to maxsize -int list_vector_px3 (Object *ob, int i, double** px, IvocVect** vv) { +int list_vector_px3 (Object *ob, int i, double** px, void** vv) { Object* obv; int sz; obv = ivoc_list_item(ob, i); - if (! ISVEC(obv)) return -1; - sz = vector_buffer_size((IvocVect*)obv->u.this_pointer); - *px = vector_vec((IvocVect*)obv->u.this_pointer); - *vv = (IvocVect*) obv->u.this_pointer; + sz = vector_buffer_size(obv->u.this_pointer); + *px = vector_vec(obv->u.this_pointer); + *vv = (void*) obv->u.this_pointer; vector_resize(*vv,sz); return sz; } -//* list_vector_px4(LIST,ITEM#,DOUBLE PTR ADDRESS,desired size) -// does resizing and returns true -int list_vector_px4 (Object *ob, int i, double** px, unsigned int n) { - Object* obv; - void* vv; - int sz; - obv = ivoc_list_item(ob, i); - if (! ISVEC(obv)) return -1; - sz = vector_buffer_size((IvocVect*)obv->u.this_pointer); - *px = vector_vec((IvocVect*)obv->u.this_pointer); - vv = (void*) obv->u.this_pointer; - if (n>sz) { - printf("List vector WARNING: unable to resize to %d requested (%d)\n",n,sz); - vector_resize((IvocVect*)vv,sz); - return 0; - } else vector_resize((IvocVect*)vv,n); - return 1; -} - //* list_vector_resize(LIST,ITEM#,NEW SIZE) -double *list_vector_resize (Object *ob, int i, int sz) { +int list_vector_resize (Object *ob, int i, int sz) { Object* obv; + int maxsz; obv = ivoc_list_item(ob, i); - if (! ISVEC(obv)) return 0x0; - vector_resize((IvocVect*)obv->u.this_pointer,sz); - return vector_vec((IvocVect*)obv->u.this_pointer); + maxsz = vector_buffer_size(obv->u.this_pointer); + if (sz>maxsz) { + printf("max:%d request:%d ",maxsz,sz); + hoc_execerror("Can't grow vector in list_vector_resize ", 0); + return -1; + } + vector_resize(obv->u.this_pointer,sz); + return sz; } ENDVERBATIM :* v1.ismono([arg]) asks whether is monotonically increasing, with arg==-1 - decreasing : with arg==0:all same; 2:no consec ==; 3: incrementing by 1 VERBATIM -double ismono1 (double *x, int n, int flag) { - int i; double last, gap, ret; - last=x[0]; ret=1.; // default return value +static int ismono1 (double *x, int n, int flag) { + int i; double last; + last=x[0]; if (flag==1) { for (i=1; i=last; i++) last=x[i]; } else if (flag==-1) { @@ -2209,16 +1694,13 @@ double ismono1 (double *x, int n, int flag) { for (i=1; ix[mid]) lt=mid+1; else if (num0) {printf("ERROR: uniq(list,vec)\n"); hxe();} - voi[1]=vector_arg(2); - if ((nz=openvec(2,&z))==0) z=vector_newsize(voi[1],nz=100); - } - } - scrset(n); - for (i=0;i0) z[0]=1.; - for (i=1, lastx=x[scr[0]], cnt=1; ilastx+hoc_epsilon) { - if (ny) { - if (cnt>=ny) y=vector_newsize(voi[0],ny*=3); - y[cnt]=x[scr[i]]; - } - if (nz>0) { - if (cnt>=nz) z=vector_newsize(voi[1],nz*=3); - z[cnt]=1.; - } - cnt++; - lastx=x[scr[i]]; - } else if (nz>0) z[cnt-1]++; - } - if (ny) vector_resize(voi[0], cnt); - if (nz>0) vector_resize(voi[1], cnt); - if (flag) { // refill z with the unique values in proper order - z=vector_newsize(voi[1], cnt); - for (i=0;iy[mid]) lt=mid+1; else if (numlastx+hoc_epsilon) { - y[cnt]=x[scr[i]]; - cnt++; - lastx=x[scr[i]]; - } - } - for (i=0;iy[mid]) lt=mid+1; else if (num100) {printf("nqsvt ERRD only 100 cols\n"); hxe();} - proc = gargstr(1); - if (!(s=hoc_lookup(proc))) {printf("nqsvt ERRA: proc %s not found\n",proc); hxe();} - fcdo=*hoc_objgetarg(2); - vector_arg_px(3, &fcd); - vl=*hoc_objgetarg(4); - if (ifarg(5)) {vector_arg_px(5, &ind); flag=1;} else flag=0; - n=list_vector_px(vl,(int)col[0],&vvo[0]); - for (i=1; i20 || x<-20) { printf("EXP(%g) called with OOB value [-20,20]\n",x) - EXP=ERR - } else { - Expo(x) - EXP=RES - } -} - -FUNCTION SUMEXP () { - VERBATIM - double i,min,max,step,sum; - if (ifarg(2)) { min=*getarg(1); max=*getarg(2); step=ifarg(3)?*getarg(3):1.; - } else { max=*getarg(1); min=0.; step=1.; } - if (max>20. || min<-20.) { - printf("SUMEXP() called with OOB value: %g %g [-20,20]\n",min,max); - sum=ERR; - } else for (i=min,sum=0;i<=max+hoc_epsilon;i+=step) { - Expo(i); - sum+=RES; - } - _lSUMEXP=sum; - ENDVERBATIM +FUNCTION AAA (x) { + Expo(x) + AAA = RES } :* dest.smgs(src,low,high,step,var) : rewrite of v.sumgauss() in nrn5.3::ivoc/ivocvect.cpp:1078 -: NEEDS DEBUGGING -- see drline.hoc:smgs() VERBATIM static double smgs (void* vv) { int i, j, nx, xv, nsum, points, maxsz; @@ -2523,9 +1771,9 @@ static double smgs (void* vv) { points = (int)((high-low)/step+hoc_epsilon); if (nsum!=points) { - maxsz=vector_buffer_size((IvocVect*)vv); + maxsz=vector_buffer_size(vv); if (points<=maxsz) { - nsum=points; vector_resize((IvocVect*)vv, nsum); + nsum=points; vector_resize(vv, nsum); } else { printf("%d > %d :: ",points,maxsz); hoc_execerror("Vector max capacity too small in smgs ", 0); @@ -2559,23 +1807,23 @@ VERBATIM static double smsy (void* vv) { int i, j, k, nx, nc, nsum, points, maxsz; double *x, *sum, *c; - double del,tstop,mdt; + double del,tstop,dtt; if (! ifarg(1)) { printf("dest.smsy(tvec,CVLV_VEC,tstop[,dt,del])\n"); return -1.; } - del=0.; mdt=0.2; + del=0.; dtt=0.2; nsum = vector_instance_px(vv, &sum); nx = vector_arg_px(1,&x); nc = vector_arg_px(2,&c); tstop = *getarg(3); - if (ifarg(4)) mdt = *getarg(4); + if (ifarg(4)) dtt = *getarg(4); if (ifarg(5)) del = *getarg(5); - points=(int)(tstop/mdt+hoc_epsilon); + points=(int)(tstop/dtt+hoc_epsilon); if (nsum!=points) { - maxsz=vector_buffer_size((IvocVect*)vv); + maxsz=vector_buffer_size(vv); if (points<=maxsz) { - vector_resize((IvocVect*)vv, points); points=nsum; + vector_resize(vv, points); points=nsum; } else { printf("%d > %d :: ",points,maxsz); hoc_execerror("Dest vector too small in smsy ", 0); @@ -2583,7 +1831,7 @@ static double smsy (void* vv) { } // don't zero out dest vec - for (i=0;iscrsz*sizeof(int)) { - if (scrsz>0) { free(scr); scr=(unsigned int *)NULL; } - scr=(unsigned int *)ecalloc(1, sz); - scrsz=sz/sizeof(int); // number of chars available - } - - xc=(char *)scr; - r=fread(xc,(size_t)sz,1,f); - - if (vflag) { - ny = vector_arg_px(2, &y); - } else { - num = ivoc_list_count(ob); - if (num>10000) { printf("rdfile ERRA: can only handle 10000 vectors"); hxe();} - for (i=0;i5) { - printf("rdfile ERRB: bad size/type: %d/%d in vec# %d\n",vsz,ty,cnt); hxe();} - if (DEBUG_VECST) printf("%d:%d:%d ",i,ty,vsz); - if (vflag) { // vector - if (k+vsz>=ny) { - printf("rdfile ERRC: No more room in vec: %d %d %d %d\n",ny,k+vsz,cnt,ty); hxe(); } - } else { // a list - if (cnt>=num) { - printf("rdfile ERRD: out of vecs: %d %d %d\n",num,cnt,ty); hxe(); } - if (vsz>nv) { - printf("rdfile ERRE: No more room in vec: %d %d %d %d\n",nv,vsz,cnt,ty); hxe(); } - } - if (ty==3) { // float must be recast - xf=(float*)(xc+i); - if (vflag) { - for (j=0;j1e7) { free(scr); scr=(unsigned int *)NULL; scrsz=0; } - return (double)num; -} -ENDVERBATIM - :* PROCEDURE install_vecst() PROCEDURE install_vecst () { if (VECST_INSTALLED==1) { - printf("$Id: vecst.mod,v 1.494 2010/12/03 17:26:01 billl Exp $\n") + printf("$Id: vecst.mod,v 1.326 2006/06/24 23:38:36 billl Exp $\n") } else { VECST_INSTALLED=1 VERBATIM { - int i,j; + int i,j; install_vector_method("indset", indset); - install_vector_method("mkind", mkind); - install_vector_method("circ", circ); install_vector_method("thresh", thresh); install_vector_method("triplet", triplet); install_vector_method("onoff", onoff); install_vector_method("bpeval", bpeval); install_vector_method("w", w); - install_vector_method("whi", whi); install_vector_method("sedit", sedit); install_vector_method("xing", xing); + install_vector_method("updown", updown); install_vector_method("scxing", scxing); install_vector_method("cvlv", cvlv); install_vector_method("sccvlv", sccvlv); install_vector_method("scl", scl); - install_vector_method("revec", revec); - install_vector_method("has", has); install_vector_method("intrp", intrp); install_vector_method("xzero", xzero); - install_vector_method("peak", peak); install_vector_method("negwrap", negwrap); install_vector_method("sw", sw); install_vector_method("ismono", ismono); install_vector_method("count", count); - install_vector_method("muladd", muladd); - install_vector_method("binfind", binfind); - install_vector_method("unq", unq); - install_vector_method("uniq", uniq); install_vector_method("rnd", rnd); install_vector_method("fewind", fewind); install_vector_method("findx", findx); - install_vector_method("lma", lma); install_vector_method("sindx", sindx); install_vector_method("sindv", sindv); install_vector_method("nind", nind); @@ -2864,18 +2004,11 @@ PROCEDURE install_vecst () { install_vector_method("slct", slct); install_vector_method("slor", slor); install_vector_method("insct", insct); - install_vector_method("linsct", linsct); install_vector_method("cull", cull); install_vector_method("redundout", redundout); install_vector_method("mredundout", mredundout); install_vector_method("d2v", d2v); install_vector_method("v2d", v2d); - install_vector_method("v2p", v2p); - install_vector_method("l2p", l2p); - install_vector_method("fetch", fetch); - install_vector_method("covar", covar); - install_vector_method("ixsort", ixsort); - install_vector_method("vlxpose", vlxpose); install_vector_method("b2v", b2v); install_vector_method("iwr", iwr); install_vector_method("ird", ird); @@ -2895,20 +2028,13 @@ PROCEDURE install_vecst () { install_vector_method("slone", slone); install_vector_method("pop", pop); install_vector_method("rdmany", rdmany); - install_vector_method("rdfile", rdfile); - install_vector_method("samp", samp); - install_vector_method("nearest", nearest); - install_vector_method("nearall", nearall); - install_vector_method("approx", approx); - install_vector_method("nqsvt", nqsvt); - install_vector_method("roton", roton); for (i=0,j=5;i<=5;i++,j--) sc[i]=pow(2,10*j); } ENDVERBATIM } } -:* isojt(OB1,EXAMPLE_OBJ) return whether OB1 is an instance of EXAMPLE_OBJ +: isojt(OB1,EXAMPLE_OBJ) return whether OB1 is an instance of EXAMPLE_OBJ FUNCTION isojt () { VERBATIM { Object *ob1, *ob2; @@ -2973,13 +2099,13 @@ FUNCTION eqojt () { :* byteswap(FILE) FUNCTION byteswap () { VERBATIM { - int n[2]; size_t r; + int n[2]; double ret; FILE* f; BYTEHEADER f = hoc_obj_file_arg(1); - r=fread(&n,sizeof(int),2,f); + fread(&n,sizeof(int),2,f); if (n[1] < 1 || n[1] > 5) { BYTESWAP_FLAG = 1; ret = 1.; @@ -3001,11 +2127,10 @@ FUNCTION mkcodf () { VERBATIM { int i; double x,a; - if (ifarg(6)) {printf("mkcodf() ERR: can only encode 5 values\n"); hxe();} for (x=0.,i=1;i<=5;i++) { - a=(ifarg(i))?*getarg(i):0.0; - if (a<0. || a>=sc[4] || floor(a+0.5)!=a) { - printf("mkcodf restricted to integers %g [0,%g]\n",a,sc[4]-1);hxe(); } + a=chkarg(i,0.,sc[4]-1); + if (floor(a+0.5)!=a) { + printf("mkcodf restricted to integers %g\n",a);hxe();} x+=a*sc[i]; } return x; @@ -3020,8 +2145,7 @@ FUNCTION uncodf () { double x,ret, *ptr; x=*getarg(1); if (hoc_is_double_arg(2)) { - i=(int)*getarg(2); - if (i<1||i>5) {printf("2nd arg must be field# 1-5 (%d)\n",i); hxe();} + i=(int)chkarg(2,1.,5.); UNCODE(x,i,ret); return ret; } else { @@ -3054,215 +2178,3 @@ FUNCTION flor () { } ENDVERBATIM } - -: ceilg(val) -FUNCTION ceilg () { - VERBATIM { - return ceil(*getarg(1)); - } - ENDVERBATIM -} - -: MINxy(val1,val2) -FUNCTION MINxy () { - VERBATIM { - return MIN(*getarg(1),*getarg(2)); - } - ENDVERBATIM -} - -: MAXxy(val1,val2) -FUNCTION MAXxy () { - VERBATIM { - return MAX(*getarg(1),*getarg(2)); - } - ENDVERBATIM -} - -:* PROCEDURE fspitchar -PROCEDURE fspitchar(c) { -VERBATIM -{ - FILE* f; - f = hoc_obj_file_arg(2); - fprintf(f, "%c", (int)_lc); -} -ENDVERBATIM -} - -:* PROCEDURE fgchar -FUNCTION fgchar() { -VERBATIM -{ - FILE* f; - f = hoc_obj_file_arg(1); - _lfgchar = (double)fgetc(f); -} -ENDVERBATIM -} - -:* FUNCTION Str2Num takes a string arg and returns the # as a double -FUNCTION Str2Num () { -VERBATIM -{ - double d; - char* c; - c = gargstr(1); - d = atof(c); - return d; -} -ENDVERBATIM -} - -:* FUNCTION vlsz() resize all the vectors in a list -FUNCTION vlsz () { -VERBATIM -{ - int i,j,c,n; double *x, sz, fill; void *vv; - ListVec* pL; Object* obl; - pL = AllocListVec(obl=*hoc_objgetarg(1)); - sz=*getarg(2); - if (ifarg(3)) fill=*getarg(3); else fill=OK; - c=pL->isz; // list length - for (i=0;ipv[i]=list_vector_resize(obl, i, (int)sz); - if (fill!=OK) for (j=0;j<(int)sz;j++) pL->pv[i][j]=fill; - } - FreeListVec(&pL); - _lvlsz = (double)sz*c; -} -ENDVERBATIM -} - -VERBATIM -void FreeListVec(ListVec** pp) { - ListVec* p = *pp; - if(p->pv){ - free(p->pv); - p->pv=0; - } - if(p->plen){ - free(p->plen); - p->plen=0; - } - free(p); - *pp=0; -} - -ListVec* AllocListVec (Object* p) { - int i, iSz; ListVec* pList; Object* obv; - if(!IsList(p)){printf("AllocListVec ERRA: arg must be list object!\n"); hxe();} - pList = (ListVec*)malloc(sizeof(ListVec)); - if(!pList) hxe(); - pList->pL=p; pList->isz=0; pList->pv=0; pList->plen=0; - iSz = pList->isz = ivoc_list_count(p); - if(iSz < 1) return pList; - pList->plen = (unsigned int*)malloc(sizeof(int)*iSz); - if(!pList->plen){printf("AllocListVec ERRB: Out of memory!\n"); hxe();} - pList->pbuflen = (unsigned int*)malloc(sizeof(int)*iSz); - pList->pv = (double**)malloc(sizeof(double*)*iSz); - if(!pList->pv){free(pList->plen); printf("AllocListVec ERRC: Out of memory!\n"); hxe();} - for(i=0;iisz;i++) { - obv = ivoc_list_item(p,i); - pList->pv[i]=vector_vec((IvocVect*)obv->u.this_pointer); - pList->plen[i]=vector_capacity((IvocVect*)obv->u.this_pointer); - pList->pbuflen[i]=vector_buffer_size((IvocVect*)obv->u.this_pointer);; - } - return pList; -} - -// Allocate a list vec that is indexed -ListVec* AllocILV (Object* p, int nx, double *x) { - int i, j, iSz, ilc; ListVec* pList; Object* obv; - if(!IsList(p)){printf("AllocILV ERRA: arg must be list object!\n"); hxe();} - pList = (ListVec*)malloc(sizeof(ListVec)); - if(!pList) hxe(); - pList->pL=p; iSz=pList->isz=nx; pList->pv=0; pList->plen=0; - ilc=ivoc_list_count(p); - if(iSz<1) return pList; - pList->plen = (unsigned int*)malloc(sizeof(int)*iSz); - if(!pList->plen){printf("AllocILV ERRB: Out of memory!\n"); hxe();} - pList->pbuflen = (unsigned int*)malloc(sizeof(int)*iSz); - pList->pv = (double**)malloc(sizeof(double*)*iSz); - if(!pList->pv){free(pList->plen); printf("AllocILV ERRC: Out of memory!\n"); hxe();} - for(i=0;i=ilc){printf("AllocILV ERRD: index OOB: %d>=%d\n",j,ilc); hxe();} - obv = ivoc_list_item(p,j); - pList->pv[i]=vector_vec((IvocVect*)obv->u.this_pointer); - pList->plen[i]=vector_capacity((IvocVect*)obv->u.this_pointer); - pList->pbuflen[i]=vector_buffer_size((IvocVect*)obv->u.this_pointer);; - } - return pList; -} - -void ListVecResize (ListVec* p,int newsz) { - int i,j; Object* obv; - for(i=0;iisz;i++){ - obv = ivoc_list_item(p->pL, i); - p->pv[i]=vector_newsize((IvocVect*)obv->u.this_pointer,newsz); - p->plen[i]=newsz; - } -} - -void FillListVec (ListVec* p,double dval) { - int i,j; - for(i=0;iisz;i++){ - for(j=0;jplen[i];j++){ - p->pv[i][j]=dval; - } - } -} - -int IsObj (Object* p,char* s){ - if(!p) return 0; - if(!s || !strlen(s)) return 0; - return !strncmp(hoc_object_name(p),s,strlen(s)); -} -int IsVector (Object* p){ return IsObj(p,"Vector"); } -int IsList (Object* p){return IsObj(p,"List"); } - -int** getint2D(int rows,int cols) { - int **pp,*pool,*curPtr; int i; - pp = (int**) malloc(sizeof(int*)*rows); - if(!pp) { printf("ERR: out of memory!\n"); return 0x0; } - pool = (int*) malloc(sizeof(int)*rows*cols); - if(!pool) { printf("ERR: out of memory!\n"); free(pp); return 0x0; } - curPtr = pool; - for(i = 0; i < rows; i++) { - pp[i] = curPtr; - curPtr += cols; - } - return pp; -} - -void freeint2D(int*** ppp,int rows) { - int** pp; - pp = *ppp; - free(pp[0]); - free(pp); - *ppp = 0; -} - -double** getdouble2D(int rows,int cols) { - double **pp,*pool,*curPtr; int i; - pp = (double**) malloc(sizeof(double*)*rows); - if(!pp) { printf("ERR: out of memory!\n"); return 0x0; } - pool = (double*) malloc(sizeof(double)*rows*cols); - if(!pool) { printf("ERR: out of memory!\n"); free(pp); return 0x0; } - curPtr = pool; - for(i = 0; i < rows; i++) { - pp[i] = curPtr; - curPtr += cols; - } - return pp; -} - -void freedouble2D(double*** ppp,int rows) { - double** pp; - pp = *ppp; - free(pp[0]); - free(pp); - *ppp = 0; -} - -ENDVERBATIM From e9011c1851083bc68f9f116b95809cd72187716d Mon Sep 17 00:00:00 2001 From: Marianne Bezaire <6456664+mbezaire@users.noreply.github.com> Date: Sun, 7 Aug 2022 01:12:32 -0400 Subject: [PATCH 5/9] python 2->3 + nrnivmodl bugs Compiles on WSL on Windows now, somewhat runs --- .gitignore | 3 +++ init.hoc | 6 +++-- network.py | 60 +++++++++++++++++++++---------------------- networkmsj.py | 70 +++++++++++++++++++++++++-------------------------- run.py | 40 ++++++++++++++--------------- stats.mod | 9 ++++--- vecst.mod | 1 + 7 files changed, 98 insertions(+), 91 deletions(-) diff --git a/.gitignore b/.gitignore index 686d33e..e91f90d 100644 --- a/.gitignore +++ b/.gitignore @@ -1,3 +1,6 @@ *.c *.o *.dll +*.tmp +x86_64/* +__pycache__/* diff --git a/init.hoc b/init.hoc index 3f43299..48ee36e 100644 --- a/init.hoc +++ b/init.hoc @@ -2,7 +2,7 @@ xwindows=1 show_panel=0 -pwd() +//pwd() - mb no longer avail load_file("xgetargs.hoc") load_file("grvec.hoc") load_file("labels.hoc") @@ -10,7 +10,9 @@ load_file("syncode.hoc") load_file("decnqs.hoc") if (! VECST_INSTALLED) install_vecst() -if (! INSTALLED_stats) install_stats() + +install_stats() +//if (! INSTALLED_stats) install_stats() //install_PLACE() //install_clust() diff --git a/network.py b/network.py index b9efe7b..5c85620 100644 --- a/network.py +++ b/network.py @@ -93,7 +93,7 @@ class Network: def __init__(self,noise=True,connections=True,DoMakeNoise=True,iseed=1234,UseNetStim=True,wseed=4321,scale=1.0,MSGain=1.0,SaveConn=False): import math - print "Setting Cells" + print("Setting Cells") self.pyr = Population(cell_type=PyrAdr,n=int(math.ceil(800*scale)), x= 0, y=0, z=0, dx=50, amp= 50e-3, dur=1e9, delay=2*h.dt) self.bas = Population(cell_type=Bwb, n=int(math.ceil(200*scale)), x=10, y=0, z=0, dx=50, amp= 0, dur= 0, delay=2*h.dt) self.olm = Population(cell_type=Ow, n=int(math.ceil(200*scale)), x=20, y=0, z=0, dx=50, amp=-25e-3, dur=1e9, delay=2*h.dt) @@ -111,7 +111,7 @@ def __init__(self,noise=True,connections=True,DoMakeNoise=True,iseed=1234,UseNet self.SaveConn = SaveConn if connections: - print "Setting Connections" + print("Setting Connections") self.set_all_conns() @@ -123,21 +123,21 @@ def set_noise_inputs(self,simdur): #simdur only used for make_all_noise self.make_all_noise(simdur,self.iseed) else: self.load_all_noise() - print "Done!" + print("Done!") def load_all_noise(self): #load noise from data files - print "Loading Noise" - print "to PYR" + print("Loading Noise") + print("to PYR") self.b_pyr_somaAMPAf=self.load_spikes("spike_noise_pyr_800_soma_AMPA_ISI_1_N_10000_noise_1.npy",self.pyr,"somaAMPAf",0.05e-3) self.b_pyr_Adend3AMPAf=self.load_spikes("spike_noise_pyr_800_Adend3_AMPA_ISI_1_N_10000_noise_1.npy",self.pyr,"Adend3AMPAf",0.05e-3) self.b_pyr_somaGABAf=self.load_spikes("spike_noise_pyr_800_soma_GABA_ISI_1_N_10000_noise_1.npy",self.pyr,"somaGABAf",0.012e-3) self.b_pyr_Adend3GABAf=self.load_spikes("spike_noise_pyr_800_Adend3_GABA_ISI_1_N_10000_noise_1.npy",self.pyr,"Adend3GABAf",0.012e-3) self.b_pyr_Adend3NMDA=self.load_spikes("spike_noise_pyr_800_Adend3_NMDA_ISI_100_N_100_noise_1.npy", self.pyr,"Adend3NMDA",6.5e-3) - print "to BAS" + print("to BAS") self.b_bas_somaAMPAf=self.load_spikes("spike_noise_bas_200_soma_AMPA_ISI_1_N_10000_noise_1.npy",self.bas,"somaAMPAf",w=0.02e-3) self.b_bas_somaGABA=self.load_spikes("spike_noise_bas_200_soma_GABA_ISI_1_N_10000_noise_1.npy",self.bas,"somaGABAf",w=0.2e-3) self.b_bas_somaGABAf=self.load_spikes("spike_noise_bas_200_soma_GABAf_ISI_150_N_65_noise_0.npy",self.bas,"somaGABAss",w=1.6e-3) - print "to OLM" + print("to OLM") self.b_olm_somaAMPAf=self.load_spikes("spike_noise_olm_200_soma_AMPA_ISI_1_N_10000_noise_1.npy",self.olm,"somaAMPAf",w=0.02e-3) self.b_olm_somaGABAf=self.load_spikes("spike_noise_olm_200_soma_GABA_ISI_1_N_10000_noise_1.npy",self.olm,"somaGABAf",w=0.2e-3) self.b_olm_somaGABAss=self.load_spikes("spike_noise_olm_200_soma_GABAss_ISI_150_N_65_noise_0.npy",self.olm,"somaGABAss",w=1.6e-3) @@ -236,25 +236,25 @@ def mkspkh(self,binsz): snq.verbose=1 def make_all_NetStims(self,simdur,rdmseed): - print "Making NetStims" + print("Making NetStims") # h.mcell_ran4_init(self.iseed) self.nsl = [] #NetStim List self.ncl = [] #NetCon List self.nrl = [] #Random List for NetStims self.nrlsead = [] #List of seeds for NetStim randoms # numpy.random.seed(rdmseed) # initialize random # generator - print "Making Noise" - print "to PYR" + print("Making Noise") + print("to PYR") rdtmp = rdmseed # starting sead value - incremented in make_NetStims rdtmp=self.make_NetStims(po=self.pyr, syn="somaAMPAf", w=0.05e-3, ISI=1, time_limit=simdur, sead=rdtmp) rdtmp=self.make_NetStims(po=self.pyr, syn="Adend3AMPAf", w=0.05e-3, ISI=1, time_limit=simdur, sead=rdtmp) rdtmp=self.make_NetStims(po=self.pyr, syn="somaGABAf", w=0.012e-3, ISI=1, time_limit=simdur, sead=rdtmp) rdtmp=self.make_NetStims(po=self.pyr, syn="Adend3GABAf", w=0.012e-3, ISI=1, time_limit=simdur, sead=rdtmp) rdtmp=self.make_NetStims(po=self.pyr, syn="Adend3NMDA", w=6.5e-3, ISI=100,time_limit=simdur, sead=rdtmp) - print "to BAS" + print("to BAS") rdtmp=self.make_NetStims(po=self.bas, syn="somaAMPAf", w=0.02e-3, ISI=1, time_limit=simdur, sead=rdtmp) rdtmp=self.make_NetStims(po=self.bas, syn="somaGABAf", w=0.2e-3, ISI=1, time_limit=simdur, sead=rdtmp) - print "to OLM" + print("to OLM") #rdtmp=self.make_NetStims(po=self.olm, syn="somaAMPAf", w=0.02e-3, ISI=1, time_limit=simdur, sead=rdtmp) rdtmp=self.make_NetStims(po=self.olm, syn="somaAMPAf", w=0.0625e-3, ISI=1, time_limit=simdur, sead=rdtmp) rdtmp=self.make_NetStims(po=self.olm, syn="somaGABAf", w=0.2e-3, ISI=1, time_limit=simdur, sead=rdtmp) @@ -278,19 +278,19 @@ def make_all_NetStims(self,simdur,rdmseed): def make_all_noise(self,simdur,rdmseed): # create noise for simdur milliseconds numpy.random.seed(rdmseed) # initialize random # generator import math - print "Making Noise" + print("Making Noise") fctr = (simdur+simdur/2) / 10000.0 - print "to PYR" + print("to PYR") self.b_pyr_somaAMPAf=self.make_spikes(self.pyr,"somaAMPAf",0.05e-3,self.pyr.n,"soma",1,math.ceil(10000*fctr),1,simdur) self.b_pyr_Adend3AMPAf=self.make_spikes(self.pyr,"Adend3AMPAf",0.05e-3,self.pyr.n,"Adend3",1,math.ceil(10000*fctr),1,simdur) self.b_pyr_somaGABAf=self.make_spikes(self.pyr,"somaGABAf",0.012e-3,self.pyr.n,"soma",1,math.ceil(10000*fctr),1,simdur) self.b_pyr_Adend3GABAf=self.make_spikes(self.pyr,"Adend3GABAf",0.012e-3,self.pyr.n,"Adend3",1,math.ceil(10000*fctr),1,simdur) self.b_pyr_Adend3NMDA=self.make_spikes(self.pyr,"Adend3NMDA",6.5e-3,self.pyr.n,"Adend3",100,math.ceil(100*fctr),1,simdur) - print "to BAS" + print("to BAS") self.b_bas_somaAMPAf=self.make_spikes(self.bas,"somaAMPAf",0.02e-3,self.bas.n,"soma",1,math.ceil(10000*fctr),1,simdur) self.b_bas_somaGABA=self.make_spikes(self.bas,"somaGABAf",0.2e-3,self.bas.n,"soma",1,math.ceil(10000*fctr),1,simdur) self.b_bas_somaGABAf=self.make_spikes(self.bas,"somaGABAss",1.6e-3,self.bas.n,"soma",150,math.ceil(65*fctr),0,simdur) - print "to OLM" + print("to OLM") self.b_olm_somaAMPAf=self.make_spikes(self.olm,"somaAMPAf",0.02e-3,self.olm.n,"soma",1,math.ceil(10000*fctr),1,simdur) self.b_olm_somaGABAf=self.make_spikes(self.olm,"somaGABAf",0.2e-3,self.olm.n,"soma",1,math.ceil(10000*fctr),1,simdur) self.b_olm_somaGABAss=self.make_spikes(self.olm,"somaGABAss",1.6e-3,self.olm.n,"soma",150,math.ceil(65*fctr),0,simdur) @@ -303,28 +303,28 @@ def make_conn(self, preN, postN, conv): def set_all_conns(self): random.seed(self.wseed) # initialize random # generator for wiring - print "PYR -> X , NMDA" # src, trg, syn, delay, weight, conv + print("PYR -> X , NMDA") # src, trg, syn, delay, weight, conv self.pyr_bas_NM=self.set_connections(self.pyr,self.bas, "somaNMDA", 2, 1.15*1.2e-3, 100) self.pyr_olm_NM=self.set_connections(self.pyr,self.olm, "somaNMDA", 2, 1.0*0.7e-3, 10) self.pyr_pyr_NM=self.set_connections(self.pyr,self.pyr, "BdendNMDA",2, 1*0.004e-3, 25) - print "PYR -> X , AMPA" + print("PYR -> X , AMPA") self.pyr_bas_AM=self.set_connections(self.pyr,self.bas, "somaAMPAf",2, 0.3*1.2e-3, 100) self.pyr_olm_AM=self.set_connections(self.pyr,self.olm, "somaAMPAf",2, 0.3*1.2e-3, 10) self.pyr_pyr_AM=self.set_connections(self.pyr,self.pyr, "BdendAMPA",2, 0.5*0.04e-3, 25) - print "BAS -> X , GABA" + print("BAS -> X , GABA") #self.bas_bas_GA=self.set_connections(self.bas,self.bas, "somaGABAf",2, 1.0e-3, 60)#orig 1 #self.bas_bas_GA=self.set_connections(self.bas,self.bas, "somaGABAf",2, 2 * 1.5*1.0e-3, 60)#new 1 self.bas_bas_GA=self.set_connections(self.bas,self.bas, "somaGABAf",2, 3 * 1.5*1.0e-3, 60)#new 2 self.bas_pyr_GA=self.set_connections(self.bas,self.pyr, "somaGABAf",2, 2 * 2*0.18e-3, 50)#new 1 - print "OLM -> PYR , GABA" + print("OLM -> PYR , GABA") #self.olm_pyr_GA=self.set_connections(self.olm,self.pyr, "Adend2GABAs",2, 3*6.0e-3, 20)#original weight value self.olm_pyr_GA=self.set_connections(self.olm,self.pyr, "Adend2GABAs",2, 4.0 * 3*6.0e-3, 20)#new weight value #pyramidal to PSR cell -- for testing only - print "PYR -> PSR, AMPA/NMDA" + print("PYR -> PSR, AMPA/NMDA") self.pyr_psr_NM=self.set_connections(self.pyr,self.psr, "BdendNMDA",2, 1*0.004e-3, 25) self.pyr_psr_AM=self.set_connections(self.pyr,self.psr, "BdendAMPA",2, 0.5*0.04e-3, 25) @@ -354,25 +354,25 @@ def set_connections(self,src,trg,syn,delay,w,conv): def load_spikes(self,fn,po,syn,w,time_limit=10000): fn = os.path.join("data",fn) events = numpy.load(fn) - print "Begin setting events...", po + print("Begin setting events..."), po print events.shape for i,ii in enumerate(events): ii=ii[ii<=time_limit] po.cell[i].__dict__[syn].append(ii) po.cell[i].__dict__[syn].syn.Vwt = w - print "End setting events" + print("End setting events") return events def make_spikes(self,po,syn,w,cellN,comp,ISI,eventN,noise,time_limit): events = numpy.random.exponential(ISI, (cellN,eventN))*noise+numpy.repeat(ISI,cellN*eventN).reshape((cellN,eventN))*(1-noise) events = numpy.cumsum(events,axis=1) - print "Begin setting events...", po + print("Begin setting events..."), po print events.shape for i,ii in enumerate(events): ii=ii[ii<=time_limit] po.cell[i].__dict__[syn].append(ii) po.cell[i].__dict__[syn].syn.Vwt = w - print "End setting events" + print("End setting events") return events def rasterplot(self,sz=2): @@ -449,9 +449,9 @@ def pravgrates(self,skipms=200): self.fnq.select("ty",ty) vf = self.fnq.getcol("freq") if vf.size() > 1: - print "ty: ", ty, " avg rate = ", vf.mean(), "+/-", vf.stderr(), " Hz" + print("ty: "), ty, " avg rate = ", vf.mean(), "+/-", vf.stderr(), " Hz" else: - print "ty: ", ty, " avg rate = ", vf.mean(), "+/-", 0.0 , " Hz" + print("ty: "), ty, " avg rate = ", vf.mean(), "+/-", 0.0 , " Hz" ty += 1 def calc_lfp(self): # lfp is modeled as a difference between voltages in distal apical and basal compartemnts @@ -480,7 +480,7 @@ def calc_psd(self,fig=3): gr = [30, 80] # Gamma frequency range t0i = int(t0/h.dt) if t0i > len(self.lfp): - print "LFP is too short! (<200 ms)" + print("LFP is too short! (<200 ms)") return 0,0,0,0,0,0 pylab.figure(fig) @@ -521,10 +521,10 @@ def get_lim_max(self, data, ind): # return the position of the maximal element fp.close() #create the network net = Network(noise=True,connections=True,DoMakeNoise=True,iseed=ISEED,UseNetStim=True,wseed=WSEED,scale=1.0,MSGain=MSG) - print "set network from rseed.txt : iseed=",ISEED,", WSEED=",WSEED,", MSG = ",MSG + print("set network from rseed.txt : iseed="),ISEED,", WSEED=",WSEED,", MSG = ",MSG except: net = Network() - print "set network from default constructor" + print("set network from default constructor") #setup some variables in hoc def sethocix(): diff --git a/networkmsj.py b/networkmsj.py index 8dbeeea..0e06895 100644 --- a/networkmsj.py +++ b/networkmsj.py @@ -93,7 +93,7 @@ class Network: def __init__(self,noise=True,connections=True,DoMakeNoise=True,iseed=1234,UseNetStim=True,wseed=4321,scale=1.0,MSGain=1.0,SaveConn=False): import math - print "Setting Cells" + print("Setting Cells") self.pyr = Population(cell_type=PyrAdr,n=int(math.ceil(800*scale)), x= 0, y=0, z=0, dx=50, amp= 50e-3, dur=1e9, delay=2*h.dt) self.bas = Population(cell_type=Bwb, n=int(math.ceil(200*scale)), x=10, y=0, z=0, dx=50, amp= 0, dur= 0, delay=2*h.dt) self.olm = Population(cell_type=Ow, n=int(math.ceil(200*scale)), x=20, y=0, z=0, dx=50, amp=-25e-3, dur=1e9, delay=2*h.dt) @@ -111,7 +111,7 @@ def __init__(self,noise=True,connections=True,DoMakeNoise=True,iseed=1234,UseNet self.SaveConn = SaveConn if connections: - print "Setting Connections" + print("Setting Connections") self.set_all_conns() @@ -123,21 +123,21 @@ def set_noise_inputs(self,simdur): #simdur only used for make_all_noise self.make_all_noise(simdur,self.iseed) else: self.load_all_noise() - print "Done!" + print("Done!") def load_all_noise(self): #load noise from data files - print "Loading Noise" - print "to PYR" + print("Loading Noise") + print("to PYR") self.b_pyr_somaAMPAf=self.load_spikes("spike_noise_pyr_800_soma_AMPA_ISI_1_N_10000_noise_1.npy",self.pyr,"somaAMPAf",0.05e-3) self.b_pyr_Adend3AMPAf=self.load_spikes("spike_noise_pyr_800_Adend3_AMPA_ISI_1_N_10000_noise_1.npy",self.pyr,"Adend3AMPAf",0.05e-3) self.b_pyr_somaGABAf=self.load_spikes("spike_noise_pyr_800_soma_GABA_ISI_1_N_10000_noise_1.npy",self.pyr,"somaGABAf",0.012e-3) self.b_pyr_Adend3GABAf=self.load_spikes("spike_noise_pyr_800_Adend3_GABA_ISI_1_N_10000_noise_1.npy",self.pyr,"Adend3GABAf",0.012e-3) self.b_pyr_Adend3NMDA=self.load_spikes("spike_noise_pyr_800_Adend3_NMDA_ISI_100_N_100_noise_1.npy", self.pyr,"Adend3NMDA",6.5e-3) - print "to BAS" + print("to BAS") self.b_bas_somaAMPAf=self.load_spikes("spike_noise_bas_200_soma_AMPA_ISI_1_N_10000_noise_1.npy",self.bas,"somaAMPAf",w=0.02e-3) self.b_bas_somaGABA=self.load_spikes("spike_noise_bas_200_soma_GABA_ISI_1_N_10000_noise_1.npy",self.bas,"somaGABAf",w=0.2e-3) self.b_bas_somaGABAf=self.load_spikes("spike_noise_bas_200_soma_GABAf_ISI_150_N_65_noise_0.npy",self.bas,"somaGABAss",w=1.6e-3) - print "to OLM" + print("to OLM") self.b_olm_somaAMPAf=self.load_spikes("spike_noise_olm_200_soma_AMPA_ISI_1_N_10000_noise_1.npy",self.olm,"somaAMPAf",w=0.02e-3) self.b_olm_somaGABAf=self.load_spikes("spike_noise_olm_200_soma_GABA_ISI_1_N_10000_noise_1.npy",self.olm,"somaGABAf",w=0.2e-3) self.b_olm_somaGABAss=self.load_spikes("spike_noise_olm_200_soma_GABAss_ISI_150_N_65_noise_0.npy",self.olm,"somaGABAss",w=1.6e-3) @@ -145,7 +145,7 @@ def load_all_noise(self): #load noise from data files #this should be called @ beginning of each sim - done in an FInitializeHandler in run.py def init_NetStims(self): # h.mcell_ran4_init(self.iseed) - for i in xrange(len(self.nrl)): + for i in range(len(self.nrl)): rds = self.nrl[i] sead = self.nrlsead[i] rds.MCellRan4(sead,sead) @@ -236,25 +236,25 @@ def mkspkh(self,binsz): snq.verbose=1 def make_all_NetStims(self,simdur,rdmseed): - print "Making NetStims" + print("Making NetStims") # h.mcell_ran4_init(self.iseed) self.nsl = [] #NetStim List self.ncl = [] #NetCon List self.nrl = [] #Random List for NetStims self.nrlsead = [] #List of seeds for NetStim randoms # numpy.random.seed(rdmseed) # initialize random # generator - print "Making Noise" - print "to PYR" + print("Making Noise") + print("to PYR") rdtmp = rdmseed # starting sead value - incremented in make_NetStims rdtmp=self.make_NetStims(po=self.pyr, syn="somaAMPAf", w=0.05e-3, ISI=1, time_limit=simdur, sead=rdtmp) rdtmp=self.make_NetStims(po=self.pyr, syn="Adend3AMPAf", w=0.05e-3, ISI=1, time_limit=simdur, sead=rdtmp) rdtmp=self.make_NetStims(po=self.pyr, syn="somaGABAf", w=0.012e-3, ISI=1, time_limit=simdur, sead=rdtmp) rdtmp=self.make_NetStims(po=self.pyr, syn="Adend3GABAf", w=0.012e-3, ISI=1, time_limit=simdur, sead=rdtmp) rdtmp=self.make_NetStims(po=self.pyr, syn="Adend3NMDA", w=6.5e-3, ISI=100,time_limit=simdur, sead=rdtmp) - print "to BAS" + print("to BAS") rdtmp=self.make_NetStims(po=self.bas, syn="somaAMPAf", w=0.02e-3, ISI=1, time_limit=simdur, sead=rdtmp) rdtmp=self.make_NetStims(po=self.bas, syn="somaGABAf", w=0.2e-3, ISI=1, time_limit=simdur, sead=rdtmp) - print "to OLM" + print("to OLM") #rdtmp=self.make_NetStims(po=self.olm, syn="somaAMPAf", w=0.02e-3, ISI=1, time_limit=simdur, sead=rdtmp) rdtmp=self.make_NetStims(po=self.olm, syn="somaAMPAf", w=0.0625e-3, ISI=1, time_limit=simdur, sead=rdtmp) rdtmp=self.make_NetStims(po=self.olm, syn="somaGABAf", w=0.2e-3, ISI=1, time_limit=simdur, sead=rdtmp) @@ -278,19 +278,19 @@ def make_all_NetStims(self,simdur,rdmseed): def make_all_noise(self,simdur,rdmseed): # create noise for simdur milliseconds numpy.random.seed(rdmseed) # initialize random # generator import math - print "Making Noise" + print("Making Noise") fctr = (simdur+simdur/2) / 10000.0 - print "to PYR" + print("to PYR") self.b_pyr_somaAMPAf=self.make_spikes(self.pyr,"somaAMPAf",0.05e-3,self.pyr.n,"soma",1,math.ceil(10000*fctr),1,simdur) self.b_pyr_Adend3AMPAf=self.make_spikes(self.pyr,"Adend3AMPAf",0.05e-3,self.pyr.n,"Adend3",1,math.ceil(10000*fctr),1,simdur) self.b_pyr_somaGABAf=self.make_spikes(self.pyr,"somaGABAf",0.012e-3,self.pyr.n,"soma",1,math.ceil(10000*fctr),1,simdur) self.b_pyr_Adend3GABAf=self.make_spikes(self.pyr,"Adend3GABAf",0.012e-3,self.pyr.n,"Adend3",1,math.ceil(10000*fctr),1,simdur) self.b_pyr_Adend3NMDA=self.make_spikes(self.pyr,"Adend3NMDA",6.5e-3,self.pyr.n,"Adend3",100,math.ceil(100*fctr),1,simdur) - print "to BAS" + print("to BAS") self.b_bas_somaAMPAf=self.make_spikes(self.bas,"somaAMPAf",0.02e-3,self.bas.n,"soma",1,math.ceil(10000*fctr),1,simdur) self.b_bas_somaGABA=self.make_spikes(self.bas,"somaGABAf",0.2e-3,self.bas.n,"soma",1,math.ceil(10000*fctr),1,simdur) self.b_bas_somaGABAf=self.make_spikes(self.bas,"somaGABAss",1.6e-3,self.bas.n,"soma",150,math.ceil(65*fctr),0,simdur) - print "to OLM" + print("to OLM") self.b_olm_somaAMPAf=self.make_spikes(self.olm,"somaAMPAf",0.02e-3,self.olm.n,"soma",1,math.ceil(10000*fctr),1,simdur) self.b_olm_somaGABAf=self.make_spikes(self.olm,"somaGABAf",0.2e-3,self.olm.n,"soma",1,math.ceil(10000*fctr),1,simdur) self.b_olm_somaGABAss=self.make_spikes(self.olm,"somaGABAss",1.6e-3,self.olm.n,"soma",150,math.ceil(65*fctr),0,simdur) @@ -303,17 +303,17 @@ def make_conn(self, preN, postN, conv): def set_all_conns(self): random.seed(self.wseed) # initialize random # generator for wiring - print "PYR -> X , NMDA" # src, trg, syn, delay, weight, conv + print("PYR -> X , NMDA") # src, trg, syn, delay, weight, conv self.pyr_bas_NM=self.set_connections(self.pyr,self.bas, "somaNMDA", 2, 1* 1.15*1.2e-3, 100) self.pyr_olm_NM=self.set_connections(self.pyr,self.olm, "somaNMDA", 2, 1* 1.0*0.7e-3, 10) self.pyr_pyr_NM=self.set_connections(self.pyr,self.pyr, "BdendNMDA",2, 1* 1*0.004e-3, 25) - print "PYR -> X , AMPA" + print("PYR -> X , AMPA") self.pyr_bas_AM=self.set_connections(self.pyr,self.bas, "somaAMPAf",2, 1* 0.3*1.2e-3, 100) self.pyr_olm_AM=self.set_connections(self.pyr,self.olm, "somaAMPAf",2, 1* 0.3*1.2e-3, 10) self.pyr_pyr_AM=self.set_connections(self.pyr,self.pyr, "BdendAMPA",2, 1* 0.5*0.04e-3, 25) - print "BAS -> X , GABA" + print("BAS -> X , GABA") #self.bas_bas_GA=self.set_connections(self.bas,self.bas, "somaGABAf",2, 1.0e-3, 60)#orig 1 #self.bas_bas_GA=self.set_connections(self.bas,self.bas, "somaGABAf",2, 2 * 1.5*1.0e-3, 60)#new 1 self.bas_bas_GA=self.set_connections(self.bas,self.bas, "somaGABAf",2, 1* 3*1.5*1.0e-3, 60)#new 2 @@ -321,14 +321,14 @@ def set_all_conns(self): self.bas_olm_GA=self.set_connections(self.bas,self.olm, "somaGABAf",2, 1* 0.05*2*2*0.18e-3, 17) #addition MSJ - print "OLM -> PYR , GABA" + print("OLM -> PYR , GABA") #self.olm_pyr_GA=self.set_connections(self.olm,self.pyr, "Adend2GABAs",2, 3*6.0e-3, 20)#original weight value self.olm_pyr_GA=self.set_connections(self.olm,self.pyr, "Adend2GABAs",2, 1* 4.0*3*6.0e-3, 20)#new weight value - self.olm_pyr_GA=self.set_connections(self.olm,self.pyr, "Adend2GABAs",2, 1* 0.08* 4.0*3*6.0e-3, 10)#new weight value + self.olm_pyr_GA=self.set_connections(self.olm,self.pyr, "Adend2GABAs",2, 1* 0.08* 4.0*3*6.0e-3, 10)#new weight value #pyramidal to PSR cell -- for testing only - print "PYR -> PSR, AMPA/NMDA" + print("PYR -> PSR, AMPA/NMDA") self.pyr_psr_NM=self.set_connections(self.pyr,self.psr, "BdendNMDA",2, 1*0.004e-3, 25) self.pyr_psr_AM=self.set_connections(self.pyr,self.psr, "BdendAMPA",2, 0.5*0.04e-3, 25) @@ -345,7 +345,7 @@ def set_connections(self,src,trg,syn,delay,w,conv): nc.append(h.NetCon(src.cell[pre_id].soma(0.5)._ref_v, trg.cell[post_id].__dict__[syn].syn, 0, delay, w, sec=src.cell[pre_id].soma)) if self.SaveConn: try: - print self.nqcon.size() + print(self.nqcon.size()) except: self.nqcon = h.NQS("id1","id2","w","syn") self.nqcon.strdec("syn") @@ -358,25 +358,25 @@ def set_connections(self,src,trg,syn,delay,w,conv): def load_spikes(self,fn,po,syn,w,time_limit=10000): fn = os.path.join("data",fn) events = numpy.load(fn) - print "Begin setting events...", po - print events.shape + print("Begin setting events..."), po + print(events.shape) for i,ii in enumerate(events): ii=ii[ii<=time_limit] po.cell[i].__dict__[syn].append(ii) po.cell[i].__dict__[syn].syn.Vwt = w - print "End setting events" + print("End setting events") return events def make_spikes(self,po,syn,w,cellN,comp,ISI,eventN,noise,time_limit): events = numpy.random.exponential(ISI, (cellN,eventN))*noise+numpy.repeat(ISI,cellN*eventN).reshape((cellN,eventN))*(1-noise) events = numpy.cumsum(events,axis=1) - print "Begin setting events...", po - print events.shape + print("Begin setting events..."), po + print(events.shape) for i,ii in enumerate(events): ii=ii[ii<=time_limit] po.cell[i].__dict__[syn].append(ii) po.cell[i].__dict__[syn].syn.Vwt = w - print "End setting events" + print("End setting events") return events def rasterplot(self,sz=2): @@ -453,9 +453,9 @@ def pravgrates(self,skipms=200): self.fnq.select("ty",ty) vf = self.fnq.getcol("freq") if vf.size() > 1: - print "ty: ", ty, " avg rate = ", vf.mean(), "+/-", vf.stderr(), " Hz" + print("ty: "), ty, " avg rate = ", vf.mean(), "+/-", vf.stderr(), " Hz" else: - print "ty: ", ty, " avg rate = ", vf.mean(), "+/-", 0.0 , " Hz" + print("ty: "), ty, " avg rate = ", vf.mean(), "+/-", 0.0 , " Hz" ty += 1 def calc_lfp(self): # lfp is modeled as a difference between voltages in distal apical and basal compartemnts @@ -484,7 +484,7 @@ def calc_psd(self,fig=3): gr = [30, 80] # Gamma frequency range t0i = int(t0/h.dt) if t0i > len(self.lfp): - print "LFP is too short! (<200 ms)" + print("LFP is too short! (<200 ms)") return 0,0,0,0,0,0 pylab.figure(fig) @@ -525,10 +525,10 @@ def get_lim_max(self, data, ind): # return the position of the maximal element fp.close() #create the network net = Network(noise=True,connections=True,DoMakeNoise=True,iseed=ISEED,UseNetStim=True,wseed=WSEED,scale=1.0,MSGain=MSG) - print "set network from rseed.txt : iseed=",ISEED,", WSEED=",WSEED,", MSG = ",MSG + print("set network from rseed.txt : iseed="),ISEED,", WSEED=",WSEED,", MSG = ",MSG except: net = Network() - print "set network from default constructor" + print("set network from default constructor") #setup some variables in hoc def sethocix(): diff --git a/run.py b/run.py index cccedf5..062ac93 100644 --- a/run.py +++ b/run.py @@ -13,7 +13,7 @@ # handler for printing out time during simulation run def fi(): for i in range(0,int(h.tstop),100): - h.cvode.event(i, "print " + str(i)) + h.cvode.event(i, "print(" + str(i) + ")") fih = h.FInitializeHandler(1, fi) @@ -41,21 +41,21 @@ def myInitNetStims(): washoutT = 0 # washout time def dowashin(): - print "washIN at ", washinT, " = ", h.t , " ", olmWash[0], basWash[0], pyrWashB[0], pyrWashA[0] + print("washIN at "), washinT, " = ", h.t , " ", olmWash[0], basWash[0], pyrWashB[0], pyrWashA[0] net.olm.set_r("somaNMDA",olmWash[0]) net.bas.set_r("somaNMDA",basWash[0]) net.pyr.set_r("BdendNMDA",pyrWashB[0]) net.pyr.set_r("Adend3NMDA",pyrWashA[0]) def dowashout(): - print "washOUT at ", washoutT, " = " , h.t, " ", olmWash[1], basWash[1], pyrWashB[1], pyrWashA[1] + print("washOUT at "), washoutT, " = " , h.t, " ", olmWash[1], basWash[1], pyrWashB[1], pyrWashA[1] net.olm.set_r("somaNMDA",olmWash[1]) net.bas.set_r("somaNMDA",basWash[1]) net.pyr.set_r("BdendNMDA",pyrWashB[1]) net.pyr.set_r("Adend3NMDA",pyrWashA[1]) def setwash(): - print "washinT ", washinT, " washoutT ", washoutT + print("washinT "), washinT, " washoutT ", washoutT h.cvode.event(washinT,"nrnpython(\"dowashin()\")") h.cvode.event(washoutT,"nrnpython(\"dowashout()\")") @@ -124,9 +124,9 @@ def r(self, n): for i4, r4 in enumerate(self.pow.x): self.net.pyr.set_r("Adend3NMDA",r4) simstr = self.getsimstr(r1,r2,r3,r4) - print "NMDA/AMPA: " + simstr + print("NMDA/AMPA: ") + simstr h.run() - print "Time: ", h.startsw() - self.pow.timer + print("Time: "), h.startsw() - self.pow.timer self.pow.arch.reset_time_stamp() @@ -171,7 +171,7 @@ def plot_fun(self, fig, var, ylabel, mode=0): for i1, r1 in enumerate(self.pow.x): for i2, r2 in enumerate(self.pow.x): for i3, r3 in enumerate(self.pow.x): - print "[:,"+str(i1)+","+str(i2)+","+str(i3)+"]" + print("[:,"+str(i1)+","+str(i2)+","+str(i3)+"]") pylab.plot(self.pow.x, self.pow.__dict__[var][:,i1,i2,i3],label="[:,"+str(i1)+","+str(i2)+","+str(i3)+"]", picker=1) #pylab.label() else: @@ -221,8 +221,8 @@ def plot_fun(self, fig, var, ylabel, mode=0): self.pow.arch.save_fig(fig,ylabel) def onpick(event): - print "REWR" - print str(event.artist.get_label())+" ("+str(event.mouseevent.xdata)+","+str(event.mouseevent.ydata)+")" + print("REWR") + print(str(event.artist.get_label())+" ("+str(event.mouseevent.xdata)+","+str(event.mouseevent.ydata)+")") return True #save vec to fn (fn is path) @@ -233,7 +233,7 @@ def mysvvec(fn,vec): vec.vwrite(fp) fp.close() else: - print "savevec ERR: couldn't open " + fn + print("savevec ERR: couldn't open ") + fn #this class is for saving output, i.e. figures and py files to backup class Archive: @@ -267,21 +267,21 @@ def minrunsv(simstr,tstop=1200,dt=0.1): h.tstop=tstop h.dt=dt h.run() - print "saving output data" + print("saving output data") net.calc_lfp() fn = "./data/"+simstr+"_lfp.vec" mysvvec(fn,net.vlfp) net.setsnq() # make NQS with spike times fn = "./data/"+simstr+"_snq.nqs" net.snq.sv(fn) - print "making and saving output figures" + print("making and saving output figures") #read a Vector from file, fn is file-path, vec is a Vector def myrdvec(fn,vec): fp=h.File() fp.ropen(fn) if not fp.isopen(): - print "myrdvec ERRA: Couldn't open " + fn + print("myrdvec ERRA: Couldn't open ") + fn return False vec.vread(fp) fp.close() @@ -303,7 +303,7 @@ def loadminrundat(simstr): try: net.snq=h.NQS(fs) except: - print "loadminrundat ERRB: couldn't read snq from " + fs + print("loadminrundat ERRB: couldn't read snq from ") + fs net.snq.verbose=0 # next, copy snq into vectors so can plot with net.rasterplot for po in net.cells: for i in xrange(len(po.lidvec)): @@ -325,7 +325,7 @@ def testrun(): arch.save_fig(2,"tmp_psr_soma_volt") net.calc_psd(3) arch.save_fig(3,"tmp_fft") - print "\a" + print("\a") def batchrun(): bat = Batch(net) @@ -349,7 +349,7 @@ def testsame(ts,v1,v2): h.run() net.calc_lfp() v2.copy(net.vlfp) - print "same = " , v1.eq(v2) + print("same = ") , v1.eq(v2) #gethilbnqs - make two NQS objects out of LFP with phase/amplitude/filered signals in theta and gamma bands def gethilbnqs(vlfp,minth=3,maxth=12,ming=30,maxg=80,usemlab=True): @@ -391,9 +391,9 @@ def getampphnq(nqtheta,nqgamma,phbins=100,skipms=200): for i in xrange(startx,sz,1): bin=int(phbins*(nqtheta.v[colp][i]-phmin)/phrng) if bin<0: - print "bin < 0!" + print("bin < 0!") if bin>=phbins+1: - print "bin >= phbins+1" + print("bin >= phbins+1") lv.o(bin).append(nqgamma.v[cola][i]) for i in xrange(0,int(vamp.size()),1): sz = lv.o(i).size() @@ -424,9 +424,9 @@ def checkbase(endt=3e3,skipms=200,justone=False): j = 0 dt = h.dt for i in xrange(1,-1,-1): - print "set olm NMDA to ", float(i) + print("set olm NMDA to "), float(i) net.olm.set_r("somaNMDA",float(i)) - print "running for " , endt , " ms " + print("running for ") , endt , " ms " h.run() net.calc_lfp() vlfp.append(net.vlfp) diff --git a/stats.mod b/stats.mod index 31b0d63..76f1150 100644 --- a/stats.mod +++ b/stats.mod @@ -433,6 +433,7 @@ ENDVERBATIM : unable to get the drand here to recognize the same fseed used in rand FUNCTION vseed () { VERBATIM + #ifdef WIN32 double seed; if (ifarg(1)) seed=*getarg(1); else { @@ -443,12 +444,12 @@ FUNCTION vseed () { set_seed(seed); return seed; #else - struct timeval tp; - struct timezone tzp; + //struct timeval tp; + //struct timezone tzp; double seed; if (ifarg(1)) seed=*getarg(1); else { - gettimeofday(&tp,&tzp); - seed=tp.tv_usec; + //gettimeofday(&tp,&tzp); // tzp unknown storage size + seed=2; //Wtp.tv_usec; } srand((unsigned)seed); // srand48 set_seed(seed); diff --git a/vecst.mod b/vecst.mod index 5bc7839..5d4d2cf 100644 --- a/vecst.mod +++ b/vecst.mod @@ -76,6 +76,7 @@ PARAMETER { EQU=-1.3470e120 EQV=-1.3469e120 : value equal to same row value in parallel vector EQW=-1.3468e120 : value found in other vector + EQX=-1.3468e120 : dunno NEQ=-1.3467e120 SEQ=-1.3466e120 RXP=-1.3465e120 From 54818f5d56b19b9c1e0a322e0809501ca6285f0a Mon Sep 17 00:00:00 2001 From: Marianne Bezaire <6456664+mbezaire@users.noreply.github.com> Date: Sun, 7 Aug 2022 01:13:59 -0400 Subject: [PATCH 6/9] xrange to range --- network.py | 16 ++++++++-------- networkmsj.py | 14 +++++++------- run.py | 10 +++++----- 3 files changed, 20 insertions(+), 20 deletions(-) diff --git a/network.py b/network.py index 5c85620..4907b46 100644 --- a/network.py +++ b/network.py @@ -145,7 +145,7 @@ def load_all_noise(self): #load noise from data files #this should be called @ beginning of each sim - done in an FInitializeHandler in run.py def init_NetStims(self): # h.mcell_ran4_init(self.iseed) - for i in xrange(len(self.nrl)): + for i in range(len(self.nrl)): rds = self.nrl[i] sead = self.nrlsead[i] rds.MCellRan4(sead,sead) @@ -199,7 +199,7 @@ def RecPYRInputs(self): self.NCV[s] = [] sidx = self.pyr.ncsidx[s] eidx = self.pyr.nceidx[s] - for i in xrange(sidx,eidx+1): + for i in range(sidx,eidx+1): self.NCV[s].append(h.Vector()) self.ncl[i].record(self.NCV[s][-1]) @@ -217,7 +217,7 @@ def setnqin(self): sidx = self.pyr.ncsidx[s] eidx = self.pyr.nceidx[s] idx = 0 - for i in xrange(0,len(self.NCV[s])): + for i in range(0,len(self.NCV[s])): nqin.append(idx,jdx,self.NCV[s][i]) idx = idx + 1 jdx = jdx + 1 @@ -227,7 +227,7 @@ def mkspkh(self,binsz): snq=self.snq snq.verbose = 0 self.spkh = h.List() - for i in xrange(0,800): + for i in range(0,800): if snq.select("id",i) > 0: vt = snq.getcol("t") self.spkh.append(vt.histogram(0,h.tstop,binsz)) @@ -383,7 +383,7 @@ def rasterplot(self,sz=2): for po in self.cells: id = h.Vector() tv = h.Vector() - for i in xrange(po.n): + for i in range(po.n): id.append(po.lidvec[i]) tv.append(po.ltimevec[i]) id.mark(h.g[0],tv,"O",sz,col[pon],1) @@ -394,7 +394,7 @@ def setrastervecs(self): self.myidvec = h.Vector() #IDs and firing times for ALL cells self.mytimevec = h.Vector() for po in self.cells: - for i in xrange(po.n): + for i in range(po.n): self.myidvec.append(po.lidvec[i]) self.mytimevec.append(po.ltimevec[i]) @@ -408,7 +408,7 @@ def setsnq(self): ty = 0 vec = h.Vector() for po in self.cells: - for i in xrange(po.n): + for i in range(po.n): self.snq.v[0].append(po.lidvec[i]) self.snq.v[1].append(po.ltimevec[i]) vec.resize(po.lidvec[i].size()) @@ -431,7 +431,7 @@ def setfnq(self,skipms=200): tf = h.tstop - skipms ty = 0 for po in self.cells: - for i in xrange(po.n): + for i in range(po.n): id = po.cell[i].id n = float( self.snq.select("t",">",skipms,"id",id) ) self.fnq.append(id, n*1e3/tf, ty) diff --git a/networkmsj.py b/networkmsj.py index 0e06895..8e55fb5 100644 --- a/networkmsj.py +++ b/networkmsj.py @@ -199,7 +199,7 @@ def RecPYRInputs(self): self.NCV[s] = [] sidx = self.pyr.ncsidx[s] eidx = self.pyr.nceidx[s] - for i in xrange(sidx,eidx+1): + for i in range(sidx,eidx+1): self.NCV[s].append(h.Vector()) self.ncl[i].record(self.NCV[s][-1]) @@ -217,7 +217,7 @@ def setnqin(self): sidx = self.pyr.ncsidx[s] eidx = self.pyr.nceidx[s] idx = 0 - for i in xrange(0,len(self.NCV[s])): + for i in range(0,len(self.NCV[s])): nqin.append(idx,jdx,self.NCV[s][i]) idx = idx + 1 jdx = jdx + 1 @@ -227,7 +227,7 @@ def mkspkh(self,binsz): snq=self.snq snq.verbose = 0 self.spkh = h.List() - for i in xrange(0,800): + for i in range(0,800): if snq.select("id",i) > 0: vt = snq.getcol("t") self.spkh.append(vt.histogram(0,h.tstop,binsz)) @@ -387,7 +387,7 @@ def rasterplot(self,sz=2): for po in self.cells: id = h.Vector() tv = h.Vector() - for i in xrange(po.n): + for i in range(po.n): id.append(po.lidvec[i]) tv.append(po.ltimevec[i]) id.mark(h.g[0],tv,"O",sz,col[pon],1) @@ -398,7 +398,7 @@ def setrastervecs(self): self.myidvec = h.Vector() #IDs and firing times for ALL cells self.mytimevec = h.Vector() for po in self.cells: - for i in xrange(po.n): + for i in range(po.n): self.myidvec.append(po.lidvec[i]) self.mytimevec.append(po.ltimevec[i]) @@ -412,7 +412,7 @@ def setsnq(self): ty = 0 vec = h.Vector() for po in self.cells: - for i in xrange(po.n): + for i in range(po.n): self.snq.v[0].append(po.lidvec[i]) self.snq.v[1].append(po.ltimevec[i]) vec.resize(po.lidvec[i].size()) @@ -435,7 +435,7 @@ def setfnq(self,skipms=200): tf = h.tstop - skipms ty = 0 for po in self.cells: - for i in xrange(po.n): + for i in range(po.n): id = po.cell[i].id n = float( self.snq.select("t",">",skipms,"id",id) ) self.fnq.append(id, n*1e3/tf, ty) diff --git a/run.py b/run.py index 062ac93..5a2d0b0 100644 --- a/run.py +++ b/run.py @@ -306,7 +306,7 @@ def loadminrundat(simstr): print("loadminrundat ERRB: couldn't read snq from ") + fs net.snq.verbose=0 # next, copy snq into vectors so can plot with net.rasterplot for po in net.cells: - for i in xrange(len(po.lidvec)): + for i in range(len(po.lidvec)): id = po.cell[i].id po.lidvec[i].resize(0) po.ltimevec[i].resize(0) @@ -384,18 +384,18 @@ def getampphnq(nqtheta,nqgamma,phbins=100,skipms=200): vamp.fill(0) vn.fill(0) # init counts to 0 lv = h.List() # list to keep amplitude samples - for i in xrange(int(vph.size())): + for i in range(int(vph.size())): lv.append(h.Vector()) sz=int(nqgamma.v[0].size()) startx=int(skipms/h.dt) - for i in xrange(startx,sz,1): + for i in range(startx,sz,1): bin=int(phbins*(nqtheta.v[colp][i]-phmin)/phrng) if bin<0: print("bin < 0!") if bin>=phbins+1: print("bin >= phbins+1") lv.o(bin).append(nqgamma.v[cola][i]) - for i in xrange(0,int(vamp.size()),1): + for i in range(0,int(vamp.size()),1): sz = lv.o(i).size() if sz > 0: # if no samples, skip av = lv.o(i).mean() @@ -423,7 +423,7 @@ def checkbase(endt=3e3,skipms=200,justone=False): h.tstop=endt j = 0 dt = h.dt - for i in xrange(1,-1,-1): + for i in range(1,-1,-1): print("set olm NMDA to "), float(i) net.olm.set_r("somaNMDA",float(i)) print("running for ") , endt , " ms " From 49a0d38f07bc9581f919c6e12bf8f28e72c3ab87 Mon Sep 17 00:00:00 2001 From: Marianne Bezaire <6456664+mbezaire@users.noreply.github.com> Date: Sun, 7 Aug 2022 01:32:18 -0400 Subject: [PATCH 7/9] runs without error --- vecst.mod | 1 + 1 file changed, 1 insertion(+) diff --git a/vecst.mod b/vecst.mod index 5d4d2cf..c4ba8c1 100644 --- a/vecst.mod +++ b/vecst.mod @@ -77,6 +77,7 @@ PARAMETER { EQV=-1.3469e120 : value equal to same row value in parallel vector EQW=-1.3468e120 : value found in other vector EQX=-1.3468e120 : dunno + EQY=-1.3468e120 : dunno NEQ=-1.3467e120 SEQ=-1.3466e120 RXP=-1.3465e120 From 4a12daff3bb7a77d4fa7936679c7df49400ae04a Mon Sep 17 00:00:00 2001 From: Marianne Bezaire <6456664+mbezaire@users.noreply.github.com> Date: Sun, 7 Aug 2022 01:33:37 -0400 Subject: [PATCH 8/9] reinstated misc.mod for wsl usage --- misc.mignoreod => misc.mod | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename misc.mignoreod => misc.mod (100%) diff --git a/misc.mignoreod b/misc.mod similarity index 100% rename from misc.mignoreod rename to misc.mod From f7946e462e7172b0f3cbe7b3036d3efaa0136f1b Mon Sep 17 00:00:00 2001 From: Marianne Bezaire <6456664+mbezaire@users.noreply.github.com> Date: Mon, 8 Aug 2022 07:41:25 -0400 Subject: [PATCH 9/9] plot in spyder --- graphstuff.py | 19 ++++++++++++++++ mosinit.py | 22 +++++++++++++------ savedata.py | 60 +++++++++++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 94 insertions(+), 7 deletions(-) create mode 100644 graphstuff.py create mode 100644 savedata.py diff --git a/graphstuff.py b/graphstuff.py new file mode 100644 index 0000000..8a13482 --- /dev/null +++ b/graphstuff.py @@ -0,0 +1,19 @@ +# -*- coding: utf-8 -*- +""" +Created on Mon Aug 8 06:23:41 2022 + +@author: maria +""" +from neuron import h +import pickle +with open('netresults.pickle', 'rb') as f: + allvars = pickle.load(f) +data = allvars[0] +lfp = allvars[1] +dt = allvars[2] +import savedata as sd +sd.graphrasterplot(data) + +sd.graphlfp(lfp,dt) + + diff --git a/mosinit.py b/mosinit.py index a60848f..9248a3d 100644 --- a/mosinit.py +++ b/mosinit.py @@ -42,12 +42,20 @@ Run.fiwash = h.FInitializeHandler(1,Run.setwash) """ - h.tstop = 5e3 + h.tstop = 100 #5e3 h.run() - net.rasterplot() + net.calc_lfp() - net.pravgrates() - myg = h.Graph() - net.vlfp.plot(myg,h.dt) - myg.exec_menu("View = plot") - myg.exec_menu("New Axis") + import savedata as sd + import pickle + + data = sd.data4rasterplot(net) + with open('netresults.pickle','wb') as f: + pickle.dump([data, net.lfp, h.dt],f,pickle.HIGHEST_PROTOCOL) + + print("Back in Spyder/windows, run graphstuff.py") + #net.pravgrates() + # myg = h.Graph() + # net.vlfp.plot(myg,h.dt) + # myg.exec_menu("View = plot") + # myg.exec_menu("New Axis") diff --git a/savedata.py b/savedata.py new file mode 100644 index 0000000..919fe33 --- /dev/null +++ b/savedata.py @@ -0,0 +1,60 @@ +# -*- coding: utf-8 -*- +""" +Created on Mon Aug 8 06:00:21 2022 + +@author: maria +""" +from neuron import h + +def data4rasterplot(net,sz=2): + + listocells = [] + for po in net.cells: + id = h.Vector() + tv = h.Vector() + for i in range(po.n): + id.append(po.lidvec[i]) + tv.append(po.ltimevec[i]) + data = {"id":id.to_python(), "tv":tv.to_python()} + listocells.append(data) + return listocells + + +def graphrasterplot(rasterdata,sz=2): + import matplotlib.pyplot as plt + + fig = plt.figure() + pon = 0 + # if h.g[0] == None: + col = ['r','g','b','o','k'] + for po in rasterdata: + id = h.Vector() + tv = h.Vector() + id.from_python(po['id']) + tv.from_python(po['tv']) + if len(tv) > 0: + plt.plot(tv,id, linestyle="None",marker='o',markersize=sz, color=col[pon]) + pon += 1 + +def data4pravgrates(self,skipms=200): + try: + self.fnq.tog("DB") + except: + self.setfnq(skipms) + ty = 0 + tf = float( h.tstop - skipms ) + for po in self.cells: + self.fnq.select("ty",ty) + vf = self.fnq.getcol("freq") + if vf.size() > 1: + print("ty: "), ty, " avg rate = ", vf.mean(), "+/-", vf.stderr(), " Hz" + else: + print("ty: "), ty, " avg rate = ", vf.mean(), "+/-", 0.0 , " Hz" + ty += 1 + +def graphlfp(lfp,dt): + import matplotlib.pyplot as plt + import numpy as np + tvec = np.arange(0,dt*len(lfp),dt) + figure = plt.figure() + plt.plot(tvec,lfp) \ No newline at end of file