: \$Id: stats.mod,v 1.74 2007/12/01 00:24:14 billl 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 } PARAMETER { : BVBASE = 0. : defined in vecst.mod INSTALLED=0 kmeasure=0 verbose=0 } ASSIGNED { seed } VERBATIM #include #include #include /* contains LONG_MAX */ #include #include extern double BVBASE; extern double* hoc_pgetarg(); Symbol *hoc_get_symbol(); char ** hoc_pgargstr(); extern Point_process* ob2pntproc(Object*); 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 double *vector_newsize(); extern int vector_instance_px(); extern void* vector_arg(); extern double* vector_vec(); extern double hoc_epsilon; extern void mcell_ran4_init(unsigned int *idum); extern double mcell_ran4(unsigned int* idum,double* ran_vec,unsigned int n,double range); extern void set_seed(); extern unsigned int *scrset(); extern int ivoc_list_count(Object*); extern Object* ivoc_list_item(Object*, int); extern int list_vector_px2(); extern int hoc_is_double_arg(int narg); extern Objectdata *hoc_objectdata; extern int openvec(int, double **); extern char* hoc_object_name(Object*); extern int nrn_mlh_gsort(); extern int cmpdfn(); int list_vector_px(); int list_vector_resize(); static void hxe() { hoc_execerror("",0); } unsigned int valseed; typedef struct BVEC { int size; int bufsize; short *x; Object* o; } bvec; union dblint { int i[2]; double d; }; #define VRRY 50 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; 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=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(voi[0], maxsz); if (lfl) vector_resize(voi[1], maxsz); return (double)maxsz; } ENDVERBATIM :* PROCEDURE install_stats() PROCEDURE install () { if (INSTALLED==1) { printf("\$Id: stats.mod,v 1.74 2007/12/01 00:24:14 billl Exp \$") } else { INSTALLED=1 VERBATIM valseed=1; install_vector_method("slope", slope); install_vector_method("vslope", vslope); install_vector_method("stats", stats); install_vector_method("pcorrel", pcorrel); 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("distance", distance); install_vector_method("ndprd", ndprd); install_vector_method("hash", hash); install_vector_method("unnan", unnan); 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 (ntop 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 #if defined(t) return *(hoc_objectdata[sym->u.oboff]._pval); #else return *(hoc_objectdata[sym->u.oboff].pval); #endif } ENDVERBATIM }