Gridap.Arrays

This module provides:

• An extension of the AbstractArray interface in order to properly deal with mutable caches.
• A mechanism to generate lazy arrays resulting from operations between arrays.
• A collection of concrete implementations of AbstractArray.

The exported names in this module are:

Exported names

AddEntriesMap, AppendedArray, ArrayBlock, ArrayBlockView, AutoDiffMap, Broadcasting, CachedArray, CachedMatrix, CachedVector, CompressedArray, ConfigMap, DualizeMap, FilterMap, IdentityVector, InverseMap, KeyToValMap, LazyArray, LocalPenaltySolveMap, LocalSolveMap, Map, MatrixBlock, MatrixBlockView, MulAddMap, Operation, PosNegPartition, PosNegReindex, Reindex, Table, TouchEntriesMap, TreeNode, UNSET, VectorBlock, VectorBlockView, VectorWithEntryInserted, VectorWithEntryRemoved, append_ptrs, append_ptrs!, append_tables_globally, append_tables_locally, array_cache, autodiff_array_gradient, autodiff_array_hessian, autodiff_array_jacobian, collect1d, dataiterator, datarange, dataview, empty_table, evaluate, evaluate!, find_inverse_index_map, find_inverse_index_map!, find_local_index, flatten_partition, generate_data_and_ptrs, get_array, get_data_eltype, get_local_item, get_ptrs_eltype, getindex!, identity_table, invalidate_cache!, inverse_map, lazy_append, lazy_map, lazy_split, length_to_ptrs!, pair_arrays, print_op_tree, return_cache, return_type, return_value, rewind_ptrs!, setsize!, similar_tree_node, test_array, test_map, testargs, testitem, testvalue, unpair_arrays, ∑,

Abstract type representing a function (mapping) that provides a cache and an in-place evaluation for performance. This is the type to be used in the lazy_map function.

Derived types must implement the following method:

• evaluate!(cache,k,x...)

and optionally these ones:

• return_cache(k,x...)
• return_type(k,x...)

The mapping interface can be tested with the test_map function.

Note that most of the functionality implemented in terms of this interface relies in duck typing. That is, it is not strictly needed to work with types that inherit from Map. This is specially useful in order to accommodate existing types into this framework without the need to implement a wrapper type that inherits from Map. For instance, a default implementation is available for Function objects. However, we recommend that new types inherit from Map.

Example

using Gridap.Arrays
using LinearAlgebra

# Define a custom Map for k(a,b,c) = α*a*b + β*c
struct MyMulAdd{T} <: Map
  α::T
  β::T
end

# Optional: provide a cache to avoid allocations in evaluate!
function Gridap.Arrays.return_cache(k::MyMulAdd, a, b, c)
  d = a*b + c # Initial allocation
  CachedArray(d)
end

# Mandatory: in-place evaluation
function Gridap.Arrays.evaluate!(cache, k::MyMulAdd, a, b, c)
  setsize!(cache, size(c))
  d = cache.array
  copyto!(d, c)
  mul!(d, a, b, k.α, k.β)
  d
end

k = MyMulAdd(2.0, 1.0)
a = rand(3,3); b = rand(3,3); c = rand(3,3)
res = evaluate(k, a, b, c)

return_cache(f,x...)

Returns the cache needed to lazy_map mapping f with arguments of the same type as the objects in x. This function returns nothing by default, i.e., no cache.

evaluate!(
  transient_space::FESpace,
  space::TransientTrialFESpace, t::Real
) -> FESpace

Replace the Dirichlet values of the space by those at time t.

Implementation of return_cache for a array of Field.

If the field vector has length nf and it is evaluated in one point, it returns an nf vector with the result. If the same array is applied to a vector of np points, it returns a matrix np x nf.

evaluate!(cache,f,x...)

Applies the mapping f at the arguments x... using the scratch data provided in the given cache object. The cache object is built with the return_cache function using arguments of the same type as in x. In general, the returned value y can share some part of its state with the cache object. If the result of two or more calls to this function need to be accessed simultaneously (e.g., in multi-threading), create and use several cache objects (e.g., one cache per thread).

evaluate(space::TransientTrialFESpace, t::Real) -> FESpace

Allocate a transient space and evaluate the Dirichlet values at time t.

evaluate(space::TransientTrialFESpace, t::Nothing) -> FESpace

Evaluating at nothing means that the Dirichlet values are not important.

evaluate(f,x...)

evaluates the mapping f at the arguments in x by creating a temporary cache internally. This functions is equivalent to

cache = return_cache(f,x...)
evaluate!(cache,f,x...)

test_map(y,f,x...;cmp=(==))

Function used to test if the mapping f has been implemented correctly. f is a Map sub-type, x is a tuple in the domain of the mapping and y is the expected result. Function cmp is used to compare the computed result with the expected one. The checks are done with the @test macro.

return_type(f,x...)

Returns the type of the result of calling mapping f with arguments of the types of the objects x.

return_value(f,x...)

Return a variable of the type of the image fx=f(x...) (possibly fx itself).

testargs(f,x...)

The default implementation of this function is testargs(f,x...) = x. One can overload it in order to use lazy_map with 0-length array and maps with non-trivial domains.

Broadcasting(f)

Returns a mapping that represents the "broadcasted" version of the function f.

Example

using Gridap.Arrays

a = [3,2]
b = [2,1]

bm = Broadcasting(+)

c = evaluate(bm,a,b)

println(c)

# output
[5, 3]

Operation(op)

Returns the map that results after applying an operation f over a set of map(s) args. That is Operation(f)(args)(x...) is formally defined as f(map(a->a(x...),args)...).

Example

using Gridap.Arrays

fa(x) = x.*x
fb(x) = sqrt.(x)

x = collect(0:5)

fab = Operation(fa)(fb)
c = evaluate(fab,x)

println(c)

# output
[0.0, 1.0, 2.0, 3.0, 4.0, 5.0]

struct InverseMap{F} <: Map

Map for the inverse of the Function or Map F.

Reindex(values) -> Map

struct representing a binary partition of a range of indices

Using this allows one to do a number of important optimizations when working with PosNegReindex

PosNegReindex(values_pos,values_neg)

struct AddEntriesMap{F} <: Map

Map for add_entries! where F is the combine function.

LocalPenaltySolveMap(; factorize! = lu!, pivot = RowMaximum())

A map for solving local constrained linear systems, relying on a factorization method.

Given a left-hand-side 2x2 block matrix matrixmat and a set of 2xN right-hand-side arrays lhs, returns an N-Tuple of arrays containing the solutions to the linear systems.

Each system is given by A*[x_i; λ_i] = b_i, where A = [App, Aλp; Apλ, 0] is the augmented matrix, and b_i = [Bp; Bλ] is the right-hand side vector. The solution is computed using a penalty method, as x_i = ldiv!(factorize!(C,pivot),d_i) with C = App + μT * Apλ * Aλp and d_i = Bp + μT * Apλ * Bλ, where μT is a penalty parameter. The penalty parameter μT is heuristically chosen as μT = norm(App)/norm(Apλ*Aλp).

LocalSolveMap(; factorize! = lu!, pivot = RowMaximum())

A map for solving local linear systems, relying on a factorization method.

Given a left-hand-side matrix mat and a set of N right-hand-side arrays lhs, returns an N-Tuple of arrays containing the solutions to the linear systems defined by

Each system is given by A*x_i = b_i, and the solution is computed as x_i = ldiv!(factorize!(A,pivot),b_i)

struct MulAddMap{T} <: Map
  α::T
  β::T
end

Map for the operation a,b,c -> LinearAlgepra.mul!(copy(c),a,b,α,β)

struct TouchEntriesMap <: Map

Dummy AddEntriesMap used to keep track of sparsity patterns in symbolic assemblies, see the symbolic methods of SparseMatrixAssembler.

struct KeyToValMap{T<:Function} <: Map

Map for lazily filling a Dict the outputs of the function T over an array of inputs.

array_cache(a::AbstractArray)

Returns a cache object to be used in the getindex! function. It defaults to

array_cache(a::T) where T = nothing

for types T such that uses_hash(T) == Val(false), and

function array_cache(a::T) where T
  hash = Dict{UInt,Any}()
  array_cache(hash,a)
end

for types T such that uses_hash(T) == Val(true), see the uses_hash function. In the later case, the type T should implement the following signature:

array_cache(hash::Dict,a::AbstractArray)

where we pass a dictionary (i.e., a hash table) in the first argument. This hash table can be used to test if the object a has already built a cache and re-use it as follows

id = objectid(a)
if haskey(hash,id)
  cache = hash[id] # Reuse cache
else
  cache = ... # Build a new cache depending on your needs
  hash[id] = cache # Register the cache in the hash table
end

This mechanism is needed, e.g., to re-use intermediate results in complex lazy operation trees. In multi-threading computations, a different hash table per thread has to be used in order to avoid race conditions.

getindex!(cache, a::AbstractArray, i...)

Returns the item of the array a associated with index i by (possibly) using the scratch data passed in the cache object. The cache is only used for single element access (i::Union{Integer,CartesianIndex,end,begin}...), otherwise getindex is used.

getindex! defaults to

getindex!(cache, a::AbstractArray, i...) = a[i...]

As for standard Julia arrays, the user needs to implement only one of the following signatures depending on the IndexStyle of the array.

getindex!(cache,a::AbstractArray,i::Integer)
getindex!(cache,a::AbstractArray{T,N},i::Vararg{Integer,N}) where {T,N}

See also array_cache, uses_hash, invalidate_cache!.

Examples

Iterating over an array using the getindex! function

using Gridap.Arrays

a = collect(10:15)

cache = array_cache(a)
for i in eachindex(a)
  ai = getindex!(cache,a,i)
  println("$i -> $ai")
end

# output
1 -> 10
2 -> 11
3 -> 12
4 -> 13
5 -> 14
6 -> 15

testitem(a::AbstractArray{T}) -> Any

Returns an arbitrary instance of eltype(a). The default returned value is the first entry in the array if length(a)>0 and testvalue(eltype(a)) if length(a)==0 See the testvalue function.

Examples

using Gridap.Arrays

a = collect(3:10)
ai = testitem(a)

b = Int[]
bi = testitem(b)

(ai, bi)

# output
(3, 0)

testvalue(::Type{T}) where T

Returns an arbitrary instance of type T. It defaults to zero(T) for non-array types and to an empty array for array types. It can be overloaded for new types T if zero(T) does not makes sense. This function is used to compute testitem for 0-length arrays.

test_array(
  a::AbstractArray{T,N}, b::AbstractArray{S,N},cmp=(==)) where {T,S,N}

Checks if the entries in a and b are equal using the comparison function cmp. It also stresses the new methods added to the AbstractArray interface.

get_array(a::AbstractArray)

Returns a.

get_array(a::CellDatum)

Get the raw array of a datas defined in the physical space.

uses_hash(::Type{<:AbstractArray})

This function is used to specify if the type T uses the hash-based mechanism to reuse caches. It should return either Val(true) or Val(false). It defaults to

uses_hash(::Type{<:AbstractArray}) = Val(false)

Once this function is defined for the type T it can also be called on instances of T.

invalidate_cache!(cache) -> nothing

Prevent possible memoized value(s) to be re-used, in order to force any future evaluation to be computed (this is aimed to be used on lazy array caches).

Alias for Base.sum.

struct LazyArray{G,T,N,F} <: AbstractArray{T,N}

Subtype of AbstractArray which is the result of lazy_map. It represents the result of lazy_mapping an (array of) Map to a set of arrays that contain the mapping arguments. This struct makes use of the cache provided by the mapping in order to compute its indices (thus allowing to prevent allocation). The array is lazy, i.e., the values are only computed on demand. It extends the AbstractArray API with three methods:

• array_cache(a::AbstractArray)
• getindex!(cache,a::AbstractArray,i...)
• invalidate_cache!(cache)

G is the type of the array of map (a Fill array if a single map is given to lazy_array), F that of the collection of argument arrays. The size of the LazyArray is size(G).

lazy_map(f,::Type{T},a::AbstractArray...) where T

Like lazy_map(f,a::AbstractArray...), but the user provides the element type of the resulting array in order to circumvent type inference.

lazy_map(f,a::AbstractArray...) -> AbstractArray

Applies the Map (or Function) f to the entries of the arrays in a (see the definition of Map).

The resulting array r is such that r[i] equals to evaluate(f,ai...) where ai is the tuple containing the i-th entry of the arrays in a (see function evaluate for more details). In other words, the resulting array is numerically equivalent to:

map( (x...)->evaluate(f,x...), a...)

Examples

Using a function as mapping

using Gridap.Arrays

a = collect(0:5)
b = collect(10:15)

c = lazy_map(+,a,b)

println(c)

# output
[10, 12, 14, 16, 18, 20]

Using a user-defined mapping

using Gridap.Arrays
import Gridap.Arrays: evaluate!

a = collect(0:5)
b = collect(10:15)

struct MySum <: Map end

evaluate!(cache,::MySum,x,y) = x + y

k = MySum()

c = lazy_map(k,a,b)

println(c)

# output
[10, 12, 14, 16, 18, 20]

struct CompressedArray{T,N,A,P} <: AbstractArray{T,N}
  values::A
  ptrs::P
end

Type representing an array with a reduced set of values. The array is represented by a short array of values, namely the field values, and a large array of indices, namely the field ptrs. The i-th component of the resulting array is defined as values[ptrs[i]].

Example

using Gridap.Arrays
values = [10, 20]
ptrs = [1, 2, 1, 1]
a = CompressedArray(values, ptrs)
# a == [10, 20, 10, 10]

CompressedArray(values::AbstractArray,ptrs::AbstractArray)

Creates a CompressedArray object by the given arrays of values and ptrs.

struct Table{T,Vd<:AbstractVector{T},Vp<:AbstractVector} <: AbstractVector{Vector{T}}
  data::Vd
  ptrs::Vp
end

Type representing a list of lists (i.e., a table) in compressed format. This is particularly efficient for representing connectivity in meshes.

Example

using Gridap.Arrays
vv = [[1,2,3], [10,20], [5]]
table = Table(vv)
table[1] # returns [1, 2, 3]
table.data # returns [1, 2, 3, 10, 20, 5]
table.ptrs # returns [1, 4, 6, 7]

Table(a::AbstractArray{<:AbstractArray})

Build a table from a vector of vectors. If the inputs are multidimensional arrays instead of vectors, they are flattened.

const UNSET = 0

append_ptrs!(pa,pb)

Similar to append_ptrs, but appends pb at the end of pa, in place in pa.

append_ptrs(pa,pb)

Concatenate two vectors of pointers in a new vector.

append_tables_globally(tables::Table...)

Concatenate multiple tables by appending their rows.

Example

t1 = Table([[1,2], [3]])
t2 = Table([[10], [20,30]])
t3 = append_tables_globally(t1, t2)
# returns Table([[1,2], [3], [10], [20,30]])

append_tables_locally(tables::Table...)
append_tables_locally(offsets, tables)

Concatenate tables of the same length by merging their rows. The optional offsets can be used to shift the data values.

Example

t1 = Table([[1,2], [3]])
t2 = Table([[10], [20]])
t3 = append_tables_locally(t1, t2)
# returns Table([[1,2,10], [3,20]])

t4 = append_tables_locally((0,10), (t1, t2))
# returns Table([[1,2,20], [3,30]])

block_identity_array(::Type{T},ptrs) where T

Given a vector of pointers of length n+1, returns a vector of length ptrs[end]-1 where the entries are the index of the block to which each entry belongs.

Example

julia> block_identity_array([1,3,7])

6-element Vector{Int64}:
 1
 1
 2
 2
 2
 2

collect1d(a)

Equivalent to

[a[i] for in 1:length(a)]

dataiterator(a::Table)

Iterate over the entries of a returning the triplets (i,j,v) where

• i is the outer index,
• j is the local inner index, and
• v is the value a[i][j].

Example

julia> t = Table([[4.,7.],[8.],[9.,2.,1.]])

julia> x = collect(dataiterator(t))

julia> x
6-element Vector{Tuple{Int64, Int64, Float64}}:
 (1, 1, 4.)
 (1, 2, 7.)
 (2, 1, 8.)
 (3, 1, 9.)
 (3, 2, 2.)
 (3, 3, 1.)

datarange(a::Table,i::Integer)
datarange(a::Table,ids::UnitRange{<:Integer})

Given a Table and an index or a range of indices, return the corresponding range of indices in the underlying data array. Similar to nzrange for sparse matrices, it allows for convenient iteration over a table:

t = Table([[4,7],[8],[9,20,1]])
for i in eachindex(t)
  for k in datarange(t,i)
    val = t.data[k]
    # stuff ... 
  end
end

dataview(a::Table, i::Integer)
dataview(a::Table, ids::UnitRange{<:Integer})

Given a Table and an index or a range of indices, return a view of the corresponding entries in the underlying data array.

Equivalent to view(a.data, datarange(a,i)).

empty_table(::Type{T},::Type{P}, l::Integer) where {T,P}
empty_table(l::Integer)

find_inverse_index_map(a_to_b, nb=maximum(a_to_b))

Given a vector of indices a_to_b, returns the inverse index map b_to_a.

find_inverse_index_map!(b_to_a, a_to_b)

In place find_inverse_index_map.

find_local_index(a_to_b, b_to_la_to_a) -> a_to_la
find_local_index(c_to_a, c_to_b, b_to_la_to_a) -> c_to_lc

find_local_nbor_index(a_to_b, a_to_lb_to_b) -> a_to_lb
find_local_nbor_index(a_to_b, a_to_c, c_to_lb_to_b) -> a_to_lb

flatten_partition(a_to_bs::Table [, nb=maximum(a_to_bs.data)])

Flattens a partition a_to_bs (map from a to multiple bs) into an array b_to_a (map from b to its unique a).

Example

a_to_bs = Table([[1,2,3],[7,8],[4,5,6]])
b_to_a = flatten_partition(a_to_bs)
# b_to_a == [1, 1, 1, 3, 3, 3, 2, 2]

data, ptrs = generate_data_and_ptrs(vv)

Given a vector of vectors, compress it and return the corresponding data and and ptrs

get_local_item(a::Table,li::Integer)

View in the liᵗʰ column of a (the liᵗʰ items in each list/row of a).

inverse_table(a_to_lb_to_b::Table [, nb=maximum(a_to_lb_to_b.data)])
inverse_table(a_to_b::AbstractVector [, nb=maximum(a_to_b)])

Returns the inverse of the input Table or non-injective array of integers, as a Table.

local_identity_array(::Type{T}, ptrs) where T

Given a vector of pointers of length n+1, returns a vector of length ptrs[end]-1 where the entries are the local index of the entry within the block it belongs to.

Example

julia> local_identity_array([1,3,7])

6-element Vector{Int64}:
 1
 2
 1
 2
 3
 4

merge_entries(a_to_lb_to_b, c_to_la_to_a) -> c_to_lb_to_b

Merge the entries of a_to_lb_to_b, grouping them by c_to_la_to_a. Returns the merged table c_to_lb_to_b.

Accepts the following keyword arguments:

• acc: Accumulator for the entries of a_to_lb_to_b. Default to a Set, ensuring that the resulting entries are unique.
• post: Postprocessing function to apply to the accumulator before storing the resulting entries. Defaults to the identity, but can be used to perform local sorts or filters, for example.

remove_empty_entries!(table::Table)

Given a Table, remove the entries that are empty by modifying its ptrs in-place.

mutable struct CachedArray{T, N, A<:AbstractArray{T, N}} <: AbstractArray{T, N}

Type providing a re-sizable array.

The size of a CachedArray is changed via the setsize! function.

A CachedArray can be build with the constructors

• CachedArray(a::AbstractArray)
• CachedArray(T,N)

using Gridap.Arrays
# Create an empty CachedArray
a = CachedArray(Float64,2)
# Resize to new shape (2,3)
setsize!(a,(2,3))
size(a)
# output
(2, 3)

CachedArray(T,N)

Constructs an empty CachedArray of element type T and N dimensions.

CachedArray(a::AbstractArray)

Constructs a CachedArray from a given array.

const CachedMatrix{T,A} = CachedArray{T,2,A}

CachedMatrix(a)

CachedMatrix(T)

const CachedVector{T,A} = CachedArray{T,1,A}

CachedVector(a)

CachedVector(T)

setsize!(a, s)

Changes the size of the CachedArray a to the size described the the tuple s. After calling setsize!, the array can store uninitialized values.

setsize_op!(op, c, args)

Sets the size of the CachedArray to accomodate the result of the operation given by op(args...).

struct AppendedArray{T,A,B} <: AbstractVector{T}

Type for a lazily appended array.

struct VectorWithEntryRemoved{T,A} <: AbstractVector{T}

mutable struct VectorWithEntryInserted{T,A} <: AbstractVector{T}

struct ArrayBlock{A,N}
  array::Array{A,N}
  touched::Array{Bool,N}
end

Block-wise storage of arrays, where touched indicates which blocks are active. Accessing a non-touched block returns nothing.

ArrayBlock is mostly used for multi-field and skeleton computations and assembly, where each block corresponds to a field or plus/minus respectively. Blocks might be nested in the case of multi-field computations on skeletons.

struct ArrayBlockView{A,N,M}

Low level container implementing view on ArrayBlock.

struct BlockBroadcasting{F} <: Map
  f::F
end

A BlockBroadcasting broadcasts a map f block-wise over the input arrays.

struct BlockMap{N} <: Map
    size::NTuple{N,Int}
    indices::Vector{CartesianIndex{N}}
end

A BlockMap maps M = length(indices) arrays to a single array of N-dimensinal blocks, where only the blocks indexed by indices are touched and contain the corresponding entry of the input arrays.

Constructors:

BlockMap(l::Integer,i::Integer) ≡ BlockMap((l,),[CartesianIndex((i,))])
BlockMap(l::Integer,inds::Vector{<:Integer}) ≡ BlockMap((l,),[CartesianIndex((i,)) for i in inds])
BlockMap(s::NTuple,i::Integer) ≡ BlockMap(s,[CartesianIndex(i)])
BlockMap(s::NTuple,inds::Vector{<:NTuple}) ≡ BlockMap(s,[CartesianIndex(i) for i in inds])

Usage:

  lazy_map(BlockMap(2,[1,2]),a,b)

const MatrixBlock = ArrayBlock{A,2}

Alias, see also ArrayBlock.

const MatrixBlockView{A} = ArrayBlockView{A,2,2}

Alias, see also ArrayBlockView.

struct MergeBlockMap{N,M} <: Map
    size::NTuple{N,Int}
    indices::Vector{Vector{Tuple{CartesianIndex{M},CartesianIndex{N}}}}
end

A MergeBlockMap create a single array of N-dimensional blocks from L=length(indices) input arrays of M-dimensional blocks.

For the l-th input array a_l, the vector of tuples in indices[l] contains the mapping between the indices of the blocks in a_l and the indices of the blocks in the output array, i.e a_out[P[2]] = a_l[P[1]] ∀ P ∈ indices[l], ∀ l.

const VectorBlock = ArrayBlock{A,1}

Alias, see also ArrayBlock.

const VectorBlockView{A} = ArrayBlockView{A,1,1}

Alias, see also ArrayBlockView.

struct AutoDiffMap <: Map

struct ConfigMap{F,T} <: Map

Map for ForwardDiff.[F]Config(T,...) where T is tag function and F is either gradient or jacobian.

struct DualizeMap <: Map

print_op_tree(a,args...;kwargs...)
print_op_tree(io::IO,a,args...;showid=false,kwargs...)

Print the operation tree of a lazy operation/array, in stdout if io is not given.

similar_tree_node(n::TreeNode,node)