R Under development (unstable) (2024-08-21 r87038 ucrt) -- "Unsuffered Consequences" Copyright (C) 2024 The R Foundation for Statistical Computing Platform: x86_64-w64-mingw32/x64 R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details. R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications. Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for an HTML browser interface to help. Type 'q()' to quit R. > library(testthat) > Sys.setenv('OMP_THREAD_LIMIT'=2) > library(rlibkriging) Attaching package: 'rlibkriging' The following objects are masked from 'package:base': load, save > > ##library(rlibkriging, lib.loc="bindings/R/Rlibs") > ##library(testthat) > default_rcond_checked = rlibkriging:::linalg_chol_rcond_checked() > default_num_nugget = rlibkriging:::linalg_get_num_nugget() > > context("LinearAlgebra") > > # pack=list.files(file.path("bindings","R"),pattern = ".tar.gz",full.names = T) > # install.packages(pack,repos=NULL) > # #library(rlibkriging) > ##library(rlibkriging, lib.loc="bindings/R/Rlibs") > ##library(testthat) > > ####################################################################### > > set.seed(123) > n = 100 > X = matrix(rnorm(n*n),n,n) > > X1 = t(X) %*% X - .001*diag(n) # will be pos def > > L1 = t(chol(X1)) > > L1_lk = rlibkriging:::linalg_chol_safe(X1) > > test_that(desc="Chol R/lK is almost the same", expect_equal(L1, L1_lk, tol=1e-10)) Test passed 😸 > > X2 = t(X) %*% X - .01*diag(n) # will not be pos def > > L2 = NULL > try(L2 <- t(chol(X2))) Error in chol.default(X2) : the leading minor of order 100 is not positive > > test_that(desc="Chol without nugget does not pass", expect_null(L2)) Test passed 🎊 > > rlibkriging:::linalg_set_num_nugget(0.001) > rlibkriging:::linalg_set_chol_warning(TRUE) > L2_lk = NULL > L2_lk = rlibkriging:::linalg_chol_safe(X2) [WARNING] Added 6 numerical nugget to force Cholesky decomposition > rlibkriging:::linalg_set_num_nugget(1e-15) # default setup > > test_that(desc="Chol with nugget passes", expect_false(is.null(L2_lk))) Test passed 😸 > > test_that(desc="Chol with nugget is close to be good", + expect_equal(L2_lk %*% t(L2_lk) - 6*0.001*diag(n), X2, tol=1e-10)) Test passed 🥳 > > ####################################################################### > > rlibkriging:::linalg_set_num_nugget(0) # force NA when chol fails > > rcond_chol = function(X) { + rlibkriging:::linalg_rcond_chol(chol(X)) + } > rcond_approx_chol = function(X) { + rlibkriging:::linalg_rcond_approx_chol(chol(X)) + } > rcond = function(X) { + 1/kappa(X, norm='1') + } > > set.seed(123) > r = array(NA, 100) > r_chol = array(NA,length(r)) > ra_chol = array(NA,length(r)) > for (i in 1:length(r)) { + n = 100 + X = matrix(rnorm(n*n),n,n) + X = t(X) %*% X + n*rnorm(1,0,1)*diag(n) + try(r[i] <- rcond((X))) + try(r_chol[i] <- rcond_chol((X))) + try(ra_chol[i] <- rcond_approx_chol((X))) + } Error in chol.default(X) : the leading minor of order 30 is not positive Error in chol.default(X) : the leading minor of order 30 is not positive Error in chol.default(X) : the leading minor of order 10 is not positive Error in chol.default(X) : the leading minor of order 10 is not positive Error in chol.default(X) : the leading minor of order 55 is not positive Error in chol.default(X) : the leading minor of order 55 is not positive Error in chol.default(X) : the leading minor of order 45 is not positive Error in chol.default(X) : the leading minor of order 45 is not positive Error in chol.default(X) : the leading minor of order 2 is not positive Error in chol.default(X) : the leading minor of order 2 is not positive Error in chol.default(X) : the leading minor of order 1 is not positive Error in chol.default(X) : the leading minor of order 1 is not positive Error in chol.default(X) : the leading minor of order 1 is not positive Error in chol.default(X) : the leading minor of order 1 is not positive Error in chol.default(X) : the leading minor of order 72 is not positive Error in chol.default(X) : the leading minor of order 72 is not positive Error in chol.default(X) : the leading minor of order 2 is not positive Error in chol.default(X) : the leading minor of order 2 is not positive Error in chol.default(X) : the leading minor of order 1 is not positive Error in chol.default(X) : the leading minor of order 1 is not positive Error in chol.default(X) : the leading minor of order 14 is not positive Error in chol.default(X) : the leading minor of order 14 is not positive Error in chol.default(X) : the leading minor of order 52 is not positive Error in chol.default(X) : the leading minor of order 52 is not positive Error in chol.default(X) : the leading minor of order 1 is not positive Error in chol.default(X) : the leading minor of order 1 is not positive Error in chol.default(X) : the leading minor of order 1 is not positive Error in chol.default(X) : the leading minor of order 1 is not positive Error in chol.default(X) : the leading minor of order 11 is not positive Error in chol.default(X) : the leading minor of order 11 is not positive Error in chol.default(X) : the leading minor of order 1 is not positive Error in chol.default(X) : the leading minor of order 1 is not positive Error in chol.default(X) : the leading minor of order 1 is not positive Error in chol.default(X) : the leading minor of order 1 is not positive Error in chol.default(X) : the leading minor of order 3 is not positive Error in chol.default(X) : the leading minor of order 3 is not positive Error in chol.default(X) : the leading minor of order 1 is not positive Error in chol.default(X) : the leading minor of order 1 is not positive Error in chol.default(X) : the leading minor of order 33 is not positive Error in chol.default(X) : the leading minor of order 33 is not positive Error in chol.default(X) : the leading minor of order 19 is not positive Error in chol.default(X) : the leading minor of order 19 is not positive Error in chol.default(X) : the leading minor of order 2 is not positive Error in chol.default(X) : the leading minor of order 2 is not positive Error in chol.default(X) : the leading minor of order 75 is not positive Error in chol.default(X) : the leading minor of order 75 is not positive Error in chol.default(X) : the leading minor of order 1 is not positive Error in chol.default(X) : the leading minor of order 1 is not positive Error in chol.default(X) : the leading minor of order 42 is not positive Error in chol.default(X) : the leading minor of order 42 is not positive Error in chol.default(X) : the leading minor of order 48 is not positive Error in chol.default(X) : the leading minor of order 48 is not positive Error in chol.default(X) : the leading minor of order 44 is not positive Error in chol.default(X) : the leading minor of order 44 is not positive Error in chol.default(X) : the leading minor of order 50 is not positive Error in chol.default(X) : the leading minor of order 50 is not positive Error in chol.default(X) : the leading minor of order 1 is not positive Error in chol.default(X) : the leading minor of order 1 is not positive Error in chol.default(X) : the leading minor of order 1 is not positive Error in chol.default(X) : the leading minor of order 1 is not positive Error in chol.default(X) : the leading minor of order 1 is not positive Error in chol.default(X) : the leading minor of order 1 is not positive Error in chol.default(X) : the leading minor of order 44 is not positive Error in chol.default(X) : the leading minor of order 44 is not positive Error in chol.default(X) : the leading minor of order 1 is not positive Error in chol.default(X) : the leading minor of order 1 is not positive Error in chol.default(X) : the leading minor of order 5 is not positive Error in chol.default(X) : the leading minor of order 5 is not positive Error in chol.default(X) : the leading minor of order 54 is not positive Error in chol.default(X) : the leading minor of order 54 is not positive Error in chol.default(X) : the leading minor of order 1 is not positive Error in chol.default(X) : the leading minor of order 1 is not positive Error in chol.default(X) : the leading minor of order 11 is not positive Error in chol.default(X) : the leading minor of order 11 is not positive Error in chol.default(X) : the leading minor of order 16 is not positive Error in chol.default(X) : the leading minor of order 16 is not positive Error in chol.default(X) : the leading minor of order 1 is not positive Error in chol.default(X) : the leading minor of order 1 is not positive Error in chol.default(X) : the leading minor of order 28 is not positive Error in chol.default(X) : the leading minor of order 28 is not positive Error in chol.default(X) : the leading minor of order 5 is not positive Error in chol.default(X) : the leading minor of order 5 is not positive Error in chol.default(X) : the leading minor of order 75 is not positive Error in chol.default(X) : the leading minor of order 75 is not positive Error in chol.default(X) : the leading minor of order 1 is not positive Error in chol.default(X) : the leading minor of order 1 is not positive Error in chol.default(X) : the leading minor of order 47 is not positive Error in chol.default(X) : the leading minor of order 47 is not positive Error in chol.default(X) : the leading minor of order 1 is not positive Error in chol.default(X) : the leading minor of order 1 is not positive Error in chol.default(X) : the leading minor of order 7 is not positive Error in chol.default(X) : the leading minor of order 7 is not positive Error in chol.default(X) : the leading minor of order 16 is not positive Error in chol.default(X) : the leading minor of order 16 is not positive Error in chol.default(X) : the leading minor of order 1 is not positive Error in chol.default(X) : the leading minor of order 1 is not positive Error in chol.default(X) : the leading minor of order 6 is not positive Error in chol.default(X) : the leading minor of order 6 is not positive Error in chol.default(X) : the leading minor of order 1 is not positive Error in chol.default(X) : the leading minor of order 1 is not positive Error in chol.default(X) : the leading minor of order 2 is not positive Error in chol.default(X) : the leading minor of order 2 is not positive Error in chol.default(X) : the leading minor of order 1 is not positive Error in chol.default(X) : the leading minor of order 1 is not positive Error in chol.default(X) : the leading minor of order 3 is not positive Error in chol.default(X) : the leading minor of order 3 is not positive Error in chol.default(X) : the leading minor of order 2 is not positive Error in chol.default(X) : the leading minor of order 2 is not positive > > plot(log10(r),log10(r_chol),col='red',ylim=range(c(range(log10(r_chol),na.rm=T),range(log10(ra_chol),na.rm=T)))) > points(log10(r),log10(ra_chol),col='orange') > abline(a=0,b=1) > > rlibkriging:::linalg_set_num_nugget(1e-15) # default setup > > ####################################################################### > > n = 1000 > R = matrix(rnorm(n*n),n,n) > R = t(R) %*% R > > ## we just want a covar (sym pos def) matrix > no = n-30 > Roo = R[1:no,1:no] > Too = t(chol(Roo)) > > test_that(desc="Full chol equals incremented chol by block", + expect_equal(t(chol(R)), rlibkriging:::linalg_chol_block(R,Too), tol=1e-9)) Test passed 🥳 > > rlibkriging:::linalg_check_chol_rcond(default_rcond_checked) > rlibkriging:::linalg_set_num_nugget(default_num_nugget) > > proc.time() user system elapsed 3.48 0.35 3.90