: \$Id: stats.mod,v 1.8 2006/06/30 19:03:36 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(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 { SUFFIX nothing GLOBAL STATS_INSTALLED } PARAMETER { : BVBASE = -1. : defined in vecst.mod STATS_INSTALLED=0 } ASSIGNED { } VERBATIM #include #include #include /* contains MAXLONG */ #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); extern int list_vector_px2(); int list_vector_px(); int list_vector_resize(); 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 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=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.hamming(v2[,v3]) compares v1 and v2 for matches, v3 gives diff vector VERBATIM 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 sum; } ENDVERBATIM :* v1.flipbits(scratch,num) flips num bits : uses scratch vector of same size as v1 to make sure doesn't flip same bit twice VERBATIM static double flipbits(void* vv) { int i, nx, ny, flip, ii; double* x, *y; nx = vector_instance_px(vv, &x); ny = vector_arg_px(1, &y); flip = (int)*getarg(2); if (nx != ny) { hoc_execerror("Scratch vector must be same size", 0); } 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) place counts for each interval VERBATIM static double bin(void* vv) { int i, j, nx, ny, next, maxsz; double* x, *y, invl, loc; nx = vector_instance_px(vv, &x); 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; jloc) { loc+=invl; j++; } y[j]++; } } vector_resize(vv, j); return (double)j; } ENDVERBATIM :* PROCEDURE install_stats() PROCEDURE install_stats () { STATS_INSTALLED=1 VERBATIM install_vector_method("slope", slope); install_vector_method("vslope", vslope); install_vector_method("stats", stats); 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); 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 }