: \$Id: stats.mod,v 1.151 2008/12/30 03:56:37 samn Exp \$ :* COMMENT COMMENT randwd randomly chooses n bits to set to 1 hamming v.hamming(v1) is hamming distance between 2 vecs flipbits v.flipbits(scratch,num) flips num rand chosen bits flipbalbits v.flipbalbits(scratch,num) balanced flipping vpr v.vpr prints out vector as 1 (x[i]>0) or 0 (x[i]<=0) 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 ENDCOMMENT NEURON { SUFFIX stats GLOBAL INSTALLED,seed,kmeasure,verbose,self_ok_combi,hretval } PARAMETER { : BVBASE = 0. : defined in vecst.mod INSTALLED=0 kmeasure=0 verbose=0 self_ok_combi=0 hretval=0 } ASSIGNED { seed } VERBATIM #include #include #include "misc.h" #if !defined(t) #define _pval pval #endif unsigned int valseed; static double *x1x, *y1y, *z1z; union dblint { int i[2]; double d; }; #define VRRY 50 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; } ENDVERBATIM :* v1.slope(num) does a linear regression to find the slope, assuming num=timestep of vector VERBATIM static double slope(void* vv) { int i, n; double *x, *y; double timestep, sigxy, sigx, sigy, sigx2; /* how to get the instance data */ n = vector_instance_px(vv, &y); if(ifarg(1)) { timestep = *getarg(1); } else { printf("You must supply a timestep\n"); return 0; } sigxy= sigx= sigy= sigx2=0; // initialize these x = (double *) malloc(sizeof(double)*n); for(i=0; i*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 timestep, sigxy, sigx, sigy, sigx2, sigy2; double r, m, b, dmeansqerr,dmaxsqerr; /* how to get the instance data */ n = vector_instance_px(vv, &y); if(ifarg(1)) { timestep = *getarg(1); } else { printf("You must supply a timestep\n"); return 0; } sigxy= sigx= sigy= sigx2=sigy2= 0; // initialize these x = (double *) malloc(sizeof(double)*n); for(i=0; iX[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*/ 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; // 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(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("scorrelsERRA: %d %d\n",n,i); hxe();} return spearman(n,x,y); } ENDVERBATIM :* vec.vnan() will reset nans and infs to selected values -- default 0,0 VERBATIM static double unnan(void *vv) { int i,nx,cnt; double newnan,newinf; union dblint xx; double *x; newnan=newinf=0; nx = vector_instance_px(vv, &x); if (ifarg(1)) newinf=newnan=*getarg(1); if (ifarg(2)) newinf=*getarg(2); for (i=0,cnt=0;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()); if (x[ii]==BVBASE) { x[ii] = 1.; if (flag) { y[jj] = ii; jj++; } } else { i--; } } return flip; } ENDVERBATIM :* v1.hash(veclist) VERBATIM static double hash (void* vv) { int i, j, nx, nv[VRRY], num, vfl; union dblint xx; Object* ob; double *x, *vvo[VRRY], big, y, prod; nx = vector_instance_px(vv, &x); if (ifarg(1)) { vfl=0; ob=*hoc_objgetarg(1); num = ivoc_list_count(ob); if (num>VRRY) {printf("vecst:hash ERR: can only handle %d vecs: %d\n",VRRY,num); hxe();} for (i=0;i= randvar and mul by step : v1.setrand(4.5,min,max[,seed]) // [min,max) : v1.setrand(5[,n,seed]) -- integers [0,100) or [0,n] : v1.setrand(5[,min,seed,max]) -- integers [min,max] -- if seed=0 it's not reset : v1.setrand(5,ind[,seed]) -- random values from ind : v1.setrand(6) -- unique integers as follows: : v1.setrand(6,min,max[,seed]) -- unique in [min,max] : v1.setrand(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; void* vhso; //vector * for changing size at end Symbol* pHocVecFunc,*pHocCompFunc; //hoc function pointers dret=-1.0; g0t=g1t=0x0; pm=pids=0x0; //init arrays to null pcombids=0x0; 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]; } } } } 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(vv,nx*augfac); for (i=1;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; } return (double)nx; } ENDVERBATIM :* v1.distance(v2) euclidean distance VERBATIM static double distance (void* vv) { int i, nx, ny; double* x, *y, sum; sum = 0.; nx = vector_instance_px(vv, &x); ny = vector_arg_px(1, &y); if (nx!=ny) {printf("Vectors must be same size %d %d\n",nx,ny); hxe();} for (i=0; i=10?"|":"",(int)x[i],x[i]>=10?"|":""); } else if (flag==-1) { for (i=0; i=10) printf("+"); else printf("%d",(int)x[i]); } else if (flag==16) { // hex for (i=0; i=16) printf("+"); else printf("%x",(int)x[i]); } else if (flag==64) { // base 64 for (i=0; i=64) { printf("+"); } else if (x[i]<16) { printf("%x",(int)x[i]); // 0-f 0-15 } else if (x[i]<36) { printf("%c",(int)x[i]+87); // g-z 16-35 } else if (x[i]<62) { printf("%c",(int)x[i]+29); // A-Z 36-61 } else if (x[i]<63) { printf("@"); // @ 62 } else if (x[i]<64) { printf("="); // = 63 } else printf("ERROR"); } } if (!ifarg(2)) printf("\n"); else printf(" "); } else { 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=0; iBVBASE) { printf("%d",1); } else { printf("%d",0); } } printf("\n"); } 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 VERBATIM static double bin (void* vv) { int i, j, nx, maxsz, lfl; double* x, *y, *ix, invl, min, max, maxf, jj; Object* ob; void* voi[2]; 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); 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();} y[j]++; } if (lfl) ix[j]=jj+min; } maxsz=(max==1e9)?(int)(maxf/invl+1):(int)((max-min)/invl+1); vector_resize((IvocVect*)voi[0], maxsz); if (lfl) vector_resize((IvocVect*)voi[1], maxsz); return (double)maxsz; } ENDVERBATIM :* PROCEDURE install_stats() PROCEDURE install () { if (INSTALLED==1) { printf("\$Id: stats.mod,v 1.151 2008/12/30 03:56:37 samn Exp \$") } else { 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("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("bin", bin); install_vector_method("setrnd", setrnd); install_vector_method("rantran", rantran); install_vector_method("distance", distance); install_vector_method("ndprd", ndprd); install_vector_method("hash", hash); 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); ENDVERBATIM } } PROCEDURE prhash (x) { VERBATIM { union dblint xx; xx.d=_lx; printf("%08x%08x\n",xx.i[0],xx.i[1]); } ENDVERBATIM } :* fac (n) : from numerical recipes p.214 FUNCTION fac (n) { VERBATIM { static int ntop=4; static double a[101]={1.,1.,2.,6.,24.}; static double cof[6]={76.18009173,-86.50532033,24.01409822, -1.231739516,0.120858003e-2,-0.536382e-5}; int j,n; n = (int)_ln; if (n<0) { hoc_execerror("No negative numbers ", 0); } if (n>100) { /* gamma function */ double x,tmp,ser; x = _ln; tmp=x+5.5; tmp -= (x+0.5)*log(tmp); ser=1.0; for (j=0;j<=5;j++) { x += 1.0; ser += cof[j]/x; } return exp(-tmp+log(2.50662827465*ser)); } else { while (ntop100) { /* gamma function */ double x,tmp,ser; x = _ln; tmp=x+5.5; tmp -= (x+0.5)*log(tmp); ser=1.0; for (j=0;j<=5;j++) { x += 1.0; ser += cof[j]/x; } return (-tmp+log(2.50662827465*ser)); } else { while (ntopUINT_MAX) y/=1e9; // UINT_MAX is 4.294967e+09 valseed*=(unsigned int)y; // keep multiplying these out } return valseed; // global } ENDVERBATIM FUNCTION hashseed () { VERBATIM int i,na,nb,sf=0; double *x=0; if (hoc_is_double_arg(2)) { for (na=1;ifarg(na);na++); na--; if (na==1) nb=2; else nb=na; x = (double *) malloc(sizeof(double)*nb); for (i=1;i<=na;i++) x[i-1]=*getarg(i); if (na==1) x[1]=x[0]*x[0]+13; // so can get something sf=1;//should free } else { nb=(int)*getarg(1); x=hoc_pgetarg(2); } hashseed2(nb,x); if(sf) free(x); //only free if malloc'd _lhashseed=(double)valseed; ENDVERBATIM } : unable to get the drand here to recognize the same fseed used in rand FUNCTION vseed () { VERBATIM #ifdef WIN32 if (ifarg(1)) seed=*getarg(1); else { printf("TIME ACCESS NOT PRESENT IN WINDOWS\n"); hxe(); } srand48((unsigned)seed); set_seed(seed); return seed; #else struct timeval tp; struct timezone tzp; if (ifarg(1)) seed=*getarg(1); else { gettimeofday(&tp,&tzp); seed=tp.tv_usec; } srand48((unsigned)seed); 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 { double x,tmp,ser; static double cof[6]={76.18009173,-86.50532033,24.01409822,-1.231739516,0.120858003e-2,-0.536382e-5}; int j; x=_lxx-1.0; tmp=x+5.5; tmp -= (x+0.5)*log(tmp); ser=1.0; for (j=0;j<=5;j++) { x += 1.0; ser += cof[j]/x; } return -tmp+log(2.50662827465*ser); } ENDVERBATIM } FUNCTION betai(a,b,x) { VERBATIM { double bt; double gammln(double),betacf(double,double,double); 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(_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(_la,_lb,_lx)/_la; else return 1.0-bt*betacf(_lb,_la,1.0-_lx)/_lb; } ENDVERBATIM } VERBATIM #define ITMAX 100 #define EPS 3.0e-7 ENDVERBATIM FUNCTION betacf(a,b,x) { VERBATIM { double qap,qam,qab,em,tem,d; double bz,bm=1.0,bp,bpp; double az=1.0,am=1.0,ap,app,aold; int m; void nrerror(); qab=_la+_lb; qap=_la+1.0; qam=_la-1.0; bz=1.0-qab*_lx/qap; for (m=1;m<=ITMAX;m++) { em=(double) m; tem=em+em; d=em*(_lb-em)*_lx/((qam+tem)*(_la+tem)); ap=az+d*am; bp=bz+d*bm; d = -(_la+em)*(qab+em)*_lx/((qap+tem)*(_la+tem)); app=ap+d*az; bpp=bp+d*bz; aold=az; am=ap/bpp; bm=bp/bpp; az=app/bpp; bz=1.0; if (fabs(az-aold) < (EPS*fabs(az))) return az; } printf("a or b too big, or ITMAX too small in BETACF"); return -1.; } 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 gammln(double); double x = *getarg(1); double dof = *getarg(2); double res = (gammln( (dof+1.0) / 2.0 ) / gammln( 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 }