# jmatrixsc

#### Juan Domingo

# Please, read this

This is a copy of the vignette of package `jmatrix`. It is included here since `jmatrix` is underlying this package and you will need to know how to prepare your data to be processed by `scellpam`. But you must NOT load package `jmatrix`; all the functions detailed below are already included into `scellpam` and are available just by loading it with `library(scellpam)`.

    library(scellpam)

# Purpose

The package `jmatrix` (Domingo (2022a)) was originally conceived as a tool for other packages, namely `parallelpam` (Domingo (2022b)) and `scellpam` (Domingo (2022c)) which needed to deal with very big matrices which might not fit in the memory of the computer, particularly if their elements are of type `double` (in most modern machines, 8 bytes per element) whereas they could fit if they were matrices of other data types, in particular of floats (4 bytes per element).

Unfortunately, R is not a strongly typed language. Double is the default type in R and it is not easy to work with other data types. Trials like the package float (Schmidt (2022)) have been done, but to use them you have to coerce a matrix already loaded in R memory to a float matrix, and then you can delete it. But, what happens if you computer has not memory enough to hold the matrix in the first place?. This is the problem this package tries to address.

Our idea is to use the disk as temporarily storage of the matrix in a file with a internal binary format (`jmatrix` format). This format has a header of 128 bytes with information like type of matrix (full, sparse or symmetric), data type of each element (char, short, int, long, float, double or long double), number of rows and columns and endianness; then comes the content as binary data (in sparse matrices zeros are not stored; in symmetric matrices only the lower-diagonal is stored) and finally the metadata (currently, names for rows/columns if needed and an optional comment).

Such files are created and loaded by functions written in C++ which are accessible from `R` with Rcpp (Eddelbuettel and François (2011)). The file, once loaded, uses strictly the needed memory for its data type and can be processed by other C++ functions (like the PAM algorithm or any other numeric library written in C++) also from inside `R`.

The matrix contained in a binary data file in `jmatrix` format cannot be loaded directly in `R` memory as a `R` matrix (that would be impossible, anyway, since precisely this package is done for the cases in which such matrix would NOT fit into the available RAM). Nevertheless, limited access through some functions is provided to read one or more rows or one or more columns as `R` vectors or matrices (obviously, coerced to double).

The package `jmatrix` must not be considered as a final, finished software. Currently is mostly an instrumental solution to address our needs and we make available as a separate package just in case it could be useful for anyone else.

# Workflow

## Debug messages

First of all, the package can show quite informative (but sometimes verbose) messages in the console. To turn on/off such messages you can use.

    # Initially, state of debug is FALSE. Turn it on exclusively for
    # the jmatrix part of scellpam with
    ScellpamSetDebug(FALSE,debparpam=FALSE,debjmat=TRUE)
    #> Debugging for jmatrix inside scellpam package set to ON.

## Data storage

As stated before, the binary matrix files should normally be created from C++ getting the data from an external source like a data file in a format used in bioinformatics or a .csv file. These files should be read by chunks. As an example, look at function `CsvToJMat` in package `scellpam` (Domingo (2022c)).

As a convenience and only for testing purposes (to be used in this vignette), we provide the function `JWriteBin` to write a R matrix as a `jmatrix` file.

    # Create a 6x8 matrix of random values
    Rf <- matrix(runif(48),nrow=6)
    # Set row and column names for it
    rownames(Rf) <- c("A","B","C","D","E","F")
    colnames(Rf) <- c("a","b","c","d","e","f","g","h")
    # Let's see the matrix
    Rf
    #>           a         b          c          d          e         f          g
    #> A 0.9727312 0.1426776 0.20763667 0.71883028 0.96708605 0.8791411 0.77231302
    #> B 0.2177657 0.1512234 0.66483899 0.71715484 0.98111583 0.8965128 0.43341403
    #> C 0.5399618 0.9992887 0.57262550 0.32554859 0.07954279 0.8596215 0.31397562
    #> D 0.2914114 0.1869364 0.67900596 0.03732338 0.23105631 0.9887593 0.01935518
    #> E 0.3118885 0.7125586 0.06934446 0.21246234 0.10425043 0.2309819 0.04643620
    #> F 0.2968369 0.9738831 0.18602511 0.57473660 0.92499299 0.6678535 0.16290861
    #>            h
    #> A 0.68539205
    #> B 0.58463073
    #> C 0.22709614
    #> D 0.02380487
    #> E 0.12273420
    #> F 0.83060626
    # and write it as the binary file Rfullfloat.bin
    JWriteBin(Rf,"Rfullfloat.bin",dtype="float",dmtype="full",
              comment="Full matrix of floats")
    #> The passed matrix has row names for the 6 rows and they will be used.
    #> The passed matrix has column names for the 8 columns and they will be used.
    #> Writing binary matrix Rfullfloat.bin of (6x8)
    #> End of block of binary data at offset 320
    #>    Writing row names (6 strings written, from A to F).
    #>    Writing column names (8 strings written, from a to h).
    #>    Writing comment: Full matrix of floats
    # Also, you can write it with double data type:
    JWriteBin(Rf,"Rfulldouble.bin",dtype="double",dmtype="full",
              comment="Full matrix of doubles")
    #> The passed matrix has row names for the 6 rows and they will be used.
    #> The passed matrix has column names for the 8 columns and they will be used.
    #> Writing binary matrix Rfulldouble.bin of (6x8)
    #> End of block of binary data at offset 512
    #>    Writing row names (6 strings written, from A to F).
    #>    Writing column names (8 strings written, from a to h).
    #>    Writing comment: Full matrix of doubles

To get information about the stored file the function `JMatInfo` is provided. Of course, this funcion does not read the complete file in memory but just the header.

    # Information about the float binary file
    JMatInfo("Rfullfloat.bin")
    #> File:               Rfullfloat.bin
    #> Matrix type:        FullMatrix
    #> Number of elements: 48
    #> Data type:          float
    #> Endianness:         little endian (same as this machine)
    #> Number of rows:     6
    #> Number of columns:  8
    #> Metadata:           Stored names of rows and columns.
    #> Metadata comment:  "Full matrix of floats"
    # Same information about the double binary file
    JMatInfo("Rfulldouble.bin")
    #> File:               Rfulldouble.bin
    #> Matrix type:        FullMatrix
    #> Number of elements: 48
    #> Data type:          double
    #> Endianness:         little endian (same as this machine)
    #> Number of rows:     6
    #> Number of columns:  8
    #> Metadata:           Stored names of rows and columns.
    #> Metadata comment:  "Full matrix of doubles"

## Data load

As stated before, no function is provided to read the whole matrix in memory which would contradict the philosophy of this package, but you can get rows or columns from a file.

    # Reads row 1 into vector vf. Float values inside the file are
    # promoted to double.
    (vf<-GetJRow("Rfullfloat.bin",1))
    #>         a         b         c         d         e         f         g         h
    #> 0.9727312 0.1426776 0.2076367 0.7188303 0.9670861 0.8791412 0.7723130 0.6853920

Obviously, storage in float provokes a loosing of precision. We have observed this not to be relevant for `PAM` (partitioning around medoids) algorihm but it can be important in other cases. It is the price to pay for halving the needed space.

    # Checks the precision lost
    max(abs(Rf[1,]-vf))
    #> [1] 2.747402e-08

Nevertheless, storing as double obviously keeps the data intact.

    vd<-GetJRow("Rfulldouble.bin",1)
    max(abs(Rf[1,]-vd))
    #> [1] 0

Now, let us see examples of some functions to read rows or columns by number or by name, or to read several rows/columns as a R matrix. In all examples numbers for rows and columns are in R-convention (i.e. starting at 1)

    # Read column number 3
    (vf<-GetJCol("Rfullfloat.bin",3))
    #>          A          B          C          D          E          F
    #> 0.20763667 0.66483897 0.57262552 0.67900598 0.06934446 0.18602511
    # Test precision
    max(abs(Rf[,3]-vf))
    #> [1] 2.421439e-08
    # Read row with name C
    (vf<-GetJRowByName("Rfullfloat.bin","C"))
    #>          a          b          c          d          e          f          g
    #> 0.53996181 0.99928868 0.57262552 0.32554859 0.07954279 0.85962147 0.31397563
    #>          h
    #> 0.22709614
    # Read column with name c
    (vf<-GetJColByName("Rfullfloat.bin","c"))
    #>          A          B          C          D          E          F
    #> 0.20763667 0.66483897 0.57262552 0.67900598 0.06934446 0.18602511
    # Get the names of all rows or columns as vectors of R strings
    (rn<-GetJRowNames("Rfullfloat.bin"))
    #> [1] "A" "B" "C" "D" "E" "F"
    (cn<-GetJColNames("Rfullfloat.bin"))
    #> [1] "a" "b" "c" "d" "e" "f" "g" "h"
    # Get the names of rows and columns simultaneosuly as a list of two elements
    (l<-GetJNames("Rfullfloat.bin"))
    #> $rownames
    #> [1] "A" "B" "C" "D" "E" "F"
    #>
    #> $colnames
    #> [1] "a" "b" "c" "d" "e" "f" "g" "h"
    # Get several rows at once. The returned matrix has the rows in the
    # same order as the passed list,
    # and this list can contain even repeated values
    (vm<-GetJManyRows("Rfullfloat.bin",c(1,4)))
    #>           a         b         c          d         e         f          g
    #> A 0.9727312 0.1426776 0.2076367 0.71883029 0.9670861 0.8791412 0.77231300
    #> D 0.2914114 0.1869364 0.6790060 0.03732338 0.2310563 0.9887593 0.01935518
    #>            h
    #> A 0.68539202
    #> D 0.02380487

    # Of course, columns can be extrated equally
    (vc<-GetJManyCols("Rfulldouble.bin",c(1,4)))
    #>           a          d
    #> A 0.9727312 0.71883028
    #> B 0.2177657 0.71715484
    #> C 0.5399618 0.32554859
    #> D 0.2914114 0.03732338
    #> E 0.3118885 0.21246234
    #> F 0.2968369 0.57473660
    # and similar functions are provided for extracting by names:
    (vm<-GetJManyRowsByNames("Rfulldouble.bin",c("A","D")))
    #>           a         b         c          d         e         f          g
    #> A 0.9727312 0.1426776 0.2076367 0.71883028 0.9670861 0.8791411 0.77231302
    #> D 0.2914114 0.1869364 0.6790060 0.03732338 0.2310563 0.9887593 0.01935518
    #>            h
    #> A 0.68539205
    #> D 0.02380487
    (vc<-GetJManyColsByNames("Rfulldouble.bin",c("a","d")))
    #>           a          d
    #> A 0.9727312 0.71883028
    #> B 0.2177657 0.71715484
    #> C 0.5399618 0.32554859
    #> D 0.2914114 0.03732338
    #> E 0.3118885 0.21246234
    #> F 0.2968369 0.57473660

The package can manage and store sparse and symmetric matrices, too.

    # Generation of a 6x8 sparse matrix
    Rsp <- matrix(rep(0,48),nrow=6)
    sparsity <- 0.1
    nnz <- round(48*sparsity)
    where <- floor(47*runif(nnz))
    val <- runif(nnz)
    for (i in 1:nnz)
    {
     Rsp[floor(where[i]/8)+1,(where[i]%%8)+1] <- val[i]
    }
    rownames(Rsp) <- c("A","B","C","D","E","F")
    colnames(Rsp) <- c("a","b","c","d","e","f","g","h")
    # Let's see the matrix
    Rsp
    #>           a b         c d         e f          g h
    #> A 0.5836277 0 0.0000000 0 0.0000000 0 0.00000000 0
    #> B 0.0000000 0 0.0000000 0 0.1568182 0 0.00000000 0
    #> C 0.0000000 0 0.9733851 0 0.0000000 0 0.00000000 0
    #> D 0.0000000 0 0.0000000 0 0.0000000 0 0.63194867 0
    #> E 0.0000000 0 0.0000000 0 0.0000000 0 0.07297231 0
    #> F 0.0000000 0 0.0000000 0 0.0000000 0 0.00000000 0
    # Write the matrix as sparse with type float
    JWriteBin(Rsp,"Rspafloat.bin",dtype="float",dmtype="sparse",
              comment="Sparse matrix of floats")
    #> The passed matrix has row names for the 6 rows and they will be used.
    #> The passed matrix has column names for the 8 columns and they will be used.
    #> Writing binary matrix Rspafloat.bin of (6x8)
    #> End of block of binary data at offset 192
    #>    Writing row names (6 strings written, from A to F).
    #>    Writing column names (8 strings written, from a to h).
    #>    Writing comment: Sparse matrix of floats

Notice that the condition of being a sparse matrix and the storage space used can be known with the matrix info.

    JMatInfo("Rspafloat.bin")
    #> File:               Rspafloat.bin
    #> Matrix type:        SparseMatrix
    #> Number of elements: 48
    #> Data type:          float
    #> Endianness:         little endian (same as this machine)
    #> Number of rows:     6
    #> Number of columns:  8
    #> Metadata:           Stored names of rows and columns.
    #> Metadata comment:  "Sparse matrix of floats"
    #> Binary data size:   64 bytes, which is 33.3333% of the full matrix size.

Be careful: trying to store as sparse a matrix which is not (it has not a majority of 0-entries) works, but produces a matrix larger than the corresponding full matrix.

With respect to symmetric matrices, `JWriteBin` works the same way. Let us generate a \(7 \times 7\) symmetric matrix.

    Rns <- matrix(runif(49),nrow=7)
    Rsym <- 0.5*(Rns+t(Rns))
    rownames(Rsym) <- c("A","B","C","D","E","F","G")
    colnames(Rsym) <- c("a","b","c","d","e","f","g")
    # Let's see the matrix
    Rsym
    #>            a         b         c          d         e          f         g
    #> A 0.28295469 0.3894980 0.3792033 0.09565392 0.3079891 0.40672747 0.4723853
    #> B 0.38949803 0.9787828 0.7560180 0.45880852 0.2931062 0.58432786 0.3476512
    #> C 0.37920326 0.7560180 0.9645777 0.46229634 0.5114315 0.55892538 0.5522704
    #> D 0.09565392 0.4588085 0.4622963 0.22508186 0.1883446 0.56806378 0.5361935
    #> E 0.30798907 0.2931062 0.5114315 0.18834461 0.9611357 0.68329620 0.2144955
    #> F 0.40672747 0.5843279 0.5589254 0.56806378 0.6832962 0.05507042 0.4192074
    #> G 0.47238535 0.3476512 0.5522704 0.53619345 0.2144955 0.41920743 0.9873061
    # Write the matrix as symmetric with type float
    JWriteBin(Rsym,"Rsymfloat.bin",dtype="float",dmtype="symmetric",
              comment="Symmetric matrix of floats")
    #> The passed matrix has row names for the 7 rows and they will be used.
    #> Writing binary matrix Rsymfloat.bin
    #> End of block of binary data at offset 240
    #>    Writing row names (7 strings written, from A to G).
    #>    Writing comment: Symmetric matrix of floats
    # Get the information
    JMatInfo("Rsymfloat.bin")
    #> File:               Rsymfloat.bin
    #> Matrix type:        SymmetricMatrix
    #> Number of elements: 49 (28 really stored)
    #> Data type:          float
    #> Endianness:         little endian (same as this machine)
    #> Number of rows:     7
    #> Number of columns:  7
    #> Metadata:           Stored only names of rows.
    #> Metadata comment:  "Symmetric matrix of floats"

Notice that if you store a R matrix which is NOT symmetric as a symmetric `jmatrix`, only the lower triangular part (including the main diagonal) will be saved. The upper-triangular part will be lost.

The functions to read rows/colums stated before works equally independently of the matrix character (full, sparse or symmetric) so you can play with them using the `Rspafloat.bin` and `Rsymfloat.bin` file to check they work.

Finally, there is a function to get the lower-diagonal part of a symmetric matrix stored in a `jmatrix` file as a R vector, ordered by columns (i.e.: column 1 from M(2,1) to M(n,1), followed by column 2 from M(3,2) to M(n,2) and so on, up to M(n,n-1) so d is a vector of \(n(n-1)/2\) components). This is to get the data of a dissimilarity/distance matrix in the format used by the Partitioning around medoids (PAM) algorithm implemented in the `cluster` (Maechler et al. (2022)) package. The purpose is to compare such implementation with our own implementation of PAM. If you are interested, see the documentation of the `parallelpam` (Domingo (2022b)) package.

    (d<-GetSubdiag("Rsymfloat.bin"))
    #>  [1] 0.38949803 0.37920326 0.09565391 0.30798906 0.40672746 0.47238535
    #>  [7] 0.75601798 0.45880851 0.29310617 0.58432788 0.34765115 0.46229634
    #> [13] 0.51143146 0.55892539 0.55227035 0.18834461 0.56806380 0.53619343
    #> [19] 0.68329620 0.21449548 0.41920742

Domingo, Juan. 2022a. Jmatrix: Read from/Write in Disks Matrices with Any Data Type in a Binary Format.
———. 2022b. Parallelpam: Applies the Partitioning-Around-Medoids (PAM) Clustering Algortihm to Big Sets of Data Using Parallel Implementation, If Several Cores Are Available.
———. 2022c. Scellpam: Applying Partitioning Around Medoids to Single Cell with High Number of Cells.
Eddelbuettel, Dirk, and Romain François. 2011. “Rcpp: Seamless R and C++ Integration.” Journal of Statistical Software 40 (8): 1–18. https://doi.org/10.18637/jss.v040.i08.
Maechler, Martin, Peter Rousseeuw, Anja Struyf, Mia Hubert, and Kurt Hornik. 2022. Cluster: Cluster Analysis Basics and Extensions. https://CRAN.R-project.org/package=cluster.
Schmidt, Drew. 2022. “float: 32-Bit Floats.” https://cran.r-project.org/package=float.