CatchAll: Species Richness Estimation with CatchAll

CalculateAnalysisVariables <- function(part1, part2, 
                                       numberParameters, r, fitsCount, 
                                       fitsExtended, 
                                       s, modelNumber, frequency, observedCount) {
  # is it better to pass this in?
  BigChiSq <- 1000000000.0
  maxGoodnessOfFit <- 10
  return_variable <- list()
  return_variable$AIC <- 2 * numberParameters - 2*(part1 + part2)

  if (s[r] - numberParameters - 1 > 0) {
    return_variable$AICc <- return_variable$AIC + (2*numberParameters*(numberParameters+1)/(s[r]-numberParameters-1))
    return_variable$AICcFlag <- 1
  }
  
  ## calculate ChiSq, no binning
  
  # removed observedCountNoF0 ? not sure what that is
  chiSqAll <- ChiSqFunction(r, fitsCount, modelNumber, frequency, observedCount, s)
  return_variable$chiSq <- chiSqAll
  
  ## calculate Goodness of Fit

  df <- frequency[r] - numberParameters 
  test <- (chiSqAll - df)/sqrt(2*df)
  message(paste("test: ", test, sep = " "))
  message(paste("chiSqAll: ", chiSqAll, sep = " "))
  GOF0 <- list()
  #bigChisq is a global variable
  if (test < maxGoodnessOfFit  & chiSqAll < BigChiSq){
    GOF0 <- GoodnessOfFit(chiSqAll, df)
    return_variable$GOF0Check <- GOF0$flag
    return_variable$GOF0 <- GOF0$gof
  } else{
    return_variable$GOF0Check <- 0
  }
  
  ## calculate ChiSq, bin 5
  chiSq5 <- ChiSqBin(r, fitsExtended, 5, df, numberParameters, frequency, s, observedCount)
  
  df <- chiSq5$df
  GOF5Check <- chiSq5$flag
  chiSq5 <- chiSq5$chiSq
 
  ## calculate goodness of fit
  test <- (chiSq5 - df)/sqrt(2*df)
  if (test < maxGoodnessOfFit & GOF5Check == 1 & chiSq5 < BigChiSq) {
    GOF5 <- GoodnessOfFit(chiSq5, df)
    return_variable$GOF5Check <- GOF5$flag
    return_variable$GOF5 <- GOF5$gof
  }
  return_variable
}


ChiSqFunction <- function(r, fitsCount, modelNumber,
                  frequency, observedCount, s) {
  message("IN CHI SQ FUNCTION")
  message(paste("r: ", r, sep = " "))
  message("fitsCount")
  message(fitsCount)
  message(paste("modelNumber: ", modelNumber, sep = " "))
  message(frequency)
  message("observedCount")
  message(observedCount)
  message("s")
  message(s)
  #why is this diff if original is same
  chiSqTemporary <- 0
  sumFit <- 0
  rr <- 1
  if (frequency[1] == 0) {
    stop("first frequency is 0?")
  }
  
  message(paste("freq[r]: ", frequency[r]))
  # this bizarre looking flow adjusts for non-contiguous frequencies
  for (t in (1:frequency[r])) {
    if (t == frequency[rr]) {
      # message(paste("fitsCount[t]: ", fitsCount[t + 1], sep = " "))
      # message(paste("observedCount[rr]: ", observedCount[rr], sep = " "))
      # message(paste("observedCount[rr] - fitsCount[t]: ",observedCount[rr] - fitsCount[t+1], sep = " " ))
      # message(paste("(observedCount[rr] - fitsCount[t])^2): ", (observedCount[rr] - fitsCount[t+1])^2, sep = "  "))
      # message(paste("((observedCount[rr] - fitsCount[t])^2)/fitsCount[t])", ((observedCount[rr] - fitsCount[t+1])^2)/fitsCount[t + 1]), sep = " ")
      if (modelNumber < 6) {
        chiSqTemporary <- chiSqTemporary + ((observedCount)[rr] - fitsCount[t])^2/fitsCount[t]
      } else {
        chiSqTemporary <- chiSqTemporary + (((observedCount[rr] - fitsCount[t + 1])^2)/fitsCount[t + 1])
      }
      rr <- rr+1
    } else {
      chiSqTemporary <- chiSqTemporary + fitsCount[t]
    }
    message(paste("chiSqTemporary in loop: ", chiSqTemporary, sep = " "))
  }
  sumFit <- sum(fitsCount)
  message(paste("fitsCount[1]: ", fitsCount[1], sep = " "))
  #is this just a temp fix for LOGTWLR? do we even need this...
  #sumFit <- sumFit - fitsCount[1]
  message(paste("sumFit: ", sumFit, sep = " "))
  message(paste("chiSqTemporary before if check: ", sep = " "))
  message(paste("s[r]: ", s[r], sep = " "))
  if(modelNumber<6) {
    message("hi there")
    chiSqTemporary <- chiSqTemporary + s[r] - sumFit
    message(paste("chiSqTemporary: ", chiSqTemporary, sep = " "))
  }
  chiSqTemporary
}

ChiSqBin <- function(r, fitsExtended, bin, 
                     df, numberParameters, 
                     frequency, s, observedCount) {
  extendedTau <- frequency[r] * 4
  
  ## find terminal indices of binned cells
  check <- rep(NA, extendedTau) 
  accumulatedFit <- 0
  df <- 0
  stop <- 0
  t <- 1

# 
  # message("extendedTau")
  # message(extendedTau)
  # message("check")
  # message(check)
  # message("r")
  # message(r)
  # message("frequency")
  # message(frequency)
  # message("bin")
  # message(bin)
  # message("s")
  # message(s)
  # message("fitsExtended")
  # message(fitsExtended)
 
  while(t <= extendedTau  &  accumulatedFit < bin & (s[r]-accumulatedFit) >= bin) {
    check[t] <- 0
    accumulatedFit  <- accumulatedFit + fitsExtended[t]
    
    # message("t")
    # message(t)
    # message("fitsExtended[t]")
    # message(fitsExtended[t])
    # message("accumulatedFit")
    # message(accumulatedFit)
    # message("extended tau")
    # message(extendedTau)
    # message("+------------------------------+")
    if (accumulatedFit >= bin  & (s[r] - accumulatedFit) >= bin) {
      check[t] <- 1
      df <- df + 1
      stop <- t
      accumulatedFit <- 0
    }
    t <- t + 1
  }
  
  ## todo: fix
 # check[t-1]<-0
  

  ## check for enough data for bininng and positive df
  chiSqTemporary <- 0
  df <- df - numberParameters
  if (stop > 0 & df > 0) {
    cellObservation <- 0
    cellFit <- 0
    rr <- 1

    for (t in 1:extendedTau) {
      if (t <= frequency[r] & t == frequency[rr]) {
        cellObservation <- cellObservation + observedCount[rr]
        rr <- rr + 1
      } 
      cellFit <- cellFit + fitsExtended[t]
      if (check[t] == 1) {
        chiSqTemporary <- chiSqTemporary + (cellObservation-cellFit)^2/cellFit
        cellObservation <- 0
        cellFit <- 0
      }
    }
    observedTail <- cellObservation
    fitTail <- cellFit
    chiSqTemporary <- chiSqTemporary + (observedTail - fitTail)^2/fitTail
    flag <- 1
  } else {
    flag <- 0
  }
  list("chiSq"=chiSqTemporary, "flag"=flag, "df"=df)
}

GoodnessOfFit <- function(chiSqAll, df) {
  v <- df/2 + 1
  x <- chiSqAll/2
  g <- 1
  p <- 1

  while (v <= 2) {
    g <- g*x
    p <- p*v + g
    v <- v + 1
  }
  
  j <- floor(2.5*(3 + abs(x))) - 1
  f <- 1/(j + v - x)
  
  while (j >= 0) {
    f <- (f*x+1)/(j+v)
    j <- j - 1
  }
  
  p <- p + (f*g*x)
  g <- (1-(2/(7*v^2))*(1-(2/(3*v^2))))/(30*v^2)
  g <- ((g-1)/(12*v)) - (v*(log(v)-1))
  f <- p*exp(g-x)*sqrt(v/(2*pi))
  
  flag <- 1
  if (is.nan(f)) flag <- 0
  gof  <- 1-f*(chiSqAll/2)^(df/2)

  if (is.nan(gof)) flag <- 0
  
  list("gof"=gof, "flag"=flag)
}


GetConfidenceBounds <- function(r, se, sHatSubset, s, observationMaximum) {
  answer <- list()
  if (sHatSubset != s[r]) {
    dtemp <- exp(1.96*sqrt(log(1+se^2/(sHatSubset-s[r])^2)))
    answer$lcb <- s[observationMaximum] + (sHatSubset-s[r])/dtemp
    answer$ucb <- s[observationMaximum] + (sHatSubset-s[r])*dtemp
    answer$test <- 1
  } else{
    answer$test  <- 0
  }
  answer
  # should return a lower and upper bound
}

CheckOutput  <- function(x) {
  ifelse(is.null(x), NA, x)
}

BracetRoot <- function(poissonConstant, momentsInit) {
  
  ## initial guess range
  factor <- 2.0
  x1 <- momentsInit / factor
  x2 <- momentsInit * factor
  
  f1 <- (x1 / (1.0 - exp(-x1))) - poissonConstant
  f2 <- (x2 / (1.0 - exp(-x2))) - poissonConstant
  
  conclusion <- 0
  i <- 0
  while (conclusion == 0 & i < 20) {
    ## one estimate is negative and the other positive
    if ((f1 * f2) < 0.0) {
      conclusion <- 1
    }
    ## move the appropriate bound
    if (abs(f1) < abs(f2)) {
      x1 <- x1 + factor * (x1- x2)
      f1 <- (x1 / (1.0 - exp(-x1))) - poissonConstant
    }
    else
    {
      x2 <- x2 + factor * (x2 - x1)
      f2 <- (x2 / (1.0 - exp(-x2))) - poissonConstant
    }
    i <- i + 1
  }
  result <- list()
  result$conclusion <- conclusion
  result$x1 <- x1
  result$x2 <- x2
  result$f1 <- f1
  result$f2 <- f2
  result
}


pow <- function(a, b) {
  a^b
}

#really dumb, delete later, still not understanding the apply family
matrix_apply <- function(m, f) {
  m2 <- m
  for (r in seq(nrow(m2)))
    for (c in seq(ncol(m2)))
      m2[[r, c]] <- f(m2[[r,c]], c)
    return(m2)
}

# DOESN'T WORK FOR LNFACTORIAL 
logFactorial <- function(x) {
  cumsum(log(x))
}

# precision
MatrixInversion <- function(sHat, a00, a0, A) {
  result <- list()
  # complete the symmetric matrix
  A <- A + t(A)
  # message("new A")
  # message(A)
  # after dividing diag / 2, we get the same answer given by C#
  
  diag(A) <- diag(A)/2
  # message("A diag")
  # message(A)
  # 
  #tol = 1e-17
  #5.04722e-22
  aInverse <- try(solve(A), silent = TRUE)
  # message("AInverse after multiplication")
  # message(aInverse)

  if (class(aInverse) != "try-error") {
    answer <- a0 %*% aInverse %*% a0
    if (a00>answer) {
      result$answer <- sqrt(a00-answer)
      result$se <- sqrt(sHat)/answer
      result$flag <- 1
    } else {
      result$flag <- 0
    }
    
  } else {
    result$se <- 0
    result$flag <- 0
  }
  result
}

GetConfidenceBoundsDiscounted <- function(r, SE,
                                            sHatSubset, cStar,excess, LCB,
                                          UCB) {
  ## Confidence Bounds
  if (sHatSubset != s[r])
  {
    dTemp <- exp(1.96 * (sqrt(log(1.0 + (SE * SE / (((sHatSubset - cStar) * (sHatSubset - cStar))))))))
    LCB <- (cStar + excess) + ((sHatSubset - cStar) / dTemp)
    UCB <- (cStar + excess) + ((sHatSubset - cStar) * dTemp)
  }
  
  list("cStar"=cStar, "SE"=SE, "LCB"=LCB, "UCB"=ucb)
}