# tryCchol.r
 #
 # demonstration of tryCatch with Cholesky factorization
 #
 # three test matrices
 # first is 3x3, positive definite
  A1 <- matrix(c(9,2,-2,2,1,0,-2,0,4),3,3)
  A1
 # second is 4x4, cov mx of multinomial, positive semi-definite
  p <- c(.1,.2,.3,.4)
  A2 <- diag(p) - outer(p,p)
  A2
 # third is 3x3, pos semi-def, but only rank 1
  A3 <- matrix(1,3,3)
  A3
 #
 # usual function
  chol(A1)
 # expressions below not executed since they will stop execution
 #  of batch job -- try these out in interactive session
 #
 #
  #v2 <- chol(A2)
  #v2
  #v3 <- chol(A3)
  #v3
 #
 #
 # first fix -- just throw a warning flag
  cholw <- function(A) {
    tryCatch(chol(A),error=function(e){print("not good")})  }
 #           expression    what to do with error
 #
 # notice response to the problem cases
 # 
 #
 # no problem
  cholw(A1)
 # usually gives error and stops
  cholw(A2)
 # also gives error and stops
  cholw(A3)
 # fancy fix -- try to fit one rank smaller
  chol1 <- function(A) {
    tryCatch(chol(A),error=function(e){ # fix if last row/col dependent
                      print(e)          # usual error message
                      print("trying one smaller")
                      lA1 <- nrow(A)-1            # one smaller
                      R1 <- chol(A[1:lA1,1:lA1])  # Cholesky factor
                      cbind(R1,forwardsolve(t(R1),A[1:lA1,lA1+1]))})  }
 # use it on the three cases
 #
 # first one should give no change
   chol1(A1)
 # second one is the interesting one
   v2 <- chol1(A2)
   v2
   t(v2)%*%v2
 # third one will still fail -- deficient in rank by 2
   chol1(A3)
 # done
  rm(list=ls())
  q()