.packageName <- "SHARE"
#line 1 "/private/tmp/Rtmpl8ol26/R.INSTALL10d63af1/SHARE/R/S4-class.R"
## Code:        S4-class.R
## Purpose:     Define general methods
## Revision
##      2008-12-03
##              - define new class: genoSet & haploSet
##      2008-12-12
##              - remove pheno data from haplo-class
##      2008-12-23
##              - rownames() instead of rowname()
##      2009-03-10
##              - combine "finalsubset" class into "share" class
##              - move class-related definitions to separate code


##---Generic Function of General Methods ---##

##--- for SNP (loci) info ---##
## nSNP
setGeneric("nSNP",
           function(.Object){
               standardGeneric("nSNP")
           }
           )

## nameSNP
setGeneric("nameSNP",
           function(.Object){
               standardGeneric("nameSNP")
           }
           )


##--- for Sequence (genotype of haplotype) info ---##
## nSeq
setGeneric("nSeq",
           function(.Object){
               standardGeneric("nSeq")
           }
           )

## nameSeq
setGeneric("nameSeq",
           function(.Object){
               standardGeneric("nameSeq")
           }
           )

#line 1 "/private/tmp/Rtmpl8ol26/R.INSTALL10d63af1/SHARE/R/genoSet-class.R"
## Code:        genoSet-class.R
## Purpose:     definition of genoSet class and associated methods
##      
## Revision
##      2009-03-10
##              - move from S4-class.R & haplo.R

setClass("genoSet",
         representation(genoSeq="data.frame",       ## k genotype sequences (0-1-2 format) * n SNPs
                        phenoData="data.frame"  ## k genotype sequences * p variables (phenotype)
                        )
         )

## show
setMethod("show", signature=signature(object="genoSet"),
          function(object){
              cat(paste("An object of class \"genoSet\"",
                        "\n# of SNP:", nSNP(object),
                        "\n# of genotype sequence: ", nSeq(object),
                        "\nPhenotype:\n", paste(colnames(object@phenoData), sep=", ", collapse=""),
                        "\n"
                        )
                  )
          }
          )

## accessor
setGeneric("genoSeq",
           function(.Object){
               standardGeneric("genoSeq")
           }
           )
setMethod("genoSeq", signature=signature(.Object="genoSet"),
          function(.Object){
              .Object@genoSeq
          }
          )

setGeneric("phenoData",
           function(.Object){
               standardGeneric("phenoData")
           }
           )
setMethod("phenoData", signature=signature(.Object="genoSet"),
          function(.Object){
              .Object@phenoData
          }
          )

## general methods
setMethod("nSNP", signature=signature(.Object="genoSet"),
          function(.Object){
              ncol(.Object@genoSeq)
          }
          )

setMethod("nameSNP", signature=signature(.Object="genoSet"),
          function(.Object){
              colnames(.Object@genoSeq)
          }
          )

setMethod("nSeq", signature=signature(.Object="genoSet"),
          function(.Object){
              nrow(.Object@genoSeq)
          }
          )

setMethod("nameSeq", signature=signature(.Object="genoSet"),
          function(.Object){
              rownames(.Object@genoSeq)
          }
          )

## specific methods
setGeneric("haplo",
           function(.Object, ...){
               standardGeneric("haplo")
           }
           )

setMethod("haplo", signature=signature(.Object="genoSet"),
          function(.Object, wgt=NULL, ...){
              #### check if wgt is a vector of length nrow(dat)
              if(is.null(wgt)){
                  wgt <- rep(1, nrow(.Object@genoSeq))
              }else{
                  wgtTry <- try(is.vector(wgt))
                  if(inherits(wgtTry, "try-error")){
                      stop("The argument wgt must be vector or NULL.")
                  }
                  if(length(wgt)!=nrow(.Object@genoSeq)){
                      stop("The length of argument wgt must equal to nrow(dat).")
                  }
              }

              ## randomize the order
              set.seed(37273)
              newOrder <- sample(1:nrow(.Object@genoSeq))
              .Object@genoSeq <- .Object@genoSeq[newOrder, ]
              if(ncol(.Object@phenoData)==1){
                  newPheno <- data.frame(.Object@phenoData[newOrder, ])
                  colnames(newPheno) <- colnames(.Object@phenoData)
                  .Object@phenoData <- newPheno
              }else{
                  .Object@phenoData <- data.frame(.Object@phenoData)[newOrder, ]
              }
              
              require(haplo.stats)
              ## construct genotype data.frame
              geno <- geno1to2(genoSeq(.Object), locus.label=nameSNP(.Object))
              
              ## haplotype reconstruction & frequency estimation
              em <- haplo.em(geno=geno, locus.label=nameSNP(.Object), 
                             miss.val=NA, weight=wgt,
                             control=haplo.em.control(loci.insert.order=NULL, insert.batch.size=3, min.posterior=0.01),
                             ...
                             )

              ## return the required values to a haplo-class object
              return(new("haplo",
                         haploSeq=em$haplotype,
                         haploFreq=em$hap.prob,
                         hap1=em$hap1code,
                         hap2=em$hap2code,
                         poolHapPair=em$subj.id,
                         nPosHapPair=em$nreps,
                         post=em$post,
                         pheno=.Object@phenoData
                         )
                     )              

          }
          )

#line 1 "/private/tmp/Rtmpl8ol26/R.INSTALL10d63af1/SHARE/R/haploSet-class.R"
setClass("haploSet",
         representation(haploSeq="data.frame",      ## m haplotype sequences * n SNPs
                        haploFreq="vector"      ## m haplotype sequences * 1
                        )
         )

## show
setMethod("show", signature=signature(object="haploSet"),
          function(object){
              cat(paste("An object of class \"haploSet\"",
                        "\n# of SNP:", nSNP(object),
                        "\n# of haplotype sequence: ", nSeq(object),
                        "\n"
                        )
                  )
          }
          )

## accessor
setGeneric("haploSeq",
           function(.Object){
               standardGeneric("haploSeq")
           }
           )
setMethod("haploSeq", signature=signature(.Object="haploSet"),
          function(.Object){
              .Object@haploSeq
          }
          )

setGeneric("haploFreq",
           function(.Object){
               standardGeneric("haploFreq")
           }
           )
setMethod("haploFreq", signature=signature(.Object="haploSet"),
          function(.Object){
              .Object@haploFreq
          }
          )

## general methods
setMethod("nSNP", signature=signature(.Object="haploSet"),
          function(.Object){
              ncol(.Object@haploSeq)
          }
          )

setMethod("nameSNP", signature=signature(.Object="haploSet"),
          function(.Object){
              colnames(.Object@haploSeq)
          }
          )

setMethod("nSeq", signature=signature(.Object="haploSet"),
          function(.Object){
              nrow(.Object@haploSeq)
          }
          )

setMethod("nameSeq", signature=signature(.Object="haploSet"),
          function(.Object){
              rownames(.Object@haploSeq)
          }
          )

#line 1 "/private/tmp/Rtmpl8ol26/R.INSTALL10d63af1/SHARE/R/haplo-class.R"
setClass("haplo",
         contains="haploSet",         
         representation(##hap == haploSeq @ haploSet class
                        ##hap="data.frame",       ## unique haplotype X loci (haplo.em$haplotype)
                        ## hapProb == haploFreq @ haploSet class
                        ##hapProb="vector",       ## haplotype frequencies for each unique haplotype (haplo.em$hap.prob)
                        hap1="vector",          ## index of first haplotype of each subject (haplo.em$hap1code)
                        hap2="vector",          ## index of second haplotype of each subject (haplo.em$hap2code)

                        poolHapPair="vector",   ## pooled haplotype ID (haplo.em$subj.id)
                        nPosHapPair="table",    ## number of possible haplotype pairs (haplo.em$nreps)

                        post="vector",          ## posterior probability for the pooled haplotype pairs
                        
                        pheno="data.frame"      ## pheno type data after re-ordering
                        )
         )
## show
setMethod("show", signature=signature(object="haplo"),
          function(object){
              cat(paste("An object of class \"haplo\" (extended from class \"haploSet\")",
                        "\n# of SNP:", nSNP(object),
                        "\n# of haplotype sequence: ", nrow(object@haploSeq),
                        "\n"
                        )
                  )
          }
          )
#line 1 "/private/tmp/Rtmpl8ol26/R.INSTALL10d63af1/SHARE/R/share-class.R"
## Revision
##      2009-03-30
##              - add new slot "inherit" to store the value of "Minherit"

setClass("share",
         representation(uhap="vector",
                        weight="vector",
                        nFold="numeric",
                        maxSNP="numeric",
                        deviance="vector",
                        bestsize="numeric",
                        finalHap="data.frame",
                        finalHapFreq="vector",
                        finalHapTest="data.frame",
                        globalP="numeric",
                        modelmethod="character",
                        inherit="character"
                        ),
         contains="haplo"
         )


##FIXME:  works, but mask graphics::plot
##        need to use namespace to prevent the conflict with graphics::plot
##
##setGeneric("plot",
##           function(.Object, ...){
##               standardGeneric("plot")
##               }
##           )
##
##setMethod("plot", signature("share"),
##          function(.Object, ...){
##              graphics::plot(.Object@deviance,
##                   type="b", xaxt="n",
##                   xlab="Model Size", ylab="Deviance")
##              axis(side=1, at=1:(length(.Object@deviance)),
##                   labels=as.character(c(0, 1:.Object@maxSNP, (.Object@maxSNP-1):1 )))
##              }
##          )

## show
setMethod("show", signature=signature(object="share"),
          function(object){

              if (object@modelmethod=="Cross-Val") {
                  cat(paste("An object of class \"share\"",
                            "\nGiven nFold=", object@nFold, " cross-validation", " and maxSNP=", object@maxSNP, 
                            ", the best number of SNP to select is ", object@bestsize, ".",
                            "\n",
                            sep=""
                            )
                      )

                  if (object@bestsize>0) {
                      testDF <- object@finalHapTest
                      for(col in colnames(testDF)){
                          testDF[, col] <- round(testDF[, col], 6)
                      }
                      ##cat("An object of class \"finalsubset\"")
                      cat(paste("\n# of selected SNP:", ncol(object@finalHap),
                                "\n# of haplotype sequence: ", nrow(object@finalHap),
                                "\n"
                                )
                          )
                      cat("\nHaplotype sequences in the final subset:\n")
                      print(cbind(object@finalHap, Frequency=round(object@finalHapFreq, 6)))
                      cat("\nTesting:\n")
                      print(testDF)
                      cat(paste("\nGlobal p-value:", round(object@globalP, 6), "\n"))
                  }
              } ## end of if (object@modelmethod=="Cross-Val") 
              else if (object@modelmethod=="BIC") {
                  cat(paste("An object of class \"share\"",
                            "\nGiven BIC criterion", " and maxSNP=", object@maxSNP, 
                            ", the best number of SNP to select is ", object@bestsize, ".",
                            "\n",
                            sep=""
                            )
                      )

                  if (object@bestsize>0) {
                      testDF <- object@finalHapTest
                      for(col in colnames(testDF)){
                          testDF[, col] <- round(testDF[, col], 6)
                      }
                      ##cat("An object of class \"finalsubset\"")
                      cat(paste("\n# of selected SNP:", ncol(object@finalHap),
                                "\n# of haplotype sequence: ", nrow(object@finalHap),
                                "\n"
                                )
                          )
                      cat("\nHaplotype sequences in the final subset:\n")
                      print(cbind(object@finalHap, Frequency=round(object@finalHapFreq, 6)))
                      cat("\nTesting:\n")
                      print(testDF)
                      cat(paste("\nGlobal p-value:", round(object@globalP, 6), "\n"))
                  }
              } ## end of if (object@modelmethod=="BIC")
          }
          )       

setGeneric("dplot",
           function(shareObj, ...){
               standardGeneric("dplot")
           })

setMethod("dplot", signature=c("share"),
          function(shareObj, ...){
              if (shareObj@modelmethod=="Cross-Val") {   
                  plot(shareObj@deviance,
                       type="b", xaxt="n",
                       xlab="Model Size", ylab="Prediction Deviance", ...)
                  axis(side=1, at=1:(length(shareObj@deviance)),
                       labels=as.character(c(0, 1:shareObj@maxSNP, (shareObj@maxSNP-1):1 )))
              } else if (shareObj@modelmethod=="BIC") {
                  plot(shareObj@deviance,
                       type="b", xaxt="n",
                       xlab="Model Size", ylab="BIC", ...)
                  axis(side=1, at=1:(length(shareObj@deviance)),
                       labels=as.character(c(0, 1:shareObj@maxSNP, (shareObj@maxSNP-1):1 )))
              }
          }
          )

## FIXME:  re-declare summary will lost the old definition
##setGeneric("summary",
##           function(.Object, ...)
##           standardGeneric("summary")
##           )
##
##
##setMethod("summary", "finalsubset",
##          function(.Object){
##              glmout <- data.frame(round(.Object@finalHapCoef, 6),
##                                   c(round(.Object@expPValue, 6), rep("", length(.Object@finalHapCoef)-1))
##                                   )
##              colnames(glmout) <- c("Coefficient", "p-Value")
##              rownames(glmout) <- rownames(.Object@finalHap)
##              return(glmout)
##              }
##          )
#line 1 "/private/tmp/Rtmpl8ol26/R.INSTALL10d63af1/SHARE/R/finalsubset.R"
## Revision
##      2008-12-03
##              - calculate the z-score = coeff / s.e., where s.e. is at the var-covar matrix
##      2008-12-12
##              - conditionally output the result by the value of basesize
##      2008-12-23
##              - not using stop(); otherwise it will faile the creation of vignette
##              - add more info in the warning message
##              - add 'PACKAGE="SHARE"' into .C()
##      2009-03-30
##              - change "STE" to "Standard Error"
##              - add (Intercept) to hap1

finalsubset <- function(shareObj,
                        phenoData,
                        Minherit,
                        tol=1e-8,
                        verbose=0,
                        phase=0  
                        ){

    ## if not significant, report NULL for most of slots
    if (shareObj@bestsize==0) {
        ##pval<-0.051
        cat("Warning:\nIt is unable to determine the final subset if the best size is 0.\nNo result was outputted.\n")
    }else{ ## otherwise, find the final subset
        if (shareObj@bestsize > shareObj@maxSNP){
            outputsnp<-rep(0, 2 * shareObj@maxSNP - shareObj@bestsize)
        }else{
            outputsnp<-rep(0, shareObj@bestsize)
        }

        ## FIXME:  need to define the 1000 in the argument
        varstore <- rep(0,1000)
        coef <- rep(0,100)
        outnhaps=0
        outhapvec <- rep(0,1000)
        outhapfreq<- rep(0,100)

        temp1 <- .C("finalsubset",
                    id=as.integer(shareObj@poolHapPair),
                    nSubj=as.integer(length(shareObj@nPosHapPair)),
                    nObs=as.integer(length(shareObj@poolHapPair)),
                    nPair=as.integer(shareObj@nPosHapPair),
                    ccStatus=as.integer(rep(phenoData, shareObj@nPosHapPair)),
                    nLoci=as.integer(ncol(shareObj@haploSeq)),
                    nHap=as.integer(length(shareObj@haploFreq)),
                    hap1=as.integer(shareObj@hap1),
                    hap2=as.integer(shareObj@hap2),

                    isUniHap=as.integer(shareObj@uhap),
                    hapFreq=as.numeric(shareObj@haploFreq),
                    postHap=as.numeric(shareObj@post),

                    maxSNP=as.integer(shareObj@maxSNP),
                    
                    as.integer(shareObj@bestsize),
                    outputsnp=as.integer(outputsnp),
                    outhapvec=as.integer(outhapvec),
                    outhapfreq=as.numeric(outhapfreq),

                    as.numeric(tol),
                    as.integer(phase),
                    varstore=as.double(varstore),
                    coef=as.double(coef),
                    outnhaps=as.integer(outnhaps),
                    Minherit=as.integer(Minherit),
                    as.integer(verbose),
                    PACKAGE="SHARE"
                    )

        ## var-covar matrix
        amat <- matrix(temp1$varstore[1:(temp1$outnhaps ^2)], nrow=temp1$outnhaps,
                       ncol=temp1$outnhaps, byrow=TRUE)
        bmat <- matrix(temp1$varstore[((temp1$outnhaps ^2)+1):(2*(temp1$outnhaps ^2))],
                       nrow=temp1$outnhaps, ncol=temp1$outnhaps, byrow=TRUE)

        require(MASS) ## for ginv funciton
        var<-ginv(amat)%*% bmat %*% ginv(amat)
        
        pval <- drop(1-pchisq(t(temp1$coef[2:temp1$outnhaps]) %*% ginv(var[2:temp1$outnhaps,2:temp1$outnhaps]) %*% temp1$coef[2:temp1$outnhaps],temp1$outnhaps-1))

        ## name of haplotypes & SNPs
        snpNames <- colnames(shareObj@haploSeq)[sort(temp1$outputsnp)]
        hapName <- paste("hap", 1:temp1$outnhaps, sep="")
        ## add (Intercept) to hap1
        ## hapName[1] <- paste(hapName[1], "(Intercept)")
        
        ## data frame for the output haplotypes
        hapDF <- data.frame(matrix(temp1$outhapvec[1:(length(temp1$outputsnp)*temp1$outnhaps)],
                                   nrow=temp1$outnhaps, ncol=length(temp1$outputsnp), byrow=TRUE)
                            )
        colnames(hapDF) <- snpNames
        rownames(hapDF) <- hapName

        ## haplotype frequency
        hapFreq <- temp1$outhapfreq[1:temp1$outnhaps]
        names(hapFreq) <- hapName

        ## coefficient
        coef <- temp1$coef[1:temp1$outnhaps]
        names(coef) <- hapName

        ## s.e.
        stdErr <- sqrt(diag(var))

        ## z-score
        zScore <- coef / stdErr

        ## p-values
        pValue <- 1-pnorm(abs(zScore))

        ## GLM-like output
        haploTest <- data.frame(coef, stdErr, zScore, pValue)
        colnames(haploTest) <- c("Coefficient", "Standard Error", "z-Score", "p-Value")
        hapbase <- which.max(temp1$outhapfreq[1:temp1$outnhaps])
        hapName1 <- paste("hap",c(hapbase,(1:temp1$outnhaps)[-hapbase]))
        hapName1[1] <- paste(hapName1[1], "(Intercept)")
        rownames(haploTest) <- hapName1

        outObj <- shareObj
        class(outObj) <- "share"
        outObj@finalHap <- hapDF
        outObj@finalHapFreq <- hapFreq
        outObj@finalHapTest <- haploTest
        outObj@globalP <- pval
        outObj@modelmethod <- "Cross-Val"
        outObj@inherit <- c("additive", "dominant", "recessive")[Minherit]

        return(outObj)
    }
    
}
#line 1 "/private/tmp/Rtmpl8ol26/R.INSTALL10d63af1/SHARE/R/cshare.R"
## Revision
##      2008-12-12
##              - add pheno data into argument
##              - change the sign of deviance
##              - add bestsize to the output
##      2008-12-23
##              - add 'PACKAGE="SHARE"' into .C()
##      2009-03-03
##              - James D: add two arguments: ModSelMethod & Minherit
##      2009-03-09
##              - change the acceptable value of ModSelMethod to "Cross-Val" and "BIC"
##              - change the acceptable value of Minherit to "additive", "dominant", and "recessive"
##      2009-03-30
##              - report standard error in the testing section
##              - add (Intercept) to hap1
##              - in shareTest, remove arguments "nfold", "maxsnps", "ModSelMethod", because they are included in outObj already
##              - compare global p-values only if bestsize != 0
##      2009-05-15
##              - make the default value of verbose=FALSE
##      2009-07-02
##              - use "status" argument to determine which clinical status in phenoData for analysis 

cshare <- function(haploObj,    ## haploObj
                   status,      ## label of status in the haploObj@pheno data.frame
                   nfold=10,    ## number of fold (numeric), NOT USED IF MODEL SELECTION IS THROUGH BIC
                   maxsnps,     ## maximum SNP selected in the model (numeric) 
                   tol=1e-8,
                   verbose=FALSE,   ## TRUE:  turn on log file
	           ModSelMethod=c("Cross-Val", "BIC"),    ## default model selection method is cross-val (=1), however can use BIC (=2) for faster computation
                   Minherit=c("additive", "dominant", "recessive")          ## the default mode of inheritance is "additive",  the other options is "dominant" and  "recessive" 
                   ){

    ## checking argument
    ModSelMethod <- match.arg(ModSelMethod)
    Minherit <- match.arg(Minherit)
    
    ## Validate arguments 
    ## check if dat is a haplo object
    if(class(haploObj)!="haplo"){
        stop("The argument dat must be object of the class 'haplo'.")
    }

    ## check if nfold is a numberic number
    nfoldTry <- try(is.numeric(nfold))
    if(inherits(nfoldTry, "try-error")){
        stop("The argument nfold must be numeric.")
    }
    if(length(nfold)!=1){
        stop("The length of argument nfold must be 1.")
    }

    ## check if maxsnp is a numberic number
    maxsnpsTry <- try(is.numeric(maxsnps))
    if(inherits(maxsnpsTry, "try-error")){
        stop("The argument maxsnps must be numeric.")
    }
    if(length(maxsnps)!=1){
        stop("The length of argument maxsnps must be 1.")
    }
    if(maxsnps>(ncol(haploObj@haploSeq)-1)){
        stop("The value of argument maxsnps must be less than the total number of SNPs.")
    }

    uhap <- rep(0, length(haploObj@haploFreq)* ncol(haploObj@haploSeq)) 
    count<-1
    for (i in 1:(length(haploObj@haploFreq))){
        for (j in 1:(ncol(haploObj@haploSeq))) {      
            uhap[count]<- as.numeric(haploObj@haploSeq[i,j])
            count<-count+1
        }
    }

    deviance <- rep(0,2*maxsnps)
    phase <- 0

    if (ModSelMethod=="Cross-Val") {
        temp <- .C("xshare",
                   id=as.integer(haploObj@poolHapPair),
                   nSubj=as.integer(length(haploObj@nPosHapPair)),
                   nObs=as.integer(length(haploObj@poolHapPair)),
                   nPair=as.integer(haploObj@nPosHapPair),
                   ccStatus=as.integer(rep(unlist(haploObj@pheno[, status]), haploObj@nPosHapPair)),
                   nLoci=as.integer(ncol(haploObj@haploSeq)),
                   nFold=as.integer(nfold),
                   nHap=as.integer(length(haploObj@haploFreq)),
                   hap1=as.integer(haploObj@hap1),
                   hap2=as.integer(haploObj@hap2),               
                   uhap=as.integer(uhap),
                   hapFreq=as.numeric(haploObj@haploFreq),
                   postHap=as.numeric(haploObj@post),
                   maxSNP=as.integer(maxsnps),
                   deviance=as.numeric(deviance),
                   converge=as.numeric(tol),
                   verbose=as.integer(sum(verbose)),
                   phase=as.integer(phase),
                   Minherit=match(Minherit, c("additive", "dominant", "recessive")),
                   PACKAGE="SHARE"
                   )

        outObj <- haploObj
        class(outObj) <- "share"
        outObj@uhap=uhap
        outObj@weight=temp$postHap
        outObj@nFold=nfold
        outObj@maxSNP=maxsnps
        outObj@deviance= - temp$deviance
        outObj@bestsize=c(0, 1:maxsnps, (maxsnps-1):1 )[which.min(outObj@deviance)]

        if (outObj@bestsize==0) {
            fSubset <- outObj
            fSubset@modelmethod <- "Cross-Val"
        } else {
            fSubset <- finalsubset(outObj, phenoData=unlist(haploObj@pheno[, status]), match(Minherit, c("additive", "dominant", "recessive")))
        }
        
        return(fSubset)
	
    } ## end of (ModSelMethod=="Cross-Val")
    else if (ModSelMethod=="BIC") {
	alpha <- 2*log(length(haploObj@nPosHapPair))
	outputsnp <- rep(0,ncol(haploObj@haploSeq))
	outhapvec <- rep(0,1000)
	outhapfreq <- rep(0,100)
	varstore <- rep(0,1000)
	coef <- rep(0,100)
	outnhaps <- 0
        bestsize <- 0
	deviance <- rep(0,2*maxsnps)
        phase <- 0
	temp <- .C("stepwise_search_alpha", 
                   id=as.integer(haploObj@poolHapPair),
                   nSubj=as.integer(length(haploObj@nPosHapPair)),
                   nObs=as.integer(length(haploObj@poolHapPair)),
                   nPair=as.integer(haploObj@nPosHapPair),
                   postHap=as.numeric(haploObj@post),
                   ccStatus=as.integer(rep(unlist(haploObj@pheno[, status]), haploObj@nPosHapPair)),
                   nLoci=as.integer(ncol(haploObj@haploSeq)),
                   nHap=as.integer(length(haploObj@haploFreq)),
                   hap1=as.integer(haploObj@hap1),
                   hap2=as.integer(haploObj@hap2),
                   uhap=as.integer(uhap),
                   hapFreq=as.numeric(haploObj@haploFreq),
                   bestsize=as.integer(bestsize),
                   Minherit=match(Minherit, c("additive", "dominant", "recessive")),
                   deviance=as.numeric(deviance),
                   maxSNP=as.integer(maxsnps),
                   converge=as.numeric(tol),
                   alpha=as.numeric(alpha),
                   verbose=as.integer(sum(verbose)),
                   phase=as.integer(phase), 
                   outputsnp=as.integer(outputsnp), 
                   outhapvec=as.integer(outhapvec), 
                   outhapfreq=as.numeric(outhapfreq), 
                   varstore=as.numeric(varstore), 
                   coef=as.numeric(coef), 
                   outnhaps=as.integer(outnhaps),
                   PACKAGE="SHARE"
                   )

        outObj <- haploObj
        class(outObj) <- "share"
        outObj@uhap=uhap
        outObj@weight=temp$postHap
        outObj@nFold=nfold
        outObj@maxSNP=maxsnps
        outObj@deviance= temp$deviance
        outObj@bestsize=temp$bestsize
	
        if (outObj@bestsize==0) {
            outObj2 <- outObj
            outObj2@modelmethod <- "BIC"
            
        } else {
            
	    ## var-covar matrix
            amat <- matrix(temp$varstore[1:(temp$outnhaps ^2)], nrow=temp$outnhaps,
                           ncol=temp$outnhaps, byrow=TRUE)
            bmat <- matrix(temp$varstore[((temp$outnhaps ^2)+1):(2*(temp$outnhaps ^2))],
                           nrow=temp$outnhaps, ncol=temp$outnhaps, byrow=TRUE)

            require(MASS) ## for ginv funciton
            var<-ginv(amat)%*% bmat %*% ginv(amat)
            
            pval <- drop(1-pchisq(t(temp$coef[2:temp$outnhaps]) %*% ginv(var[2:temp$outnhaps,2:temp$outnhaps]) %*% temp$coef[2:temp$outnhaps],temp$outnhaps-1))

            ## name of haplotypes & SNPs
            snpNames <- colnames(outObj@haploSeq)[sort(temp$outputsnp)]
            hapName <- paste("hap", 1:temp$outnhaps, sep="")
            ## add (Intercept) to hap1
            ##hapName[1] <- paste(hapName[1], "(Intercept)")
            
            ## data frame for the output haplotypes
            hapDF <- data.frame(matrix(temp$outhapvec[1:(temp$bestsize*temp$outnhaps)],
                                       nrow=temp$outnhaps, ncol=temp$bestsize, byrow=TRUE)
                                )
            colnames(hapDF) <- snpNames
            rownames(hapDF) <- hapName

            ## haplotype frequency
            hapFreq <- temp$outhapfreq[1:temp$outnhaps]
            names(hapFreq) <- hapName

            ## coefficient
            coef <- temp$coef[1:temp$outnhaps]
            names(coef) <- hapName

            ## s.e.
            stdErr <- sqrt(diag(var))

            ## z-score
            zScore <- coef / stdErr

            ## p-values
            pValue <- 1-pnorm(abs(zScore))

            ## GLM-like output
            haploTest <- data.frame(coef, stdErr, zScore, pValue)
            colnames(haploTest) <- c("Coefficient", "Standard Error", "z-Score", "p-Value")
            hapbase <- which.max(temp$outhapfreq[1:temp$outnhaps])
            hapName1 <- paste("hap",c(hapbase,(1:temp$outnhaps)[-hapbase]))
            hapName1[1] <- paste(hapName1[1], "(Intercept)")
            rownames(haploTest) <- hapName1

            outObj2 <- outObj
            outObj2@modelmethod <- "BIC"
            outObj2@finalHap <- hapDF
            outObj2@finalHapFreq <- hapFreq
            outObj2@finalHapTest <- haploTest
            outObj2@globalP <- pval
            outObj2@inherit <- c("additive", "dominant", "recessive")[Minherit]
        }  
        
        return(outObj2)
    } ## end of (ModSelMethod=="BIC")
}

shareTest <- function(outObj,            ## the returned object from chare function
                      haploObj,          ## haploObj
                      status,
                      tol=1e-8,
                      verbose=FALSE,     ## TRUE: turn on log file
                      nperm=1000         ## the number of permutations 
                      ) {
    phase <- 0

    count <- 0
    for (i in 1:nperm) {
        set.seed(38329832+i)
        if (i%%10==0) cat(i, "... ")
        haploObj@pheno[, "y"] <- sample(as.vector(haploObj@pheno[, status]))
        out <- cshare(haploObj, "y", outObj@nFold, outObj@maxSNP, tol, verbose,
                      outObj@modelmethod, outObj@inherit)

        if(out@bestsize!=0){
            if (outObj@globalP > out@globalP) count <- count+1
        }
        pval <- count/(i+1)
        if (count>=50) {
            cat("Permutation test stops early due to insignificance.","\n")
            break
        }
    }
    return(pval)
}
#line 1 "/private/tmp/Rtmpl8ol26/R.INSTALL10d63af1/SHARE/R/zzz.SHARE.R"
## without NAMESPACE
##.First.lib <- function(lib, pkg) {
##   library.dynam("SHARE", pkg, lib)
##}

## with NAMESPACE
.onLoad <- function(lib, pkg) {
   library.dynam("SHARE", pkg, lib)
}
