// $Id: infot.hoc,v 1.43 2009/12/04 01:25:55 samn Exp $
print "Loading infot.hoc..."
if(!installed_infot) {install_infot()}
//* tentropsig(v1,v2,nshuf,nbins,twoway,[xpast,ypast,hval])
//significance test of transfer entropy using shuffling
//returns (te - tes) / sds , where te is transfer entropy, tes is transfer entropy
//of shuffled data, sds is std-dev of transfer entrop of shuffled data
//should only accept as significanat values > 4-6
func tentropsig () { local nshuf,i,xp,yp,hv,nbins,te localobj v1,v2,vo
v1=new Vector() v2=new Vector() vo=new Vector(1)
v1.copy($o1) v2.copy($o2) nshuf=$3 nbins=$4
if(numarg()>4) xp=$5 else xp=1
if(numarg()>5) yp=$6 else yp=2
if(numarg()>6) hv=$7 else hv=0
te=v1.tentrop(v2,nbins,xp,yp,nshuf,vo,hv)
if(1||verbose_infot>2) printf("te=%g,sig=%g\n",te,vo.x(0))
return vo.x(0)
}
//** mutinfbshufv(v1,v2,[nshuf,nbins])
//return vector with mutual information from shuffled v1,v2
//used for significance test , i.e. : ((miorig - mishufmean) / mishufstdev) > 2
obfunc mutinfbshufv () { local nshuf,nbins,i localobj v1,v2,ve
v1=new Vector() v2=new Vector() ve=new Vector()
v1.copy($o1) v2.copy($o2)
if(numarg()>2) nshuf=$3 else nshuf=20
if(numarg()>3) nbins=$4 else nbins=10
for i=0,nshuf-1 {
v1.shuffle() v2.shuffle()
ve.append(v1.mutinfb(v2,nbins))
}
return ve
}
//** mutinfbsig(v1,v2,[nshuf,nbins])
//get significance of mutual information, should be at least > 2
func mutinfbsig () { local nshuf,nbins,st localobj ve
if(numarg()>2) nshuf=$3 else nshuf=20
if(numarg()>3) nbins=$4 else nbins=10
ve=mutinfbshufv($o1,$o2,nshuf,nbins)
st=ve.stdev
if(st<=0) st=1
return ($o1.mutinfb($o2,nbins) - ve.mean) / st
}
//** tentropspksig(v1,v2,nshuffles)
//get significance of tentropspks using shuffling
//returns (TE - AvgTEShuffle) / StdDevTEShuffle
func tentropspksig () { local nshuf,i,xp,yp,hv,nbins,te,sd localobj v1,v2,ve
v1=new Vector() v2=new Vector() ve=new Vector()
v1.copy($o1) v2.copy($o2) nshuf=$3
te=$o1.tentropspks($o2)
for i=0,nshuf-1 {
v1.shuffle() ve.append(v1.tentropspks(v2))
}
if(verbose_infot>2) printf("te=%g,ve.mean=%g,ve.stdev=%g\n",te,ve.mean,ve.stdev)
if(verbose_infot>2) ve.printf
sd=ve.stdev()
if(sd<=0)sd=1
return (te-ve.mean)/sd
}
//* normte() get normalized transfer entropy using tentropspks in output vector vo
//vo.x(0)=transfer entropy of $o1->$o2
//vo.x(1)=H($o2Future|$o2Past)
//vo.x(2)=normalized transfer entropy in 0,1 range
//$3==number of shuffles
//$o1,$o2 should both have same size and non-negative values. this func is meant for time-binned spike train data
obfunc normte () { local a localobj ve,vo
a=allocvecs(ve) vo=new Vector()
nshuf=0
nshuf=$3 vrsz(3+nshuf,vo)
te=$o1.tentropspks($o2,vo,nshuf)
if(verbose_infot>2) vo.printf
if(vo.x(1)<=0 && verbose_infot>0){printf("WARNING H(X2F|X2P)==%g<=0\n",vo.x(1)) vo.x(1)=1 }
if (nshuf>0) {
ve.copy(vo,3,vo.size-1)
vo.resize(4)
if (ve.mean!=vo.x[2]) printf("normte ERRA\n")
vo.append(ve.stdev)
}
vo.x[2]=te
dealloc(a)
return vo
}
//* GetTENQ() get an nqs with useful transfer entropy info
obfunc GetTENQ () { local te01,te10,pf01,pf10 localobj nqte,vo1,vo2
if(numarg()>3) nqte=$o4
if(nqte==nil) {
nqte=new NQS("from","to","TE","NTE","HX2|X2P","prefdir","TEshufavg","TEshufstd","sig")
} else nqte.clear()
vo1=normte($o1,$o2,$3)
vo2=normte($o2,$o1,$3)
te01=vo1.x(2)
te10=vo2.x(2)
if(vo1.x(4)<=0)vo1.x(4)=1
if(vo2.x(4)<=0)vo2.x(4)=1
if(te01>0 || te10>0) {
pf01=(te01-te10)/(te01+te10)
pf10=(te10-te01)/(te01+te10)
} else {
pf01=pf10=0
}
nqte.append(0,1,vo1.x(0),te01,vo1.x(1),pf01,vo1.x(3),vo1.x(4),(vo1.x(0)-vo1.x(3))/vo1.x(4))
nqte.append(1,0,vo2.x(0),te10,vo2.x(1),pf10,vo2.x(3),vo2.x(4),(vo2.x(0)-vo2.x(3))/vo2.x(4))
return nqte
}
//** prefdte() get preferred direction of transfer entropy
//$o1=vec 1, $o2=vec 2, $3 = # of times to shuffle
func prefdte () { local nshuf,a,te01,te10,pfd localobj v1,v2,vtmp
a=allocvecs(v1,v2,vtmp)
v1.copy($o1) v2.copy($o2) nshuf=$3 vtmp.resize(3)
v1=normte($o1,$o2,nshuf)
v2=normte($o2,$o1,nshuf)
te01=v1.x(2)
te10=v2.x(2)
pfd=(te01-te10)/(te01+te10)
dealloc(a)
return pfd
}
//** mkchist() averages entries in window into disc values and returns in new output vec
//$o1=input vec,$2=win size
obfunc mkchist () { local idx,eidx,wsz localobj vin,vout
vin=$o1 wsz=$2 vout=new Vector()
vout.resize(1+vin.size/wsz) vout.resize(0)
for(idx=0;idx<=vin.size;idx+=wsz) {
eidx=idx+wsz-1
if(eidx>=vin.size)eidx=vin.size-1
if(eidx>idx) vout.append( int(vin.mean(idx,eidx)) )
}
return vout
}
//get magnitude of difference in a preferred direction - just abs of diff, but if theyre both neg, return 0
//$1 = nTE_X->Y
//$2 = nTE_Y->X
func prefdmag () { local n1,n2,s
n1=$1 n2=$2
if(n1>0 && n2<=0) return n1-n2 //n1 is relatively strong
if(n2>0 && n1<=0) return n2-n1 //n2 is relatively strong
if(n2<0 && n1<0) return 0 //both are weak
return abs(n1-n2) //both are weak positive
}
//** simple test for nte
{declare("vb","o[2]","vs","o[2]")}
for i=0,1 {
vb[i]=new Vector()
vs[i]=new Vector()
}
//mkspktrain(Random,rate,tmax) -- make a spike train with specified rate,tmax
//Random obj must be initialized
obfunc mkspktrain () { local tmax,rate,t,dt,intt localobj rdp,vs
rdp=$o1 rate=$2 tmax=$3
intt=1e3/rate
t = 0
vs=new Vector()
while(t<=tmax) {
dt = rdp.poisson(intt)
t += dt
vs.append(t)
}
return vs
}
//make random spikes with frequency $1, tmax=$2, offset for spikes=$3, alpha=$4 -- ratio of spikes from
//vs[0] that get placed in vs[1]
//spikes in vs[0] are randomly picked, spikes in vs[1] are same as in vs[0] but shifted forward by $3 offset
//so vs[0] 'drives' vs[1], or can be used to predict it, but vs[1] cant be used to predict vs[0]
proc mkspks () { local tmax,rate,t,dt,intt,off,i,alpha localobj rdp
rate=$1 tmax=$2 off=$3
if(numarg()>3)alpha=$4 else alpha=1
intt=1e3/rate
rdp=new Random()
rdp.ACG(1234)
rdp.poisson(intt)
for i=0,1 vs[i].resize(0)
vs[0]=mkspktrain(rdp,rate,tmax)
if(alpha < 1.0) {
for vtr(&t,vs[0]) if(rdp.uniform(0,1) <= alpha) vs[1].append(t+off)
} else {
vs[1].copy(vs[0])
vs[1].add(off)
}
}
//test nTE : nTE of X0 -> X1 should be much higher than nTE of X1 -> X0
//optional $1=offset == offset to shift spikes by, in ms
//optional $2=rate == rate of spikes, in Hz
//optional $3=bin size , in ms
//optional $4=alpha == ratio of spikes of X0 that get placed in X1 with offset
//optional $5=max time, in ms
func testnte () { local a,i,bisv,maxt,alpha,off,rate,binsz,dur localobj nqt,nqout
if(numarg()>0)off=$1 else off=10
if(numarg()>1)rate=$2 else rate=50
if(numarg()>2)binsz=$3 else binsz=10
if(numarg()>3)alpha=$4 else alpha=1
if(numarg()>4)dur=$5 else dur=10000
bisv=binmin_infot binmin_infot=0
print "output should be close to:\n\t0 1 0.6707 0.9975 0.6707 0.9407 0.003435 0.001672 399.1"
print "\t1 0 0.02183 0.03049 0.6685 -0.9407 0.002828 0.001442 13.17"
mkspks(rate,dur,off,alpha)
maxt=vs[1].max
printf("maxt=%g\n",maxt)
if(vs[1].max>maxt)maxt=vs[1].max
for i=0,1 vb[i].hist(vs[i],0,(maxt+binsz-1)/binsz,binsz)
nqt=GetTENQ(vb[0],vb[1],200)
nqout=new NQS("X1","X2")
nqout.odec("X1")
nqout.odec("X2")
batch_flag=1
nqout.append(vb[0],vb[1])
nqout.sv("/u/samn/bpftest/data/09dec17.func.testnte.nqs")
batch_flag=0
nqt.pr
nqsdel(nqt)
binmin_infot=bisv
return 1
}
//get kernel smoothed prob distrib in an nqs
//$o1=input vector
//$2=increment in x , smaller values mean finer resolution
//$3=bandwidth - higher means smoother output
// $4=min value in output, $5=max value in output
obfunc khist () { local min,max,inc,h,x,i,s localobj vx,vy,nq,vin
vin=$o1
if(numarg()>1)inc=$2 else inc=0.1
if(numarg()>2)h=$3 else h=vin.getbandwidth()
if(numarg()>3)min=$4 else min=vin.min()
if(numarg()>4)max=$5 else max=vin.max()
{vx=new Vector() vy=new Vector()}
vx.indgen(min,max,inc)
vy.copy(vx)
for vtr(&x,vx,&i) vy.x(i) = vin.kprob1D(h,x)
s=vy.sum
if(s!=0) vy.div(vy.sum)
nq=new NQS("x","y")
nq.v[0]=vx
nq.v[1]=vy
return nq
}
