text:
"""
Regression test for https://github.com/JuliaPolyhedra/Polyhedra.jl/issues/224
"""
function issue224test(lib::Polyhedra.Library)
    p = polyhedron(convexhull([-5012508711568619//562949953421312, 583069911833861//562949953421312, 1219842180319805//35184372088832],
                              [-3609611117314485//1125899906842624, 5282286702292177//1125899906842624, 1219842180319805//35184372088832],
                              [-7725757995938945//1125899906842624, 1462211626014367//140737488355328, 1219842180319805//35184372088832],
                              [-532368590153764480371939662909266158793799405463504814709713817//77583810724479672721865369109640386986142477600327824161374208, 3224270984442009275405360045195456034995743746806217875031525951//310335242897918690887461476438561547944569910401311296645496832, 672457679454427424790403270603729873029078800097096746747208063//19395952681119918180466342277410096746535619400081956040343552],
                              [-99439783243032230401278200429697521956561640533421926702879//31016981899445501044376765532524576292400999567207920828416, 311994187972841311191238761811528823912588633689069650103080235//66500409192411154239143785301732691570907743072093782256123904, 288196148337611788210031345503791397986095421869031213993505867//8312551149051394279892973162716586446363467884011722782015488],
                              [-690806571294091770118800876784760427996256746750368522886675821//77583810724479672721865369109639134785033939652916888662441984, 321426694532486026572630868942976798207283537518053289963575445//310335242897918690887461476438556539140135758611667554649767936, 168114419863606861724632054332712979462833232159533470697926701//4848988170279979545116585569352445924064621228307305541402624],
                              [-1492392266911428135477270336942896864090819866924399644525722049//465502864346878027654888975401162411367555092311194971098054656, 8735837263239544956921962745120391980207061989984410077199166221//1862011457387512110619555901604649645470220369244779884392218624, 2017373038363282292032368858811720337472226359496354415171746909//58187858043359753456861121925145301420944386538899371387256832],
                              [-373098066727856630393118486909715763699958036549713602976647893//116375716086719509082798053664459230332289270193787821224886272, 1091979657904943666635059112885763141035648015680316821818842795//232751432173439018165596107328918460664578540387575642449772544, 2017373038363282329633298021940319932272991184557553168606931369//58187858043359754541399026832229615166144635096893910612443136],
                              [-7070582150880849//562949953421312, 3790773064745241//562949953421312, 1219842180319805//35184372088832],
                              [-34801539891522081227113644995895817371964429166316169581661393//2770850383017131220282711035205412611141700157032485986762752, 74633028634568024242096988893935168609715564253767685190029663//11083401532068524881130844140821650444566800628129943947051008, 6004086423700245009031028586285698644710591893089440485524287//173178148938570701267669439700338288196356259814530374172672],
                              [-15819997815895236991185625966444338695356530903287555618572799//1776728489873580456771017164550118143569876505199642857701376, 964280083597461551818957674899498352631347137057929931635288993//931005728693756159348012994224261907230615288724612857435521024, 1008686519181641264266253024289396959154227260201071539037865451//29093929021679879979625406069508184600956727772644151794860032],
                              [-2923329350887855490685984645262896737347750012999614732126096039//232751432173439018165596107328919516974055261816400618912219136, 6269174405303712066292933133855537629242176139218553155719560483//931005728693756072662384429315678067896221047265602475648876544, 1008686519181641142716748220931706730473386508489464675026988381//29093929021679877270699513416114939621756907727050077364027392]),
                  lib)
    vol = volume(p)
    @test vol < 1e-12
end

max_stars_repo_name: blegat/Polyhedra

max_issues_repo_name: blegat/Polyhedra

max_forks_repo_name: blegat/Polyhedra


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
#
#   Print function of the LaTeX backend.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

# Low-level function to print the table using the LaTeX backend.
function _pt_latex(io, pinfo;
                   tf::LatexTableFormat = latex_default,
                   cell_alignment::Dict{Tuple{Int,Int},Symbol} = Dict{Tuple{Int,Int},Symbol}(),
                   formatter::Dict = Dict(),
                   highlighters::Union{LatexHighlighter,Tuple} = (),
                   hlines::AbstractVector{Int} = Int[],
                   longtable_footer::Union{Nothing,AbstractString} = nothing,
                   noheader::Bool = false,
                   nosubheader::Bool = false,
                   row_number_vline::Bool = false,
                   show_row_number::Bool = false,
                   table_type::Symbol = :tabular,
                   vlines::AbstractVector{Int} = Int[])

    @unpack_PrintInfo pinfo
    @unpack_LatexTableFormat tf

    # Let's create a `IOBuffer` to write everything and then transfer to `io`.
    buf_io = IOBuffer()
    buf    = IOContext(buf_io)

    table_type ∉ [:tabular, :longtable] &&
    error("Unknown table type $table_type. Possible values are `:tabular` or `:longtable`.")

    table_env = table_type == :tabular ? "tabular" : "longtable"

    !noheader && num_cols != header_num_cols &&
    error("The header length must be equal to the number of columns.")

    # Additional processing necessary if the user wants to print the header.
    if !noheader
        # If the user do not want to print the sub-header but wants to print the
        # header, then just force the number of rows in header to be 1.
        if nosubheader
            # Now, `header` will be a view of the first line of the matrix that
            # has the header.
            header = @view header[1:header_num_rows:end]
            header_num_rows = 1
        end
    end

    # Make sure that `highlighters` is always a tuple.
    !(highlighters isa Tuple) && (highlighters = (highlighters,))

    # Get the string which is printed when `print` is called in each element of
    # the matrix. Notice that we must create only the matrix with the printed
    # rows and columns.
    header_str = Matrix{String}(undef, header_num_rows, num_printed_cols)
    data_str   = Matrix{AbstractString}(undef, num_printed_rows, num_printed_cols)

    @inbounds for i = 1:num_printed_cols
        # Index of the i-th printed column in `data`.
        ic = id_cols[i]

        fi = haskey(formatter, i) ? formatter[i] :
                (haskey(formatter, 0) ? formatter[0] : nothing)

        if !noheader
            for j = 1:header_num_rows
                header_str[j,i] =
                    _str_latex_escaped(sprint(print, header[(ic-1)*header_num_rows + j]))
            end
        end

        for j = 1:num_printed_rows
            # Index of the j-th printed row in `data`.
            jr = id_rows[j]

            data_ij = fi != nothing ? fi(data[jr,ic], jr) : data[jr,ic]

            # Handle `nothing` and `missing`.
            if ismissing(data_ij)
                data_str_ij = "missing"
            elseif data_ij == nothing
                data_str_ij = "nothing"
            else
                data_str_ij = sprint(print, data_ij)
            end

            data_str_ij_esc = _str_latex_escaped(data_str_ij)
            data_str[j,i]   = data_str_ij_esc
        end
    end

    # Adjustments required if the user wants to print the row number.
    if show_row_number
        alignment = pushfirst!(alignment, :l)
        vlines[vlines .!= 0] .+= 1
        row_number_vline && push!(vlines, 1)
    end

    # Print LaTeX header
    # ==========================================================================

    println(buf,"""
            \\begin{$table_env}$(_latex_table_desc(alignment,
                                                   vlines,
                                                   left_vline,
                                                   mid_vline,
                                                   right_vline))
            """ * top_line)

    # Data header
    # ==========================================================================

    # Header and sub-header texts
    # --------------------------------------------------------------------------

    if !noheader
        @inbounds @views for i = 1:header_num_rows
            # The text "Row" must appear only on the first line.
            if show_row_number
                if i == 1
                    print(buf, _latex_envs("Row", header_envs))
                end

                print(buf, " & ")
            end

            for j = 1:num_printed_cols
                # Index of the j-th printed column in `data`.
                jc = id_cols[j]

                # Configure the LaTeX environments for the header and
                # sub-headers.
                if i == 1
                    envs = header_envs
                else
                    envs = subheader_envs
                end

                print(buf, _latex_envs(header_str[i,j], envs))

                j != num_printed_cols && print(buf, " & ")
            end

            if i != header_num_rows
                println(buf, " \\\\")
            else
                println(buf, " \\\\" * header_line)
            end
        end
    end

    # If we are using `longtable`, then we must mark the end of header and also
    # create the footer.
    if table_type == :longtable
        println(buf, "\\endhead")
        println(buf, bottom_line)

        # Check if the user wants a text on the footer.
        if longtable_footer != nothing
            lvline =            0 ∈ vlines ? left_vline : ""
            rvline = id_cols[end] ∈ vlines ? right_vline : ""

            env = "multicolumn{" * string(num_printed_cols) * "}" * "{r}"

            println(buf, _latex_envs(longtable_footer, env) * "\\\\")
            println(buf, bottom_line)
        end

        println(buf, "\\endfoot")
        println(buf, "\\endlastfoot")
    end

    # Data
    # ==========================================================================

    @inbounds @views for i = 1:num_printed_rows
        ir = id_rows[i]

        if show_row_number
            print(buf, string(ir) * " & ")
        end

        for j = 1:num_printed_cols
            jc = id_cols[j]

            # If we have highlighters defined, then we need to verify if this
            # data should be highlight.
            data_str_ij = data_str[i,j]

            for h in highlighters
                if h.f(data, ir, jc)
                    data_str_ij = h.fd(data, i, j, data_str[i,j])
                    break
                end
            end

            # Check the alignment of this cell.
            if haskey(cell_alignment, (i,j))
                a = _latex_alignment(cell_alignment[(i,j)])

                # Since we are using the `multicolumn`, we need to verify if the
                # column has vertical lines.
                aux_j = show_row_number ? jc + 1 : jc

                # We only need to add left vertical line if it is the first
                # column.
                lvline = (0 ∈ vlines) && (aux_j-1 == 0) ? left_vline : ""

                # For the right vertical line, we must check if it is a mid line
                # or right line.
                if aux_j ∈ vlines
                    rvline = (j == num_printed_cols) ? right_vline : mid_vline
                else
                    rvline = ""
                end

                # Wrap the data into the multicolumn environment.
                data_str_ij = _latex_envs(data_str_ij,
                                          "multicolumn{1}{$(lvline)$(a)$(rvline)}")
            end

            print(buf, data_str_ij)
            j != num_printed_cols && print(buf, " & ")
        end

        if (i != num_printed_rows) && (i in hlines)
            # Check if we must draw a horizontal line here.
            println(buf, " \\\\" * mid_line)
        elseif (i != num_printed_rows)
            println(buf, " \\\\")
        else
            println(buf, " \\\\" * bottom_line)
        end
    end

    # Print LaTeX footer
    # ==========================================================================

    println(buf, "\\end{$table_env}")

    # Print the buffer into the io.
    # ==========================================================================

    print(io, String(take!(buf_io)))

    return nothing
end

ext: jl


text:
include("dg_tree/dg.jl")
include("dg_curved/dg.jl")
include("dg_unstructured_quad/dg.jl")
include("dg_common.jl")


text:
using BasicDataLoaders
using Distributed
using Documenter
using Test

doctest(BasicDataLoaders)

@testset "input/output operations" begin
    obj = Float32[1, 2, 3]
    dir = mktempdir(cleanup = true)
    path = joinpath(dir, "test")

    T = save(path*".bson", obj)
    @test T == Array{Float32, 1}
    @test isfile(path*".bson")

    lobj = load(path*".bson")
    @test typeof(lobj) == Array{Float32, 1}
    @test all(lobj .≈ obj)

    save(path, obj)
    @test ! isfile(path)
    @test isfile(path*".bson")

    lobj = load(path*".bson")
    @test typeof(lobj) == Array{Float32, 1}
    @test all(lobj .≈ obj)
end

@testset "Data loader" begin
    obj = Float32[i for i in 1:10]

    dl = DataLoader(obj, batchsize = 1)
    @test typeof(dl) == DataLoader{Array{Float32, 1}}
    @test_throws ArgumentError DataLoader(obj, batchsize = 0)
    @test_throws ArgumentError DataLoader(obj, batchsize = -1)
    @test_throws ArgumentError DataLoader([], batchsize = 1)

    dl = DataLoader(obj, batchsize = 3)
    @test length(dl) == 4
    @test all(dl[1] .== [1, 2, 3])
    @test all(dl[2] .== [4, 5, 6])
    @test all(dl[3] .== [7, 8, 9])
    @test all(dl[4] .== [10])
    @test all(dl[1] .== dl[begin])
    @test all(dl[4] .== dl[end])

    sobj = [[1, 1], [2, 2], [3, 3]]
    v = [3, 4]
    sdl = DataLoader(obj, batchsize = 2)
    sdl2 = DataLoader(obj, batchsize = 2, preprocess_element = x -> v .* x)
    for (i, batch) in enumerate(sdl2)
        @test all((sdl2[i]) .== [v .* a for a in sdl[i]])
    end

    dl2 = DataLoader(obj, batchsize = 3, preprocess = x -> 2 .* x)
    for (i, batch) in enumerate(dl2)
        @test all((2 .* dl[i]) .== dl2[i])
    end

    N = 10
    dl = DataLoader(1:10, batchsize = 3, preprocess = x -> 2*x)
    addprocs(2)
    @everywhere using BasicDataLoaders
    res = @distributed (+) for x in dl
        sum(x)
    end
    @test res == N*(N+1)
end

max_stars_repo_name: lucasondel/DataLoaders

max_issues_repo_name: lucasondel/DataLoaders

max_forks_repo_name: lucasondel/DataLoaders


text:
module Marshalling

using JSON
using StructTypes
using XMLDict

export Leaf, Marsh
export Parse
export fields, unmarshall

"""
    Marsh

XML nodes to be mapped to a structure

Usually types that either have children or whose associated Type is comprised of other `Marsh` types
"""
abstract type Marsh end

"""
    Leaf

A `XML node` with no children

When unmarshalling, either parsed to a Type or attributes are mapped to Type properties
"""
abstract type Leaf <: Marsh end

"""
    getfield(::Type{T}, k)

Wrapper around StructTypes.julianame(tuple::Tuple, serializationname::Symbol)

Obtain property `k` of `T`, as defined under `names(::T)`, or default to `k`
"""
function getfield(::Type{T}, k::Symbol) where T<:Marsh
    return StructTypes.julianame(StructTypes.names(T), k)
end
getfield(::Type{T}, k::String) where T<:Marsh = getfield(T, Symbol(k))

"""
    Parse(::Type{T}, xml) where T<:Union{Marsh, Vector{<:Marsh}}

Apply some parsing logic or transformation to input `xml` as a whole prior to creating an instanc of `T`
"""
function Parse(::Type{T}, xml) where T<:Union{Marsh, Vector{<:Marsh}} return xml end

# """
#     Parse(::Type{T}, prop::Val{Symbol}, v)

# Apply some transformation to value `v` of a specific `prop` of T
# """
# function Parse(::Type{T}, prop::V, v) where {T<:Union{Marsh, Vector{<:Marsh}} , V<:Val}
#     return v
# end

StructTypes.construct(::Type{T}, v) where T<:Any = v
StructTypes.construct(::Type{I}, s::String) where I<:Int = parse(I, s)
StructTypes.construct(::Type{V}, s::S) where {V<:Vector, S<:AbstractString} = [string(s)]

function StructTypes.construct(::Type{Vector{T}}, xml::Vector) where T<:Marsh
    return [StructTypes.construct(T, entry) for entry in xml]
end
function StructTypes.construct(::Type{Vector{T}}, xml::XMLDict.XMLDictElement) where T<:Marsh
    return Vector{T}([StructTypes.construct(T, xml)])
end

function StructTypes.construct(::Type{T}, xml::XMLDict.XMLDictElement) where T<:Marsh
    this = T()

    for (k, v) in Parse(T, xml)
        field    = getfield(T, k)
        typemap  = Dict([name=>type for (name, type) in zip(fieldnames(T), T.types)])

        if field in fieldnames(T)
            proptype = typemap[field]
            if isleaf(k, v)
                if !(v isa proptype)
                    # attempt to change `v` to the appropriate `type`
                    if proptype <: Number
                        if (check = tryparse(proptype, v)) !== nothing
                            v = check
                        end
                    else
                        try
                            v = StructTypes.construct(proptype, v)
                        catch err
                            v = Base.convert(proptype, v)
                        end
                    end
                end
                Base.setproperty!(this, field, v)
            else
                if !isempty(v)
                    Base.setproperty!(this, field, StructTypes.construct(proptype, v))
                end
            end
        else
            return StructTypes.construct(T, v)
        end
    end
    return this
end

"""
    unmarshall(::Type{T}, xml)

Dispatcher to prevent some inconvient overwrites when defining StructTypes.construct(T, input) over Base, Core or shared Types from other libs
"""
function unmarshall(::Type{T}, xml) where T return StructTypes.construct(T, xml) end

isleaf(value) = !(value isa XMLDict.XMLDictElement)
function isleaf(k::T, v) where T<:Union{Symbol, String}
    return k isa Symbol && isleaf(v)
end
isleaf(::Type{T}) where T<:Leaf  = true
isleaf(::Type{T}) where T<:Marsh = T <: Leaf

function Base.Dict(m::T; replace=false) where T<:Marsh
    out = Dict()
    for (prop, type) in zip(fieldnames(T), T.types)
        value = Base.getproperty(m, prop)
        if type <: Vector{<:Marsh}
            value = Dict.(value)
        elseif type <: Marsh
            value = Dict(value)
        end
        if replace
            names = StructTypes.names(T)
            if (idx = findall(x->first(x) == prop, names)) |> !isempty
                prop = last(first(names[idx]))
            end
        end
        out[string(prop)] = value
    end
    return out
end
function JSON.json(m::T, args...; replace=true, kwargs...) where T<:Marsh
    return JSON.json(Dict(m; replace=replace), args...; kwargs...)
end

end #end module


text:
module Datasets
import ..Cubes.Axes: axsym, axname, CubeAxis, findAxis, CategoricalAxis, RangeAxis, caxes
import ..Cubes: YAXArray, concatenatecubes, CleanMe
using ...YAXArrays: YAXArrays, YAXDefaults
using DataStructures: OrderedDict, counter
using Dates: Day, Hour, Minute, Second, Month, Year, Date, DateTime, TimeType, AbstractDateTime
using IntervalSets: Interval, (..)
using CFTime: timedecode, timeencode, DateTimeNoLeap, DateTime360Day, DateTimeAllLeap
using YAXArrayBase
using YAXArrayBase: iscontdimval
using DiskArrayTools: CFDiskArray, ConcatDiskArray
using Glob: glob

export Dataset, Cube, open_dataset

struct Dataset
    cubes::OrderedDict{Symbol,YAXArray}
    axes::Dict{Symbol,CubeAxis}
end
function Dataset(; cubesnew...)
    axesall = Set{CubeAxis}()
    foreach(values(cubesnew)) do c
        ax = caxes(c)
        foreach(a -> push!(axesall, a), ax)
    end
    axesall = collect(axesall)
    axnameall = axsym.(axesall)
    axesnew = Dict{Symbol,CubeAxis}(axnameall[i] => axesall[i] for i = 1:length(axesall))
    Dataset(OrderedDict(cubesnew), axesnew)
end

function Base.show(io::IO, ds::Dataset)
    println(io, "YAXArray Dataset")
    println(io, "Dimensions: ")
    foreach(a -> println(io, "   ", a), values(ds.axes))
    print(io, "Variables: ")
    foreach(i -> print(io, i, " "), keys(ds.cubes))
end
function Base.propertynames(x::Dataset, private::Bool = false)
    if private
        Symbol[:cubes; :axes; collect(keys(x.cubes)); collect(keys(x.axes))]
    else
        Symbol[collect(keys(x.cubes)); collect(keys(x.axes))]
    end
end
function Base.getproperty(x::Dataset, k::Symbol)
    if k === :cubes
        return getfield(x, :cubes)
    elseif k === :axes
        return getfield(x, :axes)
    else
        x[k]
    end
end
Base.getindex(x::Dataset, i::Symbol) =
    haskey(x.cubes, i) ? x.cubes[i] :
    haskey(x.axes, i) ? x.axes[i] : throw(ArgumentError("$i not found in Dataset"))
function Base.getindex(x::Dataset, i::Vector{Symbol})
    cubesnew = [j => x.cubes[j] for j in i]
    Dataset(; cubesnew...)
end

function fuzzyfind(s::String, comp::Vector{String})
    sl = lowercase(s)
    f = findall(i -> startswith(lowercase(i), sl), comp)
    if length(f) != 1
        throw(KeyError("Name $s not found"))
    else
        f[1]
    end
end
function Base.getindex(x::Dataset, i::Vector{String})
    istr = string.(keys(x.cubes))
    ids = map(name -> fuzzyfind(name, istr), i)
    syms = map(j -> Symbol(istr[j]), ids)
    cubesnew = [Symbol(i[j]) => x.cubes[syms[j]] for j = 1:length(ids)]
    Dataset(; cubesnew...)
end
Base.getindex(x::Dataset, i::String) = getproperty(x, Symbol(i))
function subsetcube(x::Dataset; var = nothing, kwargs...)
    if var === nothing
        cc = x.cubes
        Dataset(; map(ds -> ds => subsetcube(cc[ds]; kwargs...), collect(keys(cc)))...)
    elseif isa(var, String) || isa(var, Symbol)
        subsetcube(getproperty(x, Symbol(var)); kwargs...)
    else
        cc = x[var].cubes
        Dataset(; map(ds -> ds => subsetcube(cc[ds]; kwargs...), collect(keys(cc)))...)
    end
end
function collectdims(g)
    dlist = Set{Tuple{String,Int,Int}}()
    varnames = get_varnames(g)
    foreach(varnames) do k
        d = get_var_dims(g, k)
        v = get_var_handle(g, k)
        for (len, dname) in zip(size(v), d)
            if !occursin("bnd", dname) && !occursin("bounds", dname)
                datts = get_var_attrs(g, dname)
                offs = get(datts, "_ARRAY_OFFSET", 0)
                push!(dlist, (dname, offs, len))
            end
        end
    end
    outd = Dict(d[1] => (ax = toaxis(d[1], g, d[2], d[3]), offs = d[2]) for d in dlist)
    length(outd) == length(dlist) ||
        throw(ArgumentError("All Arrays must have the same offset"))
    outd
end

function toaxis(dimname, g, offs, len)
    axname = dimname
    if !haskey(g, dimname)
        return RangeAxis(dimname, 1:len)
    end
    ar = get_var_handle(g, dimname)
    aratts = get_var_attrs(g, dimname)
    if uppercase(axname) == "TIME" && haskey(aratts, "units")
        tsteps = try
            timedecode(ar[:], aratts["units"], lowercase(get(aratts, "calendar", "standard")))
        catch
            ar[:]
        end
        RangeAxis(dimname, tsteps[offs+1:end])
    elseif haskey(aratts, "_ARRAYVALUES")
        vals = identity.(aratts["_ARRAYVALUES"])
        CategoricalAxis(axname, vals)
    else
        axdata = cleanaxiselement.(ar[offs+1:end])
        axdata = testrange(axdata)
        if eltype(axdata) <: AbstractString ||
           (!issorted(axdata) && !issorted(axdata, rev = true))
            CategoricalAxis(axname, axdata)
        else
            RangeAxis(axname, axdata)
        end
    end
end
propfromattr(attr) = Dict{String,Any}(filter(i -> i[1] != "_ARRAY_DIMENSIONS", attr))

#there are problems with saving custom string types to netcdf, so we clean this when creating the axis:
cleanaxiselement(x::AbstractString) = String(x)
cleanaxiselement(x::String) = x
cleanaxiselement(x::TimeType) = DateTime(x)
cleanaxiselement(x::Union{Date, DateTime}) = x
cleanaxiselement(x) = x

"Test if data in x can be approximated by a step range"
function testrange(x)
    r = range(first(x), last(x), length = length(x))
    all(i -> isapprox(i...), zip(x, r)) ? r : x
end

function testrange(x::AbstractArray{<:Integer})
    steps = diff(x)
    if all(isequal(steps[1]), steps) && !iszero(steps[1])
        return range(first(x), step = steps[1], length(x))
    else
        return x
    end
end

testrange(x::AbstractArray{<:AbstractString}) = x

_glob(x) = startswith(x, "/") ? glob(x[2:end], "/") : glob(x)

open_mfdataset(g::AbstractString; kwargs...) = open_mfdataset(_glob(g); kwargs...)
open_mfdataset(g::Vector{<:AbstractString}; kwargs...) =
    merge_datasets(map(i -> open_dataset(i; kwargs...), g))

function open_dataset(g; driver = :all)
    g = to_dataset(g, driver = driver)
    isempty(get_varnames(g)) && throw(ArgumentError("Group does not contain datasets."))
    dimlist = collectdims(g)
    dnames = string.(keys(dimlist))
    varlist = filter(get_varnames(g)) do vn
        upname = uppercase(vn)
        !occursin("BNDS", upname) &&
            !occursin("BOUNDS", upname) &&
            !any(i -> isequal(upname, uppercase(i)), dnames)
    end
    allcubes = OrderedDict{Symbol,YAXArray}()
    for vname in varlist
        vardims = get_var_dims(g, vname)
        iax = [dimlist[vd].ax for vd in vardims]
        offs = [dimlist[vd].offs for vd in vardims]
        subs = if all(iszero, offs)
            nothing
        else
            ntuple(i -> (offs[i]+1):(offs[i]+length(iax[i])), length(offs))
        end
        ar = get_var_handle(g, vname)
        att = get_var_attrs(g, vname)
        if subs !== nothing
            ar = view(ar, subs...)
        end
        if !haskey(att, "name")
            att["name"] = vname
        end
        allcubes[Symbol(vname)] = YAXArray(iax, ar, propfromattr(att), cleaner = CleanMe[])
    end
    sdimlist = Dict(Symbol(k) => v.ax for (k, v) in dimlist)
    Dataset(allcubes, sdimlist)
end
Base.getindex(x::Dataset; kwargs...) = subsetcube(x; kwargs...)
YAXDataset(; kwargs...) = Dataset(YAXArrays.YAXDefaults.cubedir[]; kwargs...)

function Cube(ds::Dataset; joinname = "Variable")

    dl = collect(keys(ds.axes))
    dls = string.(dl)
    length(ds.cubes) == 1 && return first(values(ds.cubes))
    #TODO This is an ugly workaround to merge cubes with different element types,
    # There should bde a more generic solution
    eltypes = map(eltype, values(ds.cubes))
    majtype = findmax(counter(eltypes))[2]
    newkeys = Symbol[]
    for k in keys(ds.cubes)
        c = ds.cubes[k]
        if all(axn -> findAxis(axn, c) !== nothing, dls) && eltype(c) == majtype
            push!(newkeys, k)
        end
    end
    if length(newkeys) == 1
        return ds.cubes[first(newkeys)]
    else
        varax = CategoricalAxis(joinname, string.(newkeys))
        cubestomerge = [ds.cubes[k] for k in newkeys]
        foreach(
            i -> haskey(i.properties, "name") && delete!(i.properties, "name"),
            cubestomerge,
        )
        return concatenatecubes(cubestomerge, varax)
    end
end

"""
function createdataset(DS::Type,axlist; kwargs...)

  Creates a new datacube with axes specified in `axlist`. Each axis must be a subtype
  of `CubeAxis`. A new empty Zarr array will be created and can serve as a sink for
  `mapCube` operations.

  ### Keyword arguments

  * `folder=tempname()` location where the new cube is stored
  * `T=Union{Float32,Missing}` data type of the target cube
  * `chunksize = ntuple(i->length(axlist[i]),length(axlist))` chunk sizes of the array
  * `chunkoffset = ntuple(i->0,length(axlist))` offsets of the chunks
  * `persist::Bool=true` shall the disk data be garbage-collected when the cube goes out of scope?
  * `overwrite::Bool=false` overwrite cube if it already exists
  * `properties=Dict{String,Any}()` additional cube properties
  * `fillvalue= T>:Missing ? defaultfillval(Base.nonmissingtype(T)) : nothing` fill value
  * `datasetaxis="Variable"` special treatment of a categorical axis that gets written into separate zarr arrays
  """
function createdataset(
    DS,
    axlist;
    path = "",
    persist = nothing,
    T = Union{Float32,Missing},
    chunksize = ntuple(i -> length(axlist[i]), length(axlist)),
    chunkoffset = ntuple(i -> 0, length(axlist)),
    overwrite::Bool = false,
    properties = Dict{String,Any}(),
    datasetaxis = "Variable",
    kwargs...,
)
    if persist === nothing
        persist = !isempty(path)
    end
    path = getsavefolder(path, persist)
    check_overwrite(path, overwrite)
    splice_generic(x::AbstractArray, i) = [x[1:(i-1)]; x[(i+1:end)]]
    splice_generic(x::Tuple, i) = (x[1:(i-1)]..., x[(i+1:end)]...)
    myar = create_empty(DS, path)
    finalperm = nothing
    idatasetax = datasetaxis === nothing ? nothing : findAxis(datasetaxis, axlist)
    if idatasetax !== nothing
        groupaxis = axlist[idatasetax]
        axlist = splice_generic(axlist, idatasetax)
        chunksize = splice_generic(chunksize, idatasetax)
        chunkoffset = splice_generic(chunkoffset, idatasetax)
        finalperm =
            ((1:idatasetax-1)..., length(axlist) + 1, (idatasetax:length(axlist))...)
    else
        groupaxis = nothing
    end
    foreach(axlist, chunkoffset) do ax, co
        arrayfromaxis(myar, ax, co)
    end
    attr = properties
    s = map(length, axlist) .+ chunkoffset
    if all(iszero, chunkoffset)
        subs = nothing
    else
        subs = ntuple(length(axlist)) do i
            (chunkoffset[i]+1):(length(axlist[i].values)+chunkoffset[i])
        end
    end
    if groupaxis === nothing
        cubenames = ["layer"]
    else
        cubenames = groupaxis.values
    end
    cleaner = CleanMe[]
    persist || push!(cleaner, CleanMe(path, false))
    allcubes = map(cubenames) do cn
        axnames = map(axname, axlist)
        axlengths = map(i -> length(get_var_handle(myar, i)), axnames)
        v = if allow_missings(myar) || !(T >: Missing)
            add_var(myar, T, cn, axlengths, axnames, attr; chunksize = chunksize, kwargs...)
        else
            S = Base.nonmissingtype(T)
            if !haskey(attr, "missing_value")
                attr["missing_value"] = YAXArrayBase.defaultfillval(S)
            end
            v = add_var(
                myar,
                S,
                cn,
                axlengths,
                axnames,
                attr;
                chunksize = chunksize,
                kwargs...,
            )
            CFDiskArray(v, attr)
        end
        if subs !== nothing
            v = view(v, subs...)
        end
        YAXArray(axlist, v, propfromattr(attr), cleaner = cleaner)
    end
    if groupaxis === nothing
        return allcubes[1], allcubes[1]
    else
        cube = concatenatecubes(allcubes, groupaxis)
        return permutedims(cube, finalperm), cube
    end
end

function getsavefolder(name, persist)
    if isempty(name)
        name = persist ? [splitpath(tempname())[end]] : splitpath(tempname())[2:end]
        joinpath(YAXDefaults.workdir[], name...)
    else
        (occursin("/", name) || occursin("\\", name)) ? name :
        joinpath(YAXDefaults.workdir[], name)
    end
end

function check_overwrite(newfolder, overwrite)
    if isdir(newfolder) || isfile(newfolder)
        if overwrite
            rm(newfolder, recursive = true)
        else
            error(
                "$(newfolder) already exists, please pick another name or use `overwrite=true`",
            )
        end
    end
end

function arrayfromaxis(p, ax::CubeAxis, offs)
    data, attr = dataattfromaxis(ax, offs)
    attr["_ARRAY_OFFSET"] = offs
    za = add_var(p, data, axname(ax), (axname(ax),), attr)
    za
end

prependrange(r::AbstractRange, n) =
    n == 0 ? r : range(first(r) - n * step(r), last(r), length = n + length(r))
function prependrange(r::AbstractVector, n)
    if n == 0
        return r
    else
        step = r[2] - r[1]
        first = r[1] - step * n
        last = r[1] - step
        radd = range(first, last, length = n)
        return [radd; r]
    end
end

defaultcal(::Type{<:TimeType}) = "standard"
defaultcal(::Type{<:DateTimeNoLeap}) = "noleap"
defaultcal(::Type{<:DateTimeAllLeap}) = "allleap"
defaultcal(::Type{<:DateTime360Day}) = "360_day"

datetodatetime(vals::AbstractArray{<:Date}) = DateTime.(vals)
datetodatetime(vals) = vals
toaxistype(x) = x
toaxistype(x::Array{<:AbstractString}) = string.(x)
toaxistype(x::Array{String}) = x

function dataattfromaxis(ax::CubeAxis, n)
    prependrange(toaxistype(ax.values), n), Dict{String,Any}()
end

# function dataattfromaxis(ax::CubeAxis,n)
#     prependrange(1:length(ax.values),n), Dict{String,Any}("_ARRAYVALUES"=>collect(ax.values))
# end
function dataattfromaxis(ax::CubeAxis{T}, n) where {T<:TimeType}
    data = timeencode(datetodatetime(ax.values), "days since 1980-01-01", defaultcal(T))
    prependrange(data, n),
    Dict{String,Any}("units" => "days since 1980-01-01", "calendar" => defaultcal(T))
end

#The good old Cube function:
Cube(s::String; kwargs...) = Cube(open_dataset(s); kwargs...)
function Cube(; kwargs...)
    if !isempty(YAXArrays.YAXDefaults.cubedir[])
        Cube(YAXArrays.YAXDefaults.cubedir[]; kwargs...)
    else
        error("A path should be specified")
    end
end

#Defining joins of Datasets
abstract type AxisJoin end
struct AllEqual <: AxisJoin
    ax::Any
end
struct SortedRanges <: AxisJoin
    axlist::Any
    perm::Any
end
blocksize(x::AllEqual) = 1
blocksize(x::SortedRanges) = length(x.axlist)
getperminds(x::AllEqual) = 1:1
getperminds(x::SortedRanges) = x.perm
wholeax(x::AllEqual) = x.ax
wholeax(x::SortedRanges) = reduce(vcat, x.axlist[x.perm])
struct NewDim <: AxisJoin
    newax::Any
end
#Test for a range of categorical axes how to concatenate them
function analyse_axjoin_ranges(dimvallist)
    firstax = first(dimvallist)
    if all(isequal(firstax), dimvallist)
        return AllEqual(firstax)
    end
    revorder = if all(issorted, dimvallist)
        false
    elseif all(i -> issorted(i, rev = true), dimvallist)
        true
    else
        error("Dimension values are not sorted")
    end
    function ltfunc(ax1, ax2)
        min1, max1 = extrema(ax1)
        min2, max2 = extrema(ax2)
        if max1 < min2
            return true
        elseif min1 > max2
            return false
        else
            error("Dimension ranges overlap")
        end
    end
    sp = sortperm(dimvallist, rev = revorder, lt = ltfunc)
    SortedRanges(dimvallist, sp)
end
using YAXArrayBase: YAXArrayBase, getdata, getattributes, yaxcreate
function create_mergedict(dimvallist)
    allmerges = Dict{Symbol,Any}()

    for (axn, dimvals) in dimvallist
        iscont = iscontdimval.(dimvals)
        if all(iscont)
            allmerges[axn] = analyse_axjoin_ranges(dimvals)
        elseif any(iscont)
            error("Mix of continous and non-continous values")
        else
            allmerges[axn] = analyse_axjoin_categorical(dimvals)
        end
    end
    allmerges
end
function merge_datasets(dslist)
    allaxnames = counter(Symbol)
    for ds in dslist, k in keys(ds.axes)
        push!(allaxnames, k)
    end
    dimvallist = Dict(ax => map(i -> i.axes[ax].values, dslist) for ax in keys(allaxnames))
    allmerges = create_mergedict(dimvallist)
    repvars = counter(Symbol)
    for ds in dslist, v in keys(ds.cubes)
        push!(repvars, v)
    end
    tomergevars = filter(i -> i[2] == length(dslist), repvars)
    mergedvars = Dict{Symbol,Any}()
    for v in keys(tomergevars)
        dn = YAXArrayBase.dimnames(first(dslist)[v])
        howmerge = getindex.(Ref(allmerges), dn)
        sizeblockar = map(blocksize, howmerge)
        perminds = map(getperminds, howmerge)
        @assert length(dslist) == prod(sizeblockar)
        vcol = map(i -> getdata(i[v]), dslist)
        allatts =
            mapreduce(i -> getattributes(i[v]), merge, dslist, init = Dict{String,Any}())
        aa = [vcol[i] for (i, _) in enumerate(Iterators.product(perminds...))]
        dvals = map(wholeax, howmerge)
        mergedvars[v] = yaxcreate(YAXArray, ConcatDiskArray(aa), dn, dvals, allatts)
    end
    Dataset(; mergedvars...)
end

end


text:
# Numerical utilities to support implementation

using LinearAlgebra.BLAS: nrm2
using LinearAlgebra.LAPACK: potrf!, potri!, potrs!

function printf_mat(x::AbstractMatrix)
    @inbounds for i = 1:size(x,1)
        for j = 1:size(x,2)
            @printf("%8.4f ", x[i,j])
        end
        println()
    end
end

function adddiag!(A::Matrix, a::Number)
    m, n = size(A)
    m == n || error("A must be square.")
    if a != 0.0
        for i = 1:m
            @inbounds A[i,i] += a
        end
    end
    return A
end

normalize1!(a) = rmul!(a, 1 / sum(a))

function normalize1_cols!(a)
    for j = 1:size(a,2)
        normalize1!(view(a, :, j))
    end
end

function projectnn!(A::AbstractArray{T}) where T
    # project back all entries to non-negative domain
    @inbounds for i = 1:length(A)
        if A[i] < zero(T)
            A[i] = zero(T)
        end
    end
end

function posneg!(A::AbstractArray{T},
                 Ap::AbstractArray{T}, An::AbstractArray{T}) where T
    # decompose A into positive part Ap and negative part An
    # s.t. A = Ap - An

    n = length(A)
    length(Ap) == length(An) == n || error("Input dimensions mismatch.")

    @inbounds for i = 1:n
        ai = A[i]
        if ai >= zero(T)
            Ap[i] = ai
            An[i] = zero(T)
        else
            Ap[i] = zero(T)
            An[i] = -ai
        end
    end
end

function pdsolve!(A, x, uplo::Char='U')
    # A must be positive definite
    # x <- inv(A) * x
    # both A and x will be overriden

    potrf!(uplo, A)
    potrs!(uplo, A, x)
end

function pdrsolve!(A, B, x, uplo::Char='U')
    # B must be positive definite
    # x <- A * inv(B)
    # both B and x will be overriden

    # inverse B in place
    potrf!(uplo, B)
    potri!(uplo, B)
    copytri!(B, uplo)

    # x <- A * B (the inversed one)
    mul!(x, A, B)
end


text:
"""
    Sigmoid(x)

Sigmoid activation function.

    Usage:
```julia-repl
Sigmoid(x)
```
Arguments:
- `x` : Input value.
"""
function Sigmoid(x::Float64)
    x₁ = 1.0 / (1.0 + exp(-x))
    return x₁
end


text:
using Test

using Distances
using Clustering

@testset "kmedoids() (k-medoids)" begin

@testset "Argument checks" begin
    Random.seed!(34568)
    @test_throws ArgumentError kmedoids(randn(2, 3), 1)
    @test_throws ArgumentError kmedoids(randn(2, 3), 4)
    costs = inv.(max.(pairwise(Euclidean(), randn(2, 3), dims=2), 0.1))
    @test kmedoids(costs, 2) isa KmedoidsResult
    @test_throws ArgumentError kmedoids(costs, 2, display=:mylog)
    for disp in keys(Clustering.DisplayLevels)
        @test kmedoids(costs, 2, display=disp) isa KmedoidsResult
    end
end

Random.seed!(34568)

d = 3
n = 200
k = 10

X = rand(d, n)
costs = pairwise(SqEuclidean(), X, dims=2)
@assert size(costs) == (n, n)

R = kmedoids(costs, k)
@test isa(R, KmedoidsResult)
@test nclusters(R) == k
@test length(R.medoids) == length(unique(R.medoids))
@test all(a -> 1 <= a <= k, R.assignments)
@test R.assignments[R.medoids] == 1:k # Every medoid should belong to its own cluster
@test sum(counts(R)) == n
@test wcounts(R) == counts(R)
@test R.acosts == costs[LinearIndices((n, n))[CartesianIndex.(R.medoids[R.assignments], 1:n)]]
@test isapprox(sum(R.acosts), R.totalcost)
@test R.converged

# this data set has three obvious groups:
# group 1: [1, 3, 4], values: [1, 2, 3]
# group 2: [2, 5, 7], values: [6, 7, 8]
# group 3: [6, 8, 9], values: [21, 20, 22]
#

X = reshape(map(Float64, [1, 6, 2, 3, 7, 21, 8, 20, 22]), 1, 9)
costs = pairwise(SqEuclidean(), X, dims=2)

R = kmedoids!(costs, [1, 2, 6])
@test isa(R, KmedoidsResult)
@test nclusters(R) == 3
@test R.medoids == [3, 5, 6]
@test R.assignments == [1, 2, 1, 1, 2, 3, 2, 3, 3]
@test counts(R) == [3, 3, 3]
@test wcounts(R) == counts(R)
@test R.acosts ≈ [1, 1, 0, 1, 0, 0, 1, 1, 1]
@test R.totalcost ≈ 6.0
@test R.converged

end


text:
function construct_services!(
    psi_container::PSIContainer,
    sys::PSY.System,
    services_template::Dict{Symbol, ServiceModel},
    devices_template::Dict{Symbol, DeviceModel},
)
    isempty(services_template) && return
    services_mapping = PSY.get_contributing_device_mapping(sys)
    for service_model in values(services_template)
        @debug "Building $(service_model.service_type) with $(service_model.formulation) formulation"
        services = PSY.get_components(service_model.service_type, sys)
        if validate_available_services(service_model.service_type, services)
            construct_service!(
                psi_container,
                services,
                services_mapping,
                service_model,
                devices_template,
            )
        end
    end
    return
end

function construct_service!(
    psi_container::PSIContainer,
    services::IS.FlattenIteratorWrapper{SR},
    services_mapping::PSY.ServiceContributingDevicesMapping,
    model::ServiceModel{SR, RangeReserve},
    devices_template::Dict{Symbol, DeviceModel},
) where {SR <: PSY.Reserve}
    time_steps = model_time_steps(psi_container)
    names = (PSY.get_name(s) for s in services)

    incompatible_device_types = Vector{DataType}()
    for model in values(devices_template)
        formulation = get_formulation(model)
        if formulation == FixedOutput
            @info "$(formulation) for $(get_device_type(model)) is not compatible with the provision of reserve services"
            push!(incompatible_device_types, get_device_type(model))
        end
    end

    if model_has_parameters(psi_container)
        container = add_param_container!(
            psi_container,
            UpdateRef{SR}("service_requirement", "get_requirement"),
            names,
            time_steps,
        )
        get_parameter_array(container)
    end

    add_cons_container!(psi_container, constraint_name(REQUIREMENT, SR), names, time_steps)

    for service in services
        contributing_devices =
            services_mapping[(
                type = typeof(service),
                name = PSY.get_name(service),
            )].contributing_devices
        if !isempty(incompatible_device_types)
            contributing_devices =
                [d for d in contributing_devices if typeof(d) ∉ incompatible_device_types]
        end

        #Variables
        activeservice_variables!(psi_container, service, contributing_devices)
        # Constraints
        service_requirement_constraint!(psi_container, service, model)
        modify_device_model!(devices_template, model, contributing_devices)
    end
    return
end


text:
# Collection of programs to test source reads #

# Move an immediate to dmem location 0
@asap4asm function move_immediate()
    MOVI(dmem[0], 100)
end

# Read the return address.
@asap4asm function move_return_address()
    MOVE(dmem[0], ret)
end

# Read from accumulator
@asap4asm function move_accumulator()
    MOVE(dmem[0], acc)
end

# -- Bypass registers -- #
@asap4asm function bypass_1()
    MOVI(null, 10)
    MOVE(dmem[0], bypass[1])
end

@asap4asm function bypass_2()
    MOVI(null, 10)
    NOP()
    MOVE(dmem[0], bypass[2])
end

@asap4asm function bypass_3()
    MOVI(null, 10)
    NOP()
    NOP()
    MOVE(dmem[0], bypass[3])
end

# Use MOVI -> bypass register 1 -> MACL to get a multiplication
# of two numbers. Capture bypass register 4 using another MOVE
# instruction.
@asap4asm function bypass_4()
    MOVI(dmem[0], 100)
    MULTL(dmem[1], bypass[1], 10)
    NOP()
    NOP()
    NOP()
    ADD(dmem[2], bypass[4], bypass[4])
end

# This is like the bypass_4 test, but we have 1 more NOP between the
# multiply and the add.
@asap4asm function bypass_5()
    MOVI(dmem[0], 100)
    MULTL(dmem[1], bypass[1], 10)
    NOP()
    NOP()
    NOP()
    NOP()
    ADD(dmem[2], bypass[5], bypass[5])
end

# -- Pointer Dereferences -- #
@asap4asm function pointers()
    MOVE(dmem[0], pointer[0])
    MOVE(dmem[1], pointer[1])
    MOVE(dmem[2], pointer[2])
    MOVE(dmem[3], pointer[3])
end

# -- Address Generator
@asap4asm function address_generators()
    # Just dereference the address generator without incrementing.
    MOVE(dmem[0], ag[0])
    # Now dereference it, but increment.
    MOVE(dmem[1], ag_pi[0])
    # Dereference and add, but this time use an ADD instruction to make
    # sure that the address generator is not incremented twice.
    ADD(dmem[2], ag_pi[0], ag_pi[0])
    # Make sure the increment worked.
    MOVE(dmem[3], ag[0])
end

# -- Pointer bypass -- #
@asap4asm function pointer_bypass()
    # Want to dereference memory 12
    MOVI(null, 12)
    NOP()
    NOP()
    # Dereference by bypass pointer and store the result to memory 0.
    #
    # Expect the result to be 10.
    MOVE(dmem[0], pointer_bypass)
end

# -- Fifo reads -- #
# TODO: When writeback to outputs gets implemented, make sure stages 5, 6,
# and 7 don't stall if one of them is writing to an output.
@asap4asm function fifo_read()
    # This instruction should stall before writeback.
    MOVI(dmem[0], 10)
    # Read ibuf and increment
    MOVE(dmem[1], ibuf[0])
    # Read again, but don't increment
    MOVE(dmem[2], ibuf_next[0])
    MOVE(dmem[3], ibuf[0])
end


text:
"""
    basin_entropy(basin; eps_x=20, eps_y=20)
This algorithm computes the basin entropy of a computed basin of attraction on a regular grid.
The function return the basin entropy and the boundary basin entropy.

[A. Daza, A. Wagemakers, B. Georgeot, D. Guéry-Odelin and M. A. F. Sanjuán, Basin entropy:
a new tool to analyze uncertainty in dynamical systems, Sci. Rep., 6, 31416, (2016).]

## Arguments
* `basin` : the matrix containing the information of the basin.

## Keyword arguments
* `eps_x`, `eps_y` : define the size of the box that will be used to sample the basin
"""
function basin_entropy(basin; eps_x=20, eps_y=20)
    r,c= size(basin)
    vals = unique(basin)
    S=Int16(length(vals))
    pn=zeros(Float64,1,S)
    Sb=0
    Nb=0
    N=0
    for x in 1:eps_x:(r-eps_x+1), y in 1:eps_y:(c-eps_y+1)
        I = CartesianIndices((x:x+eps_x-1,y:y+eps_y-1))
        box_values=[basin[k] for k in I]
        N=N+1
        for (k,v) in enumerate(vals)
            pn[k]=count(x->x==v,box_values)/length(box_values)
        end
        Nb = Nb + (length(unique(box_values))>1)
        Sb = Sb + sum(entropy.(pn))
    end
        return Sb/N, Sb/Nb
end

function entropy(p::Float64)
    if p == 0.
        h = 0.
    else
        h=p*log(1/p)
    end
    return h
end

function basin_entropy(bsn::BasinInfo; eps_x=20, eps_y=20)
    ind  = findall(iseven.(bsn.basin) .== true)
    # Not sure it is necessary to deepcopy
    basin_test = deepcopy(bsn.basin)
    # Remove attractors from the picture
    for k in ind; basin_test[k] = basin_test[k]+1; end
    basin_entropy(basin_test; eps_x=eps_x, eps_y=eps_y)
end

"""
    basin_stability(basin)
This algorithm computes the basin stability of a computed basin of attraction on a regular grid.
This function returns a vector with the relative volume of the basins.

[P. Menck, J. Heitzig, N. Marwan et al. How basin stability complements the linear-stability paradigm. Nature Phys 9, 89–92 (2013).]

## Arguments
* `basin` : the matrix containing the information of the basin.

"""
function basin_stability(basin)

N = length(unique(basin))
v = zeros(1,N)

for (k,b) in enumerate(unique(basin))
        v[k] = count(basin .== b)/length(basin)
end

return v

end

function basin_stability(bsn::BasinInfo)
    ind  = findall(iseven.(bsn.basin) .== true)
    # Not sure it is necessary to deepcopy
    basin_test = deepcopy(bsn.basin)
    # Remove attractors from the picture
    for k in ind; basin_test[k] = basin_test[k]+1; end
    return basin_stability(basin_test)
end

"""
    fractal_test(bsn::BasinInfo)

## Arguments
* `bsn` :

## Keyword arguments
*

"""
function fractal_test(bsn::BasinInfo)

# disk radius
εₕ=5

# Sanity check.
if length(bsn.xg)/εₕ < 50
    @warn "Maybe the size of the grid is not fine enough."
end

nx = length(bsn.xg)
ny = length(bsn.yg)
xi=bsn.xg[1]; xf=bsn.xg[end]; yi=bsn.yg[1]; yf=bsn.yg[end];

# Now get the values in the boxes.
while Nb < 1e4

x = rand()*(xf-xi)+xi
y = rand()*(yf-yi)+yi
ind = find_in_range(x,y,εₕ)

end

end

function find_in_range(x,y,εₕ)

end


text:
using Statistics
using NCDatasets

using ArgParse
using JSON

include("LinearRegression.jl")

correlation = (x1, x2) -> x1' * x2 / (sum(x1.^2)*sum(x2.^2)).^0.5

function parse_commandline()

    s = ArgParseSettings()
    @add_arg_table s begin

        "--data-file"
            help = "Atm data file. New variable will be appended."
            arg_type = String
            required = true

        "--domain-file"
            help = "Domain file."
            arg_type = String
            required = true

        "--output-file"
            help = "Output file."
            arg_type = String
            required = true

        "--beg-year"
            help = "Year of begin."
            arg_type = Int64
            required = true

        "--end-year"
            help = "Year of end."
            arg_type = Int64
            required = true

    end

    return parse_args(ARGS, s)
end

parsed = parse_commandline()
print(json(parsed, 4))

output_file = parsed["output-file"]

Dataset(parsed["data-file"], "r") do ds

    beg_t = (parsed["beg-year"] - 1) * 12 + 1
    end_t = (parsed["end-year"] - 1) * 12 + 12
    rng = (:,:,beg_t:end_t)

    global total_precip = replace(ds["PRECC"][rng...] + ds["PRECL"][rng...], missing=>NaN)
    global (Nx, Ny, Nt) = size(total_precip)

    if mod(Nt, 12) != 0
        ErrorException("Time record is not multiple of 12") |> throw
    end

    global nyears = Int64(Nt / 12)
end

Dataset(parsed["domain-file"], "r") do ds
    global mask = replace(ds["mask"], missing=>NaN)
end

# M   = Mean
# VAR = Variance

# MM   = Monthly Mean
# MVAR = Monthly Variance

# ZM   = Zonal Mean
# ZVAR = Zonal Variance

# ZMM   = Zonal Monthly Mean
# ZMVAR = Zonal Monthly Variance

total_precip_M     = zeros(Float64, Nx, Ny)
total_precip_VAR   = zeros(Float64, Nx, Ny)

total_precip_MM     = zeros(Float64, Nx, Ny, 12)
total_precip_MVAR   = zeros(Float64, Nx, Ny, 12)

total_precip_ZM     = zeros(Float64, Ny)
total_precip_ZVAR   = zeros(Float64, Ny)

total_precip_ZMM    = zeros(Float64, Ny, 12)
total_precip_ZMVAR  = zeros(Float64, Ny, 12)

x = collect(Float64, 1:Nt)
for i=1:Nx, j=1:Ny

    d = view(total_precip, i, j, :)

    total_precip_M[i, j]   = mean(d)
    total_precip_VAR[i, j] = std(d)

    d_wrap = reshape(d, 12, :)
    total_precip_MM[i, j, :]   = mean( d_wrap, dims=2 )[:, 1]
    total_precip_MVAR[i, j, :] = std( d_wrap, dims=2 )[:, 1]

end

for j=1:Ny

    d = mean(view(total_precip, :, j, :), dims=1)[1, :]

    total_precip_ZM[j]   = mean(d)
    total_precip_ZVAR[j] = std(d)

    d_wrap = reshape(d, 12, :)
    total_precip_ZMM[j, :]   = mean( d_wrap, dims=2 )[:, 1]
    total_precip_ZMVAR[j, :] = std(  d_wrap, dims=2 )[:, 1]

end

Dataset(output_file, "c") do ds

    defDim(ds, "months", 12)
    defDim(ds, "time", Inf)
    defDim(ds, "Nx", Nx)
    defDim(ds, "Ny", Ny)

    for (varname, vardata, vardim, attrib) in [
        ("total_precip_M",     total_precip_M,     ("Nx", "Ny",), Dict()),
        ("total_precip_VAR",   total_precip_VAR,   ("Nx", "Ny",), Dict()),
        ("total_precip_MM",    total_precip_MM,    ("Nx", "Ny", "months"), Dict()),
        ("total_precip_MVAR",  total_precip_MVAR,  ("Nx", "Ny", "months"), Dict()),
        ("total_precip_ZM",    total_precip_ZM,    ("Ny",), Dict()),
        ("total_precip_ZVAR",  total_precip_ZVAR,  ("Ny",), Dict()),
        ("total_precip_ZMM",   total_precip_ZMM,   ("Ny", "months"), Dict()),
        ("total_precip_ZMVAR", total_precip_ZMVAR, ("Ny", "months"), Dict()),
    ]

        if ! haskey(ds, varname)
            var = defVar(ds, varname, Float64, vardim)
            var.attrib["_FillValue"] = 1e20
        end

        var = ds[varname]

        for (k, v) in attrib
            var.attrib[k] = v
        end

        rng = []
        for i in 1:length(vardim)-1
            push!(rng, Colon())
        end
        push!(rng, 1:size(vardata)[end])
        var[rng...] = vardata

    end
end

max_stars_repo_name: meteorologytoday/CESM-diagnostic

max_issues_repo_name: meteorologytoday/CESM-diagnostic

max_forks_repo_name: meteorologytoday/CESM-diagnostic


text:
module BrainImagingDataStructurePaths
using FilePaths

include("ScannerHardware.jl")
using .ScannerHardware

include("SequenceSpecifics.jl")
using .SequenceSpecifics

include("InPlaneSpatielEncoding.jl")
using .InPlaneSpatielEncoding

include("TimingParameters.jl")
using .TimingParameters

abstract type BrainImagingDataStructurePath <: AbstractPath end
const BIDSPath = BrainImagingDataStructurePath

abstract type BIDSDataType <: BIDSPath end

"""
    participant(::BIDSDataType) -> "sub-<participant>"
"""
participant(bidsp::BIDSDataType) = string("sub-", bidsp.participant)

abstract type DiffusionWeightedImaging <: BIDSDataType end
const DWI = DiffusionWeightedImaging
datatypefolder(::DWI) = "dwi"

abstract type FieldInhomogeneityMapping <: BIDSDataType end
const FMAP = FieldInhomogeneityMapping
datatypefolder(::FMAP) = "fmap"

"""
    modalitylabel(x) -> AbstractString

Returns label associated with the `x`, where `x` is a either a subtype of
`NeuroModality` or `BrainImagingDataStructurePath`.
"""
function modalitylabel end

include("participant.jl")
include("session.jl")
include("anat.jl")
include("func.jl")
include("taskevents.jl")
include("dwi.jl")
include("fmap.jl")
include("behavioral.jl")

end


text:
using MultivariateMixtures: colwise_dot!, logsumexp, logsumexp!, stats!, covariances!

@testset "Utils" begin
    Random.seed!(42)

    a = [1 2 3; 4 5 6.]
    b = [4 5 6; 7 8 9.]
    r = zeros(size(a,2))

    colwise_dot!(r, a, b)
    @testset for i in 1:length(r)
        @test r[i]  == a[:,i] ⋅ b[:,i]
    end

    c = a'/10
    d = copy(c)
    logsumexp!(r, d)
    q = mapslices(r-> log(sum(exp, r)), c, dims=2)
    @testset for (v,w) in  zip(r,q)
        @test v ≈ w atol=1e-15
    end
    @testset for (v,w) in  zip(eachrow(c),eachrow(d))
        @test  exp.(v .- maximum(v)) == w
    end

    s = big(2e-20) - big(1e-40)
    d = [1e-20 log(1e-20)]
    r = zeros(1)
    logsumexp!(r, d)
    @test r[1] ≈ s atol=1e-35
    @test d ≈ [1.0 1e-20] atol=1e-35

    d = [1e-20 log(1e-20)]
    r = logsumexp(d)
    @test r[1] ≈ s atol=1e-35

    r[1] = 0.0
    logsumexp!(r, d, similar(r))
    @test r[1] ≈ s atol=1e-35
    @test d ≈ [1.0 1e-20] atol=1e-35

    d = 2
    s = 100
    k = 3
    n = k*s
    muori = [0.0 0.0 5.0; -5.0 0.0 5.0]
    Cori = [5.0; 1.0; 0.05]
    R = vcat(repeat([1.0 0.0 0.0], s), repeat([0.0 1.0 0.0], s), repeat([0.0 0.0 1.0], s))
    X = vcat(Cori[1]*randn(s,d).+muori[:,1]', Cori[2]*randn(s,d).+muori[:,2]', Cori[3]*randn(s,d).+muori[:,3]')'

    nk, mu = zeros(k), zeros(d,k)
    stats!(nk, mu, X, R)
    @testset for i in nk
        @test i ≈ s
    end
    @testset for i in 1:k
        @test mu[:, i] ≈ muori[:, i] atol=1.0
    end

    C1 = zeros(d,d,k)
    covariances!(C1, nk, mu, X, R, FullNormal, aux=similar(X))
    @testset for i in 1:k
        @test sqrt(det(C1[1,1,i])) ≈ Cori[i] atol=0.3
    end

    C2 = zeros(d,1,k)
    covariances!(C2, nk, mu, X, R, DiagNormal)
    @testset for i in 1:k
        @test sqrt(det(C2[1,1,i])) ≈ Cori[i] atol=0.3
    end

    C3 = zeros(1,1,k)
    covariances!(C3, nk, mu, X, R, IsoNormal, aux=similar(mu))
    @testset for i in 1:k
        @test sqrt(C3[1,1,i]) ≈ Cori[i] atol=0.3
    end

    A,B,C = MultivariateMixtures.new_params(FullNormal, Float64, d, k)
    @test size(A) == (k,)
    @test size(B) == (d,k)
    @test size(C) == (d,d,k)
    A,B,C = MultivariateMixtures.new_params(DiagNormal, Float64, d, k)
    @test size(A) == (k,)
    @test size(B) == (d,k)
    @test size(C) == (d,1,k)
    A,B,C = MultivariateMixtures.new_params(IsoNormal, Float64, d, k)
    @test size(A) == (k,)
    @test size(B) == (d,k)
    @test size(C) == (1,1,k)
    A,B = MultivariateMixtures.auxiliary(FullNormal, X, k)
    @test size(A) == (n,)
    @test size(B) == (d,n)
    A,B = MultivariateMixtures.auxiliary(DiagNormal, X, k)
    @test size(A) == (n,)
    @test size(B) == (d,k)
    A,B = MultivariateMixtures.auxiliary(IsoNormal, X, k)
    @test size(A) == (n,)
    @test size(B) == (d,k)

    F1 = MultivariateMixtures.precision(FullNormal, C1[:,:,1])
    @test F1 isa LowerTriangular
    logpdf!(view(R,:,1), X, mu[:,1], F1)
    r = logpdf(MvNormal(mu[:,1], C1[:,:,1]), X)
    @test sum(view(R,:,1) .- r) ≈ 0.0 atol=1e-10

    F2 = MultivariateMixtures.precision(DiagNormal, C2[:,:,2])
    @test F2 isa Diagonal
    logpdf!(view(R,:,2), X, mu[:,2], F2)
    r = logpdf(MvNormal(mu[:,2], Diagonal(vec(C2[:,:,2]))), X)
    @test sum(view(R,:,2) .- r) ≈ 0.0 atol=1e-10

    F3 = MultivariateMixtures.precision(IsoNormal, C3[:,:,3])
    @test F3 isa UniformScaling
    logpdf!(view(R,:,3), X, mu[:,3], F3)
    r = logpdf(MvNormal(mu[:,3], sqrt(C3[:,:,3][])), X)
    @test sum(view(R,:,3) .- r) ≈ 0.0 atol=1e-8

end


text:
using Robotlib
using Robotlib.Calibration
using DSP # For filtfilt
using MAT

function demo_calibforce(filename = "")
    info("Running force calibration demo. This file uses the logfile gravitylog.mat and calculates the calibration parameters for a wrist mounted force sensor.")
    h          = 0.004032;
    dh         = DHYuMi()

    if filename != "" # Use this option if you have a textbased logfile
        pathopen   = filename
        pathsave   = "log.mat"
        data       = orcalog2mat(pathopen, pathsave)
        # data       = readmat(pathsave)
        ds         = 1
        q          = getData("robot_1.*posRawAbs", data, ds)
        q̇          = getData("robot_1.*velFlt", data, ds)
        τ          = getData("robot_1.*trqRaw", data, ds)
        f          = getData("force", data, ds)
        # f[:,[2,5]] *= -1

        abb2logical!(q)
        abb2logical!(q̇)
        abb2logical!(τ)
        q = q*dh.GR'
        q̇ = q̇*dh.GR'
        τ = τ*inv(dh.GR')

        q̈ = filtfilt(ones(50),[50.],centralDiff(q̇))

        # plot(abs([q̇, q̈]))

        lowAcc = all(abs(q̈) .< 3e-4,2)[:]
        everyN = falses(size(q,1))
        everyN[1:10:end] = true
        q      = q[lowAcc & everyN,:]
        q̇      = q̇[lowAcc & everyN,:]
        τ      = τ[lowAcc & everyN,:]
        f      = f[lowAcc & everyN,:]
    else # Use this option to use the provided demofile
        data = readmat("data/gravitylog.mat")
        q = data["q"]
        q̇ = data["qd"]
        τ = data["tau"]
        f = data["f"]
    end

    N = size(q,1)
    fkine, ikine, jacobian = get_kinematic_functions("yumi")
    T = Array(Float64,4,4,N)
    for i = 1:N
        T[:,:,i]  = fkine(q[i,:])
    end

    # plot_traj(T)
    info("Calibrating force sensor")
    Rf,m,offset     = Robotlib.Calibration.calibForce(T,f,0.2205,offset=true)
    err = cat(2,[Rf*f[i,1:3] + offset - T[1:3,1:3,i]'*[0, 0, -9.82m] for i = 1:N]...)'
    plot(f[:,1:3],lab="Force")
    plot!(err,l=:dash,lab="Error")
    println("Error: ", round(rms(err),4))
    @show Rf
    @show m
    @show offset
    info("Done")

    # The code below can be used to save the result to a .mat-file
    # using MAT
    # matwrite("Rf.mat", Dict(
    #   	"Rf" => Rf,
    #     "mf" => m,
    #     "offset" => offset
    #   ))

    Rf,m,offset
end


text:
using DashBootstrapComponents

row_content = [
    dbc_col(html_div("One of two columns"), width = 4),
    dbc_col(html_div("One of two columns"), width = 4),
];

row = html_div([
    dbc_row(row_content, justify = "start"),
    dbc_row(row_content, justify = "center"),
    dbc_row(row_content, justify = "end"),
    dbc_row(row_content, justify = "between"),
    dbc_row(row_content, justify = "around"),
    dbc_row(row_content, justify = "evenly"),
]);

max_issues_repo_name: sthagen/dash-bootstrap-components

max_forks_repo_name: sthagen/dash-bootstrap-components


text:
using BinaryProvider
using BinaryBuilder
using BinaryBuilder: preferred_runner
using Random, LibGit2, Libdl, Test, ObjectFile, SHA

# The platform we're running on
const platform = platform_key_abi()

# On windows, the `.exe` extension is very important
const exe_ext = Sys.iswindows() ? ".exe" : ""

# We are going to build/install libfoo a lot, so here's our function to make sure the
# library is working properly
function check_foo(fooifier_path = "fooifier$(exe_ext)",
                   libfoo_path = "libfoo.$(Libdl.dlext)")
    # We know that foo(a, b) returns 2*a^2 - b
    result = 2*2.2^2 - 1.1

    # Test that we can invoke fooifier
    @test !success(`$fooifier_path`)
    @test success(`$fooifier_path 1.5 2.0`)
    @test parse(Float64,readchomp(`$fooifier_path 2.2 1.1`)) ≈ result

    # Test that we can dlopen() libfoo and invoke it directly
    libfoo = Libdl.dlopen_e(libfoo_path)
    @test libfoo != C_NULL
    foo = Libdl.dlsym_e(libfoo, :foo)
    @test foo != C_NULL
    @test ccall(foo, Cdouble, (Cdouble, Cdouble), 2.2, 1.1) ≈ result
    Libdl.dlclose(libfoo)
end

libfoo_products(prefix) = [
    LibraryProduct(prefix, "libfoo", :libfoo)
    ExecutableProduct(prefix, "fooifier", :fooifier)
]
libfoo_make_script = """
make clean
make install
"""
libfoo_cmake_script = raw"""
mkdir build
cd build
cmake -DCMAKE_INSTALL_PREFIX=${prefix} -DCMAKE_TOOLCHAIN_FILE=/opt/${target}/${target}.toolchain ..
make install
"""

# Helper function to try something and panic if it doesn't work
do_try(f) = try
    f()
catch e
    bt = catch_backtrace()
    Base.display_error(stderr, e, bt)

    # If a do_try fails, panic
    Test.@test false
end

# Run all our tests
include("basic.jl")
include("building.jl")
include("auditing.jl")
include("wizard.jl")

# These are broken for now
#include("package_tests.jl")