R Under development (unstable) (2024-10-16 r87241 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.

> ## for R_DEFAULT_PACKAGES=NULL :
> library(stats)
> 
> library(Matrix)
> 
> ## Matrix Exponential
> source(system.file("test-tools.R", package = "Matrix"))
Loading required package: tools
> 
> ## e ^ 0 = 1  - for matrices:
> assert.EQ.mat(expm(Matrix(0, 3,3)), diag(3), tol = 0)# exactly
> ## e ^ diag(.) = diag(e ^ .):
> assert.EQ.mat(expm(as(diag(-1:4), "generalMatrix")), diag(exp(-1:4)))
> set.seed(1)
> rE <- replicate(100,
+             { x <- rlnorm(12)
+               relErr(as(expm(as(diag(x), "generalMatrix")),
+                         "matrix"),
+                      diag(exp(x))) })
> stopifnot(mean(rE) < 1e-15,
+           max(rE)  < 1e-14)
> summary(rE)
     Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
1.369e-17 1.931e-16 4.233e-16 5.338e-16 7.604e-16 2.247e-15 
> 
> ## Some small matrices
> 
> m1 <- Matrix(c(1,0,1,1), ncol = 2)
> e1 <- expm(m1)
> assert.EQ.mat(e1, cbind(c(exp(1),0), exp(1)))
> (p1 <- pack(m1))
2 x 2 Matrix of class "dtpMatrix"
     [,1] [,2]
[1,]    1    1
[2,]    .    1
> stopifnot(exprs = { 
+     is(p1, "dtpMatrix")
+     all.equal(pack(e1), expm(p1), tolerance = 2e-15)# failed in Matrix 1.6.1
+ })
> m2 <- Matrix(c(-49, -64, 24, 31), ncol = 2)
> e2 <- expm(m2)
> ## The true matrix exponential is 'te2':
> e_1 <-  exp(-1)
> e_17 <- exp(-17)
> te2 <- rbind(c(3*e_17 - 2*e_1, -3/2*e_17 + 3/2*e_1),
+              c(4*e_17 - 4*e_1, -2  *e_17 + 3  *e_1))
> assert.EQ.mat(e2, te2, tol = 1e-13)
> ## See the (average relative) difference:
> all.equal(as(e2,"matrix"), te2, tolerance = 0) # 2.22e-14 {was 1.48e-14} on "lynne"
[1] "Mean relative difference: 2.221749e-14"
> 
> (dsp <- pack(crossprod(matrix(-2:3, 2,3))))
3 x 3 Matrix of class "dspMatrix"
     [,1] [,2] [,3]
[1,]    5   -1   -7
[2,]   -1    1    3
[3,]   -7    3   13
> stopifnot(all(abs(expm(dsp) - expm(as.matrix(dsp))) <= 0.5)) # failed badly in Matrix 1.6.1
> 
> ## The ``surprising identity''      det(exp(A)) == exp( tr(A) )
> ## or                           log det(exp(A)) == tr(A) :
> stopifnot(all.equal(c(determinant(e2)$modulus), sum(diag(m2))))
> 
> ## a very simple nilpotent one:
> (m3 <- Matrix(cbind(0,rbind(6*diag(3),0))))# sparse
4 x 4 sparse Matrix of class "dtCMatrix"
            
[1,] . 6 . .
[2,] . . 6 .
[3,] . . . 6
[4,] . . . .
> stopifnot(all(m3 %*% m3 %*% m3 %*% m3 == 0))# <-- m3 "^" 4 == 0
> e3 <- expm(m3)
> E3 <- expm(Matrix(m3, sparse=FALSE))
> s3 <- symmpart(m3) # dsCMatrix
> es3 <- expm(s3)
> e3. <- rbind(c(1,6,18,36),
+ 	     c(0,1, 6,18),
+ 	     c(0,0, 1, 6),
+ 	     c(0,0, 0, 1))
> stopifnot(is(e3, "triangularMatrix"),
+ 	  is(es3, "symmetricMatrix"),
+ 	  identical(e3, E3),
+ 	  identical(as.mat(e3), e3.),
+ 	  all.equal(as(es3,"generalMatrix"),
+ 		    expm(as(s3,"generalMatrix")))
+ 	  )
> 
> 
> ## This used to be wrong {bug in octave-origin code}:
> M6 <- Matrix(c(0, -2, 0, 0, 0, 0,
+               10, 0, 0, 0,10,-2,
+               0,  0, 0, 0,-2, 0,
+               0, 10,-2,-2,-2,10,
+               0,  0, 0, 0, 0, 0,
+               10, 0, 0, 0, 0, 0), 6, 6)
> 
> exp.M6 <- expm(M6)
> as(exp.M6, "sparseMatrix")# prints a bit more nicely
6 x 6 sparse Matrix of class "dgCMatrix"
                                                      
[1,]  0.5397837 -1.4227470  .   8.2517610 .   2.647087
[2,] -0.5294175  0.5397837  . -10.0751692 .  -4.069834
[3,]  .          .          1  -0.8646647 .   .       
[4,]  .          .          .   0.1353353 .   .       
[5,] -4.0698344  2.6470874 -2 -26.3584982 1 -18.048090
[6,]  0.8139669 -0.5294175  .   9.6491573 .   4.609618
> stopifnot(all.equal(t(exp.M6),
+ 		    expm(t(M6)), tolerance = 1e-12),
+           all.equal(exp.M6[,3], c(0,0,1,0,-2,0), tolerance = 1e-12),
+           all.equal(exp.M6[,5], c(0,0,0,0, 1,0), tolerance = 1e-12),
+           all(exp.M6[3:4, c(1:2,5:6)] == 0)
+           )
> 
> cat('Time elapsed: ', proc.time(),'\n') # for ``statistical reasons''
Time elapsed:  1.57 0.15 1.71 NA NA 
> 
> proc.time()
   user  system elapsed 
   1.57    0.15    1.71