: \$Id: vecst.mod,v 1.499 2011/07/22 22:16:48 billl Exp \$ :* COMMENT COMMENT 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 whi find indices where vec equal to given values dest.xing(src,tvec,thresh) determines where vector crosses in a positive direction 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 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 ind.mredundout(veclistA[,INDFLAG,veclistB]) // remove repeats from parallel vectors v.cvlv(v1,v2) // convolve v1 with v2 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 smsy sum syn potentials off a tvec iwr write integers ird read integers ident give pointer addresses and sizes lcat concatentate all vectors from a list fread2 like vec.fread but uses flag==6 for unsigned int vfill fill vdest with multiple instances of vsrc until reach size inv return inverse of number multiplied by optional num slone(src,val) select indices where ==VAL from sorted vector SRC piva.join(pivb,veclista,veclistb) copies values from one set of vecs to other Non-vector routines 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 sumabs return sum of absolute values rnd round off to nearest integer ENDCOMMENT NEURON { SUFFIX nothing : BVBASE is bit vector base number (typically 0 or -1) GLOBAL BVBASE, RES, VECST_INSTALLED, DEBUG_VECST, VERBOSE_VECST, INTERP_VECST, LOOSE } PARAMETER { BVBASE = 0. 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 : 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 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 } 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; } double *dcrset (int nx) { if (nx>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; } ENDVERBATIM :* v1.ident() gives addresses and sizes VERBATIM static double ident (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",(unsigned int)vv,(unsigned int)x,nx,bsz); return (double)nx; } 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(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(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(vv,nx); y=vector_newsize(vector_arg(1),ny=nx); } for (i=0,theta=0; inx) {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 (scrsz>0) { free(scr); scr=(unsigned int *)NULL; } scrsz=nx+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;iVRRY) hoc_execerror("findx ****ERRB****: can only handle VRRY vectors", 0); for (i=0;iscrsz) { if (scrsz>0) { free(scr); scr=(unsigned int *)NULL; } scrsz=ni+10000; scr=(unsigned 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;j0 && 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("sindx ****ERRB****: can only handle VRRY vectors", 0); for (i=0;iscrsz) { if (scrsz>0) { free(scr); scr=(unsigned int *)NULL; } scrsz=ni+10000; scr=(unsigned 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;jVRRY) hoc_execerror("sindv ****ERRA****: can only handle VRRY vectors", 0); for (i=0;iscrsz) { if (scrsz>0) { free(scr); scr=(unsigned int *)NULL; } scrsz=ni+10000; scr=(unsigned 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;jVRRY) 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]) { 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); } } for (j=0;j=ALL) { field[j]=0; } else for (m=1;m<=5;m++) { if (key[j]<=EBE*(m+1) && key[j]>=ALL*(m+1)) { // m is is field key 1-5 key[j]/=(m+1); field[j]=m; } } if (field[j]==-1) {printf("vecst::slct ERRF %d %g\n",j,key[j]); hxe(); } } if (2*nk!=na) { printf("vecst::slct ERR3 %d %d\n",nk,na); hoc_execerror("Arg vector must be double key length",0); } for (i=0,n=0;i=EQV && key[i]<=EQY) n++; // special cases take 2 vec args if (nk+n!=num) { printf("vecst::slct ERR2 %d(keys)+%d(EQV/W)!=%d(vecs)\n",nk,n,num); hoc_execerror("Key length must be number of vecs + num of EQV/W",0); } for (j=0,k=0,m=0;j=0.) {fl=0; break;} else continue; } else if (key[n]==GTH) { if (val<=arg[m]) {fl=0; break;} else continue; } else if (key[n]==GTE) { if (val< arg[m]) {fl=0; break;} else continue; } else if (key[n]==LTH) { if (val>=arg[m]) {fl=0; break;} else continue; } else if (key[n]==LTE) { if (val> arg[m]) {fl=0; break;} else continue; } else if (key[n]==EQU) { if (val!=arg[m]) {fl=0; break;} else continue; } else if (key[n]==EQV) { if (val!=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;pvvo[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 } else if (key[n]==EBI) { if ((val<=arg[m])||(val> arg[m+1])) { fl=0; break; } else continue; // "(]" : exclude-bracket-include } else if (key[n]==IBI) { if ((val< arg[m])||(val> arg[m+1])) { fl=0; break; } else continue; // "[]" : include-bracket-include } else if (key[n]==EBE) { if ((val<=arg[m])||(val>=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 if (flag==1) break; } } vector_resize(vv, k); return (double)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[VRRY], num, fl, field[VRRY], lt, rt; Object* lob; double *ind, *key, *arg, *vvo[VRRY], val; 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>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]); hoc_execerror("index and searched vectors must all be same size: ", 0); } } for (j=0;j=ALL) { field[j]=0; } else for (m=1;m<=5;m++) { if (key[j]<=EBE*(m+1) && key[j]>=ALL*(m+1)) { // m is is field key 1-5 key[j]/=(m+1); field[j]=m; } } if (field[j]==-1) {printf("vecst::slor ERRF %g\n",key[j]); hxe(); } } if (2*nk!=na) { printf("vecst::slor ERR3 %d %d\n",nk,na); hoc_execerror("Arg vector must be double key length",0); } for (i=0,n=0;i=EQV && key[i]<=EQY) n++; // special case takes 2 vec args if (nk+n!=num) { printf("vecst::slor ERR2 %d(keys)+%d(EQV)!=%d(vecs)\n",nk,n,num); hoc_execerror("Key length must be number of vecs + num of EQV",0); } for (j=0,k=0,m=0;j=0.) continue; else {fl=1; break;} } else if (key[n]==GTH) { if (val<=arg[m]) continue; else {fl=1; break;} } else if (key[n]==GTE) { if (val< arg[m]) continue; else {fl=1; break;} } else if (key[n]==LTH) { if (val>=arg[m]) continue; else {fl=1; break;} } else if (key[n]==LTE) { if (val> arg[m]) continue; else {fl=1; break;} } else if (key[n]==EQU) { if (val!=arg[m]) continue; else {fl=1; break;} } else if (key[n]==EQV) { if (val!=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;pvvo[i+1][p]) lt=p+1; else if (val=arg[m+1])) { continue; } else {fl=1; break;} // IBE="[)" include-bracket-exclude } else if (key[n]==EBI) { if ((val<=arg[m])||(val> arg[m+1])) { continue; } else {fl=1; break;} // "(]" : exclude-bracket-include } else if (key[n]==IBI) { if ((val< arg[m])||(val> arg[m+1])) { continue; } else {fl=1; break;} // "[]" : include-bracket-include } else if (key[n]==EBE) { if ((val<=arg[m])||(val>=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 (double)k; } ENDVERBATIM :* src.whi(valvec,indvec[,&i0,&i1,...]) returns index of first one for each val 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=n+10000; scr=(unsigned 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); } } r=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(vv, n); } if (type==6 || type==16) { // unsigned ints unsigned int *xs; if (n>scrsz) { if (scrsz>0) { free(scr); scr=(unsigned int *)NULL; } scrsz=n+10000; scr=(unsigned int *)ecalloc(scrsz, sizeof(int)); } xs=(unsigned int*)scr; r=fread(xs,sizeof(int),n,fp); if (type==16) BYTESWAP_FLAG=1; for (i=0;i=nx) x=vector_newsize(vv,nx*=4); x[k]=*getarg(i); } else { ny=vector_arg_px(i, &y); if (k+ny>=nx) x=vector_newsize(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; } ENDVERBATIM :* ind.insct(v1,v2) : return v1 intersect v2 VERBATIM 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(vv); x=vector_newsize(vv, maxsz); nv1 = vector_arg_px(1, &v1); nv2 = vector_arg_px(2, &v2); for (i=0,k=0;iisz;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(vv, maxsz*=2); x[k++]=val; break; } } } vector_resize(vv, k); return (double)k; } ENDVERBATIM :* vdest.vfill(vsrc) : fill vdest with multiple instances of vsrc until reach size VERBATIM static double vfill (void* vv) { int i, nx, nv1; double *x, *v1; nx = vector_instance_px(vv, &x); nv1 = vector_arg_px(1, &v1); 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;iscrsz) { if (scrsz>0) { free(scr); scr=(unsigned int *)NULL; } scrsz=ns/4+10000; scr=(unsigned int *)ecalloc(scrsz, sizeof(int)); } for (j=0;j=scrsz){printf("mredundout****ERRD**** over scr size %d\n",k);hxe();} scr[k++]=i; // flag to get rid of this one } else for (j=0;jmaxsz){printf("mredundout****ERRE**** vec overflow %d>%d\n",k,maxsz);hxe();} for (i=0;iVRRY) hoc_execerror("join ****ERRA****: can only handle VRRY vectors", 0); if (num!=numb) hoc_execerror("join ****ERRB****: different #of vecs in lists", 0); for (i=0;imax) max=x[i]; if (x[i]nsrc) { hoc_execerror("sccvlv:Filter bigger than source ", 0); } if (nfilt!=ntmp){hoc_execerror("sccvlv:Filter (arg2) and tmp vector (arg3) diff size", 0);} for (i=0;i0&&kth) { // ? passing thresh if (f==0) { sum+=1; // just count the 0xing f=1; } } else { // below thresh if (f==1) { f=0; // just passed going down sum+=1; } } } if (sum>maxsum) maxsum=sum; x[i]=sum; } return (double)maxsum; } ENDVERBATIM :* vdest.cvlv(vsrc,vfilt) : convolution VERBATIM static double cvlv (void* vv) { int i, j, k, nx, nsrc, nfilt; double *x, *src, *filt, sum, lpad, rpad; nx = vector_instance_px(vv, &x); nsrc = vector_arg_px(1, &src); nfilt = vector_arg_px(2, &filt); if (nx!=nsrc) { hoc_execerror("Vectors not same size: ", 0); } if (nfilt>nsrc) { hoc_execerror("Filter bigger than source ", 0); } for (i=0;i0 && kVRRY) hoc_execerror("ERR: nind can only handle VRRY vectors", 0); num = i-2; // number of vectors to be picked apart for (i=0;iscrsz) { if (scrsz>0) { free(scr); scr=(unsigned int *)NULL; } scrsz=ni+10000; scr=(unsigned 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;jVRRY) hoc_execerror("ERR: keyind can only handle VRRY 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 (double)k; } 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.nearall(max,v2,tql) : tql from {tq=new NQS("diff","tmid","t1","t2") tq.listvecs(tql)} : max is maximum diff to add to the tq db VERBATIM static double nearall (void* vv) { register int lo, hi, mid; int i, j, k, kk, nx, ny, minind, nv[4]; Object *ob; void* vvl[4]; double *x, *y, *vvo[4], targ, dist, new, max, tmp; nx = vector_instance_px(vv, &x); max = *getarg(1); ny = vector_arg_px(2, &y); ob = *hoc_objgetarg(3); if ((nv[0]=ivoc_list_count(ob))!=4) {printf("VECST::nearall() ERRA: %d\n",nv[0]); hxe();} for (i=0;i<4;i++) { nv[i] = list_vector_px3(ob, i, &vvo[i], &vvl[i]); if (nv[i]!=nv[0]){printf("nearall ERRC: %d %d\n",nv[i],nv[0]); hxe();} } for (j=0,k=0;j0.) 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(vv); iNewSz = (int)dNewSz; if (iNewSz>maxsz) {printf("VECST samp ERRA: dest vec too small: %d %d\n",iNewSz,maxsz); hxe();} vector_resize(vv,iNewSz); scale = (double) iOrigSz / (double) iNewSz; 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 (double)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); f=c=0; 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;c++;}} } else {for (i=0; ith) { x[i]=val;c++;}}} } 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++;}}} } else if (!strcmp(op,"<=")) { if (f==1) {for (i=0; i=0;i--); for ( ;src[i]==val&&i>=0;i--); i++; // go back to the first one } for (j=0;src[i]==val && jth) { // ? 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 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]); } } else { if (tvf) dest[j++]=tvec[i]; else dest[j++]=i; } } 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]); } } else { if (tvf) dest2[j++]=tvec[i]; else dest2[j++]=i; } } } } vector_resize(vv, j); if (d2f) vector_resize(vector_arg(d2f),j); return (double)j; } ENDVERBATIM :* dest.snap(src,tvec,dt) : interpolate src with tvec to prior dt step, saves only highest value in each interval : an .interpolate that doesn't loose spikes VERBATIM static double snap (void* vv) { int i, j, nsrc, ndest, ntvec, f, maxsz, size; double *src, *dest, *tvec, mdt, tstop, tt, val; ndest = vector_instance_px(vv, &dest); nsrc = vector_arg_px(1, &src); ntvec = vector_arg_px(2, &tvec); mdt = *getarg(3); maxsz=vector_buffer_size(vv); tstop = tvec[nsrc-1]; size=(int)tstop/mdt; if (size>maxsz) { printf("%d > %d\n",size,maxsz); hoc_execerror("v.snap: insufficient room in dest", 0); } 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"); dest[i]=val; } } return (double)size; } ENDVERBATIM :* v1.xzero(vec[,thresh]) finds indices of zero [or thresh] crossings in vec VERBATIM static double xzero (void* vv) { int i, n, nv, up, cnt; double *x, *vc, th; 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(vv,(n+=100)); x[cnt++]=(double)i; } } x=vector_newsize(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; void* vy; n = vector_instance_px(vv, &x); // for indices nc = vector_arg_px(1, &vc); // source vectors if (n==0) x=vector_newsize(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(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]; } x=vector_newsize(vv,cnt); if (ny) y=vector_newsize(vy,cnt); return (double)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) { int i, n; double *x, cnt, sig; n = vector_instance_px(vv, &x); if (ifarg(1)) sig = (int)*getarg(1); else sig=1e9; // default: do wrap if (sig==0.) { for (i=0,cnt=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 (double)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=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(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,(unsigned int)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,(unsigned int)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(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) { int i, n, num; double *x, *fr; n = vector_instance_px(vv, &x); fr = 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=sc[4] || floor(vvo[j][i]+0.5)!=vvo[j][i]) { printf("vec.mkcode OOB %g>%g in vec[%d].x[%d]\n",vvo[j][i],sc[4],j,i); hxe(); } x[i]+=vvo[j][i]*sc[j+1]; } } return (double)i; } ENDVERBATIM :* v.uncode(val) -- take apart val and place in vector : v.uncode(VECLIST) -- take apart vector items and place in vectors in list (cf uncodf) : v.uncode(vec,field) -- take apart v entries and place requested field in vec : v.uncode(field,val) -- replace field in v with val (cf recodf) : v.uncode(field,vec) -- replace field in v with values from vec VERBATIM static double uncode (void* vv) { int i, j, n, ny, num, field; Object *ob; double *x, *y, *vvo[5], val, old; void *vvv[5]; n = vector_instance_px(vv, &x); field=0; if (!ifarg(1)) { // numarg()==0 printf("\tv.uncode(val) -- take apart val and place in vector\n\tv.uncode(VECLIST) -- take apart vector items and place in vectors in list (cf uncodf)\n\tv.uncode(vec,field) -- take apart vector items and place requested field in vector\n\tv.uncode(field,val) -- replace field in v with val (cf recodf)\n\tv.uncode(field,vec) -- replace field in v with values from vec\n"); return 0.; } else if (!ifarg(2)) { // numarg()==1 if (hoc_is_double_arg(1)) { val = *getarg(1); if (vector_buffer_size(vv)<5) { hoc_execerror("uncode ****ERRA****: vector too small to resize(5)", 0);} vector_resize(vv,5); for (i=1;i<=5;i++) UNCODE(val,i,x[i-1]) return x[0]; } else { ob = *hoc_objgetarg(1); num = ivoc_list_count(ob); if (num>5) hoc_execerror("uncode ****ERRA****: can only handle 5 vectors", 0); for (i=0;i0) { if (y[i]<0.||y[i]>=sc[4]||floor(y[i]+0.5)!=y[i]) { printf("vec.uncode ERRJ OOB %g (%g max) at %d\n",y[i],sc[4],i);hxe();} x[i] += sc[field]*(y[i]-old); } else { x[i] += sc[field]*(val -old); } } return (double)i; } else { // fill single vector with values ny = vector_arg_px(1, &y); field = (int)chkarg(2,1.,5.); if (ny!=n) hoc_execerror("uncode(vec) ****ERRI****: diff sized vecs", 0); for (i=0;iu.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, void** vv) { Object* obv; int sz; obv = ivoc_list_item(ob, i); if (! ISVEC(obv)) return -1; 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, void** vv) { Object* obv; int sz; obv = ivoc_list_item(ob, i); if (! ISVEC(obv)) return -1; 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(obv->u.this_pointer); *px = vector_vec(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(vv,sz); return 0; } else vector_resize(vv,n); return 1; } //* list_vector_resize(LIST,ITEM#,NEW SIZE) double *list_vector_resize (Object *ob, int i, int sz) { Object* obv; obv = ivoc_list_item(ob, i); if (! ISVEC(obv)) return 0x0; vector_resize(obv->u.this_pointer,sz); return vector_vec(obv->u.this_pointer); } 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 if (flag==1) { for (i=1; i=last; i++) last=x[i]; } else if (flag==-1) { for (i=1; ilast; i++) last=x[i]; } else if (flag==-2) { 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 } :* 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; double *x, *sum; double low , high , step , var , svar , scale , arg; nsum = vector_instance_px(vv, &sum); nx = vector_arg_px(1,&x); low = *getarg(2); high = *getarg(3); step = *getarg(4); var = *getarg(5); points = (int)((high-low)/step+hoc_epsilon); if (nsum!=points) { maxsz=vector_buffer_size(vv); if (points<=maxsz) { nsum=points; vector_resize(vv, nsum); } else { printf("%d > %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; j %d :: ",points,maxsz); hoc_execerror("Dest vector too small in smsy ", 0); } } // don't zero out dest vec for (i=0;i100) hoc_execerror("rdmany ERRA: can only handle 100 vectors", 0); if (num!=cnt) {printf("rdmany ERRB: %d != %d",num,cnt); hxe();} for (i=0;ibufsz) { if (scrsz>0) { free(scr); scr=(unsigned int *)NULL; } scrsz=vsz+10; scr=(unsigned int *)ecalloc(scrsz, sz); bufsz=scrsz*sz; // number of chars available } if (code==2) { unsigned short *xs; xs=(unsigned short *)scr; for (last=-1,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.499 2011/07/22 22:16:48 billl Exp \$\n") } else { VECST_INSTALLED=1 VERBATIM { 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("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); install_vector_method("keyind", keyind); 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); install_vector_method("smgs", smgs); install_vector_method("smsy", smsy); install_vector_method("ident", ident); install_vector_method("lcat", lcat); install_vector_method("snap", snap); install_vector_method("fread2", fread2); install_vector_method("vfill", vfill); install_vector_method("vrdh", vrdh); install_vector_method("mkcode", mkcode); install_vector_method("uncode", uncode); install_vector_method("sumabs", sumabs); install_vector_method("inv", inv); install_vector_method("join", join); 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 FUNCTION isojt () { VERBATIM { Object *ob1, *ob2; ob1 = *hoc_objgetarg(1); ob2 = *hoc_objgetarg(2); if (!ob1) if (!ob2) return 1; else return 0; if (!ob2 || ob1->template != ob2->template) { return 0; } return 1; } ENDVERBATIM } : isojn(OB1,NAME) return whether OB1 is an instance of EXAMPLE_OBJ FUNCTION isojn () { VERBATIM { Object *ob1; char* name; ob1 = *hoc_objgetarg(1); name = gargstr(2); if (strncmp(hoc_object_name(ob1),name,3)==0) _lisojn=1.; else _lisojn=0.; } ENDVERBATIM } : ojtnum(OBJ) returns object number : returns internal number of object, eg if vec[3] is Vector[432] returns 432 FUNCTION ojtnum () { VERBATIM { Object *ob1; char name[50]; int ii; ob1 = *hoc_objgetarg(1); if (!ob1) return -1; if (ifarg(2)) { strncpy(name, hoc_object_name(ob1),50); for (ii=strlen(name);ii>1;ii--) if (name[ii]==91) {name[ii]=0; break;} // 91 is [ hoc_assign_str(hoc_pgargstr(2),name); } return (double)ob1->index; } ENDVERBATIM } : eqojt(OB1,OB2) return whether OB1 and OB2 point to same object FUNCTION eqojt () { VERBATIM { Object *ob1, *ob2; ob1 = *hoc_objgetarg(1); ob2 = *hoc_objgetarg(2); if (ob1 && ob2 && ob1==ob2) { return 1; } return 0; } ENDVERBATIM } :* byteswap(FILE) FUNCTION byteswap () { VERBATIM { int n[2]; size_t r; double ret; FILE* f; BYTEHEADER f = hoc_obj_file_arg(1); r=fread(&n,sizeof(int),2,f); if (n[1] < 1 || n[1] > 5) { BYTESWAP_FLAG = 1; ret = 1.; } else ret = 0.; BYTESWAP(n[1],int) if (n[1] < 1 || n[1] > 5) { printf("byteswap: Something wrong with location sampled: %d\n",n[1]); ret = -1.; } fseek(f,-2*sizeof(int),SEEK_CUR); // go back to where we started return ret; } ENDVERBATIM } : mkcodf(val1,val2,val3,val4,val5) stuff 5 vals<=999 into a single double 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(); } x+=a*sc[i]; } return x; } ENDVERBATIM } : uncodf(code,i) returns field i (1-5) from code FUNCTION uncodf () { VERBATIM { int i; 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();} UNCODE(x,i,ret); return ret; } else { for (i=2;i<=6;i++) if (ifarg(i)) { ptr = hoc_pgetarg(i); UNCODE(x,i-1,*ptr); } else break; return *ptr; } } ENDVERBATIM } : recodf(i,code,new) replaces field i (1-5) from code with new FUNCTION recodf () { VERBATIM { int i; double x, y, old; i=(int)chkarg(1,1.,5.); x=*getarg(2); y=chkarg(3,0.,sc[4]-1); UNCODE(x,i,old); return x + sc[i]*(y-old); } ENDVERBATIM } : flor(val) FUNCTION flor () { VERBATIM { return floor(*getarg(1)); } 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(obv->u.this_pointer); pList->plen[i]=vector_capacity(obv->u.this_pointer); pList->pbuflen[i]=vector_buffer_size(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(obv->u.this_pointer); pList->plen[i]=vector_capacity(obv->u.this_pointer); pList->pbuflen[i]=vector_buffer_size(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(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