text:
"""
    AbstractDimStack

Abstract supertype for dimensional stacks.

These have multiple layers of data, but share dimensions.
"""
abstract type AbstractDimStack{L,D} end

data(s::AbstractDimStack) = s.data
dims(s::AbstractDimStack) = s.dims
dims(s::AbstractDimStack, key::Symbol) = dims(s, layerdims(s)[key])
refdims(s::AbstractDimStack) = s.refdims
layerdims(s::AbstractDimStack) = s.layerdims
layerdims(s::AbstractDimStack, key::Symbol) = layerdims(s)[key]
metadata(s::AbstractDimStack) = s.metadata
layermetadata(s::AbstractDimStack) = s.layermetadata
layermetadata(s::AbstractDimStack, key::Symbol) = layermetadata(s)[key]

layerdims(A::AbstractDimArray) = basedims(A)

@inline Base.keys(s::AbstractDimStack) = keys(data(s))
Base.values(s::AbstractDimStack) = values(dimarrays(s))
Base.first(s::AbstractDimStack) = s[first(keys(s))]
Base.last(s::AbstractDimStack) = s[last(keys(s))]
# Only compare data and dim - metadata and refdims can be different
Base.:(==)(s1::AbstractDimStack, s2::AbstractDimStack) =
    data(s1) == data(s2) && dims(s1) == dims(s2) && layerdims(s1) == layerdims(s2)
Base.length(s::AbstractDimStack) = length(data(s))
Base.iterate(s::AbstractDimStack, args...) = iterate(dimarrays(s), args...)

rebuild(s::AbstractDimStack, data, dims=dims(s), refdims=refdims(s),
        layerdims=layerdims(s), metadata=metadata(s), layermetadata=layermetadata(s)) =
    basetypeof(s)(data, dims, refdims, layerdims, metadata, layermetadata)

function rebuildsliced(f::Function, s::AbstractDimStack, data, I)
    layerdims = map(basedims, data)
    dims, refdims = slicedims(s, I)
    rebuild(s; data=map(parent, data), dims=dims, refdims=refdims, layerdims=layerdims)
end

function dimarrays(s::AbstractDimStack{<:NamedTuple{Keys}}) where Keys
    NamedTuple{Keys}(map(K -> s[K], Keys))
end

Adapt.adapt_structure(to, s::AbstractDimStack) = map(A -> adapt(to, A), s)

# Dipatch on Tuple of Dimension, and map
for func in (:index, :mode, :metadata, :sampling, :span, :bounds, :locus, :order)
    @eval ($func)(s::AbstractDimStack, args...) = ($func)(dims(s), args...)
end

"""
    Base.map(f, s::AbstractDimStack)

Apply functrion `f` to each layer of the stack `s`, and rebuild it.

If `f` returns `DimArray`s the result will be another `DimStack`.
Other values will be returned in a `NamedTuple`.
"""
Base.map(f, s::AbstractDimStack) = maybestack(map(f, dimarrays(s)))

maybestack(As::NamedTuple{<:Any,<:Tuple{Vararg{<:AbstractDimArray}}}) = DimStack(As)
maybestack(x::NamedTuple) = x

"""
    Base.copy!(dst::AbstractDimStack, src::AbstractDimStack, [keys=keys(dst)])

Copy all or a subset of layers from one stack to another.

## Example

Copy just the `:sea_surface_temp` and `:humidity` layers from `src` to `dst`.

```julia
copy!(dst::AbstractDimStack, src::AbstractDimStack, keys=(:sea_surface_temp, :humidity))
```
"""
function Base.copy!(dst::AbstractDimStack, src::AbstractDimStack, keys=keys(dst))
    # Check all keys first so we don't copy anything if there is any error
    for key in keys
        key in Base.keys(dst) || throw(ArgumentError("key $key not found in dest keys"))
        key in Base.keys(src) || throw(ArgumentError("key $key not found in source keys"))
    end
    for key in keys
        copy!(dst[key], src[key])
    end
end

# Array methods

# Methods with no arguments that return a DimStack
for (mod, fnames) in
    (:Base => (:inv, :adjoint, :transpose), :LinearAlgebra => (:Transpose,))
    for fname in fnames
        @eval ($mod.$fname)(s::AbstractDimStack) = map(A -> ($mod.$fname)(A), s)
    end
end

# Methods with an argument that return a DimStack
for fname in (:rotl90, :rotr90, :rot180, :PermutedDimsArray, :permutedims)
    @eval (Base.$fname)(s::AbstractDimStack, args...) =
        map(A -> (Base.$fname)(A, args...), s)
end

# Methods with keyword arguments that return a DimStack
for (mod, fnames) in
    (:Base => (:sum, :prod, :maximum, :minimum, :extrema, :dropdims),
     :Statistics => (:cor, :cov, :mean, :median, :std, :var))
    for fname in fnames
        @eval ($mod.$fname)(s::AbstractDimStack; kw...) =
            maybestack(map(A -> ($mod.$fname)(A; kw...), dimarrays(s)))
    end
end

# Methods that take a function
for (mod, fnames) in (:Base => (:reduce, :sum, :prod, :maximum, :minimum, :extrema),
                      :Statistics => (:mean,))
    for fname in fnames
        _fname = Symbol(:_, fname)
        @eval begin
            ($mod.$fname)(f::Function, s::AbstractDimStack; dims=Colon()) =
                ($_fname)(f, s, dims)
            # Colon returns a NamedTuple
            ($_fname)(f::Function, s::AbstractDimStack, dims::Colon) =
                map(A -> ($mod.$fname)(f, A), data(s))
            # Otherwise return a DimStack
            ($_fname)(f::Function, s::AbstractDimStack, dims) =
                map(A -> ($mod.$fname)(f, A; dims=dims), s)
        end
    end
end

"""
    DimStack <: AbstractDimStack

    DimStack(data::AbstractDimArray...)
    DimStack(data::Tuple{Vararg{<:AbstractDimArray}})
    DimStack(data::NamedTuple{Keys,Vararg{<:AbstractDimArray}})
    DimStack(data::NamedTuple, dims::DimTuple; metadata=NoMetadata())

DimStack holds multiple objects sharing some dimensions, in a `NamedTuple`.
Indexing operates as for [`AbstractDimArray`](@ref), except it occurs for all
data layers of the stack simulataneously. Layer objects can hold values of any type.

DimStack can be constructed from multiple `AbstractDimArray` or a `NamedTuple`
of `AbstractArray` and a matching `dims` `Tuple`. If `AbstractDimArray`s have
the same name they will be given the name `:layer1`, substitiuting the
layer number for `1`.

`getindex` with `Int` or `Dimension`s or `Selector`s that resolve to `Int` will
return a `NamedTuple` of values from each layer in the stack. This has very good
performace, and usually takes less time than the sum of indexing each array
separately.

Indexing with a `Vector` or `Colon` will return another `DimStack` where
all data layers have been sliced.  `setindex!` must pass a `Tuple` or `NamedTuple` maching
the layers.

Most `Base` and `Statistics` methods that apply to `AbstractArray` can be used on
all layers of the stack simulataneously. The result is a `DimStack`, or
a `NamedTuple` if methods like `mean` are used without `dims` arguments, and
return a single non-array value.

## Example

```jldoctest
julia> using DimensionalData

julia> A = [1.0 2.0 3.0; 4.0 5.0 6.0];

julia> dimz = (X([:a, :b]), Y(10.0:10.0:30.0))
(X{Vector{Symbol}, AutoMode{AutoOrder}, NoMetadata}([:a, :b], AutoMode{AutoOrder}(AutoOrder()), NoMetadata()), Y{StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}}, AutoMode{AutoOrder}, NoMetadata}(10.0:10.0:30.0, AutoMode{AutoOrder}(AutoOrder()), NoMetadata()))

julia> da1 = DimArray(1A, dimz, :one);

julia> da2 = DimArray(2A, dimz, :two);

julia> da3 = DimArray(3A, dimz, :three);

julia> s = DimStack(da1, da2, da3);

julia> s[:b, 10.0]
(one = 4.0, two = 8.0, three = 12.0)

julia> s[X(:a)] isa DimStack
true
```

"""
struct DimStack{L,D<:Tuple,R<:Tuple,LD<:NamedTuple,M,LM<:NamedTuple} <: AbstractDimStack{L,D}
    data::L
    dims::D
    refdims::R
    layerdims::LD
    metadata::M
    layermetadata::LM
end
DimStack(das::AbstractDimArray...; kwargs...) = DimStack(das; kwargs...)
function DimStack(das::Tuple{Vararg{<:AbstractDimArray}}; kwargs...)
    DimStack(NamedTuple{uniquekeys(das)}(das); kwargs...)
end
function DimStack(das::NamedTuple{<:Any,<:Tuple{Vararg{<:AbstractDimArray}}};
    refdims=(), metadata=NoMetadata(), layermetadata=map(DD.metadata, das)
)
    data = map(parent, das)
    dims = combinedims(das...)
    layerdims = map(basedims, das)
    DimStack(data, dims, refdims, layerdims, metadata, layermetadata)
end
# Same sized arrays
function DimStack(data::NamedTuple, dims::Tuple;
    refdims=(), metadata=NoMetadata(), layermetadata=map(_ -> NoMetadata(), data)
)
    all(map(d -> axes(d) == axes(first(data)), data)) || _stack_size_mismatch()
    layerdims = map(_ -> basedims(dims), data)
    DimStack(data, formatdims(first(data), dims), refdims, layerdims, metadata, layermetadata)
end

@noinline _stack_size_mismatch() = throw(ArgumentError("Arrays must have identical axes. For mixed dimensions, use DimArrays`"))


text:
module Ops

using DifferentialForms: Forms
using FillArrays
using LinearAlgebra
using SparseArrays

using ..Defs
using ..Funs
using ..Manifolds
using ..SparseOps

# TODO: Add type `S` from `Manifold{D,S}`?
export Op
"""
Operator acting on `Fun` functions
"""
struct Op{D,P1,R1,P2,R2,T} # <: AbstractMatrix{T}
    manifold::Manifold{D}
    values::AbstractMatrix{T}

    function Op{D,P1,R1,P2,R2,T}(manifold::Manifold{D},
                                 values::AbstractMatrix{T}) where {D,P1,R1,P2,
                                                                   R2,T}
        op = new{D,P1,R1,P2,R2,T}(manifold, dnz(values))
        #TODO @assert invariant(op)
        return op
    end
    function Op{D,P1,R1,P2,R2}(manifold::Manifold{D},
                               values::AbstractMatrix{T}) where {D,P1,R1,P2,R2,
                                                                 T}
        return Op{D,P1,R1,P2,R2,T}(manifold, values)
    end
end

# Drop non-structural zeros if we can do that without losing accuracy
const ExactTypes = Union{Union{Integer,Rational},
                         Complex{<:Union{Integer,Rational}}}
const FloatTypes = Union{AbstractFloat,Complex{<:AbstractFloat}}
dnz(A) = A
function dnz(A::AbstractSparseMatrix{T}) where {T<:ExactTypes}
    return dropzeros(A)::AbstractSparseMatrix{T}
end
function dnz(A::AbstractSparseMatrix{T}) where {T<:FloatTypes}
    return dropzeros!(chop.(A))::AbstractSparseMatrix{T}
end
dnz(A::Adjoint) = adjoint(dnz(adjoint(A)))
dnz(A::Transpose) = transpose(dnz(transpose(A)))
# dnz(A::Union{LowerTriangular, UpperTriangular}) = typeof(A)(dnz(A.data))

chop(f::ExactTypes) = f
chop(f::FloatTypes) = ifelse(abs2(f) < eps34sq(typeof(f)), zero(f), f)
function chop(f::Complex{T}) where {T<:AbstractFloat}
    return Complex{T}(chop(real(f)), chop(imag(f)))
end

eps34sq(::Type{T}) where {T} = sqrt(eps(T)^3)::T
eps34(::Type{T}) where {T} = sqrt(eps34sq(T))::T

function Base.show(io::IO, op::Op{D,P1,R1,P2,R2,T}) where {D,P1,R1,P2,R2,T}
    println(io)
    println(io, "Op{$D,$P1,$R1,$P2,$R2,$T}(")
    println(io, "    manifold=$(op.manifold.name)")
    println(io, "    values=$(op.values)")
    return print(io, ")")
end

function Defs.invariant(op::Op{D,P1,R1,P2,R2}) where {D,P1,R1,P2,R2}
    D::Int
    @assert D ≥ 0
    P1::PrimalDual
    R1::Int
    @assert 0 ≤ R1 ≤ D
    P2::PrimalDual
    R2::Int
    @assert 0 ≤ R2 ≤ D
    @assert size(op.values) ==
            (nsimplices(op.manifold, R1), nsimplices(op.manifold, R2))
    return true
end

# Comparison

function Base.:(==)(A::M, B::M) where {M<:Op}
    @assert A.manifold ≡ B.manifold
    A ≡ B && return true
    return A.values == B.values
end
function Base.:(<)(A::M, B::M) where {M<:Op}
    @assert A.manifold ≡ B.manifold
    return A.values < B.values
end
function Base.isequal(A::M, B::M) where {M<:Op}
    @assert A.manifold ≡ B.manifold
    return isequal(A.values, B.values)
end
function Base.hash(A::Op, h::UInt)
    return hash(0xfc734743, hash(A.manifold, hash(A.values, h)))
end

# Random operators
function Base.rand(::Type{Op{D,P1,R1,P2,R2,T}},
                   manifold::Manifold{D}) where {D,P1,R1,P2,R2,T}
    nrows, ncols = nsimplices(manifold, R1), nsimplices(manifold, R2)
    p = clamp(4 / min(nrows, ncols), 0, 1)
    return Op{D,P1,R1,P2,R2}(manifold, sprand(T, nrows, ncols, p))
end

# Operators are a collection

Base.eltype(::Type{<:Op{<:Any,<:Any,<:Any,<:Any,<:Any,T}}) where {T} = T
Base.eltype(A::Op) = eltype(typeof(A))
Base.isempty(A::Op) = isempty(A.values)
Base.iterate(A::Op, state...) = iterate(A.values, state...)
Base.length(A::Op) = length(A.values)

# function Base.map(op, A::Op{D,P1,R1,P2,R2},
#                   Bs::Op{D,P1,R1,P2,R2}...) where {D,P1,R1,P2,R2}
#     @assert all(A.manifold ≡ B.manifold for B ∈ Bs)
#     return Fun{D,R}(A.manifold, map(op, A.values, (B.values for B ∈ Bs)...))
# end
function Base.map(op, A::Op{D,P1,R1,P2,R2},
                  Bs::Op{D,P1,R1,P2,R2}...) where {D,P1,R1,P2,R2}
    @assert all(B.manifold ≡ A.manifold for B in Bs)
    U = typeof(op(zero(eltype(A)), (zero(eltype(B)) for B in Bs)...))
    return Op{D,P1,R1,P2,R2,U}(A.manifold,
                               map(op, A.values, (B.values for B in Bs)...))
end
function Base.reduce(op, A::Op{D,P1,R1,P2,R2}, Bs::Op{D,P1,R1,P2,R2}...;
                     kw...) where {D,P1,R1,P2,R2}
    @assert all(B.manifold ≡ A.manifold for B in Bs)
    return reduce(op, A.values, (B.values for B in Bs)...; kw...)
end

# Operators are an abstract matrix

Base.IndexStyle(::Type{<:Op}) = IndexStyle(Vector)
Base.axes(A::Op) = axes(A.values)
Base.axes(A::Op, dir) = axes(A.values, dir)
Base.eachindex(A::Op) = eachindex(A.values)
Base.@propagate_inbounds function Base.getindex(A::Op, inds...)
    return getindex(A.values, inds...)
end
Base.ndims(::Op) = 2
Base.size(A::Op) = size(A.values)
Base.size(A::Op, dims) = size(A.values, dims)
Base.stride(A::Op, k) = stride(A.values, k)
Base.strides(A::Op) = strides(A.values)

# Operators are a vector space

function Base.zero(::Type{Op{D,P1,R1,P2,R2,T}},
                   manifold::Manifold{D}) where {D,P1,R1,P2,R2,T}
    nrows = nsimplices(manifold, R1)
    ncols = nsimplices(manifold, R2)
    return Op{D,P1,R1,P2,R2}(manifold, spzeros(T, nrows, ncols))
end
function Base.zero(::Type{Op{D,P1,R1,P2,R2}},
                   manifold::Manifold{D}) where {D,P1,R1,P2,R2}
    return zero(Op{D,P1,R1,P2,R2,Bool}, manifold)
end
Base.zero(A::Op) = zero(typeof(A), A.manifold)
Base.iszero(A::Op) = iszero(A.values)

LinearAlgebra.norm(A::Op) = norm(A.values)
LinearAlgebra.norm(A::Op, p::Real) = norm(A.values, p)

function Forms.unit(::Type{Op{D,P1,R1,P2,R2,T}}, manifold::Manifold{D},
                    row::Int, col::Int) where {D,P1,R1,P2,R2,T}
    nrows = nsimplices(manifold, R1)
    ncols = nsimplices(manifold, R2)
    @assert 1 ≤ row ≤ nrows
    @assert 1 ≤ col ≤ ncols
    return Op{D,P1,R1,P2,R2}(manifold,
                             sparse([row], [col], [one(T)], nrows, ncols))
end

function Base.:+(A::Op{D,P1,R1,P2,R2}) where {D,P1,R1,P2,R2}
    return Op{D,P1,R1,P2,R2}(A.manifold, +A.values)
end

function Base.:-(A::Op{D,P1,R1,P2,R2}) where {D,P1,R1,P2,R2}
    return Op{D,P1,R1,P2,R2}(A.manifold, -A.values)
end

function Base.:+(A::Op{D,P1,R1,P2,R2},
                 B::Op{D,P1,R1,P2,R2}) where {D,P1,R1,P2,R2}
    @assert A.manifold ≡ B.manifold
    return Op{D,P1,R1,P2,R2}(A.manifold, A.values + B.values)
end

function Base.:-(A::Op{D,P1,R1,P2,R2},
                 B::Op{D,P1,R1,P2,R2}) where {D,P1,R1,P2,R2}
    @assert A.manifold ≡ B.manifold
    return Op{D,P1,R1,P2,R2}(A.manifold, A.values - B.values)
end

function Base.:*(a::Number, A::Op{D,P1,R1,P2,R2}) where {D,P1,R1,P2,R2}
    return Op{D,P1,R1,P2,R2}(A.manifold, a * A.values)
end

function Base.:\(a::Number, A::Op{D,P1,R1,P2,R2}) where {D,P1,R1,P2,R2}
    return Op{D,P1,R1,P2,R2}(A.manifold, a \ A.values)
end

function Base.:*(A::Op{D,P1,R1,P2,R2}, a::Number) where {D,P1,R1,P2,R2}
    return Op{D,P1,R1,P2,R2}(A.manifold, A.values * a)
end

function Base.:/(A::Op{D,P1,R1,P2,R2}, a::Number) where {D,P1,R1,P2,R2}
    return Op{D,P1,R1,P2,R2}(A.manifold, A.values / a)
end

# Operators are a ring

function Base.one(::Type{Op{D,P,R,P,R,T}},
                  manifold::Manifold{D}) where {D,P,R,T}
    nrows = nsimplices(manifold, R)
    return Op{D,P,R,P,R}(manifold, Diagonal(Fill(one(T), (nrows, nrows))))
end
function Base.one(::Type{Op{D,P,R,P,R}}, manifold::Manifold{D}) where {D,P,R}
    return one(Op{D,P,R,P,R,Bool}, manifold)
end
Base.one(A::Op) = one(typeof(A), A.manifold)
Base.isone(A::Op{D,P,R,P,R}) where {D,P,R} = A.values == I
Base.isone(A::Op) = false

function Base.:*(A::Op{D,P1,R1,P2,R2},
                 B::Op{D,P2,R2,P3,R3}) where {D,P1,R1,P2,R2,P3,R3}
    @assert A.manifold ≡ B.manifold
    return Op{D,P1,R1,P3,R3}(A.manifold, A.values * B.values)
end

# Operators are a groupoid ("division ring")

# Note: This works only for invertible operators, and only for some
# matrix representations
function Base.inv(A::Op{D,P1,R1,P2,R2}) where {D,P1,R1,P2,R2}
    @assert R1 == R2
    @assert size(R1, A.manifold) == size(R2, A.manifold)
    return Op{D,P2,R2,P1,R1}(A.manifold, inv(A.values))
end

function Base.:/(A::Op{D,P1,R1,P2,R2},
                 B::Op{D,P3,R3,P2,R2}) where {D,P1,R1,P2,R2,P3,R3}
    @assert A.manifold ≡ B.manifold
    return Op{D,P1,R1,P3,R3}(A.manifold, A.values / B.values)
end

function Base.:\(A::Op{D,P2,R2,P1,R1},
                 B::Op{D,P2,R2,P3,R3}) where {D,P1,R1,P2,R2,P3,R3}
    @assert A.manifold ≡ B.manifold
    return Op{D,P1,R1,P3,R3}(A.manifold, A.values \ B.values)
end

# There are adjoints

function Base.adjoint(A::Op{D,P1,R1,P2,R2}) where {D,P1,R1,P2,R2}
    return Op{D,P2,R2,P1,R1}(A.manifold, adjoint(A.values))
end

function Base.transpose(A::Op{D,P1,R1,P2,R2}) where {D,P1,R1,P2,R2}
    return Op{D,P2,R2,P1,R1}(A.manifold, transpose(A.values))
end

# TODO: shouldn't this be defined automatically by being a collection?
function Base.conj(A::Op{D,P1,R1,P2,R2}) where {D,P1,R1,P2,R2}
    return Op{D,P1,R1,P2,R2}(A.manifold, conj(A.values))
end

# Operators act on functions

function Base.:*(A::Op{D,P1,R1,P2,R2}, f::Fun{D,P2,R2}) where {D,P1,R1,P2,R2}
    @assert A.manifold ≡ f.manifold
    return Fun{D,P1,R1}(f.manifold, IDVector{R1}(A.values * f.values.vec))
end

function Base.:\(A::Op{D,P1,R1,P2,R2}, f::Fun{D,P1,R1}) where {D,P1,R1,P2,R2}
    @assert A.manifold ≡ f.manifold
    # Note: \ converts rationals to Float64
    return Fun{D,P2,R2}(f.manifold, IDVector{R2}(A.values \ f.values.vec))
end

end

ext: jl


text:
module ModifiedArrays

import Base: first, tail
using Base: @propagate_inbounds

export
    modify,
    ModifiedArray,
    # transpose.jl
    modified_transpose,
    modified_adjoint,
    TransposeMod,
    AdjointMod

## Main types

"The supertype of all array modifiers."
abstract type ArrayModifier end

"The special case of no modification at all."
struct NoModifier <: ArrayModifier end

"The main modified array type consists of a parent array and a modifier."
struct ModifiedArray{T,N,AA,M} <: AbstractArray{T,N}
    parent      ::  AA      # the array-like object being modified
    modifier    ::  M       # the modifier
end

ModifiedArray(A, mod::ArrayModifier) = ModifiedArray{mod_eltype(mod,A),mod_ndims(mod,A)}(A, mod)
ModifiedArray{T,N}(A, mod::ArrayModifier) where {T,N} =
    ModifiedArray{T,N,typeof(A),typeof(mod)}(A, mod)

modify(A, mod::ArrayModifier) = ModifiedArray(A, mod)

Base.parent(a::ModifiedArray) = a.parent
modifier(a::ModifiedArray) = a.modifier

## Array specific interface

abstract type Modifiable end

struct Interface{S} <: Modifiable end

# The two type variables
const IF_eltype = Interface{:eltype}
const IF_ndims = Interface{:ndims}

# Compulsory interface methods
const IF_size = Interface{:size}
const IF_getindex = Interface{:getindex}
const IF_setindex! = Interface{:setindex!}

# Optional methods
const IF_IndexStyle = Interface{:IndexStyle}
const IF_similar = Interface{:similar}
const IF_axes = Interface{:axes}

abstract type ModStyle end
struct ModRecursive <: ModStyle end
struct ModFinal <: ModStyle end
struct ModNothing <: ModStyle end

## eltype
mod_eltype(A::ModifiedArray) = mod_eltype(modifier(A), parent(A))
mod_eltype(mod::ArrayModifier, A) = mod_eltype(ModStyle(mod, IF_eltype()), mod, A)

mod_eltype(::ModNothing, mod, A) = eltype(A)
mod_eltype(::ModFinal, mod, A) = mod_eltype_final(mod)
mod_eltype(::ModRecursive, mod, A) = mod_eltype_post(mod, eltype(A))

# default definition
mod_eltype_post(mod::ArrayModifier, Z) = Z

## ndims
mod_ndims(A::ModifiedArray) = mod_ndims(modifier(A), parent(A))
mod_ndims(mod::ArrayModifier, A) = mod_ndims(ModStyle(mod, IF_ndims()), mod, A)

mod_ndims(::ModNothing, mod, A) = ndims(A)
mod_ndims(::ModFinal, mod, A) = mod_ndims_final(mod)
mod_ndims(::ModRecursive, mod, A) = mod_ndims_post(mod, ndims(A))

# default definition
mod_ndims_post(mod::ArrayModifier, Z) = Z

## size
Base.size(A::ModifiedArray) = mod_size(A)
mod_size(A::ModifiedArray) = mod_size(modifier(A), parent(A))
mod_size(mod::ArrayModifier, A) = mod_size(ModStyle(mod, IF_size()), mod, A)

mod_size(::ModNothing, mod, A) = size(A)
mod_size(::ModFinal, mod, A) = mod_size_final(mod)
mod_size(::ModRecursive, mod, A) = mod_size_post(mod, size(A))

# default definition
mod_size_post(mod::ArrayModifier, Z) = Z

## axes
Base.axes(A::ModifiedArray) = mod_axes(A)
mod_axes(A::ModifiedArray) = mod_axes(modifier(A), parent(A))
mod_axes(mod::ArrayModifier, A) = mod_axes(ModStyle(mod, IF_axes()), mod, A)

mod_axes(::ModNothing, mod, A) = axes(A)
mod_axes(::ModFinal, mod, A) = mod_axes_final(mod)
mod_axes(::ModRecursive, mod, A) = mod_axes_post(mod, axes(A))

# default definition
mod_axes_post(mod::ArrayModifier, Z) = Z

## IndexStyle
Base.IndexStyle(Mod::Type{<:ModifiedArray}) = mod_IndexStyle(Mod)
mod_IndexStyle(::Type{ModifiedArray{T,N,AA,M}}) where {T,N,AA,M} =
    mod_IndexStyle(M, AA)
mod_IndexStyle(M::Type{<:ArrayModifier}, AA) = mod_IndexStyle(ModStyle(M, IF_IndexStyle()), M, AA)

mod_IndexStyle(::ModNothing, M, AA) = IndexStyle(AA)
mod_IndexStyle(::ModFinal, M, AA) = mod_IndexStyle_final(M)
mod_IndexStyle(::ModRecursive, M, AA) = mod_IndexStyle_post(M, IndexStyle(AA))

# default definition
mod_IndexStyle_post(::Type{<:ArrayModifier}, Z) = Z

## getindex

# We only intercept calls with a linear index, or with a number of indices equal
# to the dimension of the array.
@propagate_inbounds Base.getindex(A::ModifiedArray{T,1}, i::Int) where {T} = mod_getindex(A, i)
@propagate_inbounds Base.getindex(A::ModifiedArray{T,N}, I::Vararg{Int,N}) where {T,N} =
    mod_getindex(A, I...)
# Linear indexing of arrays that are not vectors:
@propagate_inbounds Base.getindex(A::ModifiedArray{T,N}, i::Int) where {T,N} =
    __getindex(IndexStyle(A), A, i)
# - the array supports linear indexing: invoke mod_getindex
@propagate_inbounds __getindex(::IndexLinear, A, i) = mod_getindex(A, i)
# - the array supports cartesian indexing: translate the index to cartesian
@propagate_inbounds __getindex(::IndexCartesian, A, i) = mod_getindex(A, Base._to_subscript_indices(A, i)...)

@propagate_inbounds function mod_getindex(A::ModifiedArray, I...)
    @boundscheck checkbounds(A, I...)
    mod_getindex(modifier(A), parent(A), I...)
end
@propagate_inbounds function mod_getindex(mod::ArrayModifier, A, I...)
    mod_getindex(ModStyle(mod, IF_getindex()), mod, A, I...)
end

@propagate_inbounds mod_getindex(::ModNothing, mod, A, I...) = getindex(A, I...)
@propagate_inbounds mod_getindex(::ModFinal, mod, istyle, A, I...) = mod_getindex_final(mod, I...)

@propagate_inbounds function mod_getindex(::ModRecursive, mod, A, I...)
    mod_getindex_post(mod, A[mod_getindex_pre(mod, I...)...], I...)
end

# default definitions
mod_getindex_pre(mod::ArrayModifier, I...) = I
mod_getindex_post(mod::ArrayModifier, Z, I...) = Z

## setindex!
# Like with getindex above, we only intercept a few calls to setindex!.
@propagate_inbounds Base.setindex!(A::ModifiedArray{T,1}, val, i::Int) where {T} =
    mod_setindex!(A, val, i)
@propagate_inbounds Base.setindex!(A::ModifiedArray{T,N}, val, I::Vararg{Int,N}) where {T,N} =
    mod_setindex!(A, val, I...)
# Linear indexing of arrays that are not vectors, similar to getindex above:
@propagate_inbounds Base.setindex!(A::ModifiedArray{T,N}, val, i::Int) where {T,N} =
    __setindex!(IndexStyle(A), A, val, i)
@propagate_inbounds __setindex!(::IndexLinear, A, val, i) = mod_setindex!(A, val, i)
@propagate_inbounds __setindex!(::IndexCartesian, A, val, i) =
    mod_setindex!(A, val, Base._to_subscript_indices(A, i)...)

@propagate_inbounds function mod_setindex!(A::ModifiedArray, val, I...)
    @boundscheck checkbounds(A, I...)
    mod_setindex!(modifier(A), parent(A), val, I...)
    # The return value of setindex! is not prescribed by the array interface.
    # In most cases, the mutated object is returned, so we just return A.
    # That means there is no need for a mod_setindex!_post function.
    # See also Julia issue #31891
    A
end
@propagate_inbounds mod_setindex!(mod::ArrayModifier, A, val, I...) =
    mod_setindex!(ModStyle(mod, IF_setindex!()), mod, A, val, I...)

@propagate_inbounds mod_setindex!(::ModNothing, mod, A, val, I...) = setindex!(A, val, I...)
@propagate_inbounds mod_setindex!(::ModFinal, mod, A, val, I...) =
    mod_setindex!_final(mod, val, I...)
@propagate_inbounds mod_setindex!(::ModRecursive, mod, A, val, I...) =
    setindex!(A, mod_setindex!_pre(mod, val, I...)...)

# default definition
mod_setindex!_pre(mod::ArrayModifier, val, I...) = (val, I...)

## similar
Base.similar(A::ModifiedArray, ::Type{T}, dims::Dims) where {T} =
    mod_similar(A, T, dims)
mod_similar(A::ModifiedArray, ::Type{T}, dims::Dims) where {T} =
    mod_similar(modifier(A), parent(A), T, dims)
mod_similar(mod::ArrayModifier, A, T, dims) =
    mod_similar(ModStyle(mod, IF_similar()), mod, A, T, dims)

mod_similar(::ModNothing, mod, A, T, dims) = similar(A, T, dims)
mod_similar(::ModFinal, mod, A, T, dims) = mod_similar_final(mod, T, dims)
mod_similar(::ModRecursive, mod, A, T, dims) =
    mod_similar_post(mod, similar(A, mod_similar_pre(mod, T, dims)...), T, dims)

# default definition
mod_similar_pre(mod::ArrayModifier, T, dims) = (T, dims)
mod_similar_post(mod::ArrayModifier, Z, T, dims) = Z

include("composite.jl")
include("transpose.jl")
include("offset.jl")
include("eltype.jl")
include("mapped.jl")

end # module


text:
import DBInterface, PostgreSQL, DataFrames, LibPQ
using Test

@testset "PostgreSQL.jl" begin
    db_name = "postgres"
    user = "postgres"
    port = 5432
    host = "localhost"
    password = "pass"

    conn = DBInterface.connect(PostgreSQL.Connection, "dbname=$db_name user=$user port=$port host=$host password=$password")
    @test typeof(conn) <: PostgreSQL.Connection

    stmt = DBInterface.prepare(conn, "CREATE DATABASE test;")
    @test typeof(stmt) <: PostgreSQL.Stmt

    result = DBInterface.execute(stmt)
    @test LibPQ.status(result) == LibPQ.libpq_c.PGRES_COMMAND_OK

    DBInterface.close!(conn)
    @test !LibPQ.isopen(conn.inner)
    @test conn.inner.closed[] == true

    conn = DBInterface.connect(PostgreSQL.Connection, "dbname=test user=$user port=$port host=$host password=$password")

    sql = """
    CREATE TABLE Persons(
        id              SERIAL,
        name            TEXT NOT NULL,
        PRIMARY KEY (id)
    );
    """
    stmt = DBInterface.prepare(conn, sql)
    result = DBInterface.execute(stmt)
    @test LibPQ.status(result) == LibPQ.libpq_c.PGRES_COMMAND_OK

    stmt = DBInterface.prepare(conn, "INSERT INTO Persons(name) VALUES(?);")
    result = DBInterface.execute(stmt, ["adam"])
    @test LibPQ.status(result) == LibPQ.libpq_c.PGRES_COMMAND_OK

    stmt = DBInterface.prepare(conn, "INSERT INTO Persons(name) VALUES(:name);")
    result = DBInterface.execute(stmt, ("maria",))
    @test LibPQ.status(result) == LibPQ.libpq_c.PGRES_COMMAND_OK

    result = DBInterface.executemany(stmt, (name = ["eva", "paul"],))
    @test LibPQ.status(result) == LibPQ.libpq_c.PGRES_COMMAND_OK

    df = DBInterface.execute(conn, "SELECT * from Persons;") |> DataFrames.DataFrame
    @test df.name == ["adam", "maria", "eva", "paul"]

    DBInterface.close!(conn)
end


text:
# Part of the code of the function denseread come from
#The MatrixMarket.jl package which is licensed under the MIT Expat License:

#Copyright (c) 2013: Viral B. Shah.

#Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:

#The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.

#THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

function denseread(filename::AbstractString)
    file = open(filename, "r")
    firstline = chomp(readline(file))
    tokens = split(firstline)
    if tokens[1] != "%%MatrixMarket"
        throw(ParseError(string("Not a valid MatrixMarket header:", ll)))
    end
    (head1, rep, field, symm) = map(lowercase, tokens[2:5])
    eltype = field == "real" ? Float64 :
                  field == "complex" ? Complex128 :
                  field == "pattern" ? Bool :
                  throw(ParseError("Unsupported field $field (only real and complex are supported)"))
    # Skip all comments and empty lines
    ll   = readline(file)
    while length(chomp(ll))==0 || (length(ll) > 0 && ll[1] == '%')
        ll = readline(file)
    end
    return readdlm(file, eltype)
end


text:
# Used to load sequential data where the time index is preceeded by 0's: 0004
function numberDigits( number; maxdigits=4 )
   ndigs = 0
   for d in 0:maxdigits
       ndigs += number ÷ 10^d > 0
   end
   nzeros = ( maxdigits - ndigs );
   return ( nzeros < 0 ) ? string(number) : "0"^(nzeros)*string(number)
end

function vectorFieldSize( datasize, pp::PIVParameters )
    step = pp.interArea .- pp.overlap
    return 1 .+ floor.( Int64, ( datasize .- pp.interArea )./step )
end

function constructGrid( datasize::III, IA::III, step::III, mp::I )

    dims  = 1 .+ floor.( Int64, ( datasize .- IA )./step )
    grid  = Array{NTuple{12,Int64}}( undef, prod( dims ) )
    count = 1;

    # Each square/cube subdivision of the input images/volumes includes the image/volume
    # coordinates of the subdivision, and the coordinates of computed displacement
    # vector in the final vector field ( U, V, W ).
    rdelta = [ step[1], step[1], mp, mp ]
    cdelta = [ step[2], step[2], mp, mp ]
    ddelta = [ step[3], step[3], mp, mp ]

    @inbounds begin
    # ds = [ minDepthData, maxDepthData, minDepthUVW, maxDepthUVW ]
    ds = [ 1, IA[3], 1, mp ]
    for id in 1:dims[3]
        # cs = [ minColData, maxColData, minColUVW, maxColUVW ]
        cs = [ 1, IA[2], 1, mp ]
        for iw in 1:dims[2]
            # rs = [ minRowData, maxRowData, minRowUVW, maxRowUVW ]
            rs = [ 1, IA[1], 1, mp ]
            for ih in 1:dims[1]

                grid[count] = ( rs[1], rs[2], cs[1], cs[2], ds[1], ds[2],
                                rs[3], rs[4], cs[3], cs[4], ds[3], ds[4] )
                count += 1;
                @simd for idx in 1:4 ( rs[idx] += rdelta[idx] ) end
            end
            @simd for idx in 1:4 ( cs[idx] += cdelta[idx] ) end
        end
        @simd for idx in 1:4 ( ds[idx] += ddelta[idx] ) end
    end
    end #@inbounds

    return grid
end

# compute coordinates of search area after shifting them by previously computed displacements during multi-pass
function searchGrid( interGrid, SM, offh, offw, offd, h, w, d )
    r1 = max( interGrid[1] - SM[1] + offh, 1 );
    r2 = min( interGrid[2] + SM[1] + offh, h );
    c1 = max( interGrid[3] - SM[2] + offw, 1 );
    c2 = min( interGrid[4] + SM[2] + offw, w );
    d1 = max( interGrid[5] - SM[2] + offd, 1 );
    d2 = min( interGrid[6] + SM[2] + offd, d );
    so = max( 0, r1 - ( interGrid[1] - SM[1] ) - offh ),
         max( 0, c1 - ( interGrid[3] - SM[2] ) - offw ),
         max( 0, d1 - ( interGrid[5] - SM[3] ) - offd );

    return ( r1, r2, c1, c2, d1, d2 ), so
end

function setTo0!( A::Array{T,N} ) where {T<:Number,N}
    zer = T(0)
    @inbounds @simd for i in 1:length(A)
        A[ i ] = zer
    end
end

function setTo0!( A::Array{T,N}... ) where {T<:Number,N}
    setTo0!.( A )
end

# optimized to use SIMD, performing multiple comparisons at the same time
function maxval( a::AbstractArray{T,N} ) where {T<:Real,N}
    m1 = a[1]
    m2 = a[2]
    m3 = a[3]
    m4 = a[4]
    m5 = a[end]

    @inbounds begin
        @simd for i in 5:4:length(a)-4
            m1 = ( a[ i ] > m1 ) ? a[ i ] : m1
            m2 = ( a[i+1] > m2 ) ? a[i+1] : m2
            m3 = ( a[i+2] > m3 ) ? a[i+2] : m3
            m4 = ( a[i+3] > m4 ) ? a[i+3] : m4
        end
        for i in length(a)-4:length(a)-3
            m5 = ( a[i] > m5 ) ? a[i] : m5
        end
    end
    return maximum( [ m1, m2, m3, m4, m5 ] )
end

# optimized to perform parallel computations with simd
function firstPeak( A::Array{T,N} ) where {T<:Real,N}

    if length(A) < 12
        mv, mx = findmax( A )
        return Tuple(mx), mv
    end

    # maximum indices for parallel SIMD processing
    mi1 = 1
    mi2 = 2
    mi3 = 3
    mi4 = 4
    mi5 = length(A)

    @inbounds begin
    @simd for i in 5:4:length(A)-4
        mi1 = ( A[ i ] > A[mi1] ) ?  i  : mi1;
        mi2 = ( A[i+1] > A[mi2] ) ? i+1 : mi2;
        mi3 = ( A[i+2] > A[mi3] ) ? i+2 : mi3;
        mi4 = ( A[i+3] > A[mi4] ) ? i+3 : mi4;
    end

    for i in length(A)-4:length(A)-1
        mi5 = ( A[i] > A[mi5] ) ? i : mi5;
    end
    end #inbounds

    # Finding the max value among the "independently"/SIMD computed maximum indices
    val, idx = findmax( [  A[ mi1 ], A[ mi2 ], A[ mi3], A[ mi4 ], A[ mi5 ] ] );

    maxindex = (idx==1)*mi1 + (idx==2)*mi2 + (idx==3)*mi3 +
               (idx==4)*mi4 + (idx==5)*mi5;

    # Transforming from linear indexing to cartesian
    cidx = [ 0 for i in 1:N ];
    if N == 3
        h, w, d = size(A)
        z = ceil( Int64, maxindex/(h*w) )
        x = ceil( Int64, (maxindex - (z-1)*h*w)/h )
        y = maxindex - (x-1)*h - (z-1)*h*w;
        cidx = ( y, x, z );
    else
        h, w = size(A)
        x = ceil( Int64, maxindex/h )
        y = maxindex - (x-1)*h;
        cidx = ( y, x );
    end

    return cidx, A[ maxindex ]
end

function secondPeak( Arr::A{T,N}, rad::I, idx::NTuple{N,I}) where {T<:Real,N}

    # extract region around first peak, before setting it to 0
    mins   = max.(     1    , idx .- rad );
    maxs   = min.( size(Arr), idx .+ rad );
    ranges = [ mins[e]:maxs[e] for e in 1:N ];
    OGvals = copy( Arr[ ranges... ] );

    # setting region to 0
    Arr[ ranges... ] .= T(0);

    # computing first peak
    p2, val = firstPeak( Arr );

    # setting region back to its original value
    Arr[ ranges... ] .= OGvals

    return p2, val
end

# Used in subpixel approximation functions
function indexOnBorder( index::NTuple{N,I}, csize::NTuple{N,I} ) where {N}
    return any( index .== 1 ) || any( index .== csize )
end

neighIndices( idx::II  ) = [(idx[1]+y,idx[2]+x)   for y in -1:1, x in -1:1 ]
neighIndices( idx::III ) = [(idx[1]+y,idx[2]+x,idx[3]+z) for y in -1:1, x in -1:1, z in -1:1 ]

function getNeighbours( index::NTuple{N,I}, cmatrix::Array{T,N} ) where {T<:AbstractFloat,N}
    indices = neighIndices( index );
    neighs  = [ cmatrix[ x... ] for x in indices ]
    return neighs
end

function putWithinSearch!( search::Array{T,N}, A::Array{U,N},
                           so::NTuple{3,Int64}, coord::NTuple{6,Int64}
                         ) where {T<:Real,U<:Real,N}
    start  = [ x + 1 for x in so ];
    r1, r2 = (coord[1], coord[2])
    c1, c2 = (coord[3], coord[4])
    d1, d2 = (coord[5], coord[6])

	if N == 2
		start[3] = 1;
		d1, d2 = 1,1;
	end

    @inbounds begin
    for z in 0:(d2-d1)
        for c in 0:(c2-c1)
            @simd for r in 0:(r2-r1)
                search[start[1]+r, start[2]+c, start[3]+z] = A[r1+r, c1+c, d1+z]
            end
        end
    end
    end # inbounds
end

function putWithinPadded!( padA::Array{Complex{T},N}, A::AbstractArray{U,N}, mean,
                           so::NTuple{3,Int64}, coord::NTuple{6,Int64}
                         )  where {T<:AbstractFloat,U<:Real,N}
    off    = [ x + 1 for x in so ]
    r1, r2 = (coord[1], coord[2])
    c1, c2 = (coord[3], coord[4])
    d1, d2 = (coord[5], coord[6])

	if N == 2
		off[3] =  1;
		d1,d2 = 1,1;
	end

    Tzero = T(0)
    Tmean = convert( T, mean )

    @inbounds begin
    for z in 0:(d2-d1)
        for c in 0:(c2-c1)
            @simd for r in 0:(r2-r1)
                padA[ off[1]+r,
                      off[2]+c,
                      off[3]+z ] = complex( convert( T, A[r1+r,c1+c,d1+z] ) - Tmean, Tzero )
            end
        end
    end
    end # inbounds
end

function putWithinPadded!( padA::Array{Complex{T},N}, A::AbstractArray{U,N}, mean
                         )  where {T<:AbstractFloat,U<:Real,N}
    Tzero = T(0)
    Tmean = convert( T, mean )

    @inbounds begin
    for zet in 1:size( A, 3 )
        for col in 1:size( A, 2 )
            @simd for row in 1:size( A, 1 )
                padA[row,col,zet] = complex( convert( T, A[row,col,zet] ) - Tmean, Tzero )
            end
        end
    end
    end #@inbounds
end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule libsharp2_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("libsharp2")
JLLWrappers.@generate_main_file("libsharp2", UUID("180207a7-b08e-5162-af94-7d62a04fe081"))
end  # module libsharp2_jll


text:
###############################################################################
#
#   FPSymbol/FPGroupElem/FPGroup definition
#
###############################################################################

struct FPSymbol <: GSymbol
   id::Symbol
   pow::Int
end

FPGroupElem = GroupWord{FPSymbol}

mutable struct FPGroup <: AbstractFPGroup
   gens::Vector{FPSymbol}
   rels::Dict{FPGroupElem, FPGroupElem}

   function FPGroup(gens::Vector{T}, rels::Dict{FPGroupElem, FPGroupElem}) where {T<:GSymbol}
      G = new(gens)
      G.rels = Dict(G(k) => G(v) for (k,v) in rels)
      return G
   end
end

export FPGroupElem, FPGroup

###############################################################################
#
#   Type and parent object methods
#
###############################################################################

parent_type(::Type{FPGroupElem}) = FPGroup

elem_type(::FPGroup) = FPGroupElem

###############################################################################
#
#   FPSymbol constructors
#
###############################################################################

FPSymbol(s::Symbol) = FPSymbol(s, 1)
FPSymbol(s::String) = FPSymbol(Symbol(s))

convert(::Type{FPSymbol}, s::FreeSymbol) = FPSymbol(s.id, s.pow)

FPGroup(gens::Vector{FPSymbol}) = FPGroup(gens, Dict{FPGroupElem, FPGroupElem}())

FPGroup(a::Vector{String}) = FPGroup([FPSymbol(i) for i in a])

FPGroup(n::Int, symbol::String="f") = FPGroup(["$symbol$i" for i in 1:n])
FPGroup(H::FreeGroup) = FPGroup([FPSymbol(s) for s in H.gens])

###############################################################################
#
#   Parent object call overloads
#
###############################################################################

function (G::FPGroup)()
   id = FPGroupElem(FPSymbol[])
   id.parent = G
   return id
end

function (G::FPGroup)(w::GWord)
   if length(w) == 0
      return G()
   end

   if eltype(w.symbols) == FreeSymbol
      w = FPGroupElem(w.symbols)
   end

   if eltype(w.symbols) == FPSymbol
      for s in w.symbols
         i = findfirst(g -> g.id == s.id, G.gens)
         i == 0 && throw(DomainError(
            "Symbol $s does not belong to $G."))
         s.pow % G.gens[i].pow == 0 || throw(DomainError(
         "Symbol $s doesn't belong to $G."))
      end
   end
   w.parent = G
   return reduce!(w)
end

(G::FPGroup)(s::FPSymbol) = G(FPGroupElem(s))

###############################################################################
#
#   Basic manipulation
#
###############################################################################

hash(s::FPSymbol, h::UInt) = hash(s.id, hash(s.pow, hash(FPSymbol, h)))

change_pow(s::FPSymbol, n::Int) = FPSymbol(s.id, n)

length(s::FPSymbol) = abs(s.pow)

###############################################################################
#
#   String I/O
#
###############################################################################

function show(io::IO, G::FPGroup)
   print(io, "FPgroup on $(length(G.gens)) generators ")
   strrels = join(G.rels, ", ")
   if length(strrels) > 300
      print(io, "⟨ ", join(G.gens, ", "), " | $(length(G.rels)) relation(s) ⟩.")
   else
      print(io, "⟨ ", join(G.gens, ", "), " | ", join(G.rels, ", "), " ⟩.")
   end
end

###############################################################################
#
#   Comparison
#
###############################################################################

###############################################################################
#
#   Inversion
#
###############################################################################

inv(s::FPSymbol) = change_pow(s, -s.pow)

###############################################################################
#
#   Binary operations
#
###############################################################################

(*)(W::FPGroupElem, Z::FPGroupElem) = r_multiply(W, Z.symbols)
(*)(W::FPGroupElem, s::FPSymbol) = r_multiply(W, [s])
(*)(s::FPSymbol, W::FPGroupElem) = l_multiply(W, [s])

function reduce!(W::FPGroupElem)
    if length(W) < 2
        deleteat!(W.symbols, findall(x -> x.pow == 0, W.symbols))
    else
        reduced = false
        while !reduced
            reduced = freereduce!(W) || replace_all!(W, parent(W).rels)
        end
    end

    W.savedhash = hash(W.symbols, hash(typeof(W)))
    W.modified = false
    return W
end

###############################################################################
#
#   Misc
#
###############################################################################

function add_rels!(G::FPGroup, newrels::Dict{FPGroupElem,FPGroupElem})
   for w in keys(newrels)
      if !(w in keys(G.rels))
         G.rels[w] = G(newrels[w])
      end
   end
end

function /(G::FPGroup, newrels::Vector{FPGroupElem})
   for r in rels
      parent(r) == G || throw(DomainError(
      "Can not form quotient group: $r is not an element of $G"))
   end
   H = deepcopy(G)
   newrels = Dict(H(r) => H() for r in newrels)
   add_rels!(H, newrels)
   return H
end

function /(G::FreeGroup, rels::Vector{FreeGroupElem})
   for r in rels
      parent(r) == G || throw(DomainError(
         "Can not form quotient group: $r is not an element of $G"))
   end
   H = FPGroup(G)
   H.rels = Dict(H(rel) => H() for rel in unique(rels))
   return H
end


text:
@testset "JuMP direct model - Linear Constraints - Affine parameters" begin
    optimizer = POI.ParametricOptimizer(GLPK.Optimizer())

    model = direct_model(optimizer)

    @variable(model, x[i=1:2] >= 0)

    @variable(model, y in POI.Parameter(0))
    @variable(model, w in POI.Parameter(0))
    @variable(model, z in POI.Parameter(0))

    @constraint(model, 2*x[1] + x[2] + y <= 4)
    @constraint(model, 1*x[1] + 2*x[2] + z <= 4)

    @objective(model, Max, 4*x[1] + 3*x[2] + w)

    optimize!(model)

    @test isapprox.(value(x[1]), 4.0/3.0, atol = ATOL)
    @test isapprox.(value(x[2]), 4.0/3.0, atol = ATOL)
    @test isapprox.(value(y), 0, atol = ATOL)

    # ===== Set parameter value =====
    MOI.set(model, POI.ParameterValue(), y, 2.0)
    optimize!(model)

    @test isapprox.(value(x[1]), 0.0, atol = ATOL)
    @test isapprox.(value(x[2]), 2.0, atol = ATOL)
    @test isapprox.(value(y), 2.0, atol = ATOL)
end

@testset "JuMP direct model - Linear Constraints - Parameter x variable" begin
    optimizer = POI.ParametricOptimizer(GLPK.Optimizer())

    model = direct_model(optimizer)

    @variable(model, x[i=1:2] >= 0)

    @variable(model, y in POI.Parameter(0))
    @variable(model, w in POI.Parameter(0))
    @variable(model, z in POI.Parameter(0))

    @constraint(model, 2*x[1] + x[2] + y <= 4)
    @constraint(model, (1+y)*x[1] + 2*x[2] + z <= 4)

    @objective(model, Max, 4*x[1] + 3*x[2] + w)

    optimize!(model)

    @test isapprox.(value(x[1]), 4.0/3.0, atol = ATOL)
    @test isapprox.(value(x[2]), 4.0/3.0, atol = ATOL)
    @test isapprox.(value(y), 0, atol = ATOL)

    # ===== Set parameter value =====
    MOI.set(model, POI.ParameterValue(), y, 2.0)
    optimize!(model)

    @test isapprox.(value(x[1]), 0.0, atol = ATOL)
    @test isapprox.(value(x[2]), 2.0, atol = ATOL)
    @test isapprox.(value(y), 2.0, atol = ATOL)
end

@testset "JuMP - Linear Constraints - Affine parameters" begin
    model = Model(() -> POI.ParametricOptimizer(GLPK.Optimizer()))

    @variable(model, x[i=1:2] >= 0)

    @variable(model, y in POI.Parameter(0))
    @variable(model, w in POI.Parameter(0))
    @variable(model, z in POI.Parameter(0))

    @constraint(model, 2*x[1] + x[2] + y <= 4)
    @constraint(model, 1*x[1] + 2*x[2] + z <= 4)

    @objective(model, Max, 4*x[1] + 3*x[2] + w)

    optimize!(model)

    @test isapprox.(value(x[1]), 4.0/3.0, atol = ATOL)
    @test isapprox.(value(x[2]), 4.0/3.0, atol = ATOL)
    @test isapprox.(value(y), 0, atol = ATOL)

    # ===== Set parameter value =====
    MOI.set(model, POI.ParameterValue(), y, 2.0)
    optimize!(model)

    @test isapprox.(value(x[1]), 0.0, atol = ATOL)
    @test isapprox.(value(x[2]), 2.0, atol = ATOL)
    @test isapprox.(value(y), 2.0, atol = ATOL)
end

@testset "JuMP - Linear Constraints - Parameter x variable" begin
    model = Model(() -> POI.ParametricOptimizer(GLPK.Optimizer()))

    @variable(model, x[i=1:2] >= 0)

    @variable(model, y in POI.Parameter(0))
    @variable(model, w in POI.Parameter(0))
    @variable(model, z in POI.Parameter(0))

    @constraint(model, 2*x[1] + x[2] + y <= 4)
    @constraint(model, (1+y)*x[1] + 2*x[2] + z <= 4)

    @objective(model, Max, 4*x[1] + 3*x[2] + w)

    @test_broken optimize!(model)

    @test_broken isapprox.(value(x[1]), 4.0/3.0, atol = ATOL)
    @test_broken isapprox.(value(x[2]), 4.0/3.0, atol = ATOL)
    @test_broken isapprox.(value(y), 0, atol = ATOL)

    # ===== Set parameter value =====
    MOI.set(model, POI.ParameterValue(), y, 2.0)
    @test_broken optimize!(model)

    @test_broken isapprox.(value(x[1]), 0.0, atol = ATOL)
    @test_broken isapprox.(value(x[2]), 2.0, atol = ATOL)
    @test_broken isapprox.(value(y), 2.0, atol = ATOL)
end


text:
using Profile
using WaveToyDagger

# @time WaveToyDagger.main()

@profile WaveToyDagger.main()

open("profile.txt", "w") do fh
    return Profile.print(IOContext(fh, :displaysize => (24, 500)))
end


text:
# AbstractBandedMatrix must implement

abstract type AbstractBandedMatrix{T} <: AbstractSparseMatrix{T,Int} end

"""
    bandeddata(A)

returns a matrix containing the data of a banded matrix, in the
BLAS format.

This is required for gbmv! support
"""
bandeddata(A) = error("Override bandeddata(::$(typeof(A)))")

"""
    bandwidths(A)

Returns a tuple containing the upper and lower bandwidth of `A`.
"""
bandwidths(A::AbstractVecOrMat) = (size(A,1)-1 , size(A,2)-1)

bandwidths(A::AdjOrTrans{T,S}) where {T,S} = reverse(bandwidths(parent(A)))

"""
    bandwidth(A,i)

Returns the lower bandwidth (`i==1`) or the upper bandwidth (`i==2`).
"""
bandwidth(A::AbstractVecOrMat, k::Integer) = bandwidths(A)[k]

"""
    bandrange(A)

Returns the range `-bandwidth(A,1):bandwidth(A,2)`.
"""
bandrange(A) = -bandwidth(A,1):bandwidth(A,2)

# start/stop indices of the i-th column/row, bounded by actual matrix size
@inline colstart(A::AbstractVecOrMat, i::Integer) = max(i-bandwidth(A,2), 1)
@inline  colstop(A::AbstractVecOrMat, i::Integer) = max(min(i+bandwidth(A,1), size(A, 1)), 0)
@inline rowstart(A::AbstractVecOrMat, i::Integer) = max(i-bandwidth(A,1), 1)
@inline  rowstop(A::AbstractVecOrMat, i::Integer) = max(min(i+bandwidth(A,2), size(A, 2)), 0)

@inline colrange(A::AbstractVecOrMat, i::Integer) = colstart(A,i):colstop(A,i)
@inline rowrange(A::AbstractVecOrMat, i::Integer) = rowstart(A,i):rowstop(A,i)

# length of i-the column/row
@inline collength(A::AbstractVecOrMat, i::Integer) = max(colstop(A, i) - colstart(A, i) + 1, 0)
@inline rowlength(A::AbstractVecOrMat, i::Integer) = max(rowstop(A, i) - rowstart(A, i) + 1, 0)

"""
    isbanded(A)

returns true if a matrix implements the banded interface.
"""
isbanded(::AbstractBandedMatrix) = true
isbanded(_) = false

# override bandwidth(A,k) for each AbstractBandedMatrix
# override inbands_getindex(A,k,j)

# return id of first empty diagonal intersected along row k
function _firstdiagrow(A::AbstractMatrix, k::Int)
    a, b = rowstart(A, k), rowstop(A, k)
    c = a == 1 ? b+1 : a-1
    c-k
end

# return id of first empty diagonal intersected along column j
function _firstdiagcol(A::AbstractMatrix, j::Int)
    a, b = colstart(A, j), colstop(A, j)
    r = a == 1 ? b+1 : a-1
    j-r
end

function Base.maximum(B::AbstractBandedMatrix)
    m=zero(eltype(B))
    for j = 1:size(B,2), k = colrange(B,j)
        m=max(B[k,j],m)
    end
    m
end

# fallbacks for inbands_getindex and inbands_setindex!
@inline function inbands_getindex(x::AbstractMatrix, i::Integer, j::Integer)
    @inbounds r = getindex(x, i, j)
    r
end
@inline function inbands_setindex!(x::AbstractMatrix, v, i::Integer, j::Integer)
    @inbounds r = setindex!(x, v, i, j)
    r
end

inbands_getindex(x::Adjoint, i::Integer, j::Integer) =
    inbands_getindex(parent(x), j, i)'
inbands_getindex(x::Transpose, i::Integer, j::Integer) =
    transpose(inbands_getindex(parent(x), j, i))
inbands_setindex!(x::Adjoint, v, i::Integer, j::Integer) =
    inbands_setindex!(parent(x), v', j, i)
inbands_setindex!(x::Transpose, v, i::Integer, j::Integer) =
    inbands_setindex!(parent(x), transpose(v), j, i)
## Show

## structured matrix methods ##
Base.replace_in_print_matrix(A::AbstractBandedMatrix, i::Integer, j::Integer, s::AbstractString) =
    -bandwidth(A,1) ≤ j-i ≤ bandwidth(A,2) ? s : Base.replace_with_centered_mark(s)

## @inbands

function inbands_process_args!(e)
    if e.head == :ref
        e.head = :call
        pushfirst!(e.args, :(BandedMatrices.inbands_getindex))
    elseif e.head == :call && e.args[1] == :getindex
        e.args[1] = :(BandedMatrices.inbands_getindex)
    end
    e
end

function inbands_process_args_recursive!(expr)
    for (i,e) in enumerate(expr.args)
        if e isa Expr
            inbands_process_args_recursive!(e)
        end
    end
    inbands_process_args!(expr)
    expr
end

macro inbands(expr)
    esc(inbands_process_args!(expr))
end


text:
using DesignAutomation

include("functions.jl")

@testset "Exhaustive Search" begin
    x_range = [-2, 6]

    (x, fx, f_calls) = exhaustive_search(levy3, x_range, 4; verbose=false)

    @test x == [1.0, 1.0, 1.0, -2.0]
    @test fx + 1 ≈ 1 # answer is approx 0
    @test f_calls == 81^4
end

max_stars_repo_name: camirmas/DesignAutomation

max_issues_repo_name: camirmas/DesignAutomation

max_forks_repo_name: camirmas/DesignAutomation


text:
using Gadfly

set_default_plot_size(6inch, 16inch)

xs = [0.236033, 0.346517, 0.312707, 0.00790928, 0.488613, 0.210968, 0.951916, 0.999905, 0.251662, 0.986666]
ys = [0.555751, 0.437108, 0.424718, 0.773223, 0.28119, 0.209472, 0.251379, 0.0203749, 0.287702, 0.859512]
labels = String["mxhEmy", "9APlX1", "P17KC2", "2zNjTz", "g4mUef", "rcfI9s", "Bjd19K", "V1jjSR", "DEBHJk", "0J7d2j"]

l = plot(x=xs, y=ys, label=labels, Geom.point, Geom.label(position=:left))
r = plot(x=xs, y=ys, label=labels, Geom.point, Geom.label(position=:right))
a = plot(x=xs, y=ys, label=labels, Geom.point, Geom.label(position=:above))
b = plot(x=xs, y=ys, label=labels, Geom.point, Geom.label(position=:below))
c = plot(x=xs, y=ys, label=labels, Geom.point, Geom.label(position=:centered))

vstack(l, r, a, b, c)


text:
@testset "Closeness" begin
    g5 = SimpleDiGraph(4)
    add_edge!(g5, 1, 2); add_edge!(g5, 2, 3); add_edge!(g5, 1, 3); add_edge!(g5, 3, 4)

    for g in testdigraphs(g5)
        y = @inferred(closeness_centrality(g; normalize=false))
        z = @inferred(closeness_centrality(g))
        @test y == [0.75, 0.6666666666666666, 1.0, 0.0]
        @test z == [0.75, 0.4444444444444444, 0.3333333333333333, 0.0]
    end

    adjmx2 = [0 1 0; 1 0 1; 1 1 0] # digraph
    a2 = SimpleDiGraph(adjmx2)
    for g in testdigraphs(a2)
        distmx2 = [Inf 2.0 Inf; 3.2 Inf 4.2; 5.5 6.1 Inf]
        c2 = [0.24390243902439027, 0.27027027027027023, 0.1724137931034483]
        y = @inferred(closeness_centrality(g, distmx2; normalize=false))
        z = @inferred(closeness_centrality(g, distmx2))
        @test isapprox(y, c2)
        @test isapprox(z, c2)
    end

    g5 = SimpleGraph(5)
    add_edge!(g5, 1, 2)
    for g in testgraphs(g5)
        z = @inferred(closeness_centrality(g))
        @test z[1] == z[2] == 0.25
        @test z[3] == z[4] == z[5] == 0.0
    end

    adjmx1 = [0 1 0; 1 0 1; 0 1 0] # graph
    a1 = SimpleGraph(adjmx1)
    for g in testgraphs(a1)
        distmx1 = [Inf 2.0 Inf; 2.0 Inf 4.2; Inf 4.2 Inf]
        c1 = [0.24390243902439027, 0.3225806451612903, 0.1923076923076923]
        y = @inferred(closeness_centrality(g, distmx1; normalize=false))
        z = @inferred(closeness_centrality(g, distmx1))
        @test isapprox(y, c1)
        @test isapprox(z, c1)
    end
end


text:
using AbstractGPs
using BlockArrays
using Documenter
using FiniteDifferences
using LinearAlgebra
using KernelFunctions
using Random
using Statistics
using Stheno
using Test
using TimerOutputs
using Zygote

using Stheno: mean, cov, var, GPC, AV

const to = TimerOutput()

macro timedtestset(name, code)
    return esc(:(@timeit to $name @testset $name $code))
end

include("test_util.jl")

@testset "Stheno" begin

    println("util:")
    @timedtestset "util" begin
        include(joinpath("util", "zygote_rules.jl"))
        include(joinpath("util", "covariance_matrices.jl"))
        @testset "block_arrays" begin
            include(joinpath("util", "block_arrays", "test_util.jl"))
            include(joinpath("util", "block_arrays", "dense.jl"))
            include(joinpath("util", "block_arrays", "diagonal.jl"))
        end
        include(joinpath("util", "abstract_data_set.jl"))
    end

    println("gp:")
    @timedtestset "gp" begin
        include(joinpath("gp", "gp.jl"))
    end

    println("composite:")
    @timedtestset "composite" begin
        include(joinpath("composite", "test_util.jl"))
        include(joinpath("composite", "cross.jl"))
        include(joinpath("composite", "product.jl"))
        include(joinpath("composite", "addition.jl"))
        include(joinpath("composite", "compose.jl"))
    end

    println("abstract_gp:")
    @timedtestset "abstract_gp" begin
        include("abstract_gp.jl")
    end

    println("gaussian_process_probabilistic_programme:")
    include("gaussian_process_probabilistic_programme.jl")

    println("sparse_finite_gp:")
    include("sparse_finite_gp.jl")

    println("doctests")
    @timedtestset "doctests" begin
        DocMeta.setdocmeta!(
            Stheno,
            :DocTestSetup,
            :(using AbstractGPs, Stheno, Random, Documenter, LinearAlgebra);
            recursive=true,
        )
        doctest(Stheno)
    end
end

display(to)


text:
struct TextHistogram{P<:HistogramParallelization}
    table::AbstractDict
    subtable::AbstractVector{<:AbstractDict}
end

BinType(::TextHistogram{P}) where {P} = UnboundedWidth()
BinSearchAlgorithm(::TextHistogram{P}) where {P} = HashFunction()
HistogramParallelization(::TextHistogram{P}) where {P} = P()
eltype(::TextHistogram{P}) where {P} = String

"""
    create_fast_histogram(::UnboundedWidth, ::HashFunction, ::P) where {P<:HistogramParallelization}

Creates a histogram for 1D text data.
`P` can be any parallelization scheme.
"""
function create_fast_histogram(::UnboundedWidth, ::HashFunction, ::P) where {P<:HistogramParallelization}
    table = Dict{String,Int}()
    subtable = map(x -> Dict{String,Int}(), 1:Threads.nthreads())
    TextHistogram{P}(table, subtable)
end

get_table(h::TextHistogram) = h.table

# For PrivateThreads and SIMD
get_subtable(h::TextHistogram) = h.subtable

function zero!(h::TextHistogram)
    empty!(get_table(h))
    empty!.(get_subtable(h))
    nothing
end

counts(h::TextHistogram) = h.table


text:
module IndexedDims

using AcceleratedArrays

export IndexedDimsArray, bypass

include("autoaccelerate.jl")

struct BypassIndex{T}
    index::T
end

bypass(val::T) where {T} = BypassIndex{T}(val)

struct IndexedDimsArray{T, N, A<:AbstractArray{T, N}, Is<:Tuple{Vararg{<:AcceleratedVector,N}}} <: AbstractArray{T, N}
    data::A
    indexes::Is
end

Base.parent(arr::IndexedDimsArray) = arr.data
Base.size(arr::IndexedDimsArray) = size(parent(arr))

IndexedDimsArray(data::AbstractArray) = IndexedDimsArray(data, axes(data)...)

function IndexedDimsArray(data::AbstractArray{T, N}, dim_inds::Vararg{AcceleratedVector, N}) where {T, N}
    return IndexedDimsArray(data, dim_inds)
end

function IndexedDimsArray(data::AbstractArray{T, N}, dim_inds::Vararg{AbstractVector, N}) where {T, N}
    return IndexedDimsArray(data, map(autoaccelerate, dim_inds))
end

### simple fallback indexing (might delete later)
function Base.getindex(arr::IndexedDimsArray{T, N}, inds::Vararg{Int, N}) where {T, N}
    return Base.getindex(parent(arr), inds...)
end

function Base.setindex!(arr::IndexedDimsArray{T, N}, val, inds::Vararg{Int, N}) where {T, N}
    return Base.setindex!(parent(arr), val, inds...)
end
###

function Base.similar(arr::IndexedDimsArray{T, N, A}, inds::Tuple{Vararg{AbstractUnitRange, N}}) where {T, N, A}
    return IndexedDimsArray(similar(parent(arr), inds), map(reindex, arr.indexes, inds))
end

include("stdindexing.jl")  # a wrapper for non-value (regular) indexing (hypothetical)
include("dimindexing.jl")  # value indexing

end # module


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

"""
    SubString(s::AbstractString, i::Integer, j::Integer=lastindex(s))
    SubString(s::AbstractString, r::UnitRange{<:Integer})

Like [`getindex`](@ref), but returns a view into the parent string `s`
within range `i:j` or `r` respectively instead of making a copy.

# Examples
```jldoctest
julia> SubString("abc", 1, 2)
"ab"

julia> SubString("abc", 1:2)
"ab"

julia> SubString("abc", 2)
"bc"
```
"""
struct SubString{T<:AbstractString} <: AbstractString
    string::T
    offset::Int
    ncodeunits::Int

    function SubString{T}(s::T, i::Int, j::Int) where T<:AbstractString
        i ≤ j || return new(s, 0, 0)
        @boundscheck begin
            checkbounds(s, i:j)
            @inbounds isvalid(s, i) || string_index_err(s, i)
            @inbounds isvalid(s, j) || string_index_err(s, j)
        end
        return new(s, i-1, nextind(s,j)-i)
    end
end

SubString(s::T, i::Int, j::Int) where {T<:AbstractString} = SubString{T}(s, i, j)
SubString(s::AbstractString, i::Integer, j::Integer=lastindex(s)) = SubString(s, Int(i), Int(j))
SubString(s::AbstractString, r::UnitRange{<:Integer}) = SubString(s, first(r), last(r))

function SubString(s::SubString, i::Int, j::Int)
    @boundscheck i ≤ j && checkbounds(s, i:j)
    SubString(s.string, s.offset+i, s.offset+j)
end

SubString(s::AbstractString) = SubString(s, 1, lastindex(s))
SubString{T}(s::T) where {T<:AbstractString} = SubString{T}(s, 1, lastindex(s))

convert(::Type{SubString{S}}, s::AbstractString) where {S<:AbstractString} =
    SubString(convert(S, s))
convert(::Type{T}, s::T) where {T<:SubString} = s

String(s::SubString{String}) = unsafe_string(pointer(s.string, s.offset+1), s.ncodeunits)

ncodeunits(s::SubString) = s.ncodeunits
codeunit(s::SubString) = codeunit(s.string)
length(s::SubString) = length(s.string, s.offset+1, s.offset+s.ncodeunits)

function codeunit(s::SubString, i::Integer)
    @boundscheck checkbounds(s, i)
    @inbounds return codeunit(s.string, s.offset + i)
end

function iterate(s::SubString, i::Integer=firstindex(s))
    i == ncodeunits(s)+1 && return nothing
    @boundscheck checkbounds(s, i)
    y = iterate(s.string, s.offset + i)
    y === nothing && return nothing
    c, i = y
    return c, i - s.offset
end

function getindex(s::SubString, i::Integer)
    @boundscheck checkbounds(s, i)
    @inbounds return getindex(s.string, s.offset + i)
end

function isvalid(s::SubString, i::Integer)
    ib = true
    @boundscheck ib = checkbounds(Bool, s, i)
    @inbounds return ib && isvalid(s.string, s.offset + i)
end

thisind(s::SubString{String}, i::Int) = _thisind_str(s, i)
nextind(s::SubString{String}, i::Int) = _nextind_str(s, i)

function cmp(a::SubString{String}, b::SubString{String})
    na = sizeof(a)
    nb = sizeof(b)
    c = ccall(:memcmp, Int32, (Ptr{UInt8}, Ptr{UInt8}, UInt),
              pointer(a), pointer(b), min(na, nb))
    return c < 0 ? -1 : c > 0 ? +1 : cmp(na, nb)
end

# don't make unnecessary copies when passing substrings to C functions
cconvert(::Type{Ptr{UInt8}}, s::SubString{String}) = s
cconvert(::Type{Ptr{Int8}}, s::SubString{String}) = s

function unsafe_convert(::Type{Ptr{R}}, s::SubString{String}) where R<:Union{Int8, UInt8}
    convert(Ptr{R}, pointer(s.string)) + s.offset
end

pointer(x::SubString{String}) = pointer(x.string) + x.offset
pointer(x::SubString{String}, i::Integer) = pointer(x.string) + x.offset + (i-1)

"""
    reverse(s::AbstractString) -> AbstractString

Reverses a string. Technically, this function reverses the codepoints in a string and its
main utility is for reversed-order string processing, especially for reversed
regular-expression searches. See also [`reverseind`](@ref) to convert indices in `s` to
indices in `reverse(s)` and vice-versa, and `graphemes` from module `Unicode` to
operate on user-visible "characters" (graphemes) rather than codepoints.
See also [`Iterators.reverse`](@ref) for
reverse-order iteration without making a copy. Custom string types must implement the
`reverse` function themselves and should typically return a string with the same type
and encoding. If they return a string with a different encoding, they must also override
`reverseind` for that string type to satisfy `s[reverseind(s,i)] == reverse(s)[i]`.

# Examples
```jldoctest
julia> reverse("JuliaLang")
"gnaLailuJ"

julia> reverse("ax̂e") # combining characters can lead to surprising results
"êxa"

julia> using Unicode

julia> join(reverse(collect(graphemes("ax̂e")))) # reverses graphemes
"ex̂a"
```
"""
function reverse(s::Union{String,SubString{String}})::String
    sprint(sizehint=sizeof(s)) do io
        i, j = firstindex(s), lastindex(s)
        while i ≤ j
            c, j = s[j], prevind(s, j)
            write(io, c)
        end
    end
end

getindex(s::AbstractString, r::UnitRange{<:Integer}) = SubString(s, r)

max_stars_repo_name: ninjin/julia

max_issues_repo_name: ninjin/julia

max_forks_repo_name: ninjin/julia


text:
""" A small pooling problem. """
function pool1(gm::Bool = false)
    m = JuMP.Model(with_optimizer(CPLEX_SILENT))
    @variable(m, 0 <= x[1:7])
    lbs = Dict(x[2] => 3, x[3] => 1, x[4] => 2)
    [JuMP.set_lower_bound(key, value) for (key, value) in lbs]
    ubs = Dict(x[1] => 10, x[2] => 20, x[3] => 2, x[4] => 4, x[5] => 10, x[6] => 201, x[7] => 100)
    [JuMP.set_upper_bound(key, value) for (key, value) in ubs]
    @objective(m, Min, 6*x[1] + 16*x[2] -9*x[3] - 10*x[4])
    eqs = Dict("e1" => 12 - x[3]^2+x[4]^2,
           "e2_1" => x[1]^2-x[2]^2+x[4]^2 - 3, "e2_2" => 100 - x[1]^2-x[2]^2+x[4]^2,
           "e3_1" => x[1]+x[2]-5*x[3]+2*x[4] - 10, "e3_2" => 20 - x[1]+x[2]-5*x[3]+2*x[4],
           "e4" => -x[4] + x[6]^2 - 4*x[7]^2)
    if gm
        gm = GlobalModel(model = m, name = "pool1")
        add_nonlinear_constraint(gm, :(x -> 12 - x[3]^2+x[4]^2), name = "e1")
        add_nonlinear_constraint(gm, :(x -> x[1]^2-x[2]^2+x[4]^2 - 3), name = "e2_1")
        add_nonlinear_constraint(gm, :(x -> 100 - x[1]^2-x[2]^2+x[4]^2), name = "e2_2")
        add_nonlinear_or_compatible(gm, :(x -> x[1]+x[2]-5*x[3]+2*x[4] - 10), name = "e3_1")
        add_nonlinear_or_compatible(gm, :(x -> 20 - x[1]+x[2]-5*x[3]+2*x[4]), name = "e3_2")
        add_nonlinear_constraint(gm, :(x -> -x[4] + x[6]^2 - 4*x[7]^2), name = "e4")
        return gm
    else
        for (key, eq) in eqs
            @constraint(m, eq >= 0)
        end
        set_optimizer(m, BARON_SILENT)
        return m
    end
end


text:
using CSV
using DataFrames
using Gadfly
using Cairo
using Fontconfig

function plot(data::DataFrame, variable)
    p = Gadfly.plot(
        data,
        Geom.point,
        Geom.line,
        Scale.x_log10,
        Scale.y_log10,
        Theme(; background_color="white");
        x=:N,
        y=variable,
        color=:OS
    )

    draw(PNG(joinpath(pwd(), "output", "comparison", "$variable.png"), 15cm, 9cm; dpi=300), p)
end

# TODO: change paths

data1 = CSV.File(joinpath(pwd(), "output", "linux.csv")) |> DataFrame
data2 = CSV.File(joinpath(pwd(), "output", "windows.csv")) |> DataFrame

data = vcat(data1, data2, source=:OS => ["linux", "windows"])

plot(data, :Time)
plot(data, :Space)
plot(data, :Error)

max_stars_repo_name: saiteki-kai/direct-methods

max_issues_repo_name: saiteki-kai/direct-methods

max_forks_repo_name: saiteki-kai/direct-methods


text:
"""
    MeshTetrahedronModule

Module  for generation of meshes composed of tetrahedra.
"""
module MeshTetrahedronModule

using FinEtools.FTypesModule: FInt, FFlt, FCplxFlt, FFltVec, FIntVec, FFltMat, FIntMat, FMat, FVec, FDataDict
import FinEtools.FESetModule: count, FESetT4, FESetT10, setlabel!
import FinEtools.FENodeSetModule: FENodeSet
import FinEtools.MeshUtilModule: makecontainer, addhyperface!, findhyperface!, linearspace
import FinEtools.MeshSelectionModule: selectelem

"""
    T4block(Length::FFlt, Width::FFlt, Height::FFlt,
       nL::FInt, nW::FInt, nH::FInt, orientation::Symbol)

Generate a tetrahedral mesh  of the 3D block.

Four-node tetrahedra in a regular arrangement, with uniform spacing between
the nodes, with a given orientation of the diagonals.

The mesh is produced by splitting each logical  rectangular cell into six
tetrahedra.
Range =<0, Length> x <0, Width> x <0, Height>
Divided into elements: nL,  nW,  nH in the first,  second,  and
third direction (x, y, z).
"""
function T4block(Length::FFlt, Width::FFlt, Height::FFlt,
  nL::FInt, nW::FInt, nH::FInt, orientation::Symbol=:a)
    return T4blockx(collect(linearspace(0.0, Length, nL+1)),
                    collect(linearspace(0.0, Width, nW+1)),
                    collect(linearspace(0.0, Height, nH+1)), orientation);
end

"""
    T4blockx(xs::FFltMat, ys::FFltMat, zs::FFltMat, orientation::Symbol)

Generate a graded tetrahedral mesh  of a 3D block.

Four-node tetrahedra in a regular arrangement, with non-uniform given spacing
between the nodes, with a given orientation of the diagonals.

The mesh is produced by splitting each logical  rectangular cell into six
tetrahedra.
"""
function T4blockx(xs::FFltMat, ys::FFltMat, zs::FFltMat, orientation::Symbol)
    return T4blockx(vec(xs), vec(ys), vec(zs), orientation)
end

"""
    T4blockx(xs::FFltVec, ys::FFltVec, zs::FFltVec, orientation::Symbol)

Generate a graded tetrahedral mesh  of a 3D block.

Four-node tetrahedra in a regular arrangement, with non-uniform given spacing
between the nodes, with a given orientation of the diagonals.

The mesh is produced by splitting each logical  rectangular cell into six
tetrahedra.
"""
function T4blockx(xs::FFltVec, ys::FFltVec, zs::FFltVec, orientation::Symbol)
    nL =length(xs)-1;
    nW =length(ys)-1;
    nH =length(zs)-1;
    nnodes=(nL+1)*(nW+1)*(nH+1);
    ncells=6*(nL)*(nW)*(nH);
    xyzs=zeros(FFlt, nnodes, 3);
    conns=zeros(FInt, ncells, 4);
    if (orientation==:a)
        t4ia = [1  8  5  6; 3  4  2  7; 7  2  6  8; 4  7  8  2; 2  1  6  8; 4  8  1  2];
        t4ib = [1  8  5  6; 3  4  2  7; 7  2  6  8; 4  7  8  2; 2  1  6  8; 4  8  1  2];
    elseif (orientation==:b)
        t4ia = [2 7 5 6; 1 8 5 7; 1 3 4 8; 2 1 5 7; 1 2 3 7; 3 7 8 1];
        t4ib = [2 7 5 6; 1 8 5 7; 1 3 4 8; 2 1 5 7; 1 2 3 7; 3 7 8 1];
    elseif (orientation==:ca)
        t4ia = [8  4  7  5; 6  7  2  5; 3  4  2  7; 1  2  4  5; 7  4  2  5];
        t4ib = [7  3  6  8; 5  8  6  1; 2  3  1  6; 4  1  3  8; 6  3  1  8];
    elseif (orientation==:cb)
        t4ia = [7  3  6  8; 5  8  6  1; 2  3  1  6; 4  1  3  8; 6  3  1  8];
        t4ib = [8  4  7  5; 6  7  2  5; 3  4  2  7; 1  2  4  5; 7  4  2  5];
    else
        error("Unknown orientation")
    end
    f=1;
    for k=1:(nH+1)
        for j=1:(nW+1)
            for i=1:(nL+1)
                xyzs[f, 1] = xs[i]
                xyzs[f, 2] = ys[j]
                xyzs[f, 3] = zs[k];
                f=f+1;
            end
        end
    end

    fens=FENodeSet(xyzs);

    function node_numbers(i::FInt, j::FInt, k::FInt, nL::FInt, nW::FInt, nH::FInt)
        f=(k-1)*((nL+1)*(nW+1))+(j-1)*(nL+1)+i;
        nn=[f (f+1)  f+(nL+1)+1 f+(nL+1)];
        return [nn broadcast(+, nn, (nL+1)*(nW+1))];
    end

    gc=1;
    for i=1:nL
        for j=1:nW
            for k=1:nH
                nn=node_numbers(i, j, k, nL, nW, nH);
                if (mod(sum( [i, j, k] ), 2)==0)
                    t4i =t4ib;
                else
                    t4i =t4ia;
                end
                for r=1:size(t4i, 1)
                    for c1=1:size(t4i, 2)
                        conns[gc, c1] = nn[t4i[r, c1]];
                    end
                    gc=gc+1;
                end
            end
        end
    end
    fes = FESetT4(conns[1:gc-1, :]);

    return fens, fes
end

"""
    T4toT10(fens::FENodeSet,  fes::FESetT4)

Convert a mesh of tetrahedra of type T4 (four-node) to tetrahedra T10.
"""
function  T4toT10(fens::FENodeSet,  fes::FESetT4)
    nedges=6;
    ec = [1  2; 2  3; 3  1; 4  1; 4  2; 4  3];
    # Additional node numbers are numbered from here
    newn = count(fens)+1;
    # make a search structure for edges
    edges = makecontainer();
    for i= 1:length(fes.conn)
        for J = 1:nedges
            ev=fes.conn[i][ec[J,:]];
            newn = addhyperface!(edges,  ev,  newn);
        end
    end
    xyz1 =fens.xyz;             # Pre-existing nodes
    # Allocate for vertex nodes plus edge nodes plus face nodes
    xyz =zeros(FFlt, newn-1, 3);
    xyz[1:size(xyz1, 1), :] = xyz1; # existing nodes are copied over
    # calculate the locations of the new nodes
    # and construct the new nodes
    for i in keys(edges)
        C=edges[i];
        for J = 1:length(C)
          ix = vec([item for item in C[J].o])
          push!(ix,  i) # Add the anchor point as well
          xyz[C[J].n, :] = mean(xyz[ix, :], dims = 1);
        end
    end
    fens = FENodeSet(xyz);
    # construct new geometry cells
    nconn=zeros(FInt, length(fes.conn), 10);
    nc=1;
    for i= 1:length(fes.conn)
        econn=zeros(FInt, 1, nedges);
        for J = 1:nedges
            ev=fes.conn[i][ec[J,:]];
            h, n=findhyperface!(edges,  ev);
            econn[J]=n;
        end
        wn = 1
        for j in fes.conn[i]
            nconn[nc, wn] = j; wn = wn + 1
        end
        for j in econn
            nconn[nc, wn] = j; wn = wn + 1
        end
        # nconn[nc, :] = vcat([j for j in fes.conn[i]], vec(econn))
        nc= nc+ 1;
    end
    labels = deepcopy(fes.label)
    fes = FESetT10(nconn);
    fes = setlabel!(fes, labels)
    return fens, fes;
end

"""
    T10block(Length::FFlt, Width::FFlt, Height::FFlt,
      nL::FInt, nW::FInt, nH::FInt; orientation::Symbol=:a)

Generate a tetrahedral  mesh of T10 elements  of a rectangular block.
"""
function T10block(Length::FFlt, Width::FFlt, Height::FFlt,
    nL::FInt, nW::FInt, nH::FInt; orientation::Symbol=:a)
    fens, fes = T4block(Length, Width, Height, nL, nW, nH, orientation);
    fens, fes = T4toT10(fens, fes);
    return fens, fes
end

"""
    T10blockx(xs::FFltMat, ys::FFltMat, zs::FFltMat, orientation::Symbol = :a)

Generate a graded 10-node tetrahedral mesh  of a 3D block.

10-node tetrahedra in a regular arrangement, with non-uniform given spacing
between the nodes, with a given orientation of the diagonals.

The mesh is produced by splitting each logical  rectangular cell into six
tetrahedra.
"""
function T10blockx(xs::FFltMat, ys::FFltMat, zs::FFltMat, orientation::Symbol = :a)
    fens, fes =  T4blockx(vec(xs), vec(ys), vec(zs), orientation)
    fens, fes = T4toT10(fens, fes);
    return fens, fes
end

function T10blockx(xs::FFltVec, ys::FFltVec, zs::FFltVec, orientation::Symbol = :a)
    fens, fes =  T4blockx(vec(xs), vec(ys), vec(zs), orientation)
    fens, fes = T4toT10(fens, fes);
    return fens, fes
end

"""
    T10layeredplatex(xs::FFltVec, ys::FFltVec, ts::FFltVec, nts::FIntVec,
        orientation::Symbol = :a)

T10 mesh for a layered block (composite plate) with specified in plane coordinates.

xs,ys =Locations of the individual planes of nodes.
ts= Array of layer thicknesses,
nts= array of numbers of elements per layer

The finite elements of each layer are labeled with the layer number, starting
from 1 at the bottom.
"""
function T10layeredplatex(xs::FFltVec, ys::FFltVec, ts::FFltVec, nts::FIntVec,
    orientation::Symbol = :a)
    tolerance = minimum(abs.(ts))/maximum(nts)/10.;
    @assert length(ts) >= 1
    @assert sum(nts) >= length(ts)
    zs = collect(linearspace(0.0, ts[1], nts[1]+1))
    for layer = 2:length(ts)
        oz = collect(linearspace(sum(ts[1:layer-1]), sum(ts[1:layer]), nts[layer]+1))
        zs = vcat(zs, oz[2:end])
    end
    fens, fes = T4blockx(xs, ys, zs, orientation);
    List = selectelem(fens, fes, box = [-Inf Inf -Inf Inf 0.0 ts[1]], inflate = tolerance)
    fes.label[List] .= 1
    for layer = 2:length(ts)
        List = selectelem(fens, fes, box = [-Inf Inf -Inf Inf sum(ts[1:layer-1]) sum(ts[1:layer])],
            inflate = tolerance)
        fes.label[List] .= layer
    end
    fens, fes =  T4toT10(fens, fes)
    return fens, fes
end

"""
    tetv(X)

Compute the volume of a tetrahedron.

```
X = [0  4  3
9  2  4
6  1  7
0  1  5] # for these points the volume is 10.0
tetv(X)
```
"""
function tetv(X::FFltMat)
    local one6th = 1.0/6
    # @assert size(X, 1) == 4
    # @assert size(X, 2) == 3
    @inbounds let
        A1 = X[2,1]-X[1,1];
        A2 = X[2,2]-X[1,2];
        A3 = X[2,3]-X[1,3];
        B1 = X[3,1]-X[1,1];
        B2 = X[3,2]-X[1,2];
        B3 = X[3,3]-X[1,3];
        C1 = X[4,1]-X[1,1];
        C2 = X[4,2]-X[1,2];
        C3 = X[4,3]-X[1,3];
        return one6th * ((-A3*B2+A2*B3)*C1 +  (A3*B1-A1*B3)*C2 + (-A2*B1+A1*B2)*C3);
    end
end

function tetv(v11::FFlt, v12::FFlt, v13::FFlt, v21::FFlt, v22::FFlt, v23::FFlt, v31::FFlt, v32::FFlt, v33::FFlt, v41::FFlt, v42::FFlt, v43::FFlt)
    local one6th = 1.0/6
    @inbounds let
        A1 = v21 - v11;
        A2 = v22 - v12;
        A3 = v23 - v13;
        B1 = v31 - v11;
        B2 = v32 - v12;
        B3 = v33 - v13;
        C1 = v41 - v11;
        C2 = v42 - v12;
        C3 = v43 - v13;
        return one6th * ((-A3*B2+A2*B3)*C1 +  (A3*B1-A1*B3)*C2 + (-A2*B1+A1*B2)*C3);
    end
end

"""
    tetv1times6(v, i1, i2, i3, i4)

Compute 6 times the volume of the tetrahedron.
"""
function tetv1times6(v::FFltMat, i1::Int, i2::Int, i3::Int, i4::Int)
    # local one6th = 1.0/6
    # @assert size(X, 1) == 4
    # @assert size(X, 2) == 3
    @inbounds let
        A1 = v[i2,1]-v[i1,1];
        A2 = v[i2,2]-v[i1,2];
        A3 = v[i2,3]-v[i1,3];
        B1 = v[i3,1]-v[i1,1];
        B2 = v[i3,2]-v[i1,2];
        B3 = v[i3,3]-v[i1,3];
        C1 = v[i4,1]-v[i1,1];
        C2 = v[i4,2]-v[i1,2];
        C3 = v[i4,3]-v[i1,3];
        return ((-A3*B2+A2*B3)*C1 +  (A3*B1-A1*B3)*C2 + (-A2*B1+A1*B2)*C3);
    end
end

"""
    T4meshedges(t::Array{Int, 2})

Compute all the edges of the 4-node triangulation.
"""
function T4meshedges(t::Array{Int, 2})
    @assert size(t, 2) == 4
    ec = [  1  2
            2  3
            3  1
            4  1
            4  2
            4  3];
    e = vcat(t[:,ec[1,:]], t[:,ec[2,:]], t[:,ec[3,:]], t[:,ec[4,:]], t[:,ec[5,:]], t[:,ec[6,:]])
    e = sort(e; dims = 2);
    ix = sortperm(e[:,1]);
    e = e[ix,:];
    ue = deepcopy(e)
    i = 1;
    n=1;
    while n <= size(e,1)
        c = ue[n,1];
        m = n+1;
        while m <= size(e,1)
            if (ue[m,1] != c)
                break;
            end
            m = m+1;
        end
        us = unique(ue[n:m-1,2], dims=1);
        ls =length(us);
        e[i:i+ls-1,1] .= c;
        e[i:i+ls-1,2] = sort(us);
        i = i+ls;
        n = m;
    end
    e = e[1:i-1,:];
end

# Construct arrays to describe a hexahedron mesh created from voxel image.
#
# img = 3-D image (array),  the voxel values  are arbitrary
# voxval =range of voxel values to be included in the mesh,
# voxval =  [minimum value,  maximum value].  Minimum value == maximum value is
# allowed.
# Output:
# t = array of hexahedron connectivities,  one hexahedron per row
# v =Array of vertex locations,  one vertex per row
function T4voximggen(img::Array{DataT, 3},  voxval::Array{DataT, 1}) where {DataT<:Number}
    M=size(img,  1); N=size(img,  2); P=size(img,  3);
    t4ia = [8 4 7 5; 6 7 2 5; 3 4 2 7; 1 2 4 5; 7 4 2 5];
    t4ib = [7 3 6 8; 5 8 6 1; 2 3 1 6; 4 1 3 8; 6 3 1 8];

    function find_nonempty(minvoxval, maxvoxval)
        Nvoxval=0
        for I= 1:M
            for J= 1:N
                for K= 1:P
                    if (img[I, J, K]>=minvoxval) && (img[I, J, K]<=maxvoxval)
                        Nvoxval=Nvoxval+1
                    end
                end
            end
        end
        return Nvoxval
    end
    minvoxval= minimum(voxval)  # include voxels at or above this number
    maxvoxval= maximum(voxval)  # include voxels at or below this number
    Nvoxval =find_nonempty(minvoxval, maxvoxval) # how many "full" voxels are there?

    # Allocate output arrays:  one voxel is converted to 5 tetrahedra
    t =zeros(FInt, 5*Nvoxval, 4);
    v =zeros(FInt, (M+1)*(N+1)*(P+1), 3);
    tmid =zeros(FInt, 5*Nvoxval);

    Slice =zeros(FInt, 2, N+1, P+1); # auxiliary buffer
    function find_vertex(I, IJK)
        vidx = zeros(FInt, 1, size(IJK, 1));
        for r= 1:size(IJK, 1)
            if (Slice[IJK[r, 1], IJK[r, 2], IJK[r, 3]]==0)
                nv=nv+1;
                v[nv, :] =IJK[r, :]; v[nv, 1] += I-1;
                Slice[IJK[r, 1], IJK[r, 2], IJK[r, 3]] =nv;
            end
            vidx[r] =Slice[IJK[r, 1], IJK[r, 2], IJK[r, 3]];
        end
        return vidx
    end
    function store_elements(I, J, K)
        locs =[1 J K;1+1 J K;1+1 J+1 K;1 J+1 K;1 J K+1;1+1 J K+1;1+1 J+1 K+1;1 J+1 K+1];
        vidx = find_vertex(I, locs);
        for r=1:5
        nt =nt +1;
        if (mod(sum( [I,J,K] ),2) == 0)
            t[nt,:] = vidx[t4ia[r,:]];
        else
            t[nt,:] = vidx[t4ib[r,:]];
        end
        tmid[nt] = convert(FInt, img[I, J, K]);
    end

    end

    nv =0;                      # number of vertices
    nt =0;                      # number of elements
    for I= 1:M
        for J= 1:N
            for K= 1:P
                if  (img[I, J, K]>=minvoxval) && (img[I, J, K]<=maxvoxval)
                    store_elements(I, J, K);
                end
            end
        end
        Slice[1, :, :] = Slice[2, :, :] ;
        Slice[2, :, :] .= 0;
    end
    # Trim output arrays
    v = v[1:nv, :];
    t = t[1:nt, :];
    tmid = tmid[1:nt];

    return t, v, tmid
end

"""
    T4voximg(img::Array{DataT, 3}, voxdims::FFltVec,
        voxval::Array{DataT, 1}) where {DataT<:Number}

Generate a tetrahedral mesh  from three-dimensional image.
"""
function T4voximg(img::Array{DataT, 3}, voxdims::FFltVec, voxval::Array{DataT, 1}) where {DataT<:Number}
    t, v, tmid = T4voximggen(img, voxval)
    xyz = zeros(FFlt, size(v, 1), 3)
    for k=1:3
        for j=1:size(v, 1)
            xyz[j, k] = v[j, k]*voxdims[k]
        end
    end
    fens  = FENodeSet(xyz);
    fes = FESetT4(t);
    setlabel!(fes, tmid)
    return fens, fes;
end

end


text:
using EchogramImages
using Test

@testset "IMAGESC" images
	a = rand(100,100)
	img = imagesc(a, nothing, nothing, ColorScheme)
	m,n = size(img)
	@test m == 480
	@test n == 640
	@test cmap == cmap.colors
	@test typeof(vmin) == AbstractFloat
	@test typeof(vmax) == AbstractFloat
	@test vmin <= @. img <= vmax

	b = [0.1 0.23 0.33; 0.22 0.16 0.45; 0.55 0.8 0.76]
	img = imagesc(b, 0.15, 0.75, MyColors)
	x,y = size(img)
	@test y == 3
	@test x == 3
	@test typeof(vmin) == AbstractFloat
	@test typeof(vmax) == AbstractFloat
	@test vmin == 0.1
	@test vmax == 0.8
	@test vmin <= @. img <= vmax

	c = [0.12 0.54 0.1; 0.23 0.16 0.55; 0.25 0.0 0.16]
	img = imagesc(c, nothing, nothing, ColorScheme, size=(4,4))
	p,q = size(img)
	@test p = 4
	@test q = 4
	@test b[0][3] = vmin
	@test b[1][3] = vmin
	@test b[2][3] = vmin
	@test b[3]]0] = vmin
	@test b[3][1] = vmin
	@test b[3][2] = vmin
	@test typeof(vmin) == AbstractFloat
	@test typeof(vmax) == AbstractFloat
	@test vmin == 0.0
	@test vmax == 0.55
	@test vmin <= @. img <= vmax
end


text:
using HTTP

function get_url(testnet::Bool=false)
    string("http://", NODE_URL, ":", DEFAULT["rpcport"][testnet])
end

init_url(url::String, testnet::Bool) = isempty(url) ? get_url(testnet) : url

"""
    get_tx(key::String; url::String="", testnet::Bool=false)
    -> Tx

Returns the bitcoin transaction given a node url with REST server enabled and
transaction hash as an hexadecimal string.
Uses mainnet by default
"""
function get_tx(key::String; url::String="", testnet::Bool=false)
    url *= init_url(url, testnet) *"/rest/tx/" * key * ".bin"
    response = HTTP.request("GET", url)
    try
        response.status == 200
    catch
        error("Unexpected status: ", response.status)
    end
    raw = response.body
    tx = parse(IOBuffer(raw), testnet)::Tx
    if tx.segwit
        computed = id(tx)
    else
        computed = bytes2hex(reverse(hash256(raw)))
    end
    if id(tx) != key
        error("not the same id : ", id(tx),
            "\n             vs : ", tx_id)
    end
    return tx
end

"""
    get_headers(key::String; amount::Integer=1, url::String="", testnet::Bool=false)
    -> BlockHeaders[]

Returns the bitcoin transaction given a node url with REST server enabled and
transaction hash as an hexadecimal string.
"""
function get_headers(key::String; amount::Integer=1, url::String="", testnet::Bool=false)
    url *= init_url(url, testnet) * "/rest/headers/" * string(amount) * "/" * key * ".bin"
    response = HTTP.request("GET", url)
    try
        response.status == 200
    catch
        error("Unexpected status: ", response.status)
    end
    io = IOBuffer(response.body)
    headers = BlockHeader[]
    while io.ptr < io.size
        push!(headers, io2blockheader(io))
    end
    return headers
end

"""
    get_blockhashbyheight(key::Integer; url::String="", testnet::Bool=false)
    -> String

Returns the hash of the block in the current best blockchain based on its
height (how many blocks it is after the Genesis Block).
"""
function get_block(key::Integer; url::String="", testnet::Bool=false)
    url *= init_url(url, testnet) *"/rest/blockhashbyheight/" * string(key) * ".bin"
    response = HTTP.request("GET", url)
    try
        response.status == 200
    catch
        error("Unexpected status: ", response.status)
    end
    return bytes2hex(reverse(response.body))
end

"""
    get_blockhashbyheight(key::Integer; url::String="", testnet::Bool=false)
    -> String

Returns the hash of the block in the current best blockchain based on its
height (how many blocks it is after the Genesis Block).
"""
function get_blockhashbyheight(key::Integer; url::String="", testnet::Bool=false)
    url *= init_url(url, testnet) *"/rest/blockhashbyheight/" * string(key) * ".bin"
    response = HTTP.request("GET", url)
    try
        response.status == 200
    catch
        error("Unexpected status: ", response.status)
    end
    return bytes2hex(reverse(response.body))
end

@deprecate txfetch(tx_id::String, testnet::Bool=false) fetch(tx_id::String, testnet::Bool)
@deprecate fetch(tx_id::String, testnet::Bool=false) get_tx(tx_id, testnet=testnet)