rm(list=ls())

#Running this code will save Rdata files summarizing all of the methods.

#Read in data
setwd("INSERT DIRECTORY CONTAINING FILES")
source("FunctionLibrary.R")
source("FARgSCAD_functions.R")

#Movement: "H"= hand, "W"=wrist
movement="H"

#Choose subset of analysis methods
Analysis=c("SAFE","ag","LAD","FAR")

#number of blocks, SAFE stages, and delta value
nblocks<-5
nstage<-5
delta<-40
thin<-20

#sd2=1 for velocity, sd2=2 for acceleration
sd2=1

#These are norder and lambda values chosen for the analysis
if(sd2==1){norder=6;lambdaFDA<-c(rep(-15,3),rep(-14,3))}
if(sd2==2){norder=6;lambdaFDA=c(rep(-20.5,3),rep(-19,3))}


#L and M, set to 5 for faster computation (L=M=10 used for paper)
L=5
M=5

#L.FAR - should be small since no smoothing penalty done in FAR functions
L.FAR=5



#P is between 1 and 6 for different files, 1:3 fixed, 4:6 random
for(P in 1){
  data.prep=data.process(movement=movement,
                           set=P,
                           delta=delta,
                           sd=FALSE,
                           thin=thin,
                           L=L,
                           M=M,sd1=0,sd2=sd2,
                           nbasis=500,norder=norder,lambdaFDA=lambdaFDA[P])  
  attach(data.prep)
  
  #######################
  ##### SAFE gLASSO #####
  #######################
  if("SAFE" %in% Analysis){
  #Adapative Weights
  phi1<-exp(seq(-10,10,by=5))
  phi2<-phi1
  
  #Initialize weights to be 1
  WT<-list(); WT[[1]]=rep(1,16); WT[[2]]=rep(1,16); WT[[3]]=rep(1,16)
  
  #SAFEgLASSO!
  lseq<-exp(seq(-15,-2,by=0.2))
  system.time(
    SAFE <- SAFEglasso(X=Xtilde,Y=Y,
                     Omega1=Omega1,
                     Omega2=Omega2,
                     phi1=phi1,
                     phi2=phi2,
                     WT=WT,
                     nblocks=nblocks,
                     lseq=lseq,
                     nstage=nstage)
  )
  
  eval(parse(text=paste("SAFE.",P,"<-SAFE",sep="")))
  do.call(save,args=list(as.name(paste("SAFE.",P,sep="")),
                         file=paste("SAFE_",movement,P,".RData",sep=""))
  )
  }
  
  #######################
  ##### agLASSO #########
  #######################
  
  #Model includes a positional effect alpha(zi) using the scaled B2 matrix
  if("ag" %in% Analysis){
  par(mfrow=c(1,1))
  #Adapative Weights
  phi1<-exp(seq(-20,10,by=2))
  
  WT<-list(); WT[[1]]=rep(1,16); WT[[2]]=WT[[1]]; WT[[3]]=rep(0,16); 
  WT[[4]]=1; WT[[5]]=1 #WT[[4]]=0 means no sparsity penalty for position
  
  #agLASSO!
  lseq<-exp(seq(-25,-2,by=0.25))
  
  system.time(
    ag <- aglasso(X=Xomega,
                       Y=Y,
                       Omega1=Omega1,
                       Omega2=Omega2,
                       phi1=phi1,
                       WT=WT,
                       nblocks=nblocks,
                       lseq=lseq,
                       nstage=nstage,
                       pos.effect=scale(B2))
  )
  
  eval(parse(text=paste("ag.",P,"<-ag",sep="")))
  do.call(save,args=list(as.name(paste("ag.",P,sep="")),
                         file=paste("ag_",movement,P,".RData",sep=""))
  )
  }
  
  #######################
  ##### LADgLASSO #####
  #######################
  
  #Model includes a positional effect alpha(zi) using the scaled B2 matrix
  
  if("LAD" %in% Analysis){
  par(mfrow=c(1,1))
  #Adapative Weights
  phi1<-exp(seq(-16,10,by=2))
  
  WT<-list(); WT[[1]]=rep(1,16); WT[[2]]=WT[[1]]; WT[[3]]=rep(0,16); WT[[4]]=1; WT[[5]]=1
  
  #LADgLASSO!
  lseq<-exp(seq(-20,-2,by=0.2))
  system.time(
    LAD <- LADglasso(X=Xomega,Y=Y,
                  Omega1=Omega1,
                  Omega2=Omega2,
                  phi1=phi1,
                  WT=WT,
                  nblocks=nblocks,
                  lseq=lseq,
                  nstage=nstage,
                  pos.effect=scale(B2)
    )
  )
  
  eval(parse(text=paste("LAD.",P,"<-LAD",sep="")))
  do.call(save,args=list(as.name(paste("LAD.",P,sep="")),
                         file=paste("LAD_",movement,P,".RData",sep=""))
  )
  }
  
  ########################
  #### FARgSCAD 
  ########################
  
  #Model includes a positional effect alpha(zi) using the scaled B2 matrix
  
  if("FAR" %in% Analysis){
    print("FAR")
  FAR=FARgSCAD(X=X,Y=Y,L=L.FAR,delta,nblocks=nblocks,alpha=1,lam.path=exp(seq(-5,3,by=0.1)),pos.effect=scale(B2),cvmult=1)
  eval(parse(text=paste("FAR.",P,"<-FAR",sep="")))
  do.call(save,args=list(as.name(paste("FAR.",P,sep="")),
                         file=paste("FAR_",movement,P,".RData",sep=""))
  )
  }
  
  detach(data.prep)

}

