#
#
#                   S - procedures used with SuperGEM(tm) 1.0
#
#                                  (1/7/00)
#
# Coded any copyright (2000) by Jerome H. Friedman.
#
#
supgem_function (x,y=x[,p],lx=rep(1,p),nbump=10,nl=2*n/3,
   header=paste("",date()),wt=rep(1,n),alpha=0.1,
   xmiss=9.0e30,nbin=5,iseed=987654321,mtf=10,icr=3,ibu=1, 
   beta0=0.001,thr=1.0e-4,sumfile="rules.sum") {
   span_0
   if(mode(x)!="numeric") stop(" x must be numeric.")
   if(is.matrix(x) || is.data.frame(x)) { n <- nrow(x); p <- ncol(x)}
   else if(is.vector(x)) { n <- length(x); p <- 1}
   else stop(" x must be a data frame, matrix or vector.")
   if(!is.vector(y,mode="numeric")) stop (" y must be a numeric vector.")
   if(length(y)!=n) stop (" x and y are of unequal length.")
   if(length(wt)!=n) stop (" x and wt of unequal length.")
   if(!is.vector(wt,mode="numeric")) stop(" wt must be a numeric vector.")
   if(any(is.na(x))) { x(is.na(x)) <- xmiss;
       warning(" x contains NA's - xmiss substituted.")
   } 
   if(any(is.na(y))) { wt(is.na(y)) <- 0;
      warning(" y contains NA's - corresponding weights set to 0.")
   } 
   if(any(is.na(wt))) { wt(is.na(wt)) <- 0;
      warning(" wt contains NA's - zeros substituted.")
   }
   if(any(wt<0)) { wt(wt<0) <- 0;
      warning(" wt contains negative numbers - zeros substituted.")
   }
   if(missing(y)) {
      if(p==1) stop("ncol(x)=1 and y is missing.")
      else lx[p] <- 0
   }
   length(lx) <- p
   if(any(lx < 0 | lx > 2)) { lx(lx < 0 | lx > 2) <- 0;
      warning(" invalid values in lx vector - zeros substituted.")
   }
   nbump <- parchk("nbump",nbump,1,1000,10)
   alpha <- parchk("alpha",alpha,0.001,0.999,0.1)
   if(!missing(xmiss))
      xmiss <- parchk("xmiss",xmiss,2*max(x[x!=xmiss]),9.9e35,9.0e30)
   nbin <- parchk("nbin",nbin,2,1000,5)
   iseed <- parchk("iseed",iseed,1,10e10,987654321)
   span <- parchk("span",span,0,1,0.2)
   mtf <- parchk("mtf",mtf,1,1000,10)
   icr <- parchk("icr",icr,1,3,3)
   ibu <- parchk("ibu",ibu,0,1,1)
   beta0 <- parchk("beta0",beta0,0,1,0.001)
   thr <- parchk("thr",thr,0,1,1.0e-4)
   cat(n,p,nl,nbump,alpha,xmiss,nbin,iseed,span,mtf,icr,ibu,beta0,thr,"\n",
      file=paste(gemhome,"parms.gem",sep="/"))
   tfile_paste(gemhome,"data.gem",sep="/")
   write(lx,file=tfile)
   write(x,file=tfile,append=T)
   write(y,file=tfile,append=T)
   write(wt,file=tfile,append=T)
   names <- dimnames(x)[[2]]
   tfile_paste(gemhome,"vnames.gem",sep="/")
   if(!is.null(names)) { cat(p,"\n",file=tfile)
      for (i in 1:p) { cat(names[i],"             ","\n",file=tfile,append=T)}
   }
   else { cat("0","\n",file=tfile)}
   if (is.character(header)) {
      cat(header,"\n",file=paste(gemhome,"header.gem",sep="/"))
   }
   else { cat("  ","\n",file=paste(gemhome,"header.gem",sep="/"))}
   if (is.character(sumfile)) {
      cat(sumfile,"\n",file=paste(gemhome,"sfile.gem",sep="/"))
   }
   else { cat("rules.sum","\n",file=paste(gemhome,"sfile.gem",sep="/"))}
   if (gemhost == "unix") {
      cmd_paste("xterm -sb -sl 200 -e ",gemhome,"/SuperGEM ",
         gemhome," &",sep="")
      unix(cmd,output=F)
   }
   else if (gemhost=="wintel") {
      if (substring(gemhome,2,2)!=":") { dsk_""; dir_gemhome}
      else { dsk_substring(gemhome,1,2); dir_substring(gemhome,3)}
      dos(paste(dsk,dir,"/SuperGEM ",dsk," ",dir,sep=""),
         output.to.S=F,multi=T,translate=T)
   }
   else { stop("gemhost must be set to either unix or wintel.")}   
   invisible()
}
parchk <- function (ax,x,lx,ux,df) {
   if(x < lx || x > ux) { df
      warning(paste(" invalid value for",ax,"- default (",df,") used."))
   }
   else { x}
}
rules <- function(sumfile="rules.sum") {
   if (!is.character(sumfile)) stop("sumfile must be of type character.")
   if (gemhost=="unix") { unix(paste("cat ",gemhome,"/",sumfile,sep=""))}
   else if (gemhost=="wintel") {
      if (substring(gemhome,2,2)!=":") { dsk_""; dir_gemhome}
      else { dsk_substring(gemhome,1,2); dir_substring(gemhome,3)}
      dos(paste("type ",dsk,dir,"/",sumfile,sep=""),
         translate=T)
   }
   else { stop("gemhost must be set to either unix or wintel.")}
}
see <- function(argrule) { cat(argrule,sep="\n"); invisible()} 
vnames <- function(x,names) {
   if(missing(x)) stop("The matrix or frame x must be specified.")
   if(missing(names)) stop("(variable) names  must be specified.")
   if(mode(x)!="numeric") stop(" x must be numeric.")
   if(is.matrix(x) || is.data.frame(x)) { p <- ncol(x)}
   else if(is.vector(x)) { p <- 1}
   else { stop(" x must be a data frame, matrix or vector.")}
   if (!is.character(names)) stop("names must be of type character.")
   if (length(names)!=p) stop("names must be specified for all the variables.")
   list(NULL,names)
}
gemhelp <- function(fun="gemhelp") {
   if (!is.character(fun)) stop("fun must be of type character.")
   if (gemhost=="unix") {
       cat(unix(paste("cat ",gemhome,"/",fun,".hlp",sep="")),sep="\n")
   }
   else if (gemhost=="wintel") {
      if (substring(gemhome,2,2)!=":") { dsk_""; dir_gemhome}
      else { dsk_substring(gemhome,1,2); dir_substring(gemhome,3)}
      cat(dos(paste("type ",dsk,dir,"/",fun,".hlp",sep=""),
         translate=T),sep="\n")
   }
   else { stop("gemhost must be set to either unix or wintel.")}
   invisible()
}
traj_function () {
   par(mfrow=c(1,1))
   top_scan(paste(gemhome,"trajtle.spl",sep="/"),what=character(0),sep=";")
   plts_matrix(scan(paste(gemhome,"traj.spl",sep="/")),ncol=2)
   plot(plts[,1],plts[,2],xlim=c(0,1),xlab='Support',ylab='Box mean',main=top)
   identify(plts[,1],plts[,2])
   invisible()
}
relfreq_function (plt=1,nc=4) {
   title_scan(paste(gemhome,"varnms.spl",sep="/"),what=character(0),sep=";")
   nv_scan(paste(gemhome,"nvals.spl",sep="/"))
   l_length(nv)
   ir_trunc(l/nc)
   if (nc*ir!=l) {ir_ir+1}
   par(mfrow=c(ir,nc))
   if (plt!=1) {y_scan(paste(gemhome,"vals.spl",sep="/"))}
   else {y_scan(paste(gemhome,"ratios.spl",sep="/"))}
   k_1
   for (i in 1:l) {
      if (nv[i] > 0) {
         m_k+nv[i]-1
         if (plt!=1) {
            m_m-1
            bot_"bin lower limit"
         }
         else { bot_"rel. freq."}
         if (i!=l) {top_title[i]}
         else { top_"output"}
         barplot(y[k:m],main=top,ylab=bot)
         k_k+nv[i]
      }
   }
   invisible()
}
sens_function (pltvar="",miss=T,span=0.3,nc=2,ln=F,qtl=0.1,minsmu=11) {
   title_scan(paste(gemhome,"varnms.spl",sep="/"),what=character(0),sep=";")
   ifo_matrix(scan(paste(gemhome,"sensinfo.spl",sep="/")),ncol=3)
   plts_matrix(scan(paste(gemhome,"sensplot.spl",sep="/")),ncol=3)
   def_scan(paste(gemhome,"boxinfo.spl",sep="/"))
   box_matrix(scan(paste(gemhome,"boxdef.spl",sep="/")),nrow=3)
   np_nrow(ifo)
   if (is.character(pltvar)==F || pltvar!="") { par(mfrow=c(1,1))}
   else {
      ir_trunc(np/nc)
      if (nc*ir!=np) {ir_ir+1}
      par(mfrow=c(ir,nc))
   }
   pctl_1.0/qtl
   ymin_9.9e30; ymax_-ymin
   l_1
   for (i in 1:np) {
      if (ifo[i,3]==1) {
         j_l+ifo[i,1]-1
         jm1_j-1
         if (j-l >= minsmu) {
            smu_loess.smooth(plts[l:jm1,1],plts[l:jm1,2],span=span)
            y_smu$y
         }
         else { y_plts[l:jm1,2]}
         if (miss==T) { y_c(y,plts[j,2])}
         ymax_max(ymax,y); ymin_min(ymin,y)
      }
      l_l+ifo[i,1]
   }
   l_1
   for (i in 1:np) {
      if (is.character(pltvar)==T) {
         if (pltvar!="") {
            tplt_pltvar
            for (ii in 1:(10-nchar(pltvar))) { tplt_paste(tplt," ",sep="")}
            if (tplt!=title[ifo[i,2]]) { l_l+ifo[i,1]; next}
         }
      }
      else if (substring(title[ifo[i,2]],2,2) == "(") {
         if (pltvar!=as.integer(substring(title[ifo[i,2]],3,5))) {
	    l_l+ifo[i,1]; next
         }
      }
      else stop("Please use character name for this input variable.")
      j_l+ifo[i,1]-1
      if (ifo[i,3]!=1) {
         cl_rep(0,ifo[i,1])
         nme_rep(" ",ifo[i,1])
         num_rep(" ",ifo[i,1])
         for (m in 1:ifo[i,1]) {
            f_plts[m+l-1,3]
            if (m > 1 && m < ifo[i,1]) { f_f-plts[m+l-2,3]}
            num[m]_as.character(trunc(100*f+0.5))
            if (m < ifo[i,1]) {
               nme[m]_as.character(plts[m+l-1,1])
            }
            else { nme[m]_"M"}  
            for (k in 1:def[2]) {
               if (box[1,k]==ifo[i,2]) {
                  if (box[2,k] == plts[m+l-1,1]) { cl[m]_-1}
               }
            }
         }
         where_barplot(plts[l:j,2],main=title[ifo[i,2]],
                  density=cl,names=nme,ylab='output')
         lw_length(where)
	 linc_(where[lw]-where[1])*0.5/length(plts[l:j,2])
         lines(c(where[1]-linc,where[lw]+linc),c(def[1],def[1]))
         axis(side=1,at=where,line=1,labels=num,ticks=F)
      }
      else {
         jm1_j-1
         if (j-l >= minsmu) {
            smu_loess.smooth(plts[l:jm1,1],plts[l:jm1,2],span=span)
            x_smu$x; y_smu$y
         }
         else { x_plts[l:jm1,1]; y_plts[l:jm1,2]}
         if (miss==T) { x_c(x,plts[j,1]); y_c(y,plts[j,2])}
         if (ln==F) {
            plot(x,y,xlab=' ',ylim=c(ymin,ymax),
	       ylab='output',main=title[ifo[i,2]],type="n")
         }
         if (ln!=F) {
            plot(x,y,xlab=' ',ylim=c(ymin,ymax),
	       ylab='output',main=title[ifo[i,2]],log="x",type="n")
         }
         l1_length(x)
         if (miss==T) { l1_l1-1; points(plts[j:j,1],plts[j:j,2],pch="M")}
         if (j-l < minsmu) points(x[1:l1],y[1:l1])
         lines(x[1:l1],y[1:l1])
         ll_F
         ul_F
         for (k in 1:def[2]) {
            if (box[1,k]==ifo[i,2]) {
               if (box[3,k] > 0) {
                  ll_T
                  xl_box[2,k]
               }
               if (box[3,k] < 0) {
                  ul_T
                  xu_box[2,k]
               }
            }
         }
         if (ll==T) { lines(c(xl,xl),c(ymin,def[1]))}
         if (ul==T) { lines(c(xu,xu),c(ymin,def[1]))}
         if (ll==T & ul==T) { lines(c(xl,xu),c(def[1],def[1]))}
         if (ll==T & ul==F) { lines(c(xl,plts[jm1,1]),c(def[1],def[1]))}
         if (ll==F & ul==T) { lines(c(plts[l,1],xu),c(def[1],def[1]))}
         tle_(1:(pctl-1))/pctl
         eps_0.01*(plts[jm1,1]-plts[l,1])
         xl_plts[l,1]; st_0;  m_l; k_1
         repeat {
            while (plts[m,3] <= tle[k]) {
               xl_plts[m,1]; m_m+1; if (m > jm1) break;
            }
            if (m > jm1) break;
            if (plts[m,1] > xl) {
               h_(plts[m,3]-tle[k])/(plts[m,3]-plts[m-1,3])
               x[k]_h*xl+(1-h)*plts[m,1]; st_0
            }
            else { x[k]_xl+st*eps; st_st+1}
            k_k+1
            if (k > pctl-1) { k_k-1; break}
         }
         points(x[1:k],rep(ymin,k),pch="|")
         if (miss == T) { g_as.character(trunc(10*plts[j,3]+0.5))
            points(plts[j:j,1],rep(ymin,1),pch=g)
         }
      }
      l_l+ifo[i,1]
   }
   invisible()
}