#################################################################
# THE MAIN PROGRAM TO RUN THE BAYESIAN SEGMENTATION ALGORITHM.
# 
# FOR DUPLICATIONS ONLY. 
#  
# NOTE THAT THE PRIORS WERE SETUP FOR THE WGTP DATA. FOR DIFFERENT 
# APPLICATION, PLEASE CHANGE THE PRIOR VALUES.
#
##############################################################


WGTP.DUP<-function(x,y,stop.value=0.1){
   x.train<-x
   y.train<-y
   invx2.prior<-set.invx2.prior(x,y)
   last.box.avg<-Inf
   beta.num=1
   peel.obj<-ch.peel(x.train,y.train,beta.num=beta.num,last.num=1)
   last.box.avg<-peel.obj$box.avg[length(peel.obj$box.avg)]
   iter<-0
   u.bg<-0
   if(last.box.avg>stop.value){
      
   
   while(last.box.avg>stop.value & length(x.train)>1 ){
         iter<-iter+1
         eta.tmp<-set.eta.a(peel.obj,y.train)
         a<-eta.tmp$a
         cratio<-1/a 
   
         ###################################################################
         # PLEASE CHANGE THE FOLLOWING PRIOR VALUES WHEN ANALYZING
         # A NEW DATA SET. 
         #
         #

         max.box25<-call.bf(peel.obj=peel.obj,
                            invx2.prior.obj=invx2.prior,
                            cratio=cratio,
                            u.box=0.25,
                            u.bg=u.bg,
                            x=x.train,
                            y=y.train)
                                
         max.box50<-call.bf(peel.obj=peel.obj,
                            invx2.prior.obj=invx2.prior,
                             cratio=cratio,
                             u.box=0.5,
                             u.bg=u.bg,
                             x=x.train,
                             y=y.train)
                  
         max.box75<-call.bf(peel.obj=peel.obj,
                             invx2.prior.obj=invx2.prior,
                             cratio=cratio,
                             u.box=0.75,
                             u.bg=u.bg,
                             x=x.train,
                             y=y.train) 
         tmp.max.box25<-bottom.up(peel.obj,
                                max.obj=max.box25,
                                invx2.prior.obj=invx2.prior,
                                cratio=cratio,
                                u.box=0.25,
                                u.bg=u.bg,
                                x.dat=x.train,
                                y.dat=y.train)
                                  
                                  
         tmp.max.box50<-bottom.up(peel.obj,
                                max.obj=max.box50,
                                invx2.prior.obj=invx2.prior,
                                cratio=cratio,
                                u.box=0.5,
                                u.bg=u.bg,
                                x.dat=x.train,
                                y.dat=y.train)
                                                
          tmp.max.box75<-bottom.up(peel.obj,
                                 max.obj=max.box75,
                                 invx2.prior.obj=invx2.prior,
                                 cratio=cratio,
                                 u.box=0.75,
                                 u.bg=u.bg,
                                 x.dat=x.train,
                                 y.dat=y.train)
                                                  
           mg.lik25<-tmp.max.box25$max.box.lik
           mg.lik50<-tmp.max.box50$max.box.lik
           mg.lik75<-tmp.max.box75$max.box.lik
           max.lik<-max(mg.lik25,mg.lik50,mg.lik75)
           if(max.lik==mg.lik25){            
                   max.box<-tmp.max.box25
           }        
           else if(max.lik==mg.lik50){
                   max.box<-tmp.max.box50   
           }
           else{
                   max.box<-tmp.max.box75    
           }
             
            
            x.up<-max.box$max.box.def[2]
            x.low<-max.box$max.box.def[1]
            in.box<-x.train>=x.low & x.train<=x.up
            x.train<-x.train[!in.box]
            y.train<-y.train[!in.box]
           
            sol.tmp<-data.frame(max.box$max.box.def[1],
                                max.box$max.box.def[2], 
                                max.box$max.box.obs,
                                max.box$max.box.avg,
                                max.box$max.box.var, 
                                max.box$max.bf,
                                eta=max.box$max.box.eta, 
                                eta0=u.bg,
                                a=a,
                                cratio=cratio)
            if(iter==1){
                   sol<-sol.tmp
            }
            else{
                   sol<-rbind(sol,sol.tmp)
            }
           if(length(x.train)>1){     
              peel.obj<-ch.peel(x.train,y.train,beta.num=beta.num,last.num=1)
              last.box.avg<-peel.obj$box.avg[length(peel.obj$box.avg)]
              
           }
           else{
              cat("----------------- EMPTY","\n")
              break
           
           }
         
           
      }
     names(sol)=c("Start","End","Num_of_obs","Average","Variance","Bayes_Factor","Eta","Eta0","a","c_ratio")
     ssol<-sol[order(sol$"Bayes_Factor",decreasing=T),]
     ssol<-round(ssol,4)
   }
   else{
      ssol<-data.frame(0,
                       0, 
                       0,
                       0,
                       0, 
                       0,
                       0, 
                       0,
                       0,
                       0)
     names(ssol)=c("Start","End","Num_of_obs","Average","Variance","Bayes_Factor","Eta","Eta0","a","c_ratio")
   
     
  }   
     
     return(ssol)
}