tests/crv.R

In curve: Dynamic Systems Estimation - Curvature Extensions

# Expecting results from 1 tests

#R --vsize=7000000 --nsize=300000  

# gcinfo(1)
if(!require("dse"))  stop("this test requires dse.")
if(!require("curve"))stop("this test requires curve.")
if(!require("numDeriv"))stop("this test requires numDeriv.")
  Sys.info()
  DSEversion()
  
  

#######################################################################

# Test routines and data for calulating curvatures in Bates and Watts.

#######################################################################
  

curvature.function.tests <- function( verbose=T, synopsis=T, fuzz.small=1e-14, 
	fuzz.large=1e-6, show.details=F, show.extra.details=F)
{# A short set of tests of curvature functions

  max.error <- 0
  if (synopsis & !verbose) cat("All curvature tests ...")

  if (verbose) cat("curvature test 1 ... ")


signif <- 0.05

puromycin <- function(th){
   x <- c(0.02,0.02,0.06,0.06,0.11,0.11,0.22,0.22,0.56,0.56,1.10,1.10)
   y <- c(76,47,97,107,123,139,159,152,191,201,207,200)
   ( (th[1] * x)/(th[2] + x) ) - y
  }

D.anal <- function(th){
  # analytic derivatives. Note numerical approximation gives a very good
  # estimate of these, but neither give D below exactly. The results are very
  # sensitive to th, so rounding error in the reported value of th could explain
  # the difference. But more likely th is correct and D has been rounded for
  # publication - and the analytic D with published th seems to work best.
  # th = c(212.70188549 ,  0.06410027) is the nls est of th for BW published D.
  x <- c(0.02,0.02,0.06,0.06,0.11,0.11,0.22,0.22,0.56,0.56,1.10,1.10)
  y <- c(76,47,97,107,123,139,159,152,191,201,207,200)
  cbind(x/(th[2]+x), -th[1]*x/(th[2]+x)^2, 0, -x/(th[2]+x)^2, 2*th[1]*x/(th[2]+x)^3)
 }

# D matrix from p235. This may be useful for rough comparisons, but rounding
# used for publication introduces substanstial errors. check D.anal1-D.BW
D.BW <- t(matrix(c(
0.237812, -601.458, 0, -2.82773, 14303.4,
0.237812, -601.458, 0, -2.82773, 14303.4,
0.483481, -828.658, 0, -3.89590, 13354.7,
0.483481, -828.658, 0, -3.89590, 13354.7,
0.631821, -771.903, 0, -3.62907, 8867.4,
0.631821, -771.903, 0, -3.62907, 8867.4,
0.774375, -579.759, 0, -2.72571, 4081.4,
0.774375, -579.759, 0, -2.72571, 4081.4,
0.897292, -305.807, 0, -1.43774,  980.0,
0.897292, -305.807, 0, -1.43774,  980.0,
0.944936, -172.655, 0, -0.81173,  296.6,
0.944936, -172.655, 0, -0.81173,  296.6),   5,12))

   D.anal <- D.anal(c(212.7000, 0.0641))
#   D.anal2 <- D.anal(c(212.683, 0.0641194))
   D.calc <- genD(puromycin,c(212.7000, 0.0641), method.args=list(d=0.01))
#  D.calc2 <- genD(puromycin,c(212.683, 0.0641194))$D #est using nls
#check if col 4 is a mult. of col 2 max(abs(D.anal1-D.calc1$D))

   if (show.details)
     {cat("model A p329,data set 3 (table A1.3, p269) Bates & Watts (Puromycin example).\n")
      cat("With show.extra.details=T the calculations  on pages 235 - 252 are shown.\n")
      cat("Note Bates & Watts calculation uses s squared = 10.93^2, which is\n")
      cat(" based on both treated and untreat data. The treated data alone\n")
      cat("gives 10.39^2. This causes a slight difference in the curvature statistics.\n") 
      cat("Note that the R22 calculation using numerical gradients is quite\n")
      cat("different from the published result, but still gives similar results.\n")
      }

   r.BW   <- c(2, 12, 0.05, 0.105,        0.045,      0.21,       0.09,        1, 1.05)
   r.test <- c(2, 12, 0.05, 0.1046988133, 0.04541991, 0.21207185, 0.091999935, 1, 1.051959666)
   r.calc <-curvature(D.calc,signif=signif, 
                        show.extra.details=show.details)$stats 
   err <- r.calc - r.test
   max.error <- max(max.error, abs(err))
   ok <- all(abs(err) < fuzz.large)
   if (show.details | !ok) 
     {cat("Using numerical calculation of gradient and acceration from data:\n")
      m <- rbind(r.BW, r.test, r.calc, err)
      dimnames(m) <- list(
                    c("Bates&Watts","test value","calculated ", "difference "),
                    c(dimnames(r.calc)[[2]]))
      print(m)
      cat("Above should compared with Bates and Watts table 7.2, p 257, model A, data set 3.\n")
     }
   # This is by passing the constructor, and probably should be re-thought. The
   # call used to be r.calc <- curvature.Darray(list...
   # but that does not work with namespaces.
   r.calc <- curvature(classed(list(p=2,D=D.anal,
                    f0=puromycin(c(212.7000, 0.0641))), "Darray"), signif=signif, 
	            show.extra.details=show.details)$stats 
   err <- r.calc - r.test
   ok2 <- all(abs(err) < fuzz.large)
   max.error <- max(max.error, abs(err))
   if (show.details | !ok2) 
     {cat("Using D matrix of gradient and acceration calcuated with formulae Bates and Watts p234:\n")
      m <- rbind(r.BW, r.test, r.calc, err)
      dimnames(m) <- list(
                    c("Bates&Watts","test value","calculated ", "difference "),
                    c(dimnames(r.calc)[[2]]))
      print(m)
     }
   ok3 <- max(abs(D.calc$D-D.anal)) < (100 * fuzz.large)
   if (!ok3)
     {cat("max diff. between analtic and numerical D:")
      cat( max(abs(D.calc$D-D.anal)), "\n")
     }
  ok <- ok & ok2 & ok3
  all.ok <- ok
  if (verbose)  {if (ok) cat("ok\n") else  cat("failed!\n") }

  if (verbose) cat("curvature test 2 ... ")

  curvature.test.bod <- function( th){
      x <- matrix(c(1,2,3,4,5,7),6,1)
      y <- c( 8.3, 10.3, 19.0, 16.0, 15.6, 19.8) 
      y - ( th[1] * (1 - exp( - th[2] * x)))
     }

   r <- curvature(genD(curvature.test.bod, c(19.1430,  0.5311),
                       method.args=list(d=0.01)),signif=signif, 
                show.extra.details=show.details)$stats
   r.BW   <- c(2,6, .05, 1.328,   0.184,      3.5,          0.49,        1.0, 1.01 )
   r.test <- c(2,6, .05, 1.32781, 0.18440417, 3.4990419358, 0.485941626, 1.0, 1.008372434 )
   err <- r -r.test
   max.error <- max(max.error, abs(err))
   ok <- all(abs(err) < fuzz.large )
   if (show.details | !ok) 
     {m <- rbind(r.BW, r.test, r, err)
      dimnames(m) <- list(
                    c("Bates&Watts","test value","calculated ", "difference "),
                    c(dimnames(r)[[2]]))
      print(m)
     cat("Above should compared with Bates and Watts table 7.2, p 257, model B, data set 11 (A1.4, p270.).\n")
     }
  all.ok <- all.ok & ok
  if (verbose)  {if (ok) cat("ok\n") else  cat("failed!\n") }


  if (verbose) cat("curvature test 3 ... ")

  curvature.test.isomerization <- function(th)
  {x <- matrix(c(205.8,404.8,209.7,401.6,224.9,402.6,212.7,406.2,133.3,470.9,300.0,301.6
,297.3,314.0,305.7,300.1,305.4,305.2,300.1,106.6,417.2,251.0,250.3,145.1
,90.9,92.9,174.9,187.2,92.7,102.2,186.9,192.6,140.8,144.2,68.3,214.6
,142.2,146.7,142.0,143.7,141.1,141.5,83.0,209.6,83.9,294.4,148.0,291.0
,37.1,36.3,49.4,44.9,116.3,128.9,134.4,134.9,87.6,86.9,81.7,101.7
,10.5,157.1,86.0,90.2,87.4,87.0,66.4,33.0,32.9,41.5,14.7,50.2),   24,3)
   y <- c(3.541,2.397,6.6694,4.722,0.593,0.268,2.797,2.451,3.196,2.021,0.896,5.084,
5.686,1.193,2.648,3.303,3.054,3.302,1.271,11.648,2.002,9.604,7.754,11.590) 

 y - ((th[1]*th[3]*(x[,2]-x[,3]/1.632))/(1+x[,1]*th[2]+x[,2]*th[3]+x[,3]*th[4]))
}

   r <- curvature(genD(curvature.test.isomerization, c(35.9200,0.0708,0.0377,0.1670),
			   method.args=list(d=0.01, r=6)), signif=signif, 
                           show.extra.details=show.details)$stats
   r.BW   <-  c(4, 24, 0.05, 48.39,    0.048, 81.92, 0.081, 0.95,1.01)
#  r.test <-  c(4, 24, 0.05, 46.00941, 0.048, 81.92, 0.081, 0.95,1.01)
   r.test <-  c(4, 24, 0.05, 46.00942, 0.045942003, 77.891671, 0.077777537,
                1.0,   1.05959779)
   err <- r - r.test
   max.error <- max(max.error, abs(err))
   ok <- all(abs(err) < fuzz.large)
   if (ok) cat("Results check ok\n")
   else    cat("Results differ:\n")
   warning("curve test 3 results do not correspond with results in Bates and Watts.")
   if (show.details | !ok) 
     {m <- rbind(r.BW, r.test, r, err)
      dimnames(m) <- list(
                    c("Bates&Watts","test value","calculated ", "difference "),
                    c(dimnames(r)[[2]]))
      print(m)
#Above should compared with Bates and Watts table 7.2, p 257, model M, 
# data set 32, Bates & Watts (data table A1.5, p271).
# Above test has never worked. The intrinsic curvature array differs and 
# may need to be calculated analytically. The typical result is:
#calculated  4 24   0.05  4.600942e+01   0.045942003 
#        77.891671     0.077777537           1.00     1.05959779 
     }

  all.ok <- all.ok & ok
  if (verbose)  {if (ok) cat("ok\n") else  cat("failed!\n") }

  if (synopsis) 
    {if (verbose) cat("All curvature tests completed")
     if (all.ok) cat(" OK\n")
     else    
       {cat(", some FAILED!")
        if((!is.na(max.error)) && (max.error > fuzz.small))
            cat(" ( max.error = ", max.error,")")
        cat("\n")
       }
    }

  if (all.ok) invisible(T)  else stop("FAILED")
}
  
  curvature.function.tests(verbose=T) # 3 fails as usual