source("gen_simdata.smoothgamma.r")
#source("opsampling.r")
#dyn.load("opsampling.so")
dyn.load("effnetroutines.so")
source("effnetsampling.r")

fileext=""


n=100 # number of samples
p=1000  # number of predictors
gammatrue=rep(0,p) # index of covariate.
true.index=c(245:260,745:760)
#true.index=c(45:60,145:160)
gammatrue[true.index]=1

b0=2
fstr.b0=paste("b0",b0,sep="_")
b=rep(0,p) 
b[gammatrue==1]= b0

sigmaX=1 # std of X
sigma=1 # noise variance.

gammatrue=as.matrix(gammatrue)
b=as.matrix(b)
simdata = gen_simdata.smoothgamma(n,p,gammatrue,b,sigmaX,sigma)
txx=t(simdata$xmat)%*%simdata$xmat
tyx=t(simdata$y)%*%simdata$xmat
txy=t(simdata$xmat)%*%simdata$y
tyy=sum(simdata$y^2)
y=as.vector(simdata$y)
txy=as.vector(txy)
# image(t(simdata$xmat)%*%simdata$xmat)  # look at cov(X)

cat("Data set generated.\n")

# ---------------------------#
#     gibbs parameters       #
# ---------------------------#

# (1) initialize gamma and setsort.

nop=as.integer(n/2) 
index=sample(2:p,nop)
index=sort(index)
gamma.ini=rep(0,p)
gamma.ini[index]=1
nsort=p    # output file size is O(p*nsort)
            # if only output marginal and not gamma, 
            # then simply set this to p.
            
setsort=c(sort(index),rep(0,(nsort-nop)))
nset=2*nsort    # nset is always twice nsort.
set=c(setsort,rep(0,nset-nsort))


############################################################
# SAVE STATE HERE.  SIMPLY LOAD THIS AND COMMENT OUT PREVIOUS
# CODE.


save.image(file = paste("simscript_smoothgamma_startfile_n",n,"_p",p,"_b",b0,".RData",sep=""))

#load(file="simscript_smoothgamma_startfile_n100_p1000_b2.RData")

############################################################


# (2) set the number of iterations, burn-in time.

niter=10000
nburnin=5000
nloop=niter-nburnin


# (3) These are parameters to study carefully: sparsity, smoothness, and bvarsqrt.

bvarsqrt=1 #bvarsqrt is the square root of bvar
r=0.025
w.attr=c(1,3,5)
NW = length(w.attr)

p.attr=(r*w.attr-r+1)/(r*w.attr+1)
q.attr=(r*w.attr)/(r*w.attr+1)

pi0=(1-q.attr[1])/(2-p.attr[1]-q.attr[1])
pi1=1-pi0

a=log(r)-2*log(r*w.attr-r+1)
#alpha=rep(a,p)
alpha=matrix(rep(a,p),nrow=p,ncol=length(w.attr),byrow=T)

#set bmati as empty.
bi=0
bii=log(w.attr)+log(r*w.attr-r+1)

nbi=1
nbii=3
nbi=as.integer(nbi)
nbii=as.integer(nbii)
bmati=matrix(0,nrow=p,ncol=nbi)
bmatii=matrix(0,nrow=p,ncol=nbii)
bmatii[,1]=2

bmatii[,2]=c(p,1:(p-1))
bmatii[,3]=c(2:p,1)

bmati=as.integer(bmati)
bmatii=as.integer(bmatii)
bmatii=matrix(bmatii,nrow=p,byrow=F)
bmati=matrix(bmati,nrow=p,byrow=F)

cat("Gibbs initialization done.\n")


# ---------------------------#
#     start Gibbs            #
# ---------------------------#

save.image(file = paste("simscript_smoothgamma2",fileext,".RData",sep=""))

# ---------------------------#
#     gibbsop                #
# ---------------------------#

outfn.attr=vector("list",NW)
bfout=vector("list",NW)

for(wind in c(1:NW)){

    w=w.attr[wind]
    
    cat("Starting gibbsop with w = ",w,".\n",sep="")
    
    outfn.attr[[wind]]=paste("smoothgamma","r",sep="_")
    outfn.attr[[wind]]=paste(outfn.attr[[wind]],r,sep="")
    outfn.attr[[wind]]=paste(outfn.attr[[wind]],"w",sep="_")
    outfn.attr[[wind]]=paste(outfn.attr[[wind]],w,sep="")
    outfn.attr[[wind]]=paste(outfn.attr[[wind]],"v",sep="_")
    outfn.attr[[wind]]=paste(outfn.attr[[wind]],bvarsqrt,sep="")
    outfn.attr[[wind]]=paste(outfn.attr[[wind]],"b",sep="_")
    outfn.attr[[wind]]=paste(outfn.attr[[wind]],b0,sep="")
    outfn.attr[[wind]]=paste(outfn.attr[[wind]],fileext,sep="")
    bfout[[wind]]=paste(outfn.attr[[wind]],"_bf",sep="")
 
    unix.time(gibbsneteff(txx,tyx,tyy,p,n,gamma.ini,alpha[,wind],bi,bii[wind],bmati,bmatii,nbi,nbii,bvarsqrt,set,nset,setsort,nsort,nop,niter,nburnin,outfn.attr[[wind]],bfout[[wind]]))

}


# ---------------------------#
#     look at results        #
# ---------------------------#

margin.attr=vector("list",NW)
modelsize.attr=vector("list",NW)

for(wind in c(1:NW)){
    margin.attr[[wind]] = scan(outfn.attr[[wind]])
    modelsize.attr[[wind]] = margin.attr[[wind]][1:nloop]
    margin.attr[[wind]] =margin.attr[[wind]][(nloop+1):(nloop+p)]
    margin.attr[[wind]] = margin.attr[[wind]]/nloop
 
}

postscript(paste("smoothgamma_r",r,"_b",b0,"_v",bvarsqrt,fileext,"_margin.ps",sep=""),horizontal=F)
#jpeg(paste("smoothgamma_r",r,"_b",b0,"_v",bvarsqrt,"_margin.png",sep=""), height=1000, width=500)
par(mfrow=c(NW,1))
for(wind in c(1:NW)){
    plot(margin.attr[[wind]],ylim=c(0,1),ylab="marginal probability",xlab="motif index")
    points(true.index,margin.attr[[wind]][true.index],col = "red",type="h")
    title(paste("b=",b0,"v=",bvarsqrt,"r=", r, "w=", w.attr[[wind]]));
}
dev.off()


postscript(paste("smoothgamma_r",r,"_b",b0,"_v",bvarsqrt,fileext,"_modelsize.ps",sep=""),horizontal=F)
#jpeg(paste("smoothgamma_r",r,"_b",b0,"_v",bvarsqrt,"_modelsize.png",sep=""), height=1000, width=500)
par(mfrow=c(NW,1))
for(wind in c(1:NW)){
    plot(modelsize.attr[[wind]],col = "black",ylab="model size",xlab="number of iteration")
    title(paste("b=",b0,"v=",bvarsqrt,"r=", r, "w=", w.attr[[wind]]));
}
dev.off()

# ---------------------------------- #
# Do ROC curve 
# ---------------------------------- #

ROC=vector("list",NW)

totaltrue = length(true.index)
totalfalse = p-totaltrue
for(wind in c(1:NW)){
    margin=margin.attr[[wind]]
    thresholds=seq(1,0,-0.01)
    TP=rep(0,length(thresholds))
    FP=rep(0,length(thresholds))
    for(threshind in c(1:length(thresholds))){
        thresh=thresholds[threshind]
        npos = sum(margin>thresh)
        ntrue = sum(margin[true.index]>thresh)
        nfalse = npos-ntrue
        TP[threshind] = ntrue/totaltrue
        FP[threshind] = nfalse/totalfalse
    }
    ROC[[wind]] =list(TP=TP,FP=FP)
 
}

postscript(paste("smoothgamma_r",r,"_b",b0,"_v",bvarsqrt,fileext,"_ROC.ps",sep=""),horizontal=F)
plot(ROC[[1]]$FP,ROC[[1]]$TP,type="l",col=1)
for(wind in c(2:NW)){
    lines(ROC[[wind]]$FP,ROC[[wind]]$TP,col=wind)
}
legend(0.6, 0.2, c("1","3","5"),lty=1,col = c(1:NW))
dev.off()