text:
"""
    is_overpass_server_availabile()

Returns `true` if server is available, `false` if not.
"""
function is_overpass_server_availabile()
    response = String(HTTP.get(OSM_URLS[:overpass_status]).body)
    return occursin("slots available now", response)
end

"""
Checks the availability of Overpass servers by sending a request to https://overpass-api.de/api/status.
"""
function check_overpass_server_availability()
    if is_overpass_server_availabile()
        @info "Overpass server is available for download"
    else
        throw(ErrorException("Overpass server is NOT available for download, please wait and try again"))
    end
end

"""
    nominatim_request(query::Dict)::String

Sends a GET request to the Nomatim API: https://nominatim.openstreetmap.org/search.

# Arguments
- `query::Dict`: HTTP GET request query aguments.

# Return
- `String`: Response body as string.
"""
function nominatim_request(query::Dict)::String
    return String(HTTP.get(OSM_URLS[:nominatim_search], query=query).body)
end

"""
    nominatim_request(data::Dict)::String

Sends a POST request to the Overpass API: http://overpass-api.de/api/interpreter.

# Arguments
- `data::String`: HTTP POST data.

# Return
- `String`: Response body as string.
"""
function overpass_request(data::String)::String
    check_overpass_server_availability()
    return String(HTTP.post(OSM_URLS[:overpass_map], body=data).body)
end

"""
    nominatim_polygon_query(place_name::String)::Dict{String,Any}

Forms query arguments required for a Nomatim search at https://nominatim.openstreetmap.org/search.

# Arguments
- `place_name::String`: Place name string used in request to the nominatim api.

# Return
- `Dict{String,Any}`: GET request query arguments.
"""
function nominatim_polygon_query(place_name::String)::Dict{String,Any}
    return Dict(
        "format" => "json",
        "limit" => 5,
        "dedupe" => 0,
        "polygon_geojson" => 1,
        "q" => place_name
    )
end

"""
    overpass_query(filters::String,
                   metadata::Bool=false,
                   download_format::Symbol=:json,
                   bbox::Union{Vector{AbstractFloat},Nothing}=nothing
                   )::String

Forms an Overpass query string. For a guide on the OSM query language, see https://wiki.openstreetmap.org/wiki/Overpass_API/Language_Guide.

# Arguments
- `filters::String`: Filters for the query, e.g. polygon filter, highways only, traffic lights only, etc.
- `metadata::Bool=false`: Set true to return metadata.
- `download_format::Symbol=:json`: Download format, either `:osm`, `:xml` or `json`.
- `bbox::Union{Vector{AbstractFloat},Nothing}=nothing`: Optional bounding box filter.

# Return
- `String`: Overpass query string.
"""
function overpass_query(filters::String,
                        metadata::Bool=false,
                        download_format::Symbol=:json,
                        bbox::Union{Vector{<:AbstractFloat},Nothing}=nothing
                        )::String
    download_format_str = """[out:$(OSM_DOWNLOAD_FORMAT[download_format])]"""
    bbox_str = bbox === nothing ? "" : """[bbox:$(replace("$bbox", r"[\[ \]]" =>  ""))]"""
    metadata_str = metadata ? "meta" : ""
    query = """$download_format_str[timeout:180]$bbox_str;($filters);out count;out $metadata_str;"""
    @debug """Making overpass query: $query"""
    return query
end

"""
    overpass_polygon_network_query(geojson_polygons::Vector{Vector{Any}},
                                   network_type::Symbol=:drive,
                                   metadata::Bool=false,
                                   download_format::Symbol=:json
                                   )::String

Forms an Overpass query string using geojosn polygon coordinates as a filter.

# Arguments
- `geojson_polygons::Vector{Vector{Any}}`: Vector of `[lat, lon, ...]` polygon coordinates.
- `network_type::Symbol=:drive`: Network type filter, pick from `:drive`, `:drive_service`, `:walk`, `:bike`, `:all`, `:all_private`, `:none`, `:rail`.
- `metadata::Bool=false`: Set true to return metadata.
- `download_format::Symbol=:json`: Download format, either `:osm`, `:xml` or `json`.

# Return
- `String`: Overpass query string.
"""
function overpass_polygon_network_query(geojson_polygons::AbstractVector{<:AbstractVector},
                                        network_type::Symbol=:drive,
                                        metadata::Bool=false,
                                        download_format::Symbol=:json
                                        )::String
    way_filter = WAY_FILTERS_QUERY[network_type]
    relation_filter = RELATION_FILTERS_QUERY[network_type]

    filters = ""
    for polygon in geojson_polygons
        polygon = map(x -> [x[2], x[1]], polygon) # switch lon-lat to lat-lon
        polygon_str = replace("$polygon", r"[\[,\]]" =>  "")
        filters *= """way$way_filter(poly:"$polygon_str");>;"""
        if !isnothing(relation_filter)
            filters *= """rel$relation_filter(poly:"$polygon_str");>;"""
        end
    end

    return overpass_query(filters, metadata, download_format)
end

"""
    polygon_from_place_name(place_name::String)::Vector{Vector{Any}}

Retrieves polygon coordinates using any place name string.

# Arguments
- `place_name::String`: Any place name string used as a search argument to the Nominatim API.

# Return
- `Vector{Vector{Any}}`: GeoJSON polygon coordiantes.
"""
function polygon_from_place_name(place_name::String)::Vector{Vector{Any}}
    query = nominatim_polygon_query(place_name)
    response = JSON.parse(nominatim_request(query))

    for item in response
        if item["geojson"]["type"] == "Polygon"
            @info "Using Polygon for $(item["display_name"])"
            return item["geojson"]["coordinates"]
        elseif item["geojson"]["type"] ==  "MultiPolygon"
            @info "Using MultiPolygon for $(item["display_name"])"
            return collect(Iterators.flatten(item["geojson"]["coordinates"]))
        end
    end

    throw(ErrorException("Could not find valid polygon for $place_name"))
end

"""
    osm_network_from_place_name(;place_name::String,
                                network_type::Symbol=:drive,
                                metadata::Bool=false,
                                download_format::Symbol=:json
                                )::String

Downloads an OpenStreetMap network using any place name string.

# Arguments
- `place_name::String`: Any place name string used as a search argument to the Nominatim API.
- `network_type::Symbol=:drive`: Network type filter, pick from `:drive`, `:drive_service`, `:walk`, `:bike`, `:all`, `:all_private`, `:none`, `:rail`.
- `metadata::Bool=false`: Set true to return metadata.
- `download_format::Symbol=:json`: Download format, either `:osm`, `:xml` or `json`.

# Return
- `String`: OpenStreetMap network data response string.
"""
function osm_network_from_place_name(;place_name::String,
                                     network_type::Symbol=:drive,
                                     metadata::Bool=false,
                                     download_format::Symbol=:json
                                     )::String
    geojson_polygons = polygon_from_place_name(place_name)
    query = overpass_polygon_network_query(geojson_polygons, network_type, metadata, download_format)
    return overpass_request(query)
end

"""
    osm_network_from_polygon(;polygon::AbstractVector,
                             network_type::Symbol=:drive,
                             metadata::Bool=false,
                             download_format::Symbol=:json
                             )::String

Downloads an OpenStreetMap network using a polygon.

# Arguments
- `polygon::AbstractVector`: Vector of longitude-latitude pairs.
- `network_type::Symbol=:drive`: Network type filter, pick from `:drive`, `:drive_service`, `:walk`, `:bike`, `:all`, `:all_private`, `:none`, `:rail`.
- `metadata::Bool=false`: Set true to return metadata.
- `download_format::Symbol=:json`: Download format, either `:osm`, `:xml` or `json`.

# Return
- `String`: OpenStreetMap network data response string.
"""
function osm_network_from_polygon(;polygon::AbstractVector{<:AbstractVector{<:Real}},
                                  network_type::Symbol=:drive,
                                  metadata::Bool=false,
                                  download_format::Symbol=:json
                                  )::String
    query = overpass_polygon_network_query([polygon], network_type, metadata, download_format)
    return overpass_request(query)
end

"""
    overpass_bbox_network_query(bbox::Vector{AbstractFloat},
                                network_type::Symbol=:drive,
                                metadata::Bool=false,
                                download_format::Symbol=:json
                                )::String

Forms an Overpass query string using a bounding box as a filter.

# Arguments
- `bbox::Vector{AbstractFloat}`: Vector of bounding box coordinates `[minlat, minlon, maxlat, maxlon]`.
- `network_type::Symbol=:drive`: Network type filter, pick from `:drive`, `:drive_service`, `:walk`, `:bike`, `:all`, `:all_private`, `:none`, `:rail`.
- `metadata::Bool=false`: Set true to return metadata.
- `download_format::Symbol=:json`: Download format, either `:osm`, `:xml` or `json`.

# Return
- `String`: Overpass query string.
"""
function overpass_bbox_network_query(bbox::Vector{<:AbstractFloat},
                                     network_type::Symbol=:drive,
                                     metadata::Bool=false,
                                     download_format::Symbol=:json
                                     )::String
    way_filter = WAY_FILTERS_QUERY[network_type]
    relation_filter = RELATION_FILTERS_QUERY[network_type]
    filters = "way$way_filter;>;"
    if !isnothing(relation_filter)
        filters *= """rel$relation_filter;>;"""
    end
    return overpass_query(filters, metadata, download_format, bbox)
end

"""
    osm_network_from_bbox(;minlat::AbstractFloat,
                          minlon::AbstractFloat,
                          maxlat::AbstractFloat,
                          maxlon::AbstractFloat,
                          network_type::Symbol=:drive,
                          metadata::Bool=false,
                          download_format::Symbol=:json
                          )::String

Downloads an OpenStreetMap network using bounding box coordinates.

# Arguments
- `minlat::AbstractFloat`: Bottom left bounding box latitude coordinate.
- `minlon::AbstractFloat`: Bottom left bounding box longitude coordinate.
- `maxlat::AbstractFloat`: Top right bounding box latitude coordinate.
- `maxlon::AbstractFloat`: Top right bounding box longitude coordinate.
- `network_type::Symbol=:drive`: Network type filter, pick from `:drive`, `:drive_service`, `:walk`, `:bike`, `:all`, `:all_private`, `:none`, `:rail`.
- `metadata::Bool=false`: Set true to return metadata.
- `download_format::Symbol=:json`: Download format, either `:osm`, `:xml` or `json`.

# Return
- `String`: OpenStreetMap network data response string.
"""
function osm_network_from_bbox(;minlat::AbstractFloat,
                               minlon::AbstractFloat,
                               maxlat::AbstractFloat,
                               maxlon::AbstractFloat,
                               network_type::Symbol=:drive,
                               metadata::Bool=false,
                               download_format::Symbol=:json
                               )::String
    query = overpass_bbox_network_query([minlat, minlon, maxlat, maxlon], network_type, metadata, download_format)
    return overpass_request(query)
end

"""
    osm_network_from_point(;point::GeoLocation,
                           radius::Number,
                           network_type::Symbol=:drive,
                           metadata::Bool=false,
                           download_format::Symbol=:json
                           )::String

Downloads an OpenStreetMap network using bounding box coordinates calculated from a centroid point and radius (km).

# Arguments
- `point::GeoLocation`: Centroid point to draw the bounding box around.
- `radius::Number`: Distance (km) from centroid point to each bounding box corner.
- `network_type::Symbol=:drive`: Network type filter, pick from `:drive`, `:drive_service`, `:walk`, `:bike`, `:all`, `:all_private`, `:none`, `:rail`.
- `metadata::Bool=false`: Set true to return metadata.
- `download_format::Symbol=:json`: Download format, either `:osm`, `:xml` or `json`.

# Return
- `String`: OpenStreetMap network data response string.
"""
function osm_network_from_point(;point::GeoLocation,
                                radius::Number,
                                network_type::Symbol=:drive,
                                metadata::Bool=false,
                                download_format::Symbol=:json
                                )::String
    bbox = bounding_box_from_point(point, radius)
    return osm_network_from_bbox(;bbox..., network_type=network_type, metadata=metadata, download_format=download_format)
end

"""
Factory method for retrieving download functions.
"""
function osm_network_downloader(download_method::Symbol)::Function
    if download_method == :place_name
        return osm_network_from_place_name
    elseif download_method == :bbox
        return osm_network_from_bbox
    elseif download_method == :point
        return osm_network_from_point
    elseif download_method == :polygon
        return osm_network_from_polygon
    else
        throw(ArgumentError("OSM network downloader $download_method does not exist"))
    end
end

"""
    download_osm_network(download_method::Symbol;
                         network_type::Symbol=:drive,
                         metadata::Bool=false,
                         download_format::Symbol=:json,
                         save_to_file_location::Union{String,Nothing}=nothing,
                         download_kwargs...
                         )::Union{XMLDocument,Dict{String,Any}}

Downloads an OpenStreetMap network by querying with a place name, bounding box, or centroid point.

# Arguments
- `download_method::Symbol`: Download method, choose from `:place_name`, `:bbox` or `:point`.
- `network_type::Symbol=:drive`: Network type filter, pick from `:drive`, `:drive_service`, `:walk`, `:bike`, `:all`, `:all_private`, `:none`, `:rail`
- `metadata::Bool=false`: Set true to return metadata.
- `download_format::Symbol=:json`: Download format, either `:osm`, `:xml` or `json`.
- `save_to_file_location::Union{String,Nothing}=nothing`: Specify a file location to save downloaded data to disk.

# Required Kwargs for each Download Method

*`download_method=:place_name`*
- `place_name::String`: Any place name string used as a search argument to the Nominatim API.

*`download_method=:bbox`*
- `minlat::AbstractFloat`: Bottom left bounding box latitude coordinate.
- `minlon::AbstractFloat`: Bottom left bounding box longitude coordinate.
- `maxlat::AbstractFloat`: Top right bounding box latitude coordinate.
- `maxlon::AbstractFloat`: Top right bounding box longitude coordinate.

*`download_method=:point`*
- `point::GeoLocation`: Centroid point to draw the bounding box around.
- `radius::Number`: Distance (km) from centroid point to each bounding box corner.

*`download_method=:polygon`*
- `polygon::AbstractVector`: Vector of longitude-latitude pairs.

# Network Types
- `:drive`: Motorways excluding private and service ways.
- `:drive_service`: Motorways including private and service ways.
- `:walk`: Walkways only.
- `:bike`: Cycleways only.
- `:all`: All motorways, walkways and cycleways excluding private ways.
- `:all_private`: All motorways, walkways and cycleways including private ways.
- `:none`: No network filters.
- `:rail`: Railways excluding proposed and platform.

# Return
- `Union{XMLDocument,Dict{String,Any}}`: OpenStreetMap network data parsed as either XML or Dictionary object depending on the download method.
"""
function download_osm_network(download_method::Symbol;
                              network_type::Symbol=:drive,
                              metadata::Bool=false,
                              download_format::Symbol=:json,
                              save_to_file_location::Union{String,Nothing}=nothing,
                              download_kwargs...
                              )::Union{XMLDocument,Dict{String,Any}}
    downloader = osm_network_downloader(download_method)
    data = downloader(network_type=network_type, metadata=metadata, download_format=download_format; download_kwargs...)
    @info "Downloaded osm network data from $(["$k: $v" for (k, v) in download_kwargs]) in $download_format format"

    if !(save_to_file_location isa Nothing)
        save_to_file_location = validate_save_location(save_to_file_location, download_format)
        write(save_to_file_location, data)
        @info "Saved osm network data to disk: $save_to_file_location"
    end

    deserializer = string_deserializer(download_format)

    return deserializer(data)
end


text:
"""
    ArraySpace(s::Space,dims...)

is used to represent array-valued expansions in a space `s`.  The
coefficients are of each entry are interlaced.

For example,
```julia
f = Fun(x->[exp(x),sin(x)],-1..1)
space(f) == ArraySpace(Chebyshev(),2)
```
"""
struct ArraySpace{S,n,DD,RR} <: DirectSumSpace{NTuple{n,S},DD,Array{RR,n}}
     spaces::Array{S,n}
end

const VectorSpace{S,DD,RR} = ArraySpace{S,1,DD,RR}
const MatrixSpace{S,DD,RR} = ArraySpace{S,2,DD,RR}

#TODO: Think through domain/domaindominsion
ArraySpace(sp::AbstractArray{SS,N}) where {SS<:Space,N} =
    ArraySpace{SS,N,domaintype(first(sp)),mapreduce(rangetype,promote_type,sp)}(sp)
ArraySpace(S::Space,n::NTuple{N,Int}) where {N} = ArraySpace(fill(S,n))
ArraySpace(S::Space,n::Integer) = ArraySpace(S,(n,))
ArraySpace(S::Space,n,m) = ArraySpace(fill(S,(n,m)))
ArraySpace(d::Domain,n...) = ArraySpace(Space(d),n...)

Space(sp::AbstractArray{<:Space}) = ArraySpace(sp)
convert(::Type{Array}, sp::ArraySpace) = sp.spaces
convert(::Type{Vector}, sp::VectorSpace) = sp.spaces
convert(::Type{Matrix}, sp::MatrixSpace) = sp.spaces
Array(sp::ArraySpace) = sp.spaces
Vector(sp::VectorSpace) = sp.spaces
Matrix(sp::MatrixSpace) = sp.spaces

BlockInterlacer(sp::ArraySpace) = BlockInterlacer(blocklengths.(tuple(sp.spaces...)))
interlacer(sp::ArraySpace) = BlockInterlacer(sp)

for OP in (:length,:firstindex,:lastindex,:size)
    @eval begin
        $OP(S::ArraySpace) = $OP(components(S))
        $OP(f::Fun{<:ArraySpace}) = $OP(space(f))
    end
end

for OP in (:getindex,:iterate,:stride,:size,:lastindex,:firstindex)
    @eval $OP(S::ArraySpace,k) = $OP(components(S),k)
end

iterate(S::ArraySpace) = iterate(components(S))
getindex(S::ArraySpace, kr::AbstractVector) = ArraySpace(components(S)[kr])

#support tuple set

stride(f::Fun{<:ArraySpace},k) = stride(space(f),k)

getindex(f::ArraySpace,k...) = Space(component(f,k...))
iterate(f::Fun{<:ArraySpace}) = iterate(components(f))

Base.reshape(AS::ArraySpace,k...) = ArraySpace(reshape(AS.spaces,k...))
dimension(AS::ArraySpace) = mapreduce(dimension,+,AS.spaces)

# TODO: union domain
domain(AS::ArraySpace) = domain(AS.spaces[1])
setdomain(A::ArraySpace,d::Domain) = ArraySpace(map(sp->setdomain(sp,d),A.spaces))

# support for Array of PiecewiseSpace

## transforms

#TODO: rework for different spaces
points(d::ArraySpace,n) = points(d.spaces[1],n)

transform(AS::ArraySpace{SS,1},vals::AbstractVector{Vector{V}}) where {SS,V} =
    transform(AS,transpose(hcat(vals...)))

function transform(AS::ArraySpace{SS,1,T},M::AbstractArray{V,2}) where {SS,T,V<:Number}
    n=length(AS)

    @assert size(M,2) == n
    plan = plan_transform(AS.spaces[1],M[:,1])
    cfs=Vector{V}[plan*M[:,k]  for k=1:size(M,2)]

    interlace(cfs,AS)
end

# transform of array is same order as vectorizing and then transforming
transform(AS::ArraySpace{SS,n},vals::AbstractVector{Array{V,n}}) where {SS,n,V} =
    transform(vec(AS),map(vec,vals))
transform(AS::VectorSpace{SS},vals::AbstractVector{AV}) where {SS,AV<:AbstractVector} =
    transform(AS,map(Vector,vals))
transform(AS::VectorSpace{SS},vals::AbstractVector{Vec{V,n}}) where {SS,n,V} =
    transform(AS,map(Vector,vals))

function itransform(AS::VectorSpace,cfs::AbstractVector)
    vf = vec(Fun(AS, cfs))
    n = maximum(ncoefficients, vf)
    vcat.(values.(pad!.(vf, n))...)
end

Base.vec(AS::ArraySpace) = ArraySpace(vec(AS.spaces))
Base.vec(f::Fun{ArraySpace{S,n,DD,RR}}) where {S,n,DD,RR} =
    [f[j] for j=1:length(f.space)]

repeat(A::ArraySpace,n,m) = ArraySpace(repeat(A.spaces,n,m))

component(A::MatrixSpace,k::Integer,j::Integer) = A.spaces[k,j]

Base.getindex(f::Fun{DSS},k::Integer) where {DSS<:ArraySpace} = component(f,k)

Base.getindex(f::Fun{MatrixSpace{S,DD,RR}},k::Integer,j::Integer) where {S,DD,RR} =
    f[k+stride(f,2)*(j-1)]

Base.getindex(f::Fun{DSS},kj::CartesianIndex{1}) where {DSS<:ArraySpace} = f[kj[1]]
Base.getindex(f::Fun{DSS},kj::CartesianIndex{2}) where {DSS<:ArraySpace} = f[kj[1],kj[2]]

function Fun(A::AbstractArray{Fun{VectorSpace{S,DD,RR},V,VV},2}) where {S,V,VV,DD,RR}
    @assert size(A,1)==1

    M = Matrix{Fun{S,V,VV}}(undef, length(space(A[1])),size(A,2))
    for k=1:size(A,2)
        M[:,k]=vec(A[k])
    end
    Fun(M)
end

# Fun{SS,n}(v::AbstractArray{Any,n},sp::ArraySpace{SS,n}) = Fun(map((f,s)->Fun(f,s),v,sp))

# convert a vector to a Fun with ArraySpace

function Fun(v::AbstractVector,sp::Space{D,R}) where {D,R<:AbstractVector}
    if size(v) ≠ size(sp)
        throw(DimensionMismatch("Cannot convert $v to a Fun in space $sp"))
    end
    Fun(map(Fun,v,components(sp)))
end

Fun(v::AbstractArray{TT,n},sp::Space{D,R}) where {D,R<:AbstractArray{SS,n}} where {TT,SS,n} =
    reshape(Fun(vec(v),vec(sp)),size(sp))

coefficients(v::AbstractArray{TT,n},sp::ArraySpace{SS,n}) where {TT,SS,n} = coefficients(Fun(v,sp))

for (OPrule,OP) in ((:conversion_rule,:conversion_type),(:maxspace_rule,:maxspace),
                        (:union_rule,:union))
    # ArraySpace doesn't allow reordering
    @eval function $OPrule(S1::ArraySpace,S2::ArraySpace)
        sps = map($OP,S1.spaces,S2.spaces)
        for s in sps
            if isa(s,NoSpace)
                return NoSpace()
            end
        end
        ArraySpace(sps)
    end
end

## routines

spacescompatible(AS::ArraySpace,BS::ArraySpace) =
    size(AS) == size(BS) && all(spacescompatible.(AS.spaces,BS.spaces))
canonicalspace(AS::ArraySpace) = ArraySpace(canonicalspace.(AS.spaces))
evaluate(f::AbstractVector,S::ArraySpace,x) = map(g->g(x),Fun(S,f))

## choosedomainspace

function choosedomainspace(A::InterlaceOperator{T,1},sp::ArraySpace) where T
    # this ensures correct dispatch for unino
    sps = Vector{Space}(
        filter(x->!isambiguous(x),map(choosedomainspace,A.ops,sp.spaces)))
    if isempty(sps)
        UnsetSpace()
    else
        union(sps...)
    end
end

Base.reshape(f::Fun{AS},k...) where {AS<:ArraySpace} = Fun(reshape(space(f),k...),f.coefficients)

Base.diff(f::Fun{AS,T},n...) where {AS<:ArraySpace,T} = Fun(diff(Array(f),n...))

## conversion

function coefficients(f::AbstractVector, a::VectorSpace, b::VectorSpace)
    if size(a) ≠ size(b)
        throw(DimensionMismatch("dimensions must match"))
    end
    interlace(map(coefficients,Fun(a,f),b),b)
end

coefficients(Q::AbstractVector{F},rs::VectorSpace) where {F<:Fun} =
    interlace(map(coefficients,Q,rs),rs)

coefficients(Q::AbstractVector, rs::VectorSpace) = coefficients(Fun.(Q), rs)

Fun(f::AbstractVector{FF},d::VectorSpace) where {FF<:Fun} = Fun(d,coefficients(f,d))
Fun(f::AbstractMatrix{FF},d::MatrixSpace) where {FF<:Fun} = Fun(d,coefficients(f,d))

## constructor

# columns are coefficients
function Fun(M::AbstractMatrix{<:Number},sp::MatrixSpace)
    if size(M) ≠ size(sp)
        throw(DimensionMismatch())
    end
    Fun(map((f,s)->Fun(f,s),M,sp.spaces))
end

Fun(M::UniformScaling,sp::MatrixSpace) = Fun(M.λ*Matrix(I,size(sp)...),sp)

ones(::Type{T},A::ArraySpace) where {T<:Number} = Fun(ones.(T,spaces(A)))
ones(A::ArraySpace) = Fun(ones.(spaces(A)))

## EuclideanSpace

const EuclideanSpace{RR} = VectorSpace{ConstantSpace{AnyDomain},AnyDomain,RR}
EuclideanSpace(n::Integer) = ArraySpace(ConstantSpace(Float64),n)

## support pieces

npieces(f::Fun{<:ArraySpace}) = npieces(f[1])
piece(f::Fun{<:ArraySpace}, k) = Fun(piece.(Array(f),k))
pieces(f::Fun{<:ArraySpace}) = [piece(f,k) for k=1:npieces(f)]

## TODO: This is a hack to get tests working

fromcanonical(d::ProductDomain, f::Fun{<:ArraySpace}) = vcat(fromcanonical.(factors(d), vec(f))...)

function coefficients(f::AbstractVector,sp::ArraySpace{<:ConstantSpace{AnyDomain}},ts::TensorSpace{SV,D,R}) where {SV,D<:Domain2d,R}
    @assert length(ts.spaces) == 2

    if ts.spaces[1] isa ArraySpace
        coefficients(f, sp, ts.spaces[1])
    elseif ts.spaces[2] isa ArraySpace
        coefficients(f, sp, ts.spaces[2])
    else
        error("Cannot convert coefficients from $sp to $ts")
    end
end

ArraySpace(sp::TensorSpace{Tuple{S1,S2}}) where {S1<:Space{D,R},S2} where {D,R<:AbstractArray} =
    ArraySpace(map(a -> a ⊗ sp.spaces[2], sp.spaces[1]))

ArraySpace(sp::TensorSpace{Tuple{S1,S2}},k...) where {S1,S2<:Space{D,R}} where {D,R<:AbstractArray} =
    ArraySpace(map(a -> sp.spaces[1] ⊗ a, sp.spaces[2]))

function coefficients(f::AbstractVector, a::VectorSpace, b::TensorSpace{Tuple{S1,S2},<:Domain2d}) where {S1<:Space{D,R},S2} where {D,R<:AbstractArray}
    if size(a) ≠ size(b)
        throw(DimensionMismatch("dimensions must match"))
    end
    interlace(map(coefficients,Fun(a,f),b),ArraySpace(b))
end

function coefficients(f::AbstractVector, a::VectorSpace, b::TensorSpace{Tuple{S1,S2},<:Domain2d}) where {S1,S2<:Space{D,R}} where {D,R<:AbstractArray}
    if size(a) ≠ size(b)
        throw(DimensionMismatch("dimensions must match"))
    end
    interlace(map(coefficients,Fun(a,f),b),ArraySpace(b))
end

ext: jl


text:
# ensure the requisite libraries are available

const noNemo = "Nemo.jl is not found:\n  Pkg.rm(\"Nemo\"); Pkg.add(\"Nemo\"); quit()\n  Pkg.rm(\"ArbFloats\");Pkg.add(\"ArbFloats\");"
const reNemo = "Nemo.jl is not as expected:\n  Pkg.rm(\"Nemo\"); Pkg.add(\"Nemo\"); quit()\n  Pkg.rm(\"ArbFloats\");Pkg.add(\"ArbFloats\");"

function package_directory(pkgName::String)
    pkgdir = Base.find_package(pkgName)
    nothing == pkgdir && throw(ErrorException(noNemo))
    return abspath(joinpath(split(pkgdir, pkgName)[1], pkgName))
end

function library_filepath(libsdir::String, filenames::Vector{String}, libname::String)
    libfile = filenames[ findfirst([startswith(x,libname) for x in filenames]) ]
    return joinpath( libsdir, libfile )
end

NemoLibsDir = abspath(joinpath( package_directory("Nemo"), "local/lib"))
libFiles = readdir(NemoLibsDir)

libarb   = library_filepath( NemoLibsDir, libFiles, "libarb"  )
libflint = library_filepath( NemoLibsDir, libFiles, "libflint")

isfile(libarb)   || throw(ErrorException(reNemo))
isfile(libflint) || throw(ErrorException(reNemo))

# prepare the libraries for use

@static if (Sys.islinux() || Sys.isbsd())
    libarb = String(split(libarb,".so")[1])
    libflint = String(split(libflint,".so")[1])
end
@static if Sys.isapple()
    libarb = String(split(libarb,".dynlib")[1])
    libflint = String(split(libflint,".dynlib")[1])
end
@static if Sys.iswindows()
    libarb = String(split(libarb,".dll")[1])
    libflint = String(split(libflint,".dll")[1])
end

macro libarb(sym)
    (:($sym), libarb)
end

macro libflint(sym)
    (:($sym), libflint)
end


text:
import Primes

function CHD{K, V}(kv;
                   max_collisions::Integer=16777216,
                   gamma::Real=1.25,
                   check_duplicates::Bool=true) where {K, V}
  !(gamma >= 1.0) && throw(ArgumentError("gamma needs to be no less than 1"))
  count = length(kv)
  n ::UInt64 = max(17, UInt64(ceil(count * gamma)))
  n < count && throw(OverflowError("count * gamma is bigger than typemax(UInt64)"))
  #@info "Finding nextprime($n)"
  n = Primes.nextprime(n)
  #@info "The prime is $n"
  m ::UInt64 = Primes.nextprime(n ÷ 2)

  slots = zeros(UInt8, n)
  keys = Vector{K}(undef, n)
  vals = Vector{V}(undef, n)
  hasher = ChdHasher(n, m)
  buckets = [Bucket{K, V}(0, [], []) for i in 1:m]
  indices = [~UInt16(0) for i in 1:m]

  seen = Set{UInt64}()

  if check_duplicates
    duplicates = Set{K}()
    for (key, val) in kv
      (key in duplicates) && throw(ArgumentError("duplicate key $(string(key))"))
      push!(duplicates, key)
    end
  end

  for (key, val) in kv
    oh = hash_index_from_key(hasher, key)
    @inbounds bucket = buckets[oh+1]
    bucket.index = oh
    push!(bucket.keys, key)
    push!(bucket.vals, val)
  end

  filter!((arg::Bucket{K,V}) -> !isempty(arg.keys), buckets)

  collisions = 0
  sort!(buckets;
        lt=(lhs::Bucket{K,V}, rhs::Bucket{K,V}) -> length(lhs.keys) < length(rhs.keys))

  # resolve each bucket
  i = 0
  while !isempty(buckets)
    i += 1
    bucket = pop!(buckets)
    next_bucket = false

    for (ri, r) in enumerate(hasher.r)
      # works with the existing r
      if try_hash(hasher, seen, slots, keys, vals, indices, bucket, UInt16(ri-1), r)
        next_bucket = true
        break
      end
    end
    next_bucket && continue

    # look for a new r
    for j in 1:max_collisions
      collisions = max(collisions, j)
      (ri, r) = generate(hasher)

      if try_hash(hasher, seen, slots, keys, vals, indices, bucket, ri, r)
        push!(hasher.r, r)
        next_bucket = true
        break
      end
    end
    next_bucket && continue

    error("Failed to find a collision-free hash function after $(max_collisions) attempts, " *
          "for bucket $i/$(length(buckets)) with $(length(bucket.keys)) entries." *
          "max bucket collisions: $(collisions), " *
          "keys: $count, " *
          "hash functions: $(length(hasher.r))")
  end

  return CHD{K, V}(slots, keys, vals, count, hasher.r, indices)
end

function CHD{K, V}(p::Pair; kwargs...) ::CHD{K, V} where {K, V}
  return CHD{K, V}([(p.first, p.second)])
end

function CHD{K, V}(ps::Pair...; kwargs...) ::CHD{K, V} where V where K
  return CHD{K, V}(ps; kwargs...)
end

function CHD(kv; kwargs...)
  try
    Base.dict_with_eltype((K, V) -> ((x...) -> CHD{K, V}(x...; kwargs...) ), kv, eltype(kv))
  catch
    if !Base.isiterable(typeof(kv)) || !all(x->isa(x,Union{Tuple,Pair}),kv)
      throw(ArgumentError("Dict(kv): kv needs to be an iterator of tuples or pairs"))
    else
      rethrow()
    end
  end
end

CHD(; kwargs...) = CHD{Any, Any}(; kwargs...)
CHD(kv::Tuple{}; kwargs...) = CHD(; kwargs...)
copy(d::CHD{K,V}) where {K, V} = CHD{K, V}(copy(d.slots), Base.copy(d.keys), copy(d.vals), d.count, copy(d.r), copy(d.indices))
CHD(ps::Pair{K, V}...; kwargs...) where {K, V} = CHD{K, V}(ps; kwargs...)
CHD(ps::Pair...; kwargs...)                    = CHD(ps; kwargs...)

max_stars_repo_name: kyungminlee/MinimalPerfectHashing

max_issues_repo_name: kyungminlee/MinimalPerfectHashing

max_forks_repo_name: kyungminlee/MinimalPerfectHashing


text:
######################################################################
# dynamics.jl
# Each trajectory manager will have a type of dynamics
######################################################################
export Dynamics, LinearDynamics, DubinsDynamics, linearize
export IntegrationScheme, ForwardEuler, SymplecticEuler
export forward_euler
export GroupDynamics
export vvf2vvvf

type LinearDynamics <: Dynamics
    n::Int
    m::Int
    A::Matrix{Float64}
    B::Matrix{Float64}

    function LinearDynamics(A, B)
        n,m = size(B)
        return new(n,m,deepcopy(A),deepcopy(B))
    end
end

type GroupDynamics <: Dynamics
	n::Int
	m::Int
	num_agents::Int
	array

	function GroupDynamics(array)
		num_agents = length(array)
		n = m = 0
		for j = 1:num_agents
			n += array[j].n
			m += array[j].m
		end
		return new(n, m, num_agents, array)
	end
end

type DubinsDynamics <: Dynamics
	n::Int
	m::Int

	v0::Float64		# constant speed
	r::Float64		# minimum turn radius

	# Constructors
	DubinsDynamics() = DubinsDynamics(1.0, 1.0)
	DubinsDynamics(v0::Real, r::Real) = new(3, 1, v0, r)
end

# Not really types,
export SingleIntegrator, DoubleIntegrator
function SingleIntegrator(n::Int, h::Float64)
	A = eye(n)
	B = h*eye(n)
	return LinearDynamics(A, B)
end
function DoubleIntegrator(n::Int, h::Float64)
	A = eye(2*n)
	for i = 1:n
		A[i, n+i] = h
	end

	B = zeros(2*n, n)
	for i = 1:n
		B[n+i, i] = h
	end

	return LinearDynamics(A, B)
end

"""
`dynamics!(tm::TrajectoryManager, d::Dynamics)`

Not only sets `tm.dynamics = d`, but also sets reward matrices `tm.Qn`, `tm.R`, and `tm.Rn` to default matrices of the correct sizes. These defaults are:

`tm.Qn = eye(d.n)`

`tm.R = 0.01 * eye(d.m)`

`tm.Rn = eye(d.m)`

"""
function dynamics!(tm::TrajectoryManager, d::Dynamics)
	tm.dynamics = d
	tm.Qn = eye(d.n)
	tm.R = 0.01 * eye(d.m)
	tm.Rn = eye(d.m)
end

######################################################################
# linearization
######################################################################

# linearizes about a trajectory
function linearize(d::Dynamics, x::VVF, u::VVF, h::Float64)
	N = length(u)
	A = VMF()
	B = VMF()
	for n = 1:N
		An, Bn = linearize(d, x[n], u[n], h)
		push!(A, An)
		push!(B, Bn)
	end
	return A,B
end

#function linearize(vtm::Vector{TrajectoryManager}, xds::VVVF, uds::VVVF)
#	num_agents = length(vtm)
#	N = length(xds[1]) - 1
#
#	# Make A and B a matrix with N rows and num_agents cols
#	A, B = linearize(vmt[1].dynamics, xds[1], uds[1], vmt[1].h)
#	for j = 2:num_agents
#		tm = vtm[j]
#		At, Bt = linearize(vtm[j].dynamics, xds[j], uds[j], vtm[j].h)
#		A = hcat(A, At)
#		B = hcat(B, Bt)
#	end
#
#	# now concatenate along each of the rows
#	dims = ones(Int, num_agents)
#	dims[1] = 1
#	for n = 1:N
#		push!(bigA, cat(dims, A[n]...))
#		push!(bigB, cat(dims, B[n]...))
#	end
#
#	return bigA, bigB
#end

function linearize(ld::LinearDynamics, x::VF, u::VF, h::Float64)
	return ld.A, ld.B
end

function linearize(ld::DubinsDynamics, x::VF, u::VF, h::Float64)
	A = eye(3)
	A[1,3] = -h * sin(x[3]) * ld.v0
	A[2,3] = h * cos(x[3]) * ld.v0

	B = zeros(3,1)
	B[3] = h/ld.r

	return A, B
end

function gsplit(gd::GroupDynamics, x::VF)
	xarr = VVF()
	xind = 1
	for j = 1:gd.num_agents
		d = gd.array[j]
		push!(xarr, x[xind + d.n - 1])
		xind += d.n
	end
	return xarr
end

function gsplit(gd::GroupDynamics, x::VF, u::VF)
	xarr = VVF()
	uarr = VVF()
	xind = uind = 1
	for j = 1:gd.num_agents
		d = gd.array[j]
		push!(xarr, x[xind:(xind + d.n - 1)])
		push!(uarr, u[uind:(uind + d.m - 1)])
		xind += d.n
		uind += d.m
	end
	return xarr, uarr
end

function linearize(gd::GroupDynamics, x::VF, u::VF, h::Float64)

	xarr, uarr = gsplit(gd, x, u)

	Aarr = VMF()
	Barr = VMF()
	for j = 1:gd.num_agents
		A, B = linearize(gd.array[j], xarr[j], uarr[j], h)
		push!(Aarr, A)
		push!(Barr, B)
	end

	dims = 2*ones(Int, gd.num_agents)
	dims[1] = 1

	return cat(dims, Aarr...), cat(dims, Barr...)
end

######################################################################
# integration
######################################################################
type ForwardEuler <: IntegrationScheme end
type SymplecticEuler <: IntegrationScheme end

function integrate(tm::TrajectoryManager, x::VF, u::VF)
	integrate(tm.int_scheme, tm.dynamics, x, u, tm.h)
end

function integrate(::ForwardEuler, d::Dynamics, x::VF, u::VF, h::Float64)
	forward_euler(d, x, u, h)
end

function integrate(::SymplecticEuler, d::Dynamics, x::VF, u::VF, h::Float64)
	symplectic_euler(d, x, u, h)
end

function integrate(tm::TrajectoryManager, ud::VVF)
	xd = Array(Vector{Float64}, tm.N+1)

	xd[1] = deepcopy(tm.x0)
	for i = 1:tm.N
		xd[i+1] = integrate(tm.int_scheme, tm.dynamics, x, ud[i], tm.h)
	end

	return xd
end

######################################################################
# forward_euler
######################################################################
function forward_euler(tm::TrajectoryManager, x::VF, u::VF)
	forward_euler(tm.dynamics, x, u, tm.h)
end

function forward_euler(ld::LinearDynamics, x::VF, u::VF, h::Float64)
	return ld.A*x + ld.B*u
end
function forward_euler(gd::GroupDynamics, x::VF, u::VF, h::Float64)
	n_ind = m_ind = 1
	xp = VVF()
	xarr, uarr = gsplit(gd, x, u)
	for j = 1:gd.num_agents
		push!(xp, forward_euler(gd.array[j], xarr[j], uarr[j], h))
	end
	return vcat(xp...)
end

function forward_euler(dd::DubinsDynamics, x::VF, u::VF, h::Float64)
	xp = deepcopy(x)
	xp[1] += cos(x[3]) * dd.v0 * h
	xp[2] += sin(x[3]) * dd.v0 * h
	u_val = u[1]
	if u_val > 1.
		u_val = 1.
	end
	if u_val < -1.
		u_val = -1.
	end
	xp[3] += u_val / dd.r * h
	return xp
end

######################################################################
# symplectic_euler
######################################################################
function symplectic_euler(d::Dynamics, x::VF, u::VF, h::Float64)
	#A,B = linearize(d, x, u, h)
	A = zeros(4,4)
	A[1,3] = 1
	A[2,4] = 1
	A[3,1] = -1
	A[4,2] = -1
	B = zeros(4,2)
	B[3,1] = 1
	B[4,2] = 1

	#n2 = round(Int, d.n/2)
	n2 = 2

	A11 = A[1:n2, 1:n2]
	A12 = A[1:n2, n2+1:d.n]
	A21 = A[n2+1:d.n, 1:n2]
	A22 = A[n2+1:d.n, n2+1:d.n]

	B1 = B[1:n2,:]
	B2 = B[n2+1:d.n,:]

	# block matrices
	Ad22 = inv(eye(n2)-h*A22)
	Ad11 = eye(n2) + h*A11 + h*h*A12*Ad22*A21
	Ad12 = h*A12*Ad22
	Ad21 = h*Ad22*A21

	Bd1 = h*h*A12*Ad22*B2 + h*B1
	Bd2 = h*Ad22*B2

	Ase = [Ad11 Ad12; Ad21 Ad22]
	Bse = [Bd1; Bd2]

	xp = Ase*x + Bse*u
	return xp
end

# split one vector per state up
function vvf2vvvf(xd::VVF, ud::VVF, vtm::Vector{TrajectoryManager})
	N = length(xd) - 1
	num_agents = length(vtm)

	x_start = u_start = 1
	xds = VVVF(num_agents)
	uds = VVVF(num_agents)
	for j = 1:num_agents
		x_end = x_start + vtm[j].dynamics.n - 1
		u_end = u_start + vtm[j].dynamics.m - 1
		xds[j] = VVF()
		uds[j] = VVF()
		for n = 1:N
			push!(xds[j], xd[n][x_start:x_end])
			push!(uds[j], ud[n][u_start:u_end])
		end
		push!(xds[j], xd[N+1][x_start:x_end])
		x_start = x_end + 1
		u_start = u_end + 1
	end
	return xds, uds
end

function vvf2vvvf(xd::VVF, vtm::Vector{TrajectoryManager})
	N = length(xd) - 1
	num_agents = length(vtm)

	x_start = 1
	xds = VVVF(num_agents)
	for j = 1:num_agents
		x_end = x_start + vtm[j].dynamics.n - 1
		xds[j] = VVF()
		for n = 1:N+1
			push!(xds[j], xd[n][x_start:x_end])
		end
		x_start = x_end + 1
	end
	return xds
end

function vvf2vvvf(xd::VVF, gd::GroupDynamics)
	N = length(xd) - 1

	x_start = 1
	xds = VVVF(gd.num_agents)
	for j = 1:gd.num_agents
		x_end = x_start + gd.array[j].n - 1
		xds[j] = VVF()
		for n = 1:(N+1)
			push!(xds[j], xd[n][x_start:x_end])
		end
		x_start = x_end + 1
	end
	return xds
end


text:
include("benchtools.jl")
module ITensorsTimers
    using ITensors
    import ..Timers

    function add_indexlabels_as_tags(tensor, labels)
        new_inds = map(x -> settags(x[2], "l$(x[1])"), zip(labels, tensor.inds))
        return ITensor(tensor.store, new_inds)
    end

    function mpo_timer(T = Float64; Vmpo, Vmps, Vphys)
        A = randomITensor(T, Vmps, Vphys, dag(Vmps))
        M = randomITensor(T, Vmpo, Vphys, dag(Vphys), dag(Vmpo))
        FL = randomITensor(T, Vmps, dag(Vmpo), dag(Vmps))
        FR = randomITensor(T, Vmps, Vmpo, dag(Vmps))

        FL1 = add_indexlabels_as_tags(FL, (4, 2, 1))
        A1 = add_indexlabels_as_tags(A, (1, 3, 6))
        M1 = add_indexlabels_as_tags(M, (2, 5, 3, 7))
        A2 = add_indexlabels_as_tags(A, (4, 5, 8))
        FR1 = add_indexlabels_as_tags(FR, (6, 7, 8))

        return Timers.Timer(A1, A2, M1, FL1, FR1) do A1, A2, M1, FL1, FR1
            C = ((((FL1*A1)*M1)*dag(A2))*FR1)
            return C[]
        end
    end

    function pepo_timer(T = Float64; Vpepo, Vpeps, Venv, Vphys)
        Vpeps′ = dag(Vpeps)
        Vpepo′ = dag(Vpepo)
        Venv′ = dag(Venv)
        Vphys′ = dag(Vphys)

        A = randomITensor(T, Vpeps, Vpeps, Vphys, Vpeps′, Vpeps′)
        P = randomITensor(T, Vpepo, Vpepo, Vphys, Vphys′, Vpepo′, Vpepo′)
        FL = randomITensor(T, Venv, Vpeps, Vpepo′, Vpeps′, Venv′)
        FD = randomITensor(T, Venv, Vpeps, Vpepo′, Vpeps′, Venv′)
        FR = randomITensor(T, Venv, Vpeps, Vpepo, Vpeps′, Venv′)
        FU = randomITensor(T, Venv, Vpeps, Vpepo, Vpeps′, Venv′)

        FL1 = add_indexlabels_as_tags(FL, (18, 7, 4, 2, 1))
        FU1 = add_indexlabels_as_tags(FU, (1, 3, 6, 9, 10))
        A1 = add_indexlabels_as_tags(A, (2, 17, 5, 3, 11))
        P1 = add_indexlabels_as_tags(P, (4, 16, 8, 5, 6, 12))
        A2 = add_indexlabels_as_tags(A, (7, 15, 8, 9, 13))
        FR1 = add_indexlabels_as_tags(FR, (10, 11, 12, 13, 14))
        FD1 = add_indexlabels_as_tags(FD, (14, 15, 16, 17, 18))

        return Timers.Timer(A1, A2, P1, FL1, FD1, FR1, FU1) do A1, A2, P1, FL1, FD1, FR1, FU1
            C = (((((FL1*FU1)*A1)*P1)*dag(A2))*FR1)*FD1
            return C[]
        end
    end

    function mera_timer(T = Float64; Vmera)
        Vmera′ = dag(Vmera);

        u = randomITensor(T, Vmera, Vmera, Vmera′, Vmera′)
        w = randomITensor(T, Vmera, Vmera, Vmera′)
        ρ = randomITensor(T, Vmera, Vmera, Vmera, Vmera′, Vmera′, Vmera′)
        h = randomITensor(T, Vmera, Vmera, Vmera, Vmera′, Vmera′, Vmera′)

        h1 = add_indexlabels_as_tags(h, (9, 3, 4, 5, 1, 2))
        u1 = add_indexlabels_as_tags(u, (1, 2, 7, 12))
        u2 = add_indexlabels_as_tags(u, (3, 4, 11, 13))
        u3 = add_indexlabels_as_tags(u, (8, 5, 15, 6))
        w1 = add_indexlabels_as_tags(w, (6, 7, 19))
        u4 = add_indexlabels_as_tags(u, (8, 9, 17, 10))
        w2 = add_indexlabels_as_tags(w, (10, 11, 22))
        w3 = add_indexlabels_as_tags(w, (12, 14, 20))
        w4 = add_indexlabels_as_tags(w, (13, 14, 23))
        ρ1 = add_indexlabels_as_tags(ρ, (18, 19, 20, 21, 22, 23))
        w5 = add_indexlabels_as_tags(w, (16, 15, 18))
        w6 = add_indexlabels_as_tags(w, (16, 17, 21))

        return Timers.Timer(u1, u2, u3, u4, w1, w2, w3, w4, w5, w6, ρ1, h1) do u1, u2, u3, u4, w1, w2, w3, w4, w5, w6, ρ1, h1
            C = (((((( (h1*u1) * dag(u2) ) *
                                (u3*w1) ) *
                                (dag(u4)*dag(w2)) ) *
                               ( (w3 * dag(w4)) * ρ1)) * w5) * dag(w6))
            return C
        end
    end
end


text:
import Pkg
Pkg.update()

Pkg.add("JuMP")
Pkg.add("ArgParse")
Pkg.add("Nullables")

Pkg.add("Ipopt")
#Pkg.add("Requests")
Pkg.add("CSV")
Pkg.add("YAML")
Pkg.add("JSON")

Pkg.add("DocOpt")
Pkg.add("DataFrames")

Pkg.add("CSVFiles")

max_stars_repo_name: OICR/mp-biopath

max_issues_repo_name: OICR/mp-biopath

max_forks_repo_name: OICR/mp-biopath


text:
# Copyright 2018 Erik Edin
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

"Thrown at @expect failures in step definitions."
struct StepAssertFailure <: Exception
    assertion::String
    evaluated::String

    StepAssertFailure(assertion::String) = new(assertion, "")
    StepAssertFailure(assertion::String, evaluated::String) = new(assertion, evaluated)
end

"""
    expect(ex)

Assert that a condition `ex` is true. Throws a StepAssertFailure if false.

Warning: The expressions in the `ex` condition may be evaluated more than once,
         so only use expression without side effects in the condition.

# Examples
```
@then "this condition should hold" begin
    @expect 1 == 1
end
```

This will fail and throw a `StepAssertFailure`
```
@then "this condition should not hold" begin
    @expect 1 == 2
end
```
"""
macro expect(ex)
    comparisonops = [:(==), :(===), :(!=), :(!==)]

    if hasproperty(ex, :head)
        iscomparison = ex.head === :call && length(ex.args) == 3 && ex.args[1] in comparisonops
    else
        iscomparison = false
    end

    if iscomparison
        comparisonop = QuoteNode(ex.args[1])
        quote
            # Evaluate each argument separately
            local value1 = $(esc(ex.args[2]))
            local value2 = $(esc(ex.args[3]))

            # Create a en expression to show
            local showexpr = Expr(:call, $comparisonop, value1, value2)
            local msg = $(string(ex))
            local evaluated = string(showexpr)

            # Check the expectation by evaluating the original expression
            if !($(esc(ex)))
                throw(StepAssertFailure(msg, evaluated))
            end
        end
    else
        quote
            local msg = $(string(ex))
            # Check the expectation.
            if !(Bool($(esc(ex))))
                throw(StepAssertFailure(msg))
            end
        end
    end
end

"""
    @fail(message)

Fail a step with a descriptive message.
"""
macro fail(message)
    quote
        throw(StepAssertFailure($message))
    end
end


text:
## Multivariate ##

"""
    TimeSeriesModelStorage

Preallocated storage for model manipulation including preplanned ffts.
"""
abstract type TimeSeriesModelStorage end

struct Acv2EIStorage{N,M,D,L,P} <: TimeSeriesModelStorage
    allocatedarray::Array{ComplexF64,M}
    hermitianarray::Base.ReinterpretArray{
        SHermitianCompact{D,ComplexF64,L},
        N,
        ComplexF64,
        Array{ComplexF64,N},
    }
    planned_fft::P
    kernel::Vector{Float64}
    function Acv2EIStorage(
        allocatedarray::Array{ComplexF64,M},
        hermitianarray::Base.ReinterpretArray{
            SHermitianCompact{D,ComplexF64,L},
            N,
            ComplexF64,
            Array{ComplexF64,N},
        },
        kernel,
        fft_flags,
    ) where {N,M,D,L}
        planned_fft = plan_fft!(allocatedarray, M, flags = fft_flags)
        new{N,M,D,L,typeof(planned_fft)}(allocatedarray, hermitianarray, planned_fft, kernel)
    end
end

struct Sdf2AcvStorage{M,Q} <: TimeSeriesModelStorage
    allocatedarray::Array{ComplexF64,M}
    planned_ifft::Q
    function Sdf2AcvStorage(allocatedarray::Array{ComplexF64,M}, fft_flags) where {M}
        planned_ifft = plan_ifft!(allocatedarray, M, flags = fft_flags)
        new{M,typeof(planned_ifft)}(allocatedarray, planned_ifft)
    end
end

# Return the type of the fftplan from Acv2EIStorage and Sdf2AcvStorage types.
fft_plan_type(::Acv2EIStorage{N,M,D,L,P}) where {N,M,D,L,P} = P
fft_plan_type(::Sdf2AcvStorage{M,Q}) where {M,Q} = Q

# Struct to store required storage for computing EI when acv is not specified.
struct Sdf2EIStorage{N,M,D,L,Q,P} <: TimeSeriesModelStorage
    sdf2acv::Sdf2AcvStorage{M,Q}
    acv2EI::Acv2EIStorage{N,M,D,L,P}
    function Sdf2EIStorage(
        sdf_array::Array{ComplexF64,M},
        acv_array::Array{ComplexF64,M},
        hermitianarray::Base.ReinterpretArray{
            SHermitianCompact{D,ComplexF64,L},
            N,
            ComplexF64,
            Array{ComplexF64,N},
        },
        kernel,
        fft_flags,
    ) where {N,M,D,L}
        sdf2acv = Sdf2AcvStorage(sdf_array, fft_flags)
        acv2EI = Acv2EIStorage(acv_array, hermitianarray, kernel, fft_flags)
        new{N,M,D,L,fft_plan_type(sdf2acv),fft_plan_type(acv2EI)}(sdf2acv, acv2EI)
    end
end

# Generate required storage for computing the EI for a general time series model.
function EIstorage_function(model::Type{<:TimeSeriesModel{D,T}}, n, kernel, fft_flags = FFTW.MEASURE) where {D,T}
    acv_array = zeros(ComplexF64, nlowertriangle_dimension(model), 2n)
    hermitianarray = reinterpret(SHermitianCompact{D,eltype(acv_array),nlowertriangle_dimension(model)}, vec(acv_array))
    return Acv2EIStorage(acv_array, hermitianarray, kernel, fft_flags)
end
function EIstorage_function(model::Type{<:UnknownAcvTimeSeriesModel{D,T}}, n, kernel, fft_flags = FFTW.MEASURE) where {D,T}
    sdf_array = zeros(ComplexF64, nlowertriangle_dimension(model), max(2n, minbins(model)))
    acv_array = zeros(ComplexF64, nlowertriangle_dimension(model), 2n)
    hermitianarray = reinterpret(SHermitianCompact{D,eltype(acv_array),nlowertriangle_dimension(model)}, vec(acv_array))
    return Sdf2EIStorage(sdf_array, acv_array, hermitianarray, kernel, fft_flags)
end

# Generate required storage for computing the gradient of EI for a general time series model.
function EIstorage_gradient(model::Type{<:TimeSeriesModel{D,T}}, n, kernel, fft_flags = FFTW.MEASURE) where {D,T}
    acv_array = zeros(ComplexF64, nlowertriangle_dimension(model), npars(model), 2n)
    hermitianarray = reinterpret(SHermitianCompact{D,eltype(acv_array),nlowertriangle_dimension(model)}, reshape(acv_array, (:, 2n)))
    return Acv2EIStorage(acv_array, hermitianarray, kernel, fft_flags)
end
function EIstorage_gradient(model::Type{<:UnknownAcvTimeSeriesModel{D,T}}, n, kernel, fft_flags = FFTW.MEASURE) where {D,T}
    sdf_array = zeros(ComplexF64, nlowertriangle_dimension(model), npars(model), max(2n, minbins(model)))
    acv_array = zeros(ComplexF64, nlowertriangle_dimension(model), npars(model), 2n)
    hermitianarray = reinterpret(SHermitianCompact{D,eltype(acv_array),nlowertriangle_dimension(model)}, reshape(acv_array, (:, 2n)))
    return Sdf2EIStorage(sdf_array, acv_array, hermitianarray, kernel, fft_flags)
end

# Generate required storage for computing the hessian of EI for a general time series model.
function EIstorage_hessian(model::Type{<:TimeSeriesModel{D,T}}, n, kernel, fft_flags = FFTW.MEASURE) where {D,T}
    acv_array = zeros(ComplexF64, nlowertriangle_dimension(model), nlowertriangle_parameter(model), 2n)
    hermitianarray = reinterpret(SHermitianCompact{D,eltype(acv_array),nlowertriangle_dimension(model)}, reshape(acv_array, (:, 2n)))
    return Acv2EIStorage(acv_array, hermitianarray, kernel, fft_flags)
end
function EIstorage_hessian(model::Type{<:UnknownAcvTimeSeriesModel{D,T}}, n, kernel, fft_flags = FFTW.MEASURE) where {D,T}
    sdf_array = zeros(ComplexF64, nlowertriangle_dimension(model), nlowertriangle_parameter(model), max(2n, minbins(model)))
    acv_array = zeros(ComplexF64, nlowertriangle_dimension(model), nlowertriangle_parameter(model), 2n)
    hermitianarray = reinterpret(SHermitianCompact{D,eltype(acv_array),nlowertriangle_dimension(model)}, reshape(acv_array, (:, 2n)))
    return Sdf2EIStorage(sdf_array, acv_array, hermitianarray, kernel, fft_flags)
end

#  Storage for useful quanties associated with the EI memory when estimating the acv, namely n, Δ and Ωₘ.
struct FreqAcvEst
    n::Int
    Δ::Float64
    Ωₘ::Frequencies{Float64}
    function FreqAcvEst(model::Type{<:UnknownAcvTimeSeriesModel},n,Δ)
        Ωₘ = fftfreq(max(2n, minbins(model)), 2π/Δ)
        new(n,Δ,Ωₘ)
    end
    function FreqAcvEst(n,Δ)
        Ωₘ = fftfreq(n, 2π/Δ)
        new(n,Δ,Ωₘ)
    end
end

# Storage for useful quanties associated with the EI memory when the acv is known, namely n, Δ and Ωₘ.
struct LagsEI
    n::Int
    Δ::Float64
    lags::Frequencies{Float64}
    function LagsEI(n,Δ)
        lags = fftfreq(2n, 2n*Δ)
        new(n,Δ,lags)
    end
end

# Function to choose additional storage depending on whether the acv is known or not.
encodetimescale(model::Type{<:TimeSeriesModel},n,Δ) = LagsEI(n,Δ)
encodetimescale(model::Type{<:UnknownAcvTimeSeriesModel},n,Δ) = FreqAcvEst(model,n,Δ)

abstract type TimeSeriesModelStorageFunction <: TimeSeriesModelStorage end
abstract type TimeSeriesModelStorageGradient <: TimeSeriesModelStorageFunction end
abstract type TimeSeriesModelStorageHessian  <: TimeSeriesModelStorageGradient end

# Structure to store memory if only the EI is required.
struct PreallocatedEI{T,F} <: TimeSeriesModelStorageFunction
    encodedtime::T
    funcmemory::F
    function PreallocatedEI(model::Type{<:TimeSeriesModel}, n, Δ, kernel; fft_flags = FFTW.MEASURE)
        time = encodetimescale(model, n, Δ)
        func = EIstorage_function(model, n, kernel, fft_flags)
        new{typeof(time),typeof(func)}(time,func)
    end
end

# Structure to store memory if the EI and its gradient is required.
struct PreallocatedEIGradient{T,F,G} <: TimeSeriesModelStorageGradient
    encodedtime::T
    funcmemory::F
    gradmemory::G
    function PreallocatedEIGradient(model::Type{<:TimeSeriesModel}, n, Δ, kernel; fft_flags = FFTW.MEASURE)
        time = encodetimescale(model, n, Δ)
        func = EIstorage_function(model, n, kernel, fft_flags)
        grad = EIstorage_gradient(model, n, kernel, fft_flags)
        new{typeof(time),typeof(func),typeof(grad)}(time, func, grad)
    end
end

# Structure to store memory if the EI, its gradient and its hessian are required.
struct PreallocatedEIHessian{T,F,G,H} <: TimeSeriesModelStorageHessian
    encodedtime::T
    funcmemory::F
    gradmemory::G
    hessmemory::H
    function PreallocatedEIHessian(model::Type{<:TimeSeriesModel}, n, Δ, kernel; fft_flags = FFTW.MEASURE)
        time = encodetimescale(model, n, Δ)
        func = EIstorage_function(model, n, kernel, fft_flags)
        grad = EIstorage_gradient(model, n, kernel, fft_flags)
        hess = EIstorage_hessian(model, n, kernel, fft_flags)
        new{typeof(time),typeof(func),typeof(grad),typeof(hess)}(time, func, grad, hess)
    end
end

# Storage for additive models.
struct AdditiveStorage{S₁<:TimeSeriesModelStorage,S₂<:TimeSeriesModelStorage} <: TimeSeriesModelStorage
    store1::S₁
    store2::S₂
    npar1::Int
end

# Helper function for extracting the function memory from additive storage (which is only stored in the left most storage).
extract_S(store::TimeSeriesModelStorageFunction) = store.funcmemory
extract_S(store::AdditiveStorage) = extract_S(store.store1) # returns the left hand store if additive storage

# tapering kernel computation
function taper2timedomainkernel(taper::AbstractVector{Float64},n)
    length(taper) == n || throw(ArgumentError("Taper is not of length n."))
    kernel = Vector{Float64}(undef,2n)
    for (i,τ) ∈ enumerate(0:n-1)
        kernel[i] = sum(taper[t]*taper[t+τ] for t ∈ 1:length(taper)-τ)
    end
    kernel[n+1] = 0
    for i in n+2:2n
        kernel[i] = kernel[2n-i+2]
    end
    return kernel
end
function taper2timedomainkernel(taper::Nothing,n)
    kernel = Vector{Float64}(undef,2n)
    for i ∈ 1:n
        kernel[i] = (1 - (i - 1) / n)
    end
    for i ∈ n+1:2n
        kernel[i] = (1 - (2n - i + 1) / n)
    end
    return kernel
end

# Interface level function for allocating the memory for EI computation.
function allocate_memory_EI_F(model::Type{<:TimeSeriesModel}, n, Δ; taper = nothing, fft_flags = FFTW.MEASURE)
    kernel = taper2timedomainkernel(taper,n)
    return _allocate_memory_EI_F(model::Type{<:TimeSeriesModel}, n, Δ, kernel; fft_flags = fft_flags)
end
function _allocate_memory_EI_F(model::Type{<:TimeSeriesModel}, n, Δ, kernel; fft_flags)
    return PreallocatedEI(model, n, Δ, kernel; fft_flags = fft_flags)
end
function _allocate_memory_EI_F(model::Type{AdditiveTimeSeriesModel{M₁,M₂,D,T}}, n, Δ, kernel; fft_flags) where {M₁,M₂,D,T}
    return AdditiveStorage( _allocate_memory_EI_F(M₁, n, Δ, kernel; fft_flags = fft_flags),
                            _allocate_memory_EI_F(M₂, n, Δ, kernel; fft_flags = fft_flags),
                            npars(M₁))
end

# Interface level function for allocating the memory for EI computation with gradient.
function allocate_memory_EI_FG(model::Type{<:TimeSeriesModel}, n, Δ; taper = nothing, fft_flags = FFTW.MEASURE)
    kernel = taper2timedomainkernel(taper,n)
    return _allocate_memory_EI_FG(model::Type{<:TimeSeriesModel}, n, Δ, kernel; fft_flags = fft_flags)
end
function _allocate_memory_EI_FG(model::Type{<:TimeSeriesModel}, n, Δ, kernel; fft_flags)
    return PreallocatedEIGradient(model, n, Δ, kernel; fft_flags = fft_flags)
end
function _allocate_memory_EI_FG(model::Type{AdditiveTimeSeriesModel{M₁,M₂,D,T}}, n, Δ, kernel; fft_flags) where {M₁,M₂,D,T}
    return AdditiveStorage( _allocate_memory_EI_FG(M₁, n, Δ, kernel; fft_flags = fft_flags),
                            _allocate_memory_EI_FG(M₂, n, Δ, kernel; fft_flags = fft_flags),
                            npars(M₁))
end

# Interface level function for allocating the memory for EI computation with gradient and hessian.
function allocate_memory_EI_FGH(model::Type{<:TimeSeriesModel}, n, Δ; taper = nothing, fft_flags = FFTW.MEASURE)
    kernel = taper2timedomainkernel(taper,n)
    return _allocate_memory_EI_FGH(model::Type{<:TimeSeriesModel}, n, Δ, kernel; fft_flags = fft_flags)
end
function _allocate_memory_EI_FGH(model::Type{<:TimeSeriesModel}, n, Δ, kernel; fft_flags)
    return PreallocatedEIHessian(model, n, Δ, kernel; fft_flags = fft_flags)
end
function _allocate_memory_EI_FGH(model::Type{AdditiveTimeSeriesModel{M₁,M₂,D,T}}, n, Δ, kernel; fft_flags) where {M₁,M₂,D,T}
    return AdditiveStorage( _allocate_memory_EI_FGH(M₁, n, Δ, kernel; fft_flags = fft_flags),
                            _allocate_memory_EI_FGH(M₂, n, Δ, kernel; fft_flags = fft_flags),
                            npars(M₁))
end

## sdf only storage

struct SdfStorage{N,M,D,L} <: TimeSeriesModelStorage
    allocatedarray::Array{ComplexF64,M}
    hermitianarray::Base.ReinterpretArray{
        SHermitianCompact{D,ComplexF64,L},
        N,
        ComplexF64,
        Array{ComplexF64,N},
    }
    function SdfStorage(allocatedarray::Array{ComplexF64,M}, hermitianarray::Base.ReinterpretArray{
            SHermitianCompact{D,ComplexF64,L},
            N,
            ComplexF64,
            Array{ComplexF64,N},
        }) where {N,M,D,L}
        new{N,M,D,L}(allocatedarray,hermitianarray)
    end
end
# Function to allocate storage for the function part of sdf calculations.
function sdfstorage_function(model::Type{<:TimeSeriesModel{D,T}}, n) where {D,T}
    allocatedarray = zeros(ComplexF64, nlowertriangle_dimension(model), n)
    hermitianarray = reinterpret(SHermitianCompact{D,eltype(allocatedarray),nlowertriangle_dimension(model)}, vec(allocatedarray))
    return SdfStorage(allocatedarray,hermitianarray)
end
# Function to allocate storage for the gradient part of sdf calculations.
function sdfstorage_gradient(model::Type{<:TimeSeriesModel{D,T}}, n) where {D,T}
    allocatedarray = zeros(ComplexF64, nlowertriangle_dimension(model), npars(model), n)
    hermitianarray = reinterpret(SHermitianCompact{D,eltype(allocatedarray),nlowertriangle_dimension(model)}, reshape(allocatedarray, (:, n)))
    return SdfStorage(allocatedarray,hermitianarray)
end
# Function to allocate storage for the hessian part of sdf calculations.
function sdfstorage_hessian(model::Type{<:TimeSeriesModel{D,T}}, n) where {D,T}
    allocatedarray = zeros(ComplexF64, nlowertriangle_dimension(model), nlowertriangle_parameter(model), n)
    hermitianarray = reinterpret(SHermitianCompact{D,eltype(allocatedarray),nlowertriangle_dimension(model)}, reshape(allocatedarray, (:, n)))
    return SdfStorage(allocatedarray,hermitianarray)
end

# Structure to store memory if the sdf and its gradient is required.
struct PreallocatedSdf{F} <: TimeSeriesModelStorageFunction
    encodedtime::FreqAcvEst
    funcmemory::F
    function PreallocatedSdf(model::Type{<:TimeSeriesModel}, n, Δ)
        time = FreqAcvEst(n, Δ)
        func = sdfstorage_function(model, n)
        new{typeof(func)}(time, func)
    end
end

# Structure to store memory if the sdf and its gradient is required.
struct PreallocatedSdfGradient{F,G} <: TimeSeriesModelStorageGradient
    encodedtime::FreqAcvEst
    funcmemory::F
    gradmemory::G
    function PreallocatedSdfGradient(model::Type{<:TimeSeriesModel}, n, Δ)
        time = FreqAcvEst(n, Δ)
        func = sdfstorage_function(model, n)
        grad = sdfstorage_gradient(model, n)
        new{typeof(func),typeof(grad)}(time, func, grad)
    end
end

# Structure to store memory if the sdf, its gradient and its hessian are required.
struct PreallocatedSdfHessian{F,G,H} <: TimeSeriesModelStorageHessian
    encodedtime::FreqAcvEst
    funcmemory::F
    gradmemory::G
    hessmemory::H
    function PreallocatedSdfHessian(model::Type{<:TimeSeriesModel}, n, Δ)
        time = FreqAcvEst(n, Δ)
        func = sdfstorage_function(model, n)
        grad = sdfstorage_gradient(model, n)
        hess = sdfstorage_hessian(model, n)
        new{typeof(func),typeof(grad),typeof(hess)}(time, func, grad, hess)
    end
end

# Interface level function for allocating the memory for EI computation.
function allocate_memory_sdf_F(model::Type{<:TimeSeriesModel}, n, Δ)
    return PreallocatedSdf(model, n, Δ)
end
function allocate_memory_sdf_F(model::Type{AdditiveTimeSeriesModel{M₁,M₂,D,T}}, n, Δ) where {M₁,M₂,D,T}
    return AdditiveStorage( allocate_memory_sdf_F(M₁, n, Δ),
                            allocate_memory_sdf_F(M₂, n, Δ),
                            npars(M₁))
end

# Interface level function for allocating the memory for sdf computation with gradient.
function allocate_memory_sdf_FG(model::Type{<:TimeSeriesModel}, n, Δ)
    return PreallocatedSdfGradient(model, n, Δ)
end
function allocate_memory_sdf_FG(model::Type{AdditiveTimeSeriesModel{M₁,M₂,D,T}}, n, Δ) where {M₁,M₂,D,T}
    return AdditiveStorage( allocate_memory_sdf_FG(M₁, n, Δ),
                            allocate_memory_sdf_FG(M₂, n, Δ),
                            npars(M₁))
end

# Interface level function for allocating the memory for sdf computation with gradient and hessian.
function allocate_memory_sdf_FGH(model::Type{<:TimeSeriesModel}, n, Δ)
    return PreallocatedSdfHessian(model, n, Δ)
end
function allocate_memory_sdf_FGH(model::Type{AdditiveTimeSeriesModel{M₁,M₂,D,T}}, n, Δ) where {M₁,M₂,D,T}
    return AdditiveStorage( allocate_memory_sdf_FGH(M₁, n, Δ),
                            allocate_memory_sdf_FGH(M₂, n, Δ),
                            npars(M₁))
end

## Univariate ##

struct Acv2EIStorageUni{M,P} <: TimeSeriesModelStorage
    allocatedarray::Array{ComplexF64,M}
    planned_fft::P
    kernel::Vector{Float64}
    function Acv2EIStorageUni(
        allocatedarray::Array{ComplexF64,M},
        kernel,
        fft_flags,
    ) where {M}
        planned_fft = plan_fft!(allocatedarray, M, flags = fft_flags)
        new{M,typeof(planned_fft)}(allocatedarray, planned_fft, kernel)
    end
end

struct Sdf2AcvStorageUni{M,Q} <: TimeSeriesModelStorage
    allocatedarray::Array{ComplexF64,M}
    planned_ifft::Q
    function Sdf2AcvStorageUni(allocatedarray::Array{ComplexF64,M}, fft_flags) where {M}
        planned_ifft = plan_ifft!(allocatedarray, M, flags = fft_flags)
        new{M,typeof(planned_ifft)}(allocatedarray, planned_ifft)
    end
end

# Return the type of the fftplan from Acv2EIStorageUni and Sdf2AcvStorageUni types.
fft_plan_type(::Acv2EIStorageUni{M,P}) where {M,P} = P
fft_plan_type(::Sdf2AcvStorageUni{M,Q}) where {M,Q} = Q

# Struct to store required storage for computing EI when acv is not specified.
struct Sdf2EIStorageUni{M,Q,P} <: TimeSeriesModelStorage
    sdf2acv::Sdf2AcvStorageUni{M,Q}
    acv2EI::Acv2EIStorageUni{M,P}
    function Sdf2EIStorageUni(
        sdf_array::Array{ComplexF64,M},
        acv_array::Array{ComplexF64,M},
        kernel,
        fft_flags,
    ) where {M}
        sdf2acv = Sdf2AcvStorageUni(sdf_array, fft_flags)
        acv2EI = Acv2EIStorageUni(acv_array, kernel, fft_flags)
        new{M,fft_plan_type(sdf2acv),fft_plan_type(acv2EI)}(sdf2acv, acv2EI)
    end
end

# Generate required storage for computing the EI for a univariate time series model.
function EIstorage_function(model::Type{<:TimeSeriesModel{1,T}}, n, kernel, fft_flags = FFTW.MEASURE) where {T}
    acv_array = zeros(ComplexF64, 2n)
    return Acv2EIStorageUni(acv_array, kernel, fft_flags)
end
function EIstorage_function(model::Type{<:UnknownAcvTimeSeriesModel{1,T}}, n, kernel, fft_flags = FFTW.MEASURE) where {T}
    sdf_array = zeros(ComplexF64, max(2n, minbins(model)))
    acv_array = zeros(ComplexF64, 2n)
    return Sdf2EIStorageUni(sdf_array, acv_array, kernel, fft_flags)
end

# Generate required storage for computing the gradient of EI for a univariate time series model.
function EIstorage_gradient(model::Type{<:TimeSeriesModel{1,T}}, n, kernel, fft_flags = FFTW.MEASURE) where {T}
    acv_array = zeros(ComplexF64, npars(model), 2n)
    return Acv2EIStorageUni(acv_array, kernel, fft_flags)
end
function EIstorage_gradient(model::Type{<:UnknownAcvTimeSeriesModel{1,T}}, n, kernel, fft_flags = FFTW.MEASURE) where {T}
    sdf_array = zeros(ComplexF64, npars(model), max(2n, minbins(model)))
    acv_array = zeros(ComplexF64, npars(model), 2n)
    return Sdf2EIStorageUni(sdf_array, acv_array, kernel, fft_flags)
end

# Generate required storage for computing the hessian of EI for a univariate time series model.
function EIstorage_hessian(model::Type{<:TimeSeriesModel{1,T}}, n, kernel, fft_flags = FFTW.MEASURE) where {T}
    acv_array = zeros(ComplexF64, nlowertriangle_parameter(model), 2n)
    return Acv2EIStorageUni(acv_array, kernel, fft_flags)
end
function EIstorage_hessian(model::Type{<:UnknownAcvTimeSeriesModel{1,T}}, n, kernel, fft_flags = FFTW.MEASURE) where {T}
    sdf_array = zeros(ComplexF64, nlowertriangle_parameter(model), max(2n, minbins(model)))
    acv_array = zeros(ComplexF64, nlowertriangle_parameter(model), 2n)
    return Sdf2EIStorageUni(sdf_array, acv_array, kernel, fft_flags)
end

## sdf only store
struct SdfStorageUni{M} <: TimeSeriesModelStorage
    allocatedarray::Array{ComplexF64,M}
    function SdfStorageUni(allocatedarray::Array{ComplexF64,M}) where {M}
        new{M}(allocatedarray)
    end
end

# Function to allocate storage for the function part of sdf calculations for univariate case.
function sdfstorage_function(model::Type{<:TimeSeriesModel{1,T}}, n) where {T}
    allocatedarray = zeros(ComplexF64, n)
    return SdfStorageUni(allocatedarray)
end
# Function to allocate storage for the gradient part of sdf calculations for univariate case.
function sdfstorage_gradient(model::Type{<:TimeSeriesModel{1,T}}, n) where {T}
    allocatedarray = zeros(ComplexF64, npars(model), n)
    return SdfStorageUni(allocatedarray)
end
# Function to allocate storage for the hessian part of sdf calculations for univariate case.
function sdfstorage_hessian(model::Type{<:TimeSeriesModel{1,T}}, n) where {T}
    allocatedarray = zeros(ComplexF64, nlowertriangle_parameter(model), n)
    return SdfStorageUni(allocatedarray)
end

## Unitily functions for pulling out desired vector from storage for Whittle
unpack(store::Sdf2EIStorage)    = store.acv2EI.hermitianarray
unpack(store::Acv2EIStorage)    = store.hermitianarray
unpack(store::SdfStorage)       = store.hermitianarray
unpack(store::Sdf2EIStorageUni) = store.acv2EI.allocatedarray
unpack(store::Acv2EIStorageUni) = store.allocatedarray
unpack(store::SdfStorageUni)    = store.allocatedarray


text:
include("Asd.jl")
include("ReadHeader.jl")
include("ReadFrame.jl")

function readAsd(filePath)
    open(filePath, "r") do io
        headerVersion = Int64(read(io, Int32))

        if headerVersion == 0
            header = readHeaderV0(io)
            frames = []
            for i in 1:header.numFrames
                push!(frames, readFrame(io, header))
            end

            return Asd(header, frames)
        else
            # TODO
            @assert false
        end
    end
end

max_stars_repo_name: matsunagalab/BioAFM

max_issues_repo_name: matsunagalab/BioAFM

max_forks_repo_name: matsunagalab/BioAFM


text:
using ConstrainedDynamics

path = "examples/examples_files/atlas_simple.urdf"
mech = Mechanism(path, floating=true, g = -.5)

simulate!(mech, save = true)
visualize!(mech)


text:
import Pkg; Pkg.add("Crayons")

using Crayons
using Crayons.Box

function clear()
   run(`cmd /C cls`)
end

function prompt(v::String)
   print("$CYAN_FG$v$RED_FG = $LIGHT_RED_FG")
   out = "$MAGENTA_FG$(parse(Int, readline(keep = false)))"
   return out
end

function printeq(a::String, b::String, c::String)
   clear()
   println(
      "$a$(DARK_GRAY_FG)x$(YELLOW_FG)² $RED_FG+ ",
      "$b$(DARK_GRAY_FG)x $RED_FG+ ",
      "$c $RED_FG= $(MAGENTA_FG)0",
      Crayon(reset = true)
   )
end

_b = "$(CYAN_FG)b"
_c = "$(CYAN_FG)c"
printeq("$(CYAN_FG)a", _b, _c)
a = prompt("a")
printeq(a, _b, _c)
b = prompt("b")
printeq(a, b, _c)
c = prompt("c")
printeq(a, b, c)

print("Even odin was better than this...")

max_stars_repo_name: Coalpha/lang-experimentation

max_issues_repo_name: Coalpha/lang-experimentation

max_forks_repo_name: Coalpha/lang-experimentation


text:
using Justitia
using Documenter

makedocs(;
    sitename = "Justitia.jl",
    authors = "Joey Huchette",
    format = Documenter.HTML(;
        prettyurls = get(ENV, "CI", nothing) == "true",
        assets = String[],
    ),
    modules = [Justitia],
    pages = [
        "Home" => "index.md",
        "User guide" => [
            "Interface" => "manual/interface.md",
            "Included implementations" => "manual/implementations.md",
            "An example" => "manual/example.md",
        ],
    ],
)

deploydocs(repo = "github.com/joehuchette/Justitia.jl.git", devbranch = "main")


text:
module Etd

import ..Clarus

function etdim(;params...)
  return Clarus.api_request("ETD","ETDIM";params...)
end

end #MODULE-END


text:
using Commutator: commutator!
using Base.Test

using MicroLogging
using BenchmarkTools

configure_logging(min_level=:info)
#=configure_logging(min_level=:debug)=#

SUITE = BenchmarkGroup()
# Note: this file *defines* benchmarks, but does not run them.
# See benchmark/benchmarks.jl.

#==============================================================================
                                 SUPPORT
==============================================================================#

"""commutator function that we can compare against for correctness"""
function commutator(alpha, A, B, beta, C)
    C = beta * C + alpha * (A * B - B * A)
    return C  # not in place!
end

"""The "naive" in-place implementation"""
function commutator_super!(
        alpha::Number, A::AbstractMatrix, B::AbstractMatrix, beta::Number,
        C::AbstractMatrix)
    C[:,:] = beta * C[:,:] + alpha * (A * B - B * A)
end

"""
Check the correctness of `commutator!(α, A, B, β, C)`. If
`enforce_noalloc=true`, also ensure that the call to `commutator!` causes zero
allocations.
"""
function check_commutator(α, A, B, β, C; enforce_noalloc=false)
    @debug "--------------------------------"
    C = copy(C)
    N = size(C, 1)
    mat_show = repr
    if N > 5
        mat_show = summary
    end
    @debug "check_commutator"
    @debug "alpha = $(repr(α))"
    @debug "A = $(mat_show(A))"
    @debug "B = $(mat_show(B))"
    @debug "beta = $(repr(β))"
    @debug "C = $(mat_show(C))"
    C_expect = commutator(α, A, B, β, C)
    w = which(commutator!,
              (typeof(α), typeof(A), typeof(B), typeof(β), typeof(C)))
    w = replace(repr(w), r"commutator!\(.*?\)", "commutator!(…)")
    @info " → $w"
    __, t, bytes, gctime, m = @timed commutator!(α, A, B, β, C)
    @info ("$t seconds (allocations: $bytes bytes, $(m.malloc)M " *
            "$(m.realloc)R $(m.poolalloc)P $(m.bigalloc)B)")
    δ = norm(C - C_expect)
    result = true
    if enforce_noalloc
        if bytes > 0
            result = false
            @warn "commutator! allocated >0 bytes"
        end
    end
    limit = 10 * eps(norm(C_expect))
    # My limit is somewhat heuristic, I'm not sure if it's more or less strict
    # than the builtin ≉. For added safety, require both.
    if C ≉ C_expect && abs(δ) > limit
        @warn "C  (out) = $(mat_show(C))"
        @warn "C_expect = $(mat_show(C_expect))"
        @warn "δ = $δ > $limit"
        result = false
    else
        @debug "δ = $δ < $limit"
    end
    @debug "--------------------------------"
    return result
end

"""Custom rand function that returns small values for integers"""
function myrand(S, dims...)
    if S<:Int
        return rand(1:1000, dims...)
    else
        return rand(S, dims...)
    end
end

SPRAND_DICT = Dict()
SPRAND_DUMP = joinpath(dirname(@__FILE__), "sprand.dump")
if isfile(SPRAND_DUMP)
    @info "Loading SPRAND_DICT from $SPRAND_DUMP"
    open(SPRAND_DUMP) do fh
        merge!(SPRAND_DICT, deserialize(fh))
    end
end

"""
Custom sprand function that returns small values for integers

Also caches returned matrices in `SPRAND_DICT` using `key`. This is so that we
can run repeated benchmarks (within the same REPL session) with the exact same
number of non-zero elements (nnz) in sparse matrices. Otherwise, since the
runtime of sparse-operations depends on nnz, we wouldn't be able to compare
benchmarks
"""
function mysprand(S, n, m, p, key)
    if key ∈ keys(SPRAND_DICT)
        result = SPRAND_DICT[key]
    else
        if S<:Int
            result = rand(1:1000, n, m) .* sprand(Bool, n, m, p)
        else
            result = sprand(S, n, m, p)
        end
        SPRAND_DICT[key] = result
    end
    return result
end

"""String representation of a number of two-digit precision (for labels)"""
reprnum(v) = repr(v)
reprnum(v::Float64) = @sprintf("%.2f", v)
reprnum(v::Complex128) =  @sprintf("%.2f%+.2fim", real(v), imag(v))

"""String representation of a type (for labels)"""
function reprtype(t)
    if t == Complex128
        return "Complex128"
    else
        return repr(t)
    end
end

#==============================================================================
                                  TESTS
==============================================================================#

@testset "commutator!" begin

SUITE["NDM"] = BenchmarkGroup(["commutator"])
@testset "numerical dense matrices; T=$T, β=$β, N=$N" for
        T ∈ (Complex128, Float64),
        β ∈ (0, 1, -1, rand(T)),
        N ∈ (4, 500)
    label = "T=$T, β=$(reprnum(β)), N=$N"
    α, A, B, C = rand(T), rand(T, N, N), rand(T, N, N), rand(T, N, N)
    t_super = (Number, AbstractMatrix, AbstractMatrix, Number, AbstractMatrix)
    t_args = (typeof(α), typeof(A), typeof(B), typeof(β), typeof(C))
    @test check_commutator(α, A, B, β, C; enforce_noalloc=true)
    SUITE["NDM"]["$label - super"] =
        @benchmarkable(commutator_super!($α, $A, $B, $β, Cb), setup=(Cb=copy($C)))
    SUITE["NDM"]["$label - direct"] =
        @benchmarkable(commutator!($α, $A, $B, $β, Cb), setup=(Cb=copy($C)))
end

SUITE["NDMMT"] = BenchmarkGroup(["commutator"])
@testset "numerical dense matrices (mixed types); S = $S, T=$T, β=$β" for
        S ∈ (Complex128, Float64, Int),
        T ∈ (Complex128, Float64, Int),
        β ∈ (myrand(T), myrand(S))
    # these generally cannot be mapped to gemm! and may require allocation of
    # temporary storage
    N = 4
    label = "S=$(reprtype(S)), T=$(reprtype(T)), β=$(reprnum(β)), N=$N"
    α, A, B, C = (
        myrand(T), myrand(S, N, N), myrand(T, N, N), myrand(Complex128, N, N))
    t_super = (Number, AbstractMatrix, AbstractMatrix, Number, AbstractMatrix)
    t_args = (typeof(α), typeof(A), typeof(B), typeof(β), typeof(C))
    @test check_commutator(α, A, B, β, C)
    SUITE["NDMMT"]["$label - direct"] =
        @benchmarkable(commutator!($α, $A, $B, $β, Cb), setup=(Cb = copy($C)))
    SUITE["NDMMT"]["$label - super"] =
        @benchmarkable(commutator_super!($α, $A, $B, $β, Cb), setup=(Cb=copy($C)))
end

SUITE["SPA"] = BenchmarkGroup(["commutator"])
@testset "sparse A; S = $S, T=$T" for
        S ∈ (Complex128, Float64, Int),
        T ∈ (Complex128, Float64, Int)
    N = 10
    if S == T
        for β ∈ (0, 1, -1, myrand(S))
            label = "S=$(reprtype(S)), T=$(reprtype(T)), β=$(reprnum(β)), N=$N"
            key = "SPA $(label)"
            α, A, B, C = (
                myrand(T), mysprand(S, N, N, 0.3, key), myrand(T, N, N),
                myrand(S, N, N))
            t_super = (Number, AbstractMatrix, AbstractMatrix, Number,
                       AbstractMatrix)
            t_args = (typeof(α), typeof(A), typeof(B), typeof(β), typeof(C))
            @test check_commutator(α, A, B, β, C; enforce_noalloc=true)
            SUITE["SPA"]["$label - direct"] =
                @benchmarkable(commutator!($α, $A, $B, $β, Cb), setup=(Cb = copy($C)))
            SUITE["SPA"]["$label - super"] =
                @benchmarkable(commutator_super!($α, $A, $B, $β, Cb), setup=(Cb=copy($C)))
        end
    else  # S ≠ T
        for β ∈ (0, 1, -1, myrand(Complex128))
            label = "S=$(reprtype(S)), T=$(reprtype(T)), β=$(reprnum(β)), N=$N"
            key = "SPA $(label)"
            α, A, B, C = (
                myrand(T), mysprand(S, N, N, 0.3, key), myrand(T, N, N),
                myrand(Complex128, N, N))
            C2 = copy(C); C3 = copy(C)
            t_super = (Number, AbstractMatrix, AbstractMatrix, Number,
                       AbstractMatrix)
            t_args = (typeof(α), typeof(A), typeof(B), typeof(β), typeof(C))
            @test check_commutator(α, A, B, β, C; enforce_noalloc=true)
            SUITE["SPA"]["$label - direct"] =
                @benchmarkable(commutator!($α, $A, $B, $β, Cb), setup=(Cb = copy($C)))
            SUITE["SPA"]["$label - super"] =
                @benchmarkable(commutator_super!($α, $A, $B, $β, Cb), setup=(Cb=copy($C)))
        end
    end
end

end

@info "Writing SPRAND_DICT to $SPRAND_DUMP"
open(SPRAND_DUMP, "w") do fh serialize(fh, SPRAND_DICT) end