R/supportRegion.old.R

In hisse: Hidden State Speciation and Extinction

######################################################################################################################################
######################################################################################################################################
### Adaptive Bootstrap -- Simulating confidence intervals for parameters estimated in HiSSE
######################################################################################################################################
######################################################################################################################################

SupportRegion.old <- function(hisse.old.obj, n.points=1000, scale.int=0.1, desired.delta=2, min.number.points=10, output.type="turnover", hidden.states=TRUE, condition.on.survival=TRUE, root.type="madfitz", root.p=NULL, verbose=TRUE){
    if(inherits(hisse.old.obj, what="hisse.null4.fit")){
        phy <- hisse.old.obj$phy
        data <- hisse.old.obj$data
        data.new<-data.frame(data[,2], data[,2], row.names=data[,1])
        data.new<-data.new[phy$tip.label,]
        f <- hisse.old.obj$f
        np <- max(hisse.old.obj$index.par)
        par <- numeric(np)
        free.parameters <- which(hisse.old.obj$index.par < max(hisse.old.obj$index.par))
        np.sequence <- 1:np
        for(i in np.sequence){
            par[i] <- hisse.old.obj$solution[which(hisse.old.obj$index.par == np.sequence[i])[1]]
        }

        lower <- exp(hisse.old.obj$lower.bounds)
        upper <- exp(hisse.old.obj$upper.bounds)

        params_and_bounds <- data.frame(lower=lower, par=par, upper=upper)
        while(any(c(par<lower, par>upper))) {
            if(any(par<lower)) {
                print("Problem with bounds")
                print(params_and_bounds)
                print(paste("At least one of your parameters is lower than your lower bound: parameter(s)", paste(which(par<lower), collapse=', ')))
                print("We're resetting the bounds automatically")
                lower <- lower*0.1
            }
            if(any(par>upper)) {
                print("Problem with bounds")
                print(params_and_bounds)
                print(paste("At least one of your parameters is higher than your upper bound: parameter(s)", paste(which(par>upper), collapse=', ')))
                print("We're resetting the bounds automatically")
                upper <- upper*10
            }
        }


        #Bad Jeremy! Hard-coded column headers...
        if(output.type == "turnover"){
            interval.names <- c("lnLik", "turn.0A", "turn.0B", "turn.0C", "turn.0D", "turn.1A", "turn.1B", "turn.1C", "turn.1D", "eps.0A", "eps.0B", "eps.0C", "eps.0D", "eps.1A", "eps.1B", "eps.1C", "eps.1D", "q0B0A", "q0C0A", "q0D0A", "q1A0A", "q0A0B", "q0C0B", "q0D0B", "q1B0B", "q0A0C", "q0B0C", "q0D0C", "q1C0C", "q0A0D", "q0B0D", "q0C0D", "q1D0D", "q0A1A", "q1B1A", "q1C1A", "q1D1A", "q0B1B", "q1A1B", "q1C1B", "q1D1B", "q0C1C", "q1A1C", "q1B1C", "q1D1C", "q0D1D", "q1A1D", "q1B1D", "q1C1D")
        }
        if(output.type == "net.div"){
            interval.names <- c("lnLik", "netdiv.0A", "netdiv.0B", "netdiv.0C", "netdiv.0D", "netdiv.1A", "netdiv.1B", "netdiv.1C", "netdiv.1D", "eps.0A", "eps.0B", "eps.0C", "eps.0D", "eps.1A", "eps.1B", "eps.1C", "eps.1D", "q0B0A", "q0C0A", "q0D0A", "q1A0A", "q0A0B", "q0C0B", "q0D0B", "q1B0B", "q0A0C", "q0B0C", "q0D0C", "q1C0C", "q0A0D", "q0B0D", "q0C0D", "q1D0D", "q0A1A", "q1B1A", "q1C1A", "q1D1A", "q0B1B", "q1A1B", "q1C1B", "q1D1B", "q0C1C", "q1A1C", "q1B1C", "q1D1C", "q0D1D", "q1A1D", "q1B1D", "q1C1D")
        }
        if(output.type == "raw"){
            interval.names <- c("lnLik", "lambda.0A", "lambda.0B", "lambda.0C", "lambda.0D", "lambda.1A", "lambda.1B", "lambda.1C", "lambda.1D", "mu.0A", "mu.0B", "mu.0C", "mu.0D", "mu.1A", "mu.1B", "mu.1C", "mu.1D", "q0B0A", "q0C0A", "q0D0A", "q1A0A", "q0A0B", "q0C0B", "q0D0B", "q1B0B", "q0A0C", "q0B0C", "q0D0C", "q1C0C", "q0A0D", "q0B0D", "q0C0D", "q1D0D", "q0A1A", "q1B1A", "q1C1A", "q1D1A", "q0B1B", "q1A1B", "q1C1B", "q1D1B", "q0C1C", "q1A1C", "q1B1C", "q1D1C", "q0D1D", "q1A1D", "q1B1D", "q1C1D")
        }

        interval.results <- AdaptiveConfidenceIntervalSampling.old(par, lower=lower, upper=upper, desired.delta = desired.delta, n.points=n.points, verbose=verbose, phy=phy, data=data.new, index.par=hisse.old.obj$index.par, f=f, hidden.states=hidden.states, condition.on.survival=condition.on.survival, root.type=root.type, root.p=root.p, scale.int=scale.int, hisse.null.four=TRUE, min.number.points=min.number.points)
        interval.results.final <- matrix(0, n.points+1, length(hisse.old.obj$index.par))
        for(i in 1:(n.points+1)){
            par.rep <- unlist(interval.results[i,-1],use.names=FALSE)
            interval.results.final[i,] <- c(par.rep,0)[hisse.old.obj$index.par]
        }
        interval.results.final <- cbind(interval.results[,1], interval.results.final)
        if(output.type == "net.div"){
            lambda.0A <- interval.results.final[,2] / (1 + interval.results.final[,10])
            lambda.0B <- interval.results.final[,3] / (1 + interval.results.final[,11])
            lambda.0C <- interval.results.final[,4] / (1 + interval.results.final[,12])
            lambda.0D <- interval.results.final[,5] / (1 + interval.results.final[,13])
            lambda.1A <- interval.results.final[,6] / (1 + interval.results.final[,14])
            lambda.1B <- interval.results.final[,7] / (1 + interval.results.final[,15])
            lambda.1C <- interval.results.final[,8] / (1 + interval.results.final[,16])
            lambda.1D <- interval.results.final[,9] / (1 + interval.results.final[,17])
            mu.0A <- (interval.results.final[,2] * interval.results.final[,10]) / (1 + interval.results.final[,10])
            mu.0B <- (interval.results.final[,3] * interval.results.final[,11]) / (1 + interval.results.final[,11])
            mu.0C <- (interval.results.final[,4] * interval.results.final[,12]) / (1 + interval.results.final[,12])
            mu.0D <- (interval.results.final[,5] * interval.results.final[,13]) / (1 + interval.results.final[,13])
            mu.1A <- (interval.results.final[,6] * interval.results.final[,14]) / (1 + interval.results.final[,14])
            mu.1B <- (interval.results.final[,7] * interval.results.final[,15]) / (1 + interval.results.final[,15])
            mu.1C <- (interval.results.final[,8] * interval.results.final[,16]) / (1 + interval.results.final[,16])
            mu.1D <- (interval.results.final[,9] * interval.results.final[,17]) / (1 + interval.results.final[,17])
            interval.results.final[,2] <- lambda.0A - mu.0A
            interval.results.final[,3] <- lambda.0B - mu.0B
            interval.results.final[,4] <- lambda.0C - mu.0C
            interval.results.final[,5] <- lambda.0D - mu.0D
            interval.results.final[,6] <- lambda.1A - mu.1A
            interval.results.final[,7] <- lambda.1B - mu.1B
            interval.results.final[,8] <- lambda.1C - mu.1C
            interval.results.final[,9] <- lambda.1D - mu.1D
        }
        if(output.type == "raw"){
            lambda.0A <- interval.results.final[,2] / (1 + interval.results.final[,10])
            lambda.0B <- interval.results.final[,3] / (1 + interval.results.final[,11])
            lambda.0C <- interval.results.final[,4] / (1 + interval.results.final[,12])
            lambda.0D <- interval.results.final[,5] / (1 + interval.results.final[,13])
            lambda.1A <- interval.results.final[,6] / (1 + interval.results.final[,14])
            lambda.1B <- interval.results.final[,7] / (1 + interval.results.final[,15])
            lambda.1C <- interval.results.final[,8] / (1 + interval.results.final[,16])
            lambda.1D <- interval.results.final[,9] / (1 + interval.results.final[,17])
            mu.0A <- (interval.results.final[,2] * interval.results.final[,10]) / (1 + interval.results.final[,10])
            mu.0B <- (interval.results.final[,3] * interval.results.final[,11]) / (1 + interval.results.final[,11])
            mu.0C <- (interval.results.final[,4] * interval.results.final[,12]) / (1 + interval.results.final[,12])
            mu.0D <- (interval.results.final[,5] * interval.results.final[,13]) / (1 + interval.results.final[,13])
            mu.1A <- (interval.results.final[,6] * interval.results.final[,14]) / (1 + interval.results.final[,14])
            mu.1B <- (interval.results.final[,7] * interval.results.final[,15]) / (1 + interval.results.final[,15])
            mu.1C <- (interval.results.final[,8] * interval.results.final[,16]) / (1 + interval.results.final[,16])
            mu.1D <- (interval.results.final[,9] * interval.results.final[,17]) / (1 + interval.results.final[,17])
            interval.results.final[,2] <- lambda.0A
            interval.results.final[,3] <- lambda.0B
            interval.results.final[,4] <- lambda.0C
            interval.results.final[,5] <- lambda.0D
            interval.results.final[,6] <- lambda.1A
            interval.results.final[,7] <- lambda.1B
            interval.results.final[,8] <- lambda.1C
            interval.results.final[,9] <- lambda.1D
            interval.results.final[,10] <- mu.0A
            interval.results.final[,11] <- mu.0B
            interval.results.final[,12] <- mu.0C
            interval.results.final[,13] <- mu.0D
            interval.results.final[,14] <- mu.1A
            interval.results.final[,15] <- mu.1B
            interval.results.final[,16] <- mu.1C
            interval.results.final[,17] <- mu.1D
        }
        interval.results.in <- interval.results.final[which(interval.results.final[,1] - min(interval.results.final[,1])<=desired.delta),]
        if(inherits(interval.results.in, what="numeric")){
            stop("Only the MLE is in the desired range. Try reducing scale.int.", call.=FALSE)
        }else{
            ci.interval = apply(interval.results.in, 2, quantile)
            colnames(interval.results.final) <- colnames(interval.results.in) <- colnames(ci.interval) <- interval.names
            obj = NULL
            obj$ci <- ci.interval
            obj$points.within.region = interval.results.in
            obj$all.points = interval.results.final
            class(obj) = "hisse.support"
            return(obj)
        }
    }else{
        phy <- hisse.old.obj$phy
        data <- hisse.old.obj$data
        data.new<-data.frame(data[,2], data[,2], row.names=data[,1])
        data.new<-data.new[phy$tip.label,]
        f <- hisse.old.obj$f
        np <- max(hisse.old.obj$index.par)-1
        par <- numeric(np)
        free.parameters <- which(hisse.old.obj$index.par < max(hisse.old.obj$index.par))
        np.sequence <- 1:np
        for(i in np.sequence){
            par[i] <- hisse.old.obj$solution[which(hisse.old.obj$index.par == np.sequence[i])[1]]
        }

        lower <- exp(hisse.old.obj$lower.bounds)
        upper <- exp(hisse.old.obj$upper.bounds)

        #Bad Jeremy! Hard-coded column headers...
        if(output.type == "turnover"){
            interval.names <- c("lnLik", "turn.0A", "turn.1A", "turn.0B", "turn.1B", "eps.0A", "eps.1A", "eps.0B", "eps.1B","q1A0A","q0B0A","q1B0A","q0A1A","q0B1A","q1B1A","q0A0B","q1A0B","q1B0B","q0A1B","q1A1B","q0B1B","turn.alpha.0A","turn.alpha.1A", "turn.alpha.0B", "turn.alpha.1B", "turn.beta.0A","turn.beta.1A", "turn.beta.0B", "turn.beta.1B", "eps.alpha.0A","eps.alpha.1A", "eps.alpha.0B", "eps.alpha.1B", "eps.beta.0A","eps.beta.1A", "eps.beta.0B", "eps.beta.1B", "turn.slice.0A","turn.slice.1A", "turn.slice.0B", "turn.slice.1B", "eps.slice.0A","eps.slice.1A", "eps.slice.0B", "eps.slice.1B", "q0A1A.slice","q1A0A.slice","q0A0B.slice","q0B0A.slice","q1A1B.slice","q1B1A.slice","q0A1B.slice","q1B0A.slice","q1A0B.slice","q0B1A.slice","q1B0B.slice","q0B1B.slice")
        }
        if(output.type == "net.div"){
            interval.names <- c("lnLik", "netdiv.0A", "netdiv.1A", "netdiv.0B", "netdiv.1B", "eps.0A", "eps.1A", "eps.0B", "eps.1B","q1A0A","q0B0A","q1B0A","q0A1A","q0B1A","q1B1A","q0A0B","q1A0B","q1B0B","q0A1B","q1A1B","q0B1B","turn.alpha.0A","turn.alpha.1A", "turn.alpha.0B", "turn.alpha.1B", "turn.beta.0A","turn.beta.1A", "turn.beta.0B", "turn.beta.1B", "eps.alpha.0A","eps.alpha.1A", "eps.alpha.0B", "eps.alpha.1B", "eps.beta.0A","eps.beta.1A", "eps.beta.0B", "eps.beta.1B", "turn.slice.0A","turn.slice.1A", "turn.slice.0B", "turn.slice.1B", "eps.slice.0A","eps.slice.1A", "eps.slice.0B", "eps.slice.1B", "q0A1A.slice","q1A0A.slice","q0A0B.slice","q0B0A.slice","q1A1B.slice","q1B1A.slice","q0A1B.slice","q1B0A.slice","q1A0B.slice","q0B1A.slice","q1B0B.slice","q0B1B.slice")
        }
        if(output.type == "raw"){
            interval.names <- c("lnLik", "lambda.0A", "lambda.1A", "lambda.0B", "lambda.1B", "mu.0A", "mu.1A", "mu.0B", "mu.1B","q1A0A","q0B0A","q1B0A","q0A1A","q0B1A","q1B1A","q0A0B","q1A0B","q1B0B","q0A1B","q1A1B","q0B1B","turn.alpha.0A","turn.alpha.1A", "turn.alpha.0B", "turn.alpha.1B", "turn.beta.0A","turn.beta.1A", "turn.beta.0B", "turn.beta.1B", "eps.alpha.0A","eps.alpha.1A", "eps.alpha.0B", "eps.alpha.1B", "eps.beta.0A","eps.beta.1A", "eps.beta.0B", "eps.beta.1B", "turn.slice.0A","turn.slice.1A", "turn.slice.0B", "turn.slice.1B", "eps.slice.0A","eps.slice.1A", "eps.slice.0B", "eps.slice.1B", "q0A1A.slice","q1A0A.slice","q0A0B.slice","q0B0A.slice","q1A1B.slice","q1B1A.slice","q0A1B.slice","q1B0A.slice","q1A0B.slice","q0B1A.slice","q1B0B.slice","q0B1B.slice")
        }

        interval.results <- AdaptiveConfidenceIntervalSampling.old(par, lower=lower, upper=upper, desired.delta = desired.delta, n.points=n.points, verbose=verbose, phy=phy, data=data.new, index.par=hisse.old.obj$index.par, f=f, hidden.states=hidden.states, condition.on.survival=condition.on.survival, root.type=root.type, root.p=root.p, scale.int=scale.int, hisse.null.four=FALSE, min.number.points=min.number.points)
        interval.results.final <- matrix(0, n.points+1, length(hisse.old.obj$index.par))
        for(i in 1:(n.points+1)){
            par.rep <- unlist(interval.results[i,-1],use.names=FALSE)
            interval.results.final[i,] <- c(par.rep,0)[hisse.old.obj$index.par]
        }
        interval.results.final[,21:56] = 1
        interval.results.final <- cbind(interval.results[,1], interval.results.final)
        if(output.type == "net.div"){
            lambda.0A <- interval.results.final[,2] / (1 + interval.results.final[,6])
            lambda.1A <- interval.results.final[,3] / (1 + interval.results.final[,7])
            lambda.0B <- interval.results.final[,4] / (1 + interval.results.final[,8])
            lambda.1B <- interval.results.final[,5] / (1 + interval.results.final[,9])
            mu.0A <- (interval.results.final[,2] * interval.results.final[,6]) / (1 + interval.results.final[,6])
            mu.1A <- (interval.results.final[,3] * interval.results.final[,7]) / (1 + interval.results.final[,7])
            mu.0B <- (interval.results.final[,4] * interval.results.final[,8]) / (1 + interval.results.final[,8])
            mu.1B <- (interval.results.final[,5] * interval.results.final[,9]) / (1 + interval.results.final[,9])
            interval.results.final[,2] <- lambda.0A - mu.0A
            interval.results.final[,3] <- lambda.1A - mu.1A
            interval.results.final[,4] <- lambda.0B - mu.0B
            interval.results.final[,5] <- lambda.1B - mu.1B
        }
        if(output.type == "raw"){
            lambda.0A <- interval.results.final[,2] / (1 + interval.results.final[,6])
            lambda.1A <- interval.results.final[,3] / (1 + interval.results.final[,7])
            lambda.0B <- interval.results.final[,4] / (1 + interval.results.final[,8])
            lambda.1B <- interval.results.final[,5] / (1 + interval.results.final[,9])
            mu.0A <- (interval.results.final[,2] * interval.results.final[,6]) / (1 + interval.results.final[,6])
            mu.1A <- (interval.results.final[,3] * interval.results.final[,7]) / (1 + interval.results.final[,7])
            mu.0B <- (interval.results.final[,4] * interval.results.final[,8]) / (1 + interval.results.final[,8])
            mu.1B <- (interval.results.final[,5] * interval.results.final[,9]) / (1 + interval.results.final[,9])
            interval.results.final[,2] <- lambda.0A
            interval.results.final[,3] <- lambda.1A
            interval.results.final[,4] <- lambda.0B
            interval.results.final[,5] <- lambda.1B
            interval.results.final[,6] <- mu.0A
            interval.results.final[,7] <- mu.1A
            interval.results.final[,8] <- mu.0B
            interval.results.final[,9] <- mu.1B
        }
        interval.results.in <- interval.results.final[which(interval.results.final[,1] - min(interval.results.final[,1])<=desired.delta),]
        if(inherits(interval.results.in, what="numeric")){
            stop("Only the MLE is in the desired range. Try reducing scale.int.", call.=FALSE)
        }else{
            ci.interval = apply(interval.results.in, 2, quantile)
            colnames(interval.results.final) <- colnames(interval.results.in) <- colnames(ci.interval) <- interval.names
            if(inherits(interval.results.in, what="numeric")){
                stop("Only the MLE is in the desired range. Try reducing scale.int.", call.=FALSE)
            }else{
                obj = NULL
                obj$ci <- ci.interval[,1:21]
                obj$points.within.region = interval.results.in[,1:21]
                obj$all.points = interval.results.final[,1:21]
                class(obj) = "hisse.support"
                return(obj)
            }
        }
    }
}


AdaptiveConfidenceIntervalSampling.old <- function(par, lower, upper, desired.delta=2, n.points=5000, verbose=TRUE, phy, data, index.par, f, hidden.states, condition.on.survival, root.type, root.p, scale.int, hisse.null.four=FALSE, min.number.points=10) {
                                        #Wrangle the data so that we can make use of DownPass easily:
    actual.params = which(index.par < max(index.par))
    model.vec <- numeric(length(index.par))
    model.vec[] <- c(par,0)[index.par]
    if(hisse.null.four == TRUE){
        cache = ParametersToPassNull(phy, data[,1], f=f, model.vec)
        phy$node.label <- NULL
        starting <- -DownPassNull(phy, cache, condition.on.survival=condition.on.survival, root.type=root.type, root.p=root.p)
    }else{
        model.vec.tmp = model.vec[21:56]
        model.vec.tmp[model.vec.tmp==0] = 1
        model.vec[21:56] = model.vec.tmp
        cache = ParametersToPass(phy, data[,1], model.vec, f=f, timeslice=NULL, hidden.states=hidden.states)
        cache$turnover.beta.factor0 = 1 / dbeta(0.1, model.vec[21], model.vec[25])
        cache$turnover.beta.factor1 = 1 / dbeta(0.1, model.vec[22], model.vec[26])
        cache$turnover.beta.factorA = 1 / dbeta(0.1, model.vec[23], model.vec[27])
        cache$turnover.beta.factorB = 1 / dbeta(0.1, model.vec[24], model.vec[28])
        cache$eps.beta.factor0 = 1 / dbeta(0.1, model.vec[29], model.vec[33])
        cache$eps.beta.factor1 = 1 / dbeta(0.1, model.vec[30], model.vec[34])
        cache$eps.beta.factorA = 1 / dbeta(0.1, model.vec[31], model.vec[35])
        cache$eps.beta.factorB = 1 / dbeta(0.1, model.vec[32], model.vec[36])
        phy$node.label <- NULL
        #############################################################
        #Now assess the likelihood at the MLE:
        starting <- -DownPass(phy, cache, hidden.states=hidden.states, condition.on.survival=condition.on.survival, root.type=root.type, root.p=root.p)
    }
    #Generate the multipliers for feeling the boundaries:
    min.multipliers <- rep(1, length(par))
    max.multipliers <- rep(1, length(par))
    results <- data.frame(data.frame(matrix(nrow=n.points+1, ncol=1+length(par))))
    results[1,] <- unname(c(starting, par))
    for (i in sequence(n.points)) {
        sim.points <- NA
        while(is.na(sim.points[1])) {
            sim.points <- GenerateValues(par, lower=lower, upper=upper, scale.int=scale.int, examined.max=max.multipliers*apply(results[which(results[,1]-min(results[,1], na.rm=TRUE)<=desired.delta),-1], 2, max, na.rm=TRUE), examined.min=min.multipliers*apply(results[which(results[,1]-min(results[,1], na.rm=TRUE)<=desired.delta),-1], 2, min, na.rm=TRUE))
        }
        par <- sim.points
        model.vec <- numeric(length(index.par))
        model.vec[] <- c(sim.points,0)[index.par]
        if(hisse.null.four == TRUE){
            cache = ParametersToPassNull(phy, data[,1], f=f, model.vec)
            second <- -DownPassNull(phy, cache, condition.on.survival=condition.on.survival, root.type=root.type, root.p=root.p)
        }else{
            model.vec.tmp = model.vec[21:56]
            model.vec.tmp[model.vec.tmp==0] = 1
            model.vec[21:56] = model.vec.tmp
            cache = ParametersToPass(phy, data[,1], model.vec, f=f, timeslice=NULL, hidden.states=hidden.states)
            cache$turnover.beta.factor0 = 1 / dbeta(0.1, model.vec[21], model.vec[25])
            cache$turnover.beta.factor1 = 1 / dbeta(0.1, model.vec[22], model.vec[26])
            cache$turnover.beta.factorA = 1 / dbeta(0.1, model.vec[23], model.vec[27])
            cache$turnover.beta.factorB = 1 / dbeta(0.1, model.vec[24], model.vec[28])
            cache$eps.beta.factor0 = 1 / dbeta(0.1, model.vec[29], model.vec[33])
            cache$eps.beta.factor1 = 1 / dbeta(0.1, model.vec[30], model.vec[34])
            cache$eps.beta.factorA = 1 / dbeta(0.1, model.vec[31], model.vec[35])
            cache$eps.beta.factorB = 1 / dbeta(0.1, model.vec[32], model.vec[36])
            second <- -DownPass(phy, cache, hidden.states=hidden.states, condition.on.survival=condition.on.survival, root.type=root.type, root.p=root.p)
        }
        results[i+1,] <- c(second, sim.points)
        if(i%%20==0) {
            for (j in sequence(length(par))) {
                returned.range <- range(results[which((results[,1]-min(results[,1], na.rm=TRUE))<desired.delta), j+1], na.rm=TRUE)
                total.range <- range(results[,j+1], na.rm=TRUE)
                width.ratio <- diff(returned.range)/diff(total.range)
                if(is.na(width.ratio)) {
                    width.ratio=1
                }
                if(width.ratio > 0.5) { #we are not sampling widely enough
                    min.multipliers[j] <- min.multipliers[j] * (1-scale.int)
                    max.multipliers[j] <- max.multipliers[j] * (1+scale.int) #expand the range
                } else {
                    if(width.ratio < 0.02 & i>100) { # we are sampling too widely
                        min.multipliers[j] <- min.multipliers[j] * (1+scale.int)
                        max.multipliers[j] <- max.multipliers[j] * (1-scale.int) #contract the range
                    } else {
                        min.multipliers[j] <- 1
                        max.multipliers[j] <- 1
                    }
                }
            }
        }
        if (verbose && i%%100==0) {
            cat(paste(i, "of", n.points, "points done"), "\n")
        }
    }
    while(length(which((results[,1]-min(results[,1], na.rm=TRUE))<desired.delta))<min.number.points) {
        warning("Did not generate enough points in the region; restarting to create additional points")
        print(paste("Now doing an additional", 2+round(n.points/4), "points to the", dim(results)[1], "ones already done because not enough points in the good enough region were sampled"))
        scale.int <- scale.int*0.5
        new.results <- AdaptiveConfidenceIntervalSampling.old(par=par, lower=lower, upper=upper, desired.delta=desired.delta, n.points=2+round(n.points/4), verbose=verbose, phy=phy, data=data, index.par=index.par, f=f, hidden.states=hidden.states, condition.on.survival=condition.on.survival, root.type=root.type, root.p=root.p, scale.int=scale.int, hisse.null.four=hisse.null.four, min.number.points=0)
        results <- rbind(results, new.results[-1,])
    }
    return(results)
}



GenerateValues <- function(par, lower, upper, scale.int, max.tries=100, expand.prob=0, examined.max, examined.min) {
    pass=FALSE
    tries=0
    while(!pass && tries<=max.tries) {
        tries <- tries+1
        pass=TRUE
        new.vals <- rep(NA, length(par))
        for(i in sequence(length(par))) {
            examined.max[i] <- max(0.001, examined.max[i])
            min.val <- min(max(lower[i], (1-scale.int)*examined.min[i]), examined.max[i]) #just in case min is greater than max
            max.val <- max(min(upper[i], (1+scale.int)*examined.max[i]), examined.min[i])
            if(isTRUE(all.equal(min.val, max.val))) {
                min.val <- min.val * (1-scale.int)
                max.val <- max.val * (1+scale.int)
            }
            new.vals[i] <- runif(1, min.val, max.val)
            if(rbinom(1,1,.1)==1) { #ten percent of the time, try something else
                new.vals[i] <- max.val+100
                ntries <- 0
                while(new.vals[i]>max.val & ntries < 10) {
                    new.vals[i] <- min.val + rexp(1, 1/((max.val-min.val)/2))
                    ntries <- ntries + 1
                }
                if(new.vals[i] > max.val) {
                    new.vals[i] <- runif(1, min.val, max.val) #we give up
                }
            }
            if(new.vals[i]<lower[i]) {
                pass=FALSE
            }
            if(new.vals[i]>upper[i]) {
                pass=FALSE
            }
        }
    }
    if(tries>max.tries) {
        return(NA)
    }
    return(new.vals)
}


print.hisse.support <- function(x,...){

    cat("\nSupport Region\n")
    print(x$ci[,-1])
    cat("\n")
}