btoolbox: random toolbox

context("spectral decompositions of adjacency matrices")

WIN_32 = .Platform$OS.type == "windows" && .Machine$sizeof.pointer != 8

test_that("results correspond to R's dense counterparts", {

    skip_if(WIN_32, "skipping large matrix test for windows 32-bit")
    
    library("Matrix")
    
    keep_largest = function(x, n) {
        
        # magnitude (for cut off)
        om = order(-abs(x$values))
        x$values = x$values[om[1:n]]
        x$vectors = x$vectors[, om[1:n]]
        
        # algebraic (for sorting)
        oa = order(-x$values)
        x$values = x$values[oa]
        x$vectors = x$vectors[, oa]
        
        return(x)
        
    }
    
    keep_smallest = function(x, n) {
        
        # magnitude (for cut off)
        om = order(abs(x$values))
        x$values = x$values[om[1:n]]
        x$vectors = x$vectors[, om[1:n]]
        
        # algebraic (for sorting)
        oa = order(x$values)
        x$values = x$values[oa]
        x$vectors = x$vectors[, oa]
        
        return(x)
        
    }
    
    
    r = 5
    k = 20
    m = 10
    
    nullmat = Matrix(nrow = k, ncol = k, data = 0L, sparse = TRUE)
    
    for (rr in 1:r) {
        
        # adjacency matrix without isolates
        rsums = 0
        while (any(rsums == 0)) {
            
            x = lapply(1:m, function(w) sample(1:k, sample(2:6, 1), replace = F))
            A = update_affiliation(nullmat, x)
            rsums = Matrix::rowSums(A)
            
        }
        
        # normalized adjacency
        dsqrt_inv = diag(1/sqrt(rsums))
        d_inv = dsqrt_inv^2
        
        A_left = d_inv %*% A
        A_sym = dsqrt_inv %*% A %*% dsqrt_inv
        
        # Laplacian
        L = diag(rsums) - A
        
        # normalized Laplacian
        L_left = diag(k) - A_left
        L_sym = diag(k) - A_sym
        
        # check
        p = 3
        r_testmat = keep_largest(eigen(A), k - 1)
        c_testmat = sp_eigen_adj(A, k - 1, "adj")
        
        for (q in 1:p) {
            expect_equal(r_testmat$values[p], c_testmat$values[p])
            expect_equal(
                abs(r_testmat$vectors[, p]), abs(c_testmat$vectors[, p])
            )
        }
        
        r_testmat = keep_largest(eigen(A_left), k - 1)
        c_testmat = sp_eigen_adj(A, k - 1, "adj_left")
        
        for (q in 1:p) {
            expect_equal(r_testmat$values[p], c_testmat$values[p])
            expect_equal(
                abs(r_testmat$vectors[, p]), abs(c_testmat$vectors[, p])
            )
        }
        
        r_testmat = keep_largest(eigen(A_sym), k - 1)
        c_testmat = sp_eigen_adj(A, k - 1, "adj_sym")
        
        for (q in 1:p) {
            expect_equal(r_testmat$values[p], c_testmat$values[p])
            expect_equal(
                abs(r_testmat$vectors[, p]), abs(c_testmat$vectors[, p])
            )
        }
        
        r_testmat = keep_smallest(eigen(L), k - 1)
        c_testmat = sp_eigen_adj(A, k - 1, "lap")
        
        for (q in 1:p) {
            expect_equal(r_testmat$values[p], c_testmat$values[p])
            expect_equal(
                abs(r_testmat$vectors[, p]), abs(c_testmat$vectors[, p])
            )
        }
        
        r_testmat = keep_smallest(eigen(L_left), k - 1)
        c_testmat = sp_eigen_adj(A, k - 1, "lap_left")
        
        for (q in 1:p) {
            expect_equal(r_testmat$values[p], c_testmat$values[p])
            expect_equal(
                abs(r_testmat$vectors[, p]), abs(c_testmat$vectors[, p])
            )
        }
        
        r_testmat = keep_smallest(eigen(L_sym), k - 1)
        c_testmat = sp_eigen_adj(A, k - 1, "lap_sym")
        
        for (q in 1:p) {
            expect_equal(r_testmat$values[p], c_testmat$values[p])
            expect_equal(
                abs(r_testmat$vectors[, p]), abs(c_testmat$vectors[, p])
            )
        }
        
    }
    
})


test_that("function throws appropriate errors", {
    
    skip_if(WIN_32, "skipping large matrix test for windows 32-bit")

    library("Matrix")
    
    r = 5
    k = 20
    m = 10
    
    nullmat = Matrix::Matrix(nrow = k, ncol = k, data = 0L, sparse = TRUE)
    
    for (rr in 1:r) {
        
        # adjacency matrix with isolates
        rsums = 1
        while (all(rsums > 0)) {
            
            x = lapply(1:m, function(w) sample(1:k, sample(2:6, 1), replace = F))
            A = update_affiliation(nullmat, x)
            rsums = Matrix::rowSums(A)
            
        }
        
        for (cc in c("adj", "adj_sym", "lap", "lap_sym")) {
            
            expect_error(sp_eigen_adj(A, 7, cc))
            expect_error(sp_eigen_adj(A, 7, cc, check = "none"), NA)
            expect_error(sp_eigen_adj(A, 7, cc, check = "symmetric"), NA)
            
        }
        
        
    }
    
    for (rr in 1:r) {
        
        # adjacency matrix without isolates
        rsums = 0
        while (any(rsums == 0)) {
            
            x = lapply(1:m, function(w) sample(1:k, sample(2:6, 1), replace = F))
            A = update_affiliation(nullmat, x)
            rsums = Matrix::rowSums(A)
            
        }
        
        # induce asymmetry
        i = sample(1:k, 1)
        j = i
        while (j == i)
            j = sample(1:k, 1)
        
        A[i, j] = A[i, j] + 1
        
        for (cc in c("adj", "adj_sym", "lap", "lap_sym")) {
            
            expect_error(sp_eigen_adj(A, 7, cc))
            expect_error(sp_eigen_adj(A, 7, cc, "none"), NA)
            expect_error(sp_eigen_adj(A, 7, cc, "isolates"), NA)
        }
        
    }
    
})


test_that("Laplacian has non-negative eigenvalues", {

    skip_if(WIN_32, "skipping large matrix test for windows 32-bit")
    
    library("Matrix")
    
    r = 5
    k = 20
    m = 10
    
    nullmat = Matrix::Matrix(nrow = k, ncol = k, data = 0L, sparse = TRUE)

    for (rr in 1:r) {
        
        # adjacency matrix without isolates
        rsums = 0
        while (any(rsums == 0)) {
            
            x = lapply(1:m, function(w) sample(1:k, sample(2:6, 1), replace = F))
            A = update_affiliation(nullmat, x)
            rsums = Matrix::rowSums(A)
            
        }
        
        evals = sp_eigen_adj(A, k - 1, "lap")$values
        evals_sym = sp_eigen_adj(A, k - 1, "lap_sym")$values
        
        pos_or_zero = function(w, tol) 
            all(ifelse(w > 0.0, T, ifelse(abs(w) < tol, T, F)))
        
        expect_true(pos_or_zero(evals, 1e-10))
        expect_true(pos_or_zero(evals_sym, 1e-10))
        
    }

})


test_that("function passes matrix by value", {
    
    skip_if(WIN_32, "skipping large matrix test for windows 32-bit")

    library("Matrix")
    
    r = 5
    k = 20
    m = 10
    
    nullmat = Matrix::Matrix(nrow = k, ncol = k, data = 0L, sparse = TRUE)
    
    for (rr in 1:r) {
        
        # adjacency matrix without isolates
        rsums = 0
        while (any(rsums == 0)) {
            
            x = lapply(1:m, function(w) sample(1:k, sample(2:6, 1), replace = F))
            A = update_affiliation(nullmat, x)
            rsums = Matrix::rowSums(A)
            
        }
        
        A_old = A
        
        for (cc in c("adj", "adj_sym", "lap", "lap_sym")) {
            
            B = sp_eigen_adj(A, 7, cc)
            expect_identical(A, A_old)
            
        }
        
    }
    
})