: \$Id: vecst.mod,v 1.178 2005/03/07 22:26:47 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) thresh turns analog vec to BVBASE,1 vec separating at thresh (scalar or vec) 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 xing determines where vector crosses in a positive direction (destructive) xzero count 0 crossings by putting in a 1 or -1 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 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 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 vdest.intrp(flag) // interpolate numbers replacing numbers given as flag v.insct(v1,v2) // return v1 intersect v2 vdest.cull(vsrc,key) // remove values found in key vdest.redundout(vsrc[,INDFLAG]) // remove repeat values, with flag return indices v.cvlv(v1,v2) // convolve v1 with v2 fac not vec related - returns factorial logfac not vec related - returns log factorial v2d copy into a double block d2v copy out of a double block NB: same as vec.from_double() smgs rewrite of sumgauss form ivovect.cpp iwr write integers ird read integers info give pointer addresses and sizes lcat concatentate all vectors from a list Non-vector routines vseed set some C level randomizer seeds isojt compare 2 objects to see if they're of the same type eqojt compare 2 object pointers to see if they point to same thing ENDCOMMENT NEURON { SUFFIX nothing GLOBAL BVBASE, RES, VECST_INSTALLED : bit vector base number (typically 0 or -1) } PARAMETER { BVBASE = -1. VECST_INSTALLED=0 : misc ERR=-1.3479e121 : 0 args ALL=-1.3479e120 NEG=-1.3478e120 POS=-1.3477e120 CHK=-1.3476e120 NOZ=-1.3475e120 : 1 arg GTH=-1.3474e120 GTE=-1.3473e120 LTH=-1.3472e120 LTE=-1.3471e120 EQU=-1.3470e120 EQV=-1.3469e120 : value equal to same row value in parallel vector EQW=-1.3468e120 : value found in other vector NEQ=-1.3467e120 SEQ=-1.3466e120 RXP=-1.3465e120 : 2 args IBE=-1.3464e120 EBI=-1.3463e120 IBI=-1.3462e120 EBE=-1.3461e120 } ASSIGNED { RES } VERBATIM #include #include #include /* contains MAXLONG */ #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 void set_seed(); extern int ivoc_list_count(Object*); extern Object* ivoc_list_item(Object*, int); static int list_vector_px(); static int list_vector_px2(); static int list_vector_resize(); 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.info() gives addresses and sizes VERBATIM static double info(void* vv) { int nx,bsz; double* x; nx = vector_instance_px(vv, &x); bsz=vector_buffer_size(vv); printf("Obj*%x Dbl*%x Size: %d Bufsize: %d\n",vv,x,nx,bsz); } ENDVERBATIM :* v1.indset(ind,x[,y]) sets indexed values to x and other values to optional y VERBATIM static double indset(void* vv) { int i, nx, ny, nz, flag; double* x, *y, *z, val, val2 ; nx = vector_instance_px(vv, &x); ny = vector_arg_px(1, &y); if (hoc_is_object_arg(2)) { flag=1; nz = vector_arg_px(2, &z); if (ny!=nz) { hoc_execerror("v.indset: Vector sizes don't match.", 0); } } else { flag=0; val=*getarg(2); } if (ifarg(3)) { val2 = *getarg(3); for (i=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; } return i; } ENDVERBATIM VERBATIM /* Maintain parallel int vector to avoid slowness of repeated casts */ static int scrsz=0; static int *scr; static double dcr[100]; // scratch area for doubles ENDVERBATIM :* tmp.fewind(ind,veclist) : picks out numbers from multiple vectors using index ind VERBATIM static double fewind (void* vv) { int i, j, nx, ni, nv[10], num; Object* ob; double *x, *ind, *vvo[10]; nx = vector_instance_px(vv, &x); ni = vector_arg_px(1, &ind); ob = *hoc_objgetarg(2); num = ivoc_list_count(ob); if (num>10) hoc_execerror("ERR: fewind can only handle 10 vectors", 0); for (i=0;iscrsz) { if (scrsz>0) { free(scr); scr=(int *)NULL; } scrsz=ni+10000; scr=(int *)ecalloc(scrsz, sizeof(int)); } 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); } } for (j=0;j11) hoc_execerror("findx ****ERRB****: can only handle 11 vectors", 0); for (i=0;iscrsz) { if (scrsz>0) { free(scr); scr=(int *)NULL; } scrsz=ni+10000; scr=(int *)ecalloc(scrsz, sizeof(int)); } for (i=0;i=nx || scr[i]<0) { printf("findx ****ERRE**** **** IND:%d SZ:%d\n",scr[i],nx); hoc_execerror("Index vector out-of-bounds", 0); } } for (j=0,i=0;j11) hoc_execerror("sindx ****ERRB****: can only handle 11 vectors", 0); for (i=0;iscrsz) { if (scrsz>0) { free(scr); scr=(int *)NULL; } scrsz=ni+10000; scr=(int *)ecalloc(scrsz, sizeof(int)); } for (i=0;i=nx || scr[i]<0) { printf("sindx ****ERRE**** IND:%d SZ:%d\n",scr[i],nx); hoc_execerror("Index vector out-of-bounds", 0); } } for (j=0,i=0;j11) hoc_execerror("sindv ****ERRA****: can only handle 11 vectors", 0); for (i=0;iscrsz) { if (scrsz>0) { free(scr); scr=(int *)NULL; } scrsz=ni+10000; scr=(int *)ecalloc(scrsz, sizeof(int)); } for (i=0;i=nx || scr[i]<0) { printf("sindv ****ERRE**** IND:%d SZ:%d\n",scr[i],nx); hoc_execerror("Index vector out-of-bounds", 0); } } for (j=0,i=0;j10) hoc_execerror("ERR: vecst::slct can only handle 10 vectors", 0); for (i=0,j=0,k=0;i=0.) {fl=0; break;} else continue; } else if (key[n]==GTH) { if (vvo[i][j]<=arg[m]) {fl=0; break;} else continue; } else if (key[n]==GTE) { if (vvo[i][j]< arg[m]) {fl=0; break;} else continue; } else if (key[n]==LTH) { if (vvo[i][j]>=arg[m]) {fl=0; break;} else continue; } else if (key[n]==LTE) { if (vvo[i][j]> arg[m]) {fl=0; break;} else continue; } else if (key[n]==EQU) { if (vvo[i][j]!=arg[m]) {fl=0; break;} else continue; } else if (key[n]==EQV) { if (vvo[i][j]!=vvo[i+1][j]) { fl=0; break;} else { i++; continue; } } else if (key[n]==EQW) { // check value against values in following vec fl=0; // assume it's not going to match for (p=0;p=arg[m+1])) { fl=0; break; } else continue; // IBE="[)" include-bracket-exclude } else if (key[n]==EBI) { if ((vvo[i][j]<=arg[m])||(vvo[i][j]> arg[m+1])) { fl=0; break; } else continue; // "(]" : exclude-bracket-include } else if (key[n]==IBI) { if ((vvo[i][j]< arg[m])||(vvo[i][j]> arg[m+1])) { fl=0; break; } else continue; // "[]" : include-bracket-include } else if (key[n]==EBE) { if ((vvo[i][j]<=arg[m])||(vvo[i][j]>=arg[m+1])) { 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 } vector_resize(vv, k); return k; } 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[10], num, fl; Object* lob; double *ind, *key, *arg, *vvo[10]; ni = vector_instance_px(vv, &ind); // vv is ind nk = vector_arg_px(1, &key); na = vector_arg_px(2, &arg); lob = *hoc_objgetarg(3); num = ivoc_list_count(lob); if (num>10) hoc_execerror("ERR: vecst::slor can only handle 10 vectors", 0); for (i=0,j=0,k=0;i=0.) continue; else {fl=1; break;} } else if (key[n]==GTH) { if (vvo[i][j]<=arg[m]) continue; else {fl=1; break;} } else if (key[n]==GTE) { if (vvo[i][j]< arg[m]) continue; else {fl=1; break;} } else if (key[n]==LTH) { if (vvo[i][j]>=arg[m]) continue; else {fl=1; break;} } else if (key[n]==LTE) { if (vvo[i][j]> arg[m]) continue; else {fl=1; break;} } else if (key[n]==EQU) { if (vvo[i][j]!=arg[m]) continue; else {fl=1; break;} } else if (key[n]==EQV) { if (vvo[i][j]!=vvo[i+1][j]) continue; else { i++; fl=1; break; } } else if (key[n]==EQW) { // check value against values in following vec fl=0; // assume it's not going to match for (p=0;p=arg[m+1])) { continue; } else {fl=1; break;} // IBE="[)" include-bracket-exclude } else if (key[n]==EBI) { if ((vvo[i][j]<=arg[m])||(vvo[i][j]> arg[m+1])) { continue; } else {fl=1; break;} // "(]" : exclude-bracket-include } else if (key[n]==IBI) { if ((vvo[i][j]< arg[m])||(vvo[i][j]> arg[m+1])) { continue; } else {fl=1; break;} // "[]" : include-bracket-include } else if (key[n]==EBE) { if ((vvo[i][j]<=arg[m])||(vvo[i][j]>=arg[m+1])) { continue; } else {fl=1; break;} // "()" : exclude-bracket-exclude } else {printf("vecst::slor ERR4 %g\n",key[n]); hoc_execerror("Unknown key",0);} } if (fl) ind[k++]=j; // all equal } vector_resize(vv, k); return k; } ENDVERBATIM :* v.iwr() VERBATIM static double iwr(void* vv) { int i, j, nx; double *x; FILE* f, *hoc_obj_file_arg(); f = hoc_obj_file_arg(1); nx = vector_instance_px(vv, &x); if (nx>scrsz) { if (scrsz>0) { free(scr); scr=(int *)NULL; } scrsz=nx+10000; scr=(int *)ecalloc(scrsz, sizeof(int)); } for (i=0;iscrsz) { if (scrsz>0) { free(scr); scr=(int *)NULL; } scrsz=n+10000; scr=(int *)ecalloc(scrsz, sizeof(int)); } if (n!=nx) { nx=vector_buffer_size(vv); if (n<=nx) { vector_resize(vv, n); nx=n; } else { printf("%d > %d :: ",n,nx); hoc_execerror("Vector max capacity too small for ird ", 0); } } fread(scr,sizeof(int),n,f); for (i=0;imaxsz) { printf("\tinsct WARNING: ran out of room: %d<%d\n",maxsz,k); } else { vector_resize(vv, k); } return (double)k; } ENDVERBATIM :* vec.cull(src,key) : remove numbers in vec that are found in the key VERBATIM static double cull(void* vv) { int i, j, k, nx, nv1, nv2, maxsz, flag; double *x, *v1, *v2; nx = vector_instance_px(vv, &x); 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("\tcull WARNING: ran out of room: %d<%d\n",maxsz,k); } else { vector_resize(vv, k); } return (double)k; } ENDVERBATIM :* dest.redundout(src[,INDFLAG]) : flag redundant numbers; must sort src first VERBATIM static double redundout(void* vv) { int i, j, nx, nv1, maxsz, indflag; double *x, *v1, val; if (ifarg(2)) indflag=1; else indflag=0; nx = vector_instance_px(vv, &x); maxsz=vector_buffer_size(vv); vector_resize(vv, maxsz); nv1 = vector_arg_px(1, &v1); val=v1[0]; x[0]=(indflag?0:val); for (j=1,i=1;insrc) { hoc_execerror("Filter bigger than source ", 0); } for (i=0;i0 && k9) hoc_execerror("ERR: nind can only handle 9 vectors", 0); num = i-2; /* number of vectors to be picked apart */ for (i=0;iscrsz) { if (scrsz>0) { free(scr); scr=(int *)NULL; } scrsz=ni+10000; scr=(int *)ecalloc(scrsz, sizeof(int)); } for (i=0,last=-1;i=nx) hoc_execerror("nind(): Index out of bounds", 0); if (scr[i]<=last) hoc_execerror("nind(): indices should mono increase", 0); last=scr[i]; } for (j=0;j10) hoc_execerror("ERR: keyind can only handle 9 vectors", 0); num = i-1; /* number of vectors to be picked apart */ for (i=0;i=0.) break; else continue; } else if (key[i]!=vvo[i][j]) break; // default } if (i==nk) ind[k++]=j; // all equal } vector_resize(vv, k); return k; } 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.thresh() threshold above and below thresh VERBATIM static double thresh(void* vv) { int i, nx, ny, cnt; double *x, *y, th; nx = vector_instance_px(vv, &x); cnt=0; if (hoc_is_object_arg(1)) { ny = vector_arg_px(1, &y); th=0; if (nx!=ny) { hoc_execerror("Vector sizes don't match in thresh.", 0); } for (i=0; i=y[i]) { x[i]= 1.; cnt++;} else { x[i]= BVBASE; } } } else { th = *getarg(1); for (i=0; i= th) { x[i]= 1.; cnt++;} else { x[i]= BVBASE; } } } return cnt; } ENDVERBATIM :* v1.triplet() return location of a triplet VERBATIM static double triplet(void* vv) { int i, nx; double *x, *y, a, b; nx = vector_instance_px(vv, &x); a = *getarg(1); b = *getarg(2); for (i=0; i0.) { st[i]=1.; continue; } /* cell must fire */ if (vol[i]>=thr[i] && obon[i]<= -refr[i]) { /* past refractory period */ st[i]=1.; obon[i]=dur[i]; num++; } else { st[i]= BVBASE; } } return num; } ENDVERBATIM :* vo.bpeval(outp,del) : service routine for back-prop: vo=outp*(1-outp)*del VERBATIM static double bpeval(void* vv) { int i, n, no, nd, flag=0; double add,div; double *vo, *outp, *del; n = vector_instance_px(vv, &vo); no = vector_arg_px(1, &outp); nd = vector_arg_px(2, &del); if (ifarg(3) && ifarg(4)) { add= *getarg(3); div= *getarg(4); flag=1;} if (n!=no||n!=nd) hoc_execerror("v.bpeval: vectors not all same size", 0); if (flag) { for (i=0;i only look at these elements\n"); printf(" 'op'=='function name' is a .apply targeted by v2 called as func(x[i],thresh,val)\n"); return -1.; } op = gargstr(1); n = vector_instance_px(vv, &x); th = *getarg(2); if (ifarg(3)) { val = *getarg(3); } else { val = 0.0; } if (ifarg(4)) {ni = vector_arg_px(4, &ind); f=1;} // just look at the spots indexed if (!strcmp(op,"==")) { if (f==1) {for (i=0; i")) { 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;}} } else {for (i=0; i=th) { x[i] = val;}}} } else if (!strcmp(op,"<=")) { if (f==1) {for (i=0; ith) { /* ? passing thresh */ if (f==0) { if (j>=maxsz) { printf("(%d) :: ",maxsz); hoc_execerror("Dest vec too small in xing ", 0); } if (i>0) { /* don't record if first value is above thresh */ dest[j++] = tvec[i-1] + (tvec[i]-tvec[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 */ } } vector_resize(vv, j); return (double)i; } ENDVERBATIM :* v1.xzero() looks for zero crossings VERBATIM 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=x[i]=1; cnt++; } } return cnt; } ENDVERBATIM :* v1.sw(FROM,TO) switchs all FROMs to TO VERBATIM static double sw(void* vv) { int i, n; double *x, fr, to; n = vector_instance_px(vv, &x); fr = *getarg(1); to = *getarg(2); for (i=0; iu.this_pointer; // doesn't check that this is actually a bvec if (to->size!=n) { hoc_execerror("Vector and bytevec sizes don't match.", 0); } for (i=0; ix[i]; return n; } ENDVERBATIM :* v.v2d(&x) copies from vector to double area -- a seg error waiting to happen VERBATIM static double v2d(void* vv) { int i, n, num; double *x, *to; n = vector_instance_px(vv, &x); to = hoc_pgetarg(1); if (ifarg(2)) { num = *getarg(2); } else { num=-1;} if (num>-1 && num!=n) { hoc_execerror("Vector size doesn't match.", 0); } for (i=0; i-1 && num!=n) { hoc_execerror("Vector size doesn't match.", 0); } for (i=0; iu.this_pointer); *px = vector_vec(obv->u.this_pointer); return sz; } //* list_vector_px(LIST,ITEM#,DOUBLE PTR ADDRESS,VEC POINTER ADDRESS) // returns the vector pointer as well as the double pointer static int list_vector_px2 (Object *ob, int i, double** px, void** vv) { Object* obv; int sz; obv = ivoc_list_item(ob, i); sz = vector_capacity(obv->u.this_pointer); *px = vector_vec(obv->u.this_pointer); *vv = (void*) obv->u.this_pointer; return sz; } //* list_vector_resize(LIST,ITEM#,NEW SIZE) static int list_vector_resize (Object *ob, int i, int sz) { Object* obv; int maxsz; obv = ivoc_list_item(ob, i); 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 VERBATIM static double ismono1(double *x, double 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) { for (i=1; i100) { /* 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 %d :: ",points,maxsz); hoc_execerror("Vector max capacity too small in smgs ", 0); } } svar = -2.*var*var/step/step; scale = 1./sqrt(2.*M_PI)/var; for (j=0; j-20;j++) { Expo(arg); sum[j] += RES; } for (j=xv-1; j>=0 && (arg=(j-xv)*(j-xv)/svar)>-20;j--) { Expo(arg); sum[j] += RES; } } for (j=0; jtemplate != ob2->template) { return 0; } return 1; ENDVERBATIM } FUNCTION eqojt () { VERBATIM Object *ob1, *ob2; ob1 = *hoc_objgetarg(1); ob2 = *hoc_objgetarg(2); if (ob1 && ob2 && ob1==ob2) { return 1; } return 0; ENDVERBATIM }