How do I reinterpret and reshape a multi-dimensional array `J` into a two-dimensional array `j` such that `typeof(j) <: Matrix` is true in Julia 0.7

General Usage
1

When upgrading code from version 0.6.4 to 0.7.0 I came across the following issue:

Julia 0.6.4

J = rand(4,10,12)
j = reinterpret(Float64,J,(4*10,12))
typeof(j) <: Matrix  # Evaluates to true

Julia 0.7.0

J = rand(4,10,12)
j = reshape(reinterpret(Float64,vec(J)),(4*10,12))
typeof(j) <: Matrix #Evaluates to false

Why is this a problem?

This became an issue in a portion of my code where I wanted to reshape a Jacobian matrix into a multi-dimensional array. The reshaping facilitated convenient indexing into the elements of the Jacobian matrix whilst I was looping over the data that is required to construct the entries of the matrix.

I am using LsqFit.jl and the required Jacobian matrix is defined to be of type Matrix, i.e. jacobian::Matrix{T} (https://github.com/JuliaNLSolvers/LsqFit.jl/blob/v0.6.0/src/curve_fit.jl). I wanted to construct my Jacobian as a multi-dimensional array, and then subsequently reshape it into a matrix (thereby avoiding reallocation/copy).

In the latest versions of Julia, the reshaped and reinterpreted array was no longer a subtype of Matrix which broke my code.

Thanks for some great suggestions on the Julia Slack channel I was able to resolve my problem. I will post the solution below.

2

Kristoffer Carlsson suggested an unsanctioned way of getting around this by using unsafe_wrap. As an example he showed:

julia> a = rand(5,5,5);

julia> unsafe_wrap(Vector{Float64}, pointer(a), (15,))
15-element Array{Float64,1}:
 0.39667876246883926
 0.22245479132465285

with the understanding that one would have to ensure that a is protected from the garbage collector.
He also linked to a related discussion on reshaping StaticArrays.

https://github.com/JuliaArrays/StaticArrays.jl/pull/496#issuecomment-421806480

Moritz Schauer suggested a clever alternative and safer approach which he attributed to Keno Fisher and dubbed Keno’s obvious trick:

Fill the reshaped view of a full matrix of the shape you need, instead of reshaping the full matrix with a view into the shape you need.

This means that I start by declaring a Jacobian matrix, e.g. J = zeros(4*10,12) and then create a reshaped view of the Jacobian matrix: Jv = reshape(reinterpret(Float64,J), 4, 10, 12).

I can then conveniently index into Jv. Since Jv is a view, I can pass J into LsqFit.jl because it contains the same data as Jv and is also a Matrix.

3

Why not use

j = reshape(J, (4*10,12))
typeof(j) <: Matrix #Evaluates to true

There’s no need for reinterpret as the elements are going to be Float64 still.

Note that a reshaped array in Julia 0.7 is similar to a view, in that it shares its data with the original array.

5 Likes