Is abind console output incorrect?

Question

I am trying to combine multiple 2d vectors to obtain a 3dvector, the 3d vector should be [rows, columns, 2dvectors]. So far I have only worked with 3d vectors made from a only one 2d vector:

x <- matrix(1:12, 3,4)

#    [,1] [,2] [,3] [,4]
#[1,]  1     4   7   10
#[2,]  2     5   8   11
#[3,]  3     6   9   12

dim(x) <- c(dim(x), 1))
dim(x)
#[1] 3 4 1

And when I select a row shows this:

x[1,,]
#[1] 1 4 7 10

Now to combine 2d vectors into a 3d vector I have done as follow:

x <- matrix(1:12, 3, 4)
y <- x + 100
z <- abind(x, y, along=3)
dim(z)
#[1] 3 4 2

As you can see the dimensions are correct [rows, columns, 2dvectors]. Therefore I expect this:

z[1,,]
#     [,1] [,2] [,3] [,4]
#[1,]  1    4    7    10
#[2,]  101 104  107  110

But instead I get this:

z[1,,]
#     [,1] [,2]
#[1,]  1   101
#[2,]  4   104
#[3,]  7   107
#[4,]  10  110

So I want to know if is the output displaying correctly at the console or if there is something I am not understanding.

Answer 1

There might be several confusing aspects that lead to the unexpected output:

1. A matrix in R is restricted to two-dimensional matrices, and in fact is.matrix(...) returns TRUE if dim(...) <= 2. Matrices with higher dimensions are actually called array. Consider

x <- matrix(1:12, 3,4)

class(x) # x is matrix, array
# [1] "matrix" "array" 

dim(x) <- c(dim(x), 1))
dim(x)
# [1] 3 4 1

class(x) # x is array
# [1] "array" 

2. Extract or replace parts of an object (that is, slicing of a matrix) with bracket-notation like x[1,,] involves a default parameter drop = TRUE.

For matrices and arrays, if drop = TRUE, then the result is coerced to the lowest possible dimension.

When drop = FALSE, function drop(...) will NOT be invoked at all, that is x[1,,,drop=FALSE] will has its dim(...) as c(1,4,1), and console output turns into:

x[1,,,drop=FALSE]
# , , 1
#
#      [,1] [,2] [,3] [,4]
# [1,]    1    4    7   10

Then the abind output appears more consistent:

y <- x + 100
z <- abind(x, y, along=3)

dim(z)
# [1] 3 4 2

z[1,,,drop=FALSE]

# , , 1
# 
#      [,1] [,2] [,3] [,4]
# [1,]    1    4    7   10
# 
# , , 2
#
#      [,1] [,2] [,3] [,4]
# [1,]  101  104  107  110

3. The way drop(...) works can be confusing too because it

Delete the dimensions of an array which have only one level.

that is, for x[1,,,drop=FALSE] with dimensions c(1, 4, 1), drop(...) will delete its first and third dimension, reducing dimensions to c(4). But will return NULL instead, because dim(...) is defined to return NULL for a vector:

drop(x[1,,,drop=FALSE])
# [1]  1  4  7 10

class(drop(x[1,,,drop=FALSE]))
# [1] "integer"

dim(drop(x[1,,,drop=FALSE]))
# NULL

For z[1,,,drop=FALSE] with dimensions c(1,4,2), drop(...) will reduce its dimensions to c(4,2), like:

drop(z[1,,,drop=FALSE])

#      [,1] [,2]
# [1,]    1  101
# [2,]    4  104
# [3,]    7  107
# [4,]   10  110

How to make drop(...) more intuitive?

An alternative way of implementing drop(...) might be more intuitive -- delete the heading dimensions of an array which have only one level. That is for an array with dimension like c(1,1,3,1,2,1), our drop2(...) will reduce into c(3,1,2,1). As a crude implementation, consider

dummy = array(1:6, dim=c(1,1,3,1,2,1))

drop2 = function(x) {
    x.dim = dim(x)
    if (is.null(x.dim)) return(x)
    array(x, dim=x.dim[cumsum(x.dim > 1) > 0])
}

dim(drop2(dummy))
# [1] 3 1 2 1

drop2(x[1,,,drop=FALSE]) # resulting dim c(4, 1)
#      [,1]
# [1,]    1
# [2,]    4
# [3,]    7
# [4,]   10

drop2(z[1,,,drop=FALSE])
#      [,1] [,2]
# [1,]    1  101
# [2,]    4  104
# [3,]    7  107
# [4,]   10  110