text:
using Base: datatype_haspadding
@testset "Discretization" begin
  @testset "FIST" begin
    𝒫 = Chain(P2[(0,0),(1,0),(1,1),(2,1),(2,2),(1,2),(0,0)])
    @test TopOptDash.ears(𝒫) == [2,4,5]

    𝒫 = Chain(P2[(0,0),(1,0),(1,1),(2,1),(1,2),(0,0)])
    @test TopOptDash.ears(𝒫) == [2,4]

    𝒫 = Chain(P2[(0,0),(1,0),(1,1),(1,2),(0,0)])
    @test TopOptDash.ears(𝒫) == [2,4]

    𝒫 = Chain(P2[(0,0),(1,1),(1,2),(0,0)])
    @test TopOptDash.ears(𝒫) == []

    𝒫 = Chain(P2[(0.443339268495331, 0.283757618605357),
                 (0.497822414616971, 0.398142813114205),
                 (0.770343126156527, 0.201815462842808),
                 (0.761236456732531, 0.330085709922366),
                 (0.985658085510286, 0.221530395507904),
                 (0.877899962498139, 0.325516131702896),
                 (0.561404274882782, 0.540334008885703),
                 (0.949459768187313, 0.396227653478068),
                 (0.594962560615951, 0.584927547374551),
                 (0.324208409133154, 0.607290684450708),
                 (0.424085089823892, 0.493532112641353),
                 (0.209843417261654, 0.590030658255966),
                 (0.27993878548962, 0.525162463476181),
                 (0.385557753911967, 0.322338556632868),
                 (0.443339268495331, 0.283757618605357)])
    @test TopOptDash.ears(𝒫) == [1,3,5,6,8,10,12,14]

    points = P2[(0,0),(1,0),(1,1),(2,1),(2,2),(1,2),(0,0)]
    connec = connect.([(4,5,6),(3,4,6),(3,6,1),(1,2,3)], Triangle)
    target = SimpleMesh(points[1:end-1], connec)
    poly = PolyArea(points)
    mesh = discretize(poly, FIST(false))
    @test mesh == target
    @test Set(vertices(poly)) == Set(vertices(mesh))
    @test nelements(mesh) == length(vertices(mesh)) - 2
  end

  @testset "Miscellaneous" begin
    for method in [FIST(), Dehn1899()]
      triangle = Triangle(P2(0,0), P2(1,0), P2(0,1))
      mesh = discretize(triangle, method)
      @test vertices(mesh) == [P2(0,0), P2(1,0), P2(0,1)]
      @test collect(elements(mesh)) == [triangle]

      quadrangle = Quadrangle(P2(0,0), P2(1,0), P2(1,1), P2(0,1))
      mesh = discretize(quadrangle, method)
      elms = collect(elements(mesh))
      @test vertices(mesh) == [P2(0,0), P2(1,0), P2(1,1), P2(0,1)]
      @test eltype(elms) <: Triangle
      @test length(elms) == 2

      q = Quadrangle(P2(0,0), P2(1,0), P2(1,1), P2(0,1))
      t = Triangle(P2(1,0), P2(2,1), P2(1,1))
      m = Multi([q, t])
      mesh = discretize(m, method)
      elms = collect(elements(mesh))
      @test vertices(mesh) == [vertices(q); vertices(t)]
      @test vertices(elms[1]) ⊆ vertices(q)
      @test vertices(elms[2]) ⊆ vertices(q)
      @test vertices(elms[3]) ⊆ vertices(t)
      @test eltype(elms) <: Triangle
      @test length(elms) == 3

      outer = P2[(0,0),(1,0),(1,1),(0,1),(0,0)]
      hole1 = P2[(0.2,0.2),(0.4,0.2),(0.4,0.4),(0.2,0.4),(0.2,0.2)]
      hole2 = P2[(0.6,0.2),(0.8,0.2),(0.8,0.4),(0.6,0.4),(0.6,0.2)]
      poly  = PolyArea(outer, [hole1, hole2])
      chain, _ = bridge(poly, width=0.01)
      mesh  = discretize(chain, method)
      @test nvertices(mesh) == 16
      @test nelements(mesh) == 14
      @test all(t -> area(t) > zero(T), mesh)
    end
  end

  @testset "Difficult examples" begin
    for method in [FIST(), Dehn1899()]
      poly = readpoly(T, joinpath(datadir, "taubin.line"))
      mesh = discretize(poly, method)
      @test Set(vertices(poly)) == Set(vertices(mesh))
      @test nelements(mesh) == length(vertices(mesh)) - 2

      poly = readpoly(T, joinpath(datadir, "poly1.line"))
      mesh = discretize(poly, method)
      @test Set(vertices(poly)) == Set(vertices(mesh))
      @test nelements(mesh) == length(vertices(mesh)) - 2

      poly = readpoly(T, joinpath(datadir, "poly2.line"))
      mesh = discretize(poly, method)
      @test Set(vertices(poly)) == Set(vertices(mesh))
      @test nelements(mesh) == length(vertices(mesh)) - 2

      poly = readpoly(T, joinpath(datadir, "poly3.line"))
      mesh = discretize(poly, method)
      @test Set(vertices(poly)) == Set(vertices(mesh))
      @test nelements(mesh) == length(vertices(mesh)) - 2

      poly = readpoly(T, joinpath(datadir, "poly4.line"))
      mesh = discretize(poly, method)
      @test Set(vertices(poly)) == Set(vertices(mesh))
      @test nelements(mesh) == length(vertices(mesh)) - 2

      poly = readpoly(T, joinpath(datadir, "poly5.line"))
      mesh = discretize(poly, method)
      @test Set(vertices(poly)) == Set(vertices(mesh))
      @test nelements(mesh) == length(vertices(mesh)) - 2

      poly = readpoly(T, joinpath(datadir, "smooth1.line"))
      mesh = discretize(poly, method)
      @test Set(vertices(poly)) == Set(vertices(mesh))
      @test nelements(mesh) == length(vertices(mesh)) - 2

      poly = readpoly(T, joinpath(datadir, "smooth2.line"))
      mesh = discretize(poly, method)
      @test Set(vertices(poly)) == Set(vertices(mesh))
      @test nelements(mesh) == length(vertices(mesh)) - 2

      poly = readpoly(T, joinpath(datadir, "smooth3.line"))
      mesh = discretize(poly, method)
      @test Set(vertices(poly)) == Set(vertices(mesh))
      @test nelements(mesh) == length(vertices(mesh)) - 2

      poly = readpoly(T, joinpath(datadir, "smooth4.line"))
      mesh = discretize(poly, method)
      @test Set(vertices(poly)) == Set(vertices(mesh))
      @test nelements(mesh) == length(vertices(mesh)) - 2

      poly = readpoly(T, joinpath(datadir, "smooth5.line"))
      mesh = discretize(poly, method)
      @test Set(vertices(poly)) == Set(vertices(mesh))
      @test nelements(mesh) == length(vertices(mesh)) - 2

      poly = readpoly(T, joinpath(datadir, "hole1.line"))
      mesh = discretize(poly, method)
      @test Set(vertices(poly)) == Set(vertices(mesh))
      @test nelements(mesh) == 32

      poly = readpoly(T, joinpath(datadir, "hole2.line"))
      mesh = discretize(poly, method)
      @test Set(vertices(poly)) == Set(vertices(mesh))
      @test nelements(mesh) == 30

      poly = readpoly(T, joinpath(datadir, "hole3.line"))
      mesh = discretize(poly, method)
      @test Set(vertices(poly)) == Set(vertices(mesh))
      @test nelements(mesh) == 32

      poly = readpoly(T, joinpath(datadir, "hole4.line"))
      mesh = discretize(poly, method)
      @test Set(vertices(poly)) == Set(vertices(mesh))
      @test nelements(mesh) == 30

      poly = readpoly(T, joinpath(datadir, "hole5.line"))
      mesh = discretize(poly, method)
      @test Set(vertices(poly)) == Set(vertices(mesh))
      @test nelements(mesh) == 32
    end
  end
end


text:
using MLUtils: obsview

@testset "BatchView" begin

    @testset "constructor" begin
        @test_throws DimensionMismatch BatchView((rand(2,10),rand(9)))
        @test_throws DimensionMismatch BatchView((rand(2,10),rand(9)))
        @test_throws DimensionMismatch BatchView((rand(2,10),rand(4,9,10),rand(9)))
        @test_throws MethodError BatchView(EmptyType())
        for var in (vars..., tuples..., Xs, ys)
            @test_throws MethodError BatchView(var...)
            @test_throws MethodError BatchView(var, 16)

            A = BatchView(var, batchsize=3)
            @test length(A) == 5
            @test batchsize(A) == 3
            @test numobs(A) == length(A)
            @test @inferred(parent(A)) === var
        end
        A = BatchView(X, batchsize=16)
        @test length(A) == 1
        @test batchsize(A) == 15
    end

    @testset "typestability" begin
        for var in (vars..., tuples..., Xs, ys)
            @test typeof(@inferred(BatchView(var))) <: BatchView
            @test typeof(@inferred(BatchView(var, batchsize=3))) <: BatchView
            @test typeof(@inferred(BatchView(var, batchsize=3, partial=true))) <: BatchView
            @test typeof(@inferred(BatchView(var, batchsize=3, partial=false))) <: BatchView
        end
        @test typeof(@inferred(BatchView(CustomType()))) <: BatchView
    end

    @testset "AbstractArray interface" begin
        for var in (vars..., tuples..., Xs, ys)
            A = BatchView(var, batchsize=5)
            @test_throws BoundsError A[-1]
            @test_throws BoundsError A[4]
            @test @inferred(numobs(A)) == 3
            @test @inferred(length(A)) == 3
            @test @inferred(batchsize(A)) == 5
            @test @inferred(size(A)) == (3,)
            @test @inferred(getobs(A[2:3])) == getobs(BatchView(ObsView(var, 6:15), batchsize=5))
            @test @inferred(getobs(A[[1,3]])) == getobs(BatchView(ObsView(var, [1:5..., 11:15...]), batchsize=5))
            @test @inferred(A[1]) == obsview(var, 1:5)
            @test @inferred(A[2]) == obsview(var, 6:10)
            @test @inferred(A[3]) == obsview(var, 11:15)
            @test A[end] == A[3]
            @test @inferred(getobs(A,1)) == getobs(var, 1:5)
            @test @inferred(getobs(A,2)) == getobs(var, 6:10)
            @test @inferred(getobs(A,3)) == getobs(var, 11:15)
            @test typeof(@inferred(collect(A))) <: Vector
        end
        for var in (vars..., tuples...)
            A = BatchView(var, batchsize=5)
            @test @inferred(getobs(A)) == var
            @test A[2:3] == obsview(var, [6:15;])
            @test A[[1,3]] == obsview(var, [1:5..., 11:15...])
        end
    end

    @testset "subsetting" begin
        # for var in (vars..., tuples..., Xs, ys)
            for var in (vars[1], )
            A = BatchView(var, batchsize=3)
            @test getobs(@inferred(ObsView(A))) == @inferred(getobs(A))
            @test_throws BoundsError ObsView(A,1:6)

            S = @inferred(ObsView(A, 1:2))
            @test typeof(S) <: ObsView
            @test @inferred(numobs(S)) == 2
            @test @inferred(length(S)) == numobs(S)
            @test @inferred(size(S)) == (length(S),)
            @test getobs(@inferred(A[1:2])) == getobs(S)
            @test @inferred(getobs(A,1:2)) == getobs(S)
            @test @inferred(getobs(S)) == getobs(BatchView(ObsView(var,1:6),batchsize=3))

            S = @inferred(ObsView(A, 1:2))
            @test typeof(S) <: ObsView
        end
        A = BatchView(X, batchsize=3)
        @test typeof(A.data) <: Array
        S = @inferred(obsview(A))
        S === A
    end

    # @testset "nesting with ObsView" begin
    #     for var in vars
    #         @test eltype(@inferred(BatchView(ObsView(var)))[1]) <: Union{SubArray,String}
    #     end
    #     for var in tuples
    #         @test eltype(@inferred(BatchView(ObsView(var)))[1]) <: Tuple
    #     end
    #     for var in (Xs, ys)
    #         @test eltype(@inferred(BatchView(ObsView(var)))[1]) <: SubArray
    #     end
    # end
end

ext: jl


text:
# include("breast_cancer.jl")
include("simulation.jl")
# include("rui_example.jl")


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

"""
Tools for collecting and manipulating stack traces. Mainly used for building errors.
"""
module StackTraces

import Base: hash, ==, show
using Base.Printf: @printf
using Base: coalesce

export StackTrace, StackFrame, stacktrace

"""
    StackFrame

Stack information representing execution context, with the following fields:

- `func::Symbol`

  The name of the function containing the execution context.

- `linfo::Union{Core.MethodInstance, CodeInfo, Nothing}`

  The MethodInstance containing the execution context (if it could be found).

- `file::Symbol`

  The path to the file containing the execution context.

- `line::Int`

  The line number in the file containing the execution context.

- `from_c::Bool`

  True if the code is from C.

- `inlined::Bool`

  True if the code is from an inlined frame.

- `pointer::UInt64`

  Representation of the pointer to the execution context as returned by `backtrace`.

"""
struct StackFrame # this type should be kept platform-agnostic so that profiles can be dumped on one machine and read on another
    "the name of the function containing the execution context"
    func::Symbol
    "the path to the file containing the execution context"
    file::Symbol
    "the line number in the file containing the execution context"
    line::Int
    "the MethodInstance or CodeInfo containing the execution context (if it could be found)"
    linfo::Union{Core.MethodInstance, Core.CodeInfo, Nothing}
    "true if the code is from C"
    from_c::Bool
    "true if the code is from an inlined frame"
    inlined::Bool
    "representation of the pointer to the execution context as returned by `backtrace`"
    pointer::UInt64  # Large enough to be read losslessly on 32- and 64-bit machines.
end

StackFrame(func, file, line) = StackFrame(Symbol(func), Symbol(file), line,
                                          nothing, false, false, 0)

"""
    StackTrace

An alias for `Vector{StackFrame}` provided for convenience; returned by calls to
`stacktrace`.
"""
const StackTrace = Vector{StackFrame}

const empty_sym = Symbol("")
const UNKNOWN = StackFrame(empty_sym, empty_sym, -1, nothing, true, false, 0) # === lookup(C_NULL)

#=
If the StackFrame has function and line information, we consider two of them the same if
they share the same function/line information.
=#
function ==(a::StackFrame, b::StackFrame)
    a.line == b.line && a.from_c == b.from_c && a.func == b.func && a.file == b.file && a.inlined == b.inlined
end

function hash(frame::StackFrame, h::UInt)
    h += 0xf4fbda67fe20ce88 % UInt
    h = hash(frame.line, h)
    h = hash(frame.file, h)
    h = hash(frame.func, h)
    h = hash(frame.from_c, h)
    h = hash(frame.inlined, h)
end

"""
    lookup(pointer::Union{Ptr{Cvoid}, UInt}) -> Vector{StackFrame}

Given a pointer to an execution context (usually generated by a call to `backtrace`), looks
up stack frame context information. Returns an array of frame information for all functions
inlined at that point, innermost function first.
"""
function lookup(pointer::Ptr{Cvoid})
    infos = ccall(:jl_lookup_code_address, Any, (Ptr{Cvoid}, Cint), pointer - 1, false)
    isempty(infos) && return [StackFrame(empty_sym, empty_sym, -1, nothing, true, false, convert(UInt64, pointer))]
    res = Vector{StackFrame}(uninitialized, length(infos))
    for i in 1:length(infos)
        info = infos[i]
        @assert(length(info) == 7)
        res[i] = StackFrame(info[1], info[2], info[3], info[4], info[5], info[6], info[7])
    end
    return res
end

lookup(pointer::UInt) = lookup(convert(Ptr{Cvoid}, pointer))

const top_level_scope_sym = Symbol("top-level scope")

using Base.Meta
is_loc_meta(expr, kind) = isexpr(expr, :meta) && length(expr.args) >= 1 && expr.args[1] === kind
function lookup(ip::Base.InterpreterIP)
    i = ip.stmt
    foundline = false
    if ip.code isa Core.MethodInstance
        codeinfo = ip.code.inferred
        func = ip.code.def.name
        file = ip.code.def.file
        line = ip.code.def.line
    elseif ip.code === nothing
        # interpreted top-level expression with no CodeInfo
        return [StackFrame(top_level_scope_sym, empty_sym, 0, nothing, false, false, 0)]
    else
        @assert ip.code isa Core.CodeInfo
        codeinfo = ip.code
        func = top_level_scope_sym
        file = empty_sym
        line = 0
    end
    while i >= 1
        expr = codeinfo.code[i]
        if isa(expr, LineNumberNode)
            if line == 0
                line = expr.line
            end
            if expr.file !== nothing
                file = expr.file
            end
            foundline = true
        elseif foundline && is_loc_meta(expr, :push_loc)
            file = expr.args[2]
            if length(expr.args) >= 3
                func = expr.args[3]
            else
                # Note: This is not quite correct. See issue #23971
                func = Symbol("macro expansion")
            end
            scopes = lookup(Base.InterpreterIP(ip.code, i-1))
            pushfirst!(scopes, StackFrame(
                func, file, line, nothing, false, true, 0
            ))
            return scopes
        elseif is_loc_meta(expr, :pop_loc)
            npops = 1
            while npops >= 1
                i -= 1
                expr = codeinfo.code[i]
                is_loc_meta(expr, :pop_loc) && (npops += 1)
                is_loc_meta(expr, :push_loc) && (npops -= 1)
            end
        end
        i -= 1
    end
    return [StackFrame(func, file, line, ip.code, false, false, 0)]
end

# allow lookup on already-looked-up data for easier handling of pre-processed frames
lookup(s::StackFrame) = StackFrame[s]
lookup(s::Tuple{StackFrame,Int}) = StackFrame[s[1]]

"""
    backtrace()

Get a backtrace object for the current program point.
"""
function Base.backtrace()
    bt, bt2 = ccall(:jl_backtrace_from_here, Any, (Int32,), false)
    if length(bt) > 2
        # remove frames for jl_backtrace_from_here and backtrace()
        if bt[2] == Ptr{Cvoid}(-1%UInt)
            # backtrace() is interpreted
            # Note: win32 is missing the top frame (see https://bugs.chromium.org/p/crashpad/issues/detail?id=53)
            @static if Base.Sys.iswindows() && Int === Int32
                deleteat!(bt, 1:2)
            else
                deleteat!(bt, 1:3)
            end
            pushfirst!(bt2)
        else
            @static if Base.Sys.iswindows() && Int === Int32
                deleteat!(bt, 1)
            else
                deleteat!(bt, 1:2)
            end
        end
    end
    return Base._reformat_bt(bt, bt2)
end

"""
    stacktrace([trace::Vector{Ptr{Cvoid}},] [c_funcs::Bool=false]) -> StackTrace

Returns a stack trace in the form of a vector of `StackFrame`s. (By default stacktrace
doesn't return C functions, but this can be enabled.) When called without specifying a
trace, `stacktrace` first calls `backtrace`.
"""
function stacktrace(trace::Vector{<:Union{Base.InterpreterIP,Ptr{Cvoid}}}, c_funcs::Bool=false)
    stack = vcat(StackTrace(), map(lookup, trace)...)::StackTrace

    # Remove frames that come from C calls.
    if !c_funcs
        filter!(frame -> !frame.from_c, stack)
    end

    # Remove frame for this function (and any functions called by this function).
    remove_frames!(stack, :stacktrace)

    # is there a better way?  the func symbol has a number suffix which changes.
    # it's possible that no test is needed and we could just popfirst! all the time.
    # this line was added to PR #16213 because otherwise stacktrace() != stacktrace(false).
    # not sure why.  possibly b/c of re-ordering of base/sysimg.jl
    !isempty(stack) && startswith(string(stack[1].func),"jlcall_stacktrace") && popfirst!(stack)
    stack
end

stacktrace(c_funcs::Bool=false) = stacktrace(backtrace(), c_funcs)

"""
    remove_frames!(stack::StackTrace, name::Symbol)

Takes a `StackTrace` (a vector of `StackFrames`) and a function name (a `Symbol`) and
removes the `StackFrame` specified by the function name from the `StackTrace` (also removing
all frames above the specified function). Primarily used to remove `StackTraces` functions
from the `StackTrace` prior to returning it.
"""
function remove_frames!(stack::StackTrace, name::Symbol)
    splice!(stack, 1:coalesce(findlast(frame -> frame.func == name, stack), 0))
    return stack
end

function remove_frames!(stack::StackTrace, names::Vector{Symbol})
    splice!(stack, 1:coalesce(findlast(frame -> frame.func in names, stack), 0))
    return stack
end

"""
    remove_frames!(stack::StackTrace, m::Module)

Returns the `StackTrace` with all `StackFrame`s from the provided `Module` removed.
"""
function remove_frames!(stack::StackTrace, m::Module)
    filter!(f -> !from(f, m), stack)
    return stack
end

is_top_level_frame(f::StackFrame) = f.linfo isa Core.CodeInfo || (f.linfo === nothing && f.func === top_level_scope_sym)

function show_spec_linfo(io::IO, frame::StackFrame)
    if frame.linfo == nothing
        if frame.func === empty_sym
            @printf(io, "ip:%#x", frame.pointer)
        elseif frame.func === top_level_scope_sym
            print(io, "top-level scope")
        else
            color = get(io, :color, false) && get(io, :backtrace, false) ?
                        Base.stackframe_function_color() :
                        :nothing
            printstyled(io, string(frame.func), color=color)
        end
    elseif frame.linfo isa Core.MethodInstance
        if isa(frame.linfo.def, Method)
            Base.show_tuple_as_call(io, frame.linfo.def.name, frame.linfo.specTypes)
        else
            Base.show(io, frame.linfo)
        end
    elseif frame.linfo isa Core.CodeInfo
        print(io, "top-level scope")
    end
end

function show(io::IO, frame::StackFrame; full_path::Bool=false)
    show_spec_linfo(io, frame)
    if frame.file !== empty_sym
        file_info = full_path ? string(frame.file) : basename(string(frame.file))
        print(io, " at ")
        Base.with_output_color(get(io, :color, false) && get(io, :backtrace, false) ? Base.stackframe_lineinfo_color() : :nothing, io) do io
            print(io, file_info, ":")
            if frame.line >= 0
                print(io, frame.line)
            else
                print(io, "?")
            end
        end
    end
    if frame.inlined
        print(io, " [inlined]")
    end
end

"""
    from(frame::StackFrame, filter_mod::Module) -> Bool

Returns whether the `frame` is from the provided `Module`
"""
function from(frame::StackFrame, m::Module)
    finfo = frame.linfo
    result = false

    if finfo isa Core.MethodInstance
        frame_m = finfo.def
        isa(frame_m, Method) && (frame_m = frame_m.module)
        result = nameof(frame_m) === nameof(m)
    end

    return result
end

end

max_stars_repo_name: blegat/julia

max_issues_repo_name: blegat/julia

max_forks_repo_name: blegat/julia


text:
abstract type CrossValGenerator end

struct LOOCV <: CrossValGenerator
    n::Int
end

length(c::LOOCV) = c.n

function iterate(c::LOOCV, s::Int=1)
    (s > c.n) && return nothing
    return (leave_one_out(c.n, s), s + 1)
end

function leave_one_out(n::Int, i::Int)
    @assert 1 <= i <= n
    x = Array{Int}(undef, n - 1)
    for j = 1:i-1
        x[j] = j
    end
    for j = i+1:n
        x[j-1] = j
    end
    return x
end

function split_database(database::Array{Int,1}, k::Int)
    n = size(database, 1)
    [database[(i-1)*(ceil(Int,n/k))+1:min(i*(ceil(Int,n/k)),n)] for i in 1:k]
end

function kfoldcrossvalidation!(
    previousresult::ModelSelection.ModelSelectionData,
    data::ModelSelection.ModelSelectionData,
    k::Int,
    s::Float64)
    kfoldcrossvalidation(previousresult, data, k, s)
end

function kfoldcrossvalidation(
    previousresult::ModelSelection.ModelSelectionData,
    data::ModelSelection.ModelSelectionData,
    k::Int,
    s::Float64)

    #db = randperm(data.nobs)
    db = collect(1:data.nobs)
    folds = split_database(db, k)

    # if data.time != nothing
    #     if data.panel != nothing
    #         # time & panel -> vemos que pasa acá
    #         folds = []
    #     else
    #         # time -> divisiones sin permutación
    #         folds = []
    #     end
    # else
    #     folds = []
    # end

    bestmodels = []

    for obs in LOOCV(k)
        dataset = collect(Iterators.flatten(folds[obs]))
        testset = setdiff(1:data.nobs, dataset)

        reduced = ModelSelection.copy_data(data)
        reduced.depvar_data = data.depvar_data[dataset]
        reduced.expvars_data = data.expvars_data[dataset, :]

        if reduced.time != nothing
            reduced.time_data = data.time_data[dataset]
        end

        if reduced.panel != nothing
            reduced.panel_data = data.panel_data[dataset]
        end

        reduced.nobs = size(dataset, 1)
        _, vars = ModelSelection.PreliminarySelection.lasso!(reduced, addextrasflag=false)

        backup = ModelSelection.copy_data(data)
        backup.expvars = data.expvars[vars]
        backup.expvars_data = data.expvars_data[:,vars]

        ModelSelection.AllSubsetRegression.ols!(backup,
            outsample = testset,
            criteria = [ :rmseout ],
            ttest = previousresult.results[1].ttest,
            residualtest = previousresult.results[1].residualtest
        )

        push!(bestmodels, Dict(
            :data => backup.results[1].bestresult_data,
            :datanames => backup.results[1].datanames
        ))
    end

    datanames = unique(Iterators.flatten(model[:datanames] for model in bestmodels))

    data = Array{Any,2}(zeros(size(bestmodels,1), size(datanames,1)))

    for (i, model) in enumerate(bestmodels)
        for (f, col) in enumerate(model[:datanames])
            pos = ModelSelection.get_column_index(col, datanames)
            data[i, pos] = model[:data][f]
        end
    end

    replace!(data, NaN => 0)

    average_data = mean(data, dims=1)
    median_data = median(data, dims=1)

    result = CrossValidationResult(k, 0, previousresult.results[1].ttest, datanames, average_data, median_data, data)

    previousresult = ModelSelection.addresult!(previousresult, result)

    addextras(previousresult, result)

    return previousresult
end

function to_string(data::ModelSelection.ModelSelectionData, result::CrossValidationResult)
    datanames_index = ModelSelection.create_datanames_index(result.datanames)

    out = ""
    out *= @sprintf("\n")
    out *= @sprintf("══════════════════════════════════════════════════════════════════════════════\n")
    out *= @sprintf("                       Cross validation average results                       \n")
    out *= @sprintf("══════════════════════════════════════════════════════════════════════════════\n")
    out *= @sprintf("                                                                              \n")
    out *= @sprintf("                                     Dependent variable: %s                   \n", data.depvar)
    out *= @sprintf("                                     ─────────────────────────────────────────\n")
    out *= @sprintf("                                                                              \n")
    out *= @sprintf(" Selected covariates                 Coef.")
    if result.ttest
        out *= @sprintf("        Std.         t-test")
    end
    out *= @sprintf("\n")
    out *= @sprintf("──────────────────────────────────────────────────────────────────────────────\n")

    for varname in data.expvars
        if Symbol(string(varname, "_b")) in keys(datanames_index)
            out *= @sprintf(" %-35s", varname)
            out *= @sprintf(" %-10f", result.average_data[datanames_index[Symbol(string(varname, "_b"))]])
            if result.ttest
                out *= @sprintf("   %-10f", result.average_data[datanames_index[Symbol(string(varname, "_bstd"))]])
                out *= @sprintf("   %-10f", result.average_data[datanames_index[Symbol(string(varname, "_t"))]])
            end
            out *= @sprintf("\n")
        end
    end

    out *= @sprintf("──────────────────────────────────────────────────────────────────────────────\n")
    out *= @sprintf(" Observations                        %-10d\n", result.average_data[datanames_index[:nobs]])
    out *= @sprintf(" RMSE OUT                            %-10f\n", result.average_data[datanames_index[:rmseout]])

    out *= @sprintf("\n")
    out *= @sprintf("\n")
    out *= @sprintf("══════════════════════════════════════════════════════════════════════════════\n")
    out *= @sprintf("                       Cross validation median results                        \n")
    out *= @sprintf("══════════════════════════════════════════════════════════════════════════════\n")
    out *= @sprintf("                                                                              \n")
    out *= @sprintf("                                     Dependent variable: %s                   \n", data.depvar)
    out *= @sprintf("                                     ─────────────────────────────────────────\n")
    out *= @sprintf("                                                                              \n")
    out *= @sprintf(" Covariates                          Coef.")
    if result.ttest
        out *= @sprintf("        Std.         t-test")
    end
    out *= @sprintf("\n")
    out *= @sprintf("──────────────────────────────────────────────────────────────────────────────\n")

    for varname in data.expvars
        if Symbol(string(varname, "_b")) in keys(datanames_index)
            out *= @sprintf(" %-35s", varname)
            out *= @sprintf(" %-10f", result.median_data[datanames_index[Symbol(string(varname, "_b"))]])
            if result.ttest
                out *= @sprintf("   %-10f", result.median_data[datanames_index[Symbol(string(varname, "_bstd"))]])
                out *= @sprintf("   %-10f", result.median_data[datanames_index[Symbol(string(varname, "_t"))]])
            end
            out *= @sprintf("\n")
        end
    end

    out *= @sprintf("\n")
    out *= @sprintf("──────────────────────────────────────────────────────────────────────────────\n")
    out *= @sprintf(" Observations                        %-10d\n", result.median_data[datanames_index[:nobs]])
    out *= @sprintf(" RMSE OUT                            %-10f\n", result.median_data[datanames_index[:rmseout]])
    out *= @sprintf("──────────────────────────────────────────────────────────────────────────────\n")

    return out
end


text:
using Revise; using DynamicalBilliards, BenchmarkTools

bt = billiard_sinai()
p = MagneticParticle(0.1,0.1, π/8, 0.5)

@btime next_collision($p, $bt)

@btime bounce!($p, $bt)

# before sincos():
julia> @btime bounce!($p, $bt)
  278.529 ns (0 allocations: 0 bytes)
(1, 0.44691013238163385, [0.707091, 0.640047], [-0.353179, 0.935556])

# after sincos():
julia> @btime bounce!($p, $bt)
  57.476 ns (0 allocations: 0 bytes)
(0, Inf, [0.415228, 1.10629], [0.957396, -0.288778])


text:
@inline static(::Val{N}) where {N} = StaticInt{N}()
@inline static(::Nothing) = nothing
@generated function static_sizeof(::Type{T}) where {T}
    st = Base.allocatedinline(T) ? sizeof(T) : sizeof(Int)
    Expr(:block, Expr(:meta, :inline), Expr(:call, Expr(:curly, :StaticInt, st)))
end

# These have versions that may allow for more optimizations, so we override base methods with a single `StaticInt` argument.
for (f,ff) ∈ [
  # (:(Base.:+),:vadd_fast), (:(Base.:-),:vsub_fast), (:(Base.:*),:vmul_fast),
  (:(Base.:+),:vadd_nsw), (:(Base.:-),:vsub_nsw), (:(Base.:*),:vmul_nsw),
  (:(Base.:+),:vadd_nuw), (:(Base.:-),:vsub_nuw), (:(Base.:*),:vmul_nuw),
  (:(Base.:+),:vadd_nw), (:(Base.:-),:vsub_nw), (:(Base.:*),:vmul_nw),
  # (:(Base.:<<),:vshl), (:(Base.:÷),:(Base.div)), (:(Base.:%), :(Base.rem)), (:(Base.:>>>),:(Base.:>>>))
]
  @eval begin
    # @inline $f(::StaticInt{M}, ::StaticInt{N}) where {M, N} = StaticInt{$f(M, N)}()
    # If `M` and `N` are known at compile time, there's no need to add nsw/nuw flags.
    @inline $ff(::StaticInt{M}, ::StaticInt{N}) where {M, N} = $f(StaticInt{M}(),StaticInt{N}())
    # @inline $f(::StaticInt{M}, x) where {M} = $ff(M, x)
    # @inline $f(x, ::StaticInt{M}) where {M} = $ff(x, M)
    @inline $ff(::StaticInt{M}, x::T) where {M,T<:IntegerTypesHW} = $ff(M%T, x)
    @inline $ff(x::T, ::StaticInt{M}) where {M,T<:IntegerTypesHW} = $ff(x, M%T)
    @inline $ff(::StaticInt{M}, x) where {M} = $ff(M, x)
    @inline $ff(x, ::StaticInt{M}) where {M} = $ff(x, M)
  end
end
for f ∈ [:vsub_nsw, :vsub_nuw, :vsub_nw]
  @eval begin
    @inline $f(::Zero, m::Number) = -m
    @inline $f(::Zero, m::IntegerTypesHW) = -m
    @inline $f(m::Number, ::Zero) =  m
    @inline $f(m::IntegerTypesHW, ::Zero) =  m
    @inline $f(::Zero, ::Zero) = Zero()
    @inline $f(::Zero, ::StaticInt{N}) where {N} = -StaticInt{N}()
    @inline $f(::StaticInt{N}, ::Zero) where {N} = StaticInt{N}()
  end
end
for f ∈ [:vadd_nsw, :vadd_nuw, :vadd_nw]
  @eval begin
    @inline $f(::StaticInt{N}, ::Zero) where {N} = StaticInt{N}()
    @inline $f(::Zero, ::StaticInt{N}) where {N} = StaticInt{N}()
    @inline $f(::Zero, ::Zero) = Zero()
    @inline $f(a::Number, ::Zero) = a
    @inline $f(a::IntegerTypesHW, ::Zero) = a
    @inline $f(::Zero, a::Number) = a
    @inline $f(::Zero, a::IntegerTypesHW) = a
  end
end

@inline Base.:(+)(m::AbstractMask{W}, ::StaticInt{N}) where {N,W} = m + vbroadcast(Val{W}(), N)
@inline Base.:(+)(::StaticInt{N}, m::AbstractMask{W}) where {N,W} = vbroadcast(Val{W}(), N) + m

for V ∈ [:AbstractSIMD, :MM, :VecUnroll, :AbstractMask]
  for (mod,m,a,s) ∈ ((:FastMath, :mul_fast, :add_fast, :sub_fast),(:Base,:(*),:(+),:(-)))
    @eval begin
      @inline $mod.$m(::Zero, x::$V) = Zero()
      @inline $mod.$m(::One, x::$V) = x
      @inline $mod.$m(x::$V, ::Zero) = Zero()
      @inline $mod.$m(x::$V, ::One) = x

      @inline $mod.$a(::Zero, x::$V) = x
      @inline $mod.$a(x::$V, ::Zero) = x

      @inline $mod.$s(::Zero, x::$V) = -x
      @inline $mod.$s(x::$V, ::Zero) =  x
    end
  end
end

@inline vload(::StaticInt{N}, args...) where {N} = StaticInt{N}()
@inline stridedpointer(::StaticInt{N}) where {N} = StaticInt{N}()
@inline zero_offsets(::StaticInt{N}) where {N} = StaticInt{N}()


text:
include("llvmtest.jl")
include("cxxstr.jl")
include("icxxstr.jl")
include("llvmgraph.jl")


text:
# TODO: Check.
function logprior(state::State, prior::Prior)
  lp = sum(logpdf.(prior.mu, state.mu))
  lp += sum(logpdf.(prior.omega, state.omega))
  lp += sum(logpdf.(prior.psi, state.psi))
  lp += sum(logpdf.(prior.nu, state.nu))
  lp += logpdf(prior.gamma, state.gammaC)
  lp += logpdf(prior.gamma, state.gammaT)
  lp += logpdf(prior.eta, state.etaC)
  lp += logpdf(prior.eta, state.etaT)

  # lambda prior
  etaT_star = state.beta ? state.etaT : state.etaC
  lp += sum(log.(state.etaC[state.lambdaC]))
  lp += sum(log.(etaT_star[state.lambdaT]))

  # v prior
  nuC = state.nu[state.lambdaC]
  nuT = state.nu[state.lambdaT]
  lp += sum(gammalogpdf.(nuC/2, 2 ./ nuC, state.vC))
  lp += sum(gammalogpdf.(nuT/2, 2 ./ nuT, state.vT))

  # zeta prior
  lp += sum(logpdf(truncated(Normal(0, 1/sqrt(state.vC[n])), 0, Inf), state.zetaC[n])
            for n in eachindex(state.zetaC))
  lp += sum(logpdf(truncated(Normal(0, 1/sqrt(state.vT[n])), 0, Inf), state.zetaT[n])
            for n in eachindex(state.zetaT))
  return lp
end

function marginal_logprior(state::State, prior::Prior)
  lp = sum(logpdf.(prior.mu, state.mu))
  lp += sum(logpdf.(prior.omega, state.omega))
  lp += sum(logpdf.(prior.psi, state.psi))
  lp += sum(logpdf.(prior.nu, state.nu))
  lp += logpdf(prior.gamma, state.gammaC)
  lp += logpdf(prior.eta, state.etaC)
  if state.beta
    lp += logpdf(prior.gamma, state.gammaT)
    lp += logpdf(prior.eta, state.etaT)
  end
  return lp
end

max_stars_repo_name: luiarthur/CytofDensityEstimation

max_issues_repo_name: luiarthur/CytofDensityEstimation

max_forks_repo_name: luiarthur/CytofDensityEstimation


text:
# # Heat equation tutorial

# In this tutorial, we'll be solving the [heat
# equation](https://en.wikipedia.org/wiki/Heat_equation):

# ``
# \frac{∂ ρcT}{∂ t} + ∇ ⋅ (-α ∇ρcT) = 0
# ``

# where
#  - `t` is time
#  - `α` is the thermal diffusivity
#  - `T` is the temperature
#  - `ρ` is the density
#  - `c` is the heat capacity
#  - `ρcT` is the thermal energy

# To put this in the form of ClimateMachine's [`BalanceLaw`](@ref
# ClimateMachine.BalanceLaws.BalanceLaw), we'll re-write the equation as:

# ``
# \frac{∂ ρcT}{∂ t} + ∇ ⋅ (F(α, ρcT, t)) = 0
# ``

# where
#  - ``F(α, ρcT, t) = -α ∇ρcT`` is the second-order flux

# with boundary conditions
#  - Fixed temperature ``T_{surface}`` at ``z_{min}`` (non-zero Dirichlet)
#  - No thermal flux at ``z_{min}`` (zero Neumann)

# Solving these equations is broken down into the following steps:
# 1) Preliminary configuration
# 2) PDEs
# 3) Space discretization
# 4) Time discretization / solver
# 5) Solver hooks / callbacks
# 6) Solve
# 7) Post-processing

# # Preliminary configuration

# ## [Loading code](@id Loading-code-heat)

# First, we'll load our pre-requisites:
#  - load external packages:
using MPI
using OrderedCollections
using Plots
using StaticArrays
using OrdinaryDiffEq
using DiffEqBase

#  - load CLIMAParameters and set up to use it:

using CLIMAParameters
struct EarthParameterSet <: AbstractEarthParameterSet end
const param_set = EarthParameterSet()

#  - load necessary ClimateMachine modules:
using ClimateMachine
using ClimateMachine.Mesh.Topologies
using ClimateMachine.Mesh.Grids
using ClimateMachine.DGMethods
using ClimateMachine.DGMethods.NumericalFluxes
using ClimateMachine.BalanceLaws:
    BalanceLaw, Prognostic, Auxiliary, Gradient, GradientFlux

using ClimateMachine.Mesh.Geometry: LocalGeometry
using ClimateMachine.MPIStateArrays
using ClimateMachine.GenericCallbacks
using ClimateMachine.ODESolvers
using ClimateMachine.VariableTemplates
using ClimateMachine.SingleStackUtils

#  - import necessary ClimateMachine modules: (`import`ing enables us to
#  provide implementations of these structs/methods)
import ClimateMachine.BalanceLaws:
    vars_state,
    source!,
    flux_second_order!,
    flux_first_order!,
    compute_gradient_argument!,
    compute_gradient_flux!,
    nodal_update_auxiliary_state!,
    nodal_init_state_auxiliary!,
    init_state_prognostic!,
    boundary_state!

import ClimateMachine.DGMethods: calculate_dt

# ## Initialization

# Define the float type (`Float64` or `Float32`)
FT = Float64;
# Initialize ClimateMachine for CPU.
ClimateMachine.init(; disable_gpu = true);

const clima_dir = dirname(dirname(pathof(ClimateMachine)));

# Load some helper functions for plotting
include(joinpath(clima_dir, "docs", "plothelpers.jl"));

# # Define the set of Partial Differential Equations (PDEs)

# ## Define the model

# Model parameters can be stored in the particular [`BalanceLaw`](@ref
# ClimateMachine.BalanceLaws.BalanceLaw), in this case, a `HeatModel`:

Base.@kwdef struct HeatModel{FT, APS} <: BalanceLaw
    "Parameters"
    param_set::APS
    "Heat capacity"
    ρc::FT = 1
    "Thermal diffusivity"
    α::FT = 0.01
    "Initial conditions for temperature"
    initialT::FT = 295.15
    "Bottom boundary value for temperature (Dirichlet boundary conditions)"
    T_bottom::FT = 300.0
    "Top flux (α∇ρcT) at top boundary (Neumann boundary conditions)"
    flux_top::FT = 0.0
end

# Create an instance of the `HeatModel`:
m = HeatModel{FT, typeof(param_set)}(; param_set = param_set);

# This model dictates the flow control, using [Dynamic Multiple
# Dispatch](https://en.wikipedia.org/wiki/Multiple_dispatch), for which
# kernels are executed.

# ## Define the variables

# All of the methods defined in this section were `import`ed in
# [Loading code](@ref Loading-code-heat) to let us provide
# implementations for our `HeatModel` as they will be used by
# the solver.

# Specify auxiliary variables for `HeatModel`
vars_state(::HeatModel, ::Auxiliary, FT) = @vars(z::FT, T::FT);

# Specify prognostic variables, the variables solved for in the PDEs, for
# `HeatModel`
vars_state(::HeatModel, ::Prognostic, FT) = @vars(ρcT::FT);

# Specify state variables whose gradients are needed for `HeatModel`
vars_state(::HeatModel, ::Gradient, FT) = @vars(ρcT::FT);

# Specify gradient variables for `HeatModel`
vars_state(::HeatModel, ::GradientFlux, FT) = @vars(α∇ρcT::SVector{3, FT});

# ## Define the compute kernels

# Specify the initial values in `aux::Vars`, which are available in
# `init_state_prognostic!`. Note that
# - this method is only called at `t=0`
# - `aux.z` and `aux.T` are available here because we've specified `z` and `T`
# in `vars_state` given `Auxiliary`
# in `vars_state`
function nodal_init_state_auxiliary!(
    m::HeatModel,
    aux::Vars,
    tmp::Vars,
    geom::LocalGeometry,
)
    aux.z = geom.coord[3]
    aux.T = m.initialT
end;

# Specify the initial values in `state::Vars`. Note that
# - this method is only called at `t=0`
# - `state.ρcT` is available here because we've specified `ρcT` in
# `vars_state` given `Prognostic`
function init_state_prognostic!(
    m::HeatModel,
    state::Vars,
    aux::Vars,
    coords,
    t::Real,
)
    state.ρcT = m.ρc * aux.T
end;

# The remaining methods, defined in this section, are called at every
# time-step in the solver by the [`BalanceLaw`](@ref
# ClimateMachine.BalanceLaws.BalanceLaw) framework.

# Compute/update all auxiliary variables at each node. Note that
# - `aux.T` is available here because we've specified `T` in
# `vars_state` given `Auxiliary`
function nodal_update_auxiliary_state!(
    m::HeatModel,
    state::Vars,
    aux::Vars,
    t::Real,
)
    aux.T = state.ρcT / m.ρc
end;

# Since we have second-order fluxes, we must tell `ClimateMachine` to compute
# the gradient of `ρcT`. Here, we specify how `ρcT` is computed. Note that
#  - `transform.ρcT` is available here because we've specified `ρcT` in
#  `vars_state` given `Gradient`
function compute_gradient_argument!(
    m::HeatModel,
    transform::Vars,
    state::Vars,
    aux::Vars,
    t::Real,
)
    transform.ρcT = state.ρcT
end;

# Specify where in `diffusive::Vars` to store the computed gradient from
# `compute_gradient_argument!`. Note that:
#  - `diffusive.α∇ρcT` is available here because we've specified `α∇ρcT` in
#  `vars_state` given `Gradient`
#  - `∇transform.ρcT` is available here because we've specified `ρcT`  in
#  `vars_state` given `Gradient`
function compute_gradient_flux!(
    m::HeatModel,
    diffusive::Vars,
    ∇transform::Grad,
    state::Vars,
    aux::Vars,
    t::Real,
)
    diffusive.α∇ρcT = -m.α * ∇transform.ρcT
end;

# We have no sources, nor non-diffusive fluxes.
function source!(m::HeatModel, _...) end;
function flux_first_order!(m::HeatModel, _...) end;

# Compute diffusive flux (``F(α, ρcT, t) = -α ∇ρcT`` in the original PDE).
# Note that:
# - `diffusive.α∇ρcT` is available here because we've specified `α∇ρcT` in
# `vars_state` given `GradientFlux`
function flux_second_order!(
    m::HeatModel,
    flux::Grad,
    state::Vars,
    diffusive::Vars,
    hyperdiffusive::Vars,
    aux::Vars,
    t::Real,
)
    flux.ρcT += diffusive.α∇ρcT
end;

# ### Boundary conditions

# Second-order terms in our equations, ``∇⋅(F)`` where ``F = -α∇ρcT``, are
# internally reformulated to first-order unknowns.
# Boundary conditions must be specified for all unknowns, both first-order and
# second-order unknowns which have been reformulated.

# The boundary conditions for `ρcT` (first order unknown)
function boundary_state!(
    nf,
    m::HeatModel,
    state⁺::Vars,
    aux⁺::Vars,
    n⁻,
    state⁻::Vars,
    aux⁻::Vars,
    bctype,
    t,
    _...,
)
    ## Apply Dirichlet BCs
    if bctype == 1 # At bottom
        state⁺.ρcT = m.ρc * m.T_bottom
    elseif bctype == 2 # At top
        nothing
    end
end;

# The boundary conditions for `ρcT` are specified here for second-order
# unknowns
function boundary_state!(
    nf,
    m::HeatModel,
    state⁺::Vars,
    diff⁺::Vars,
    aux⁺::Vars,
    n⁻,
    state⁻::Vars,
    diff⁻::Vars,
    aux⁻::Vars,
    bctype,
    t,
    _...,
)
    ## Apply Neumann BCs
    if bctype == 1 # At bottom
        nothing
    elseif bctype == 2 # At top
        diff⁺.α∇ρcT = n⁻ * m.flux_top
    end
end;

# # Spatial discretization

# Prescribe polynomial order of basis functions in finite elements
N_poly = 5;

# Specify the number of vertical elements
nelem_vert = 10;

# Specify the domain height
zmax = FT(1);

# Establish a `ClimateMachine` single stack configuration
driver_config = ClimateMachine.SingleStackConfiguration(
    "HeatEquation",
    N_poly,
    nelem_vert,
    zmax,
    param_set,
    m,
    numerical_flux_first_order = CentralNumericalFluxFirstOrder(),
);

# # Time discretization / solver

# Specify simulation time (SI units)
t0 = FT(0)
timeend = FT(40)

# In this section, we initialize the state vector and allocate memory for
# the solution in space (`dg` has the model `m`, which describes the PDEs
# as well as the function used for initialization). `SolverConfiguration`
# initializes the ODE solver, by default an explicit Low-Storage
# [Runge-Kutta](https://en.wikipedia.org/wiki/Runge%E2%80%93Kutta_methods)
# method. In this tutorial, we prescribe an option for an implicit
# `Kvaerno3` method.

# First, let's define how the time-step is computed, based on the
# [Fourier number](https://en.wikipedia.org/wiki/Fourier_number)
# (i.e., diffusive Courant number) is defined. Because
# the `HeatModel` is a custom model, we must define how both are computed.
# First, we must define our own implementation of `DGMethods.calculate_dt`,
# (which we imported):
function calculate_dt(dg, model::HeatModel, Q, Courant_number, t, direction)
    Δt = one(eltype(Q))
    CFL = DGMethods.courant(diffusive_courant, dg, model, Q, Δt, t, direction)
    return Courant_number / CFL
end

# Next, we'll define our implementation of `diffusive_courant`:
function diffusive_courant(
    m::HeatModel,
    state::Vars,
    aux::Vars,
    diffusive::Vars,
    Δx,
    Δt,
    t,
    direction,
)
    return Δt * m.α / (Δx * Δx)
end

# Finally, we initialize the state vector and solver
# configuration based on the given Fourier number.
# Note that, we can use a much larger Fourier number
# for implicit solvers as compared to explicit solvers.
use_implicit_solver = false
if use_implicit_solver
    given_Fourier = FT(30)

    solver_config = ClimateMachine.SolverConfiguration(
        t0,
        timeend,
        driver_config;
        ode_solver_type = ImplicitSolverType(OrdinaryDiffEq.Kvaerno3(
            autodiff = false,
            linsolve = LinSolveGMRES(),
        )),
        Courant_number = given_Fourier,
        CFL_direction = VerticalDirection(),
    )
else
    given_Fourier = FT(0.7)

    solver_config = ClimateMachine.SolverConfiguration(
        t0,
        timeend,
        driver_config;
        Courant_number = given_Fourier,
        CFL_direction = VerticalDirection(),
    )
end;

grid = solver_config.dg.grid;
Q = solver_config.Q;
aux = solver_config.dg.state_auxiliary;

# ## Inspect the initial conditions

# Let's export a plot of the initial state
output_dir = @__DIR__;

mkpath(output_dir);

z_scale = 100; # convert from meters to cm
z_key = "z";
z_label = "z [cm]";
z = get_z(grid, z_scale);

# Create an array to store the solution:
all_data = Dict[dict_of_nodal_states(solver_config, [z_key])]  # store initial condition at ``t=0``
time_data = FT[0]                                      # store time data

export_plot(
    z,
    all_data,
    ("ρcT",),
    joinpath(output_dir, "initial_condition.png");
    xlabel = "ρcT",
    ylabel = z_label,
    time_data = time_data,
);
# ![](initial_condition.png)

# It matches what we have in `init_state_prognostic!(m::HeatModel, ...)`, so
# let's continue.

# # Solver hooks / callbacks

# Define the number of outputs from `t0` to `timeend`
const n_outputs = 5;

# This equates to exports every ceil(Int, timeend/n_outputs) time-step:
const every_x_simulation_time = ceil(Int, timeend / n_outputs);

# The `ClimateMachine`'s time-steppers provide hooks, or callbacks, which
# allow users to inject code to be executed at specified intervals. In this
# callback, a dictionary of prognostic and auxiliary states are appended to
# `all_data` for time the callback is executed. In addition, time is collected
# and appended to `time_data`.
callback = GenericCallbacks.EveryXSimulationTime(every_x_simulation_time) do
    push!(all_data, dict_of_nodal_states(solver_config, [z_key]))
    push!(time_data, gettime(solver_config.solver))
    nothing
end;

# # Solve

# This is the main `ClimateMachine` solver invocation. While users do not have
# access to the time-stepping loop, code may be injected via `user_callbacks`,
# which is a `Tuple` of callbacks in [`GenericCallbacks`](@ref ClimateMachine.GenericCallbacks).
ClimateMachine.invoke!(solver_config; user_callbacks = (callback,));

# Append result at the end of the last time step:
push!(all_data, dict_of_nodal_states(solver_config, [z_key]));
push!(time_data, gettime(solver_config.solver));

# # Post-processing

# Our solution is stored in the array of dictionaries `all_data` whose keys are
# the output interval. The next level keys are the variable names, and the
# values are the values along the grid:

# To get `T` at ``t=0``, we can use `T_at_t_0 = all_data[1]["T"][:]`
@show keys(all_data[1])

# Let's plot the solution:

export_plot(
    z,
    all_data,
    ("ρcT",),
    joinpath(output_dir, "solution_vs_time.png");
    xlabel = "ρcT",
    ylabel = z_label,
    time_data = time_data,
);
# ![](solution_vs_time.png)

# The results look as we would expect: a fixed temperature at the bottom is
# resulting in heat flux that propagates up the domain. To run this file, and
# inspect the solution in `all_data`, include this tutorial in the Julia REPL
# with:

# ```julia
# include(joinpath("tutorials", "Land", "Heat", "heat_equation.jl"))
# ```

max_stars_repo_name: mwarusz/CLIMA

max_issues_repo_name: mwarusz/CLIMA

max_forks_repo_name: mwarusz/CLIMA


text:
@test mysignif(SpecificG(3.0, 0.307_842_258E3), TestDigits) ≈ mysignif(SpecificG_Ps(3.0, 0.5), TestDigits)
@test mysignif(SpecificG(3.0u"MPa", 0.307_842_258E3u"K"), TestDigits) ≈ mysignif(SpecificG_Ps(3.0u"MPa", 0.5u"kJ/kg/K"), TestDigits)

@test mysignif(SpecificG(80.0, 0.309_979_785E3), TestDigits) ≈ mysignif(SpecificG_Ps(80.0, 0.5), TestDigits)
@test mysignif(SpecificG(80.0u"MPa", 0.309_979_785E3u"K"), TestDigits) ≈ mysignif(SpecificG_Ps(80.0u"MPa", 0.5u"kJ/kg/K"), TestDigits)

@test mysignif(SpecificG(80.0, 0.565_899_909E3), TestDigits) ≈ mysignif(SpecificG_Ps(80.0, 3.0), TestDigits)
@test mysignif(SpecificG(80.0u"MPa", 0.565_899_909E3u"K"), TestDigits) ≈ mysignif(SpecificG_Ps(80.0u"MPa", 3.0u"kJ/kg/K"), TestDigits)

@test mysignif(SpecificF(3.0, 0.307_842_258E3), TestDigits) ≈ mysignif(SpecificF_Ps(3.0, 0.5), TestDigits)
@test mysignif(SpecificF(3.0u"MPa", 0.307_842_258E3u"K"), TestDigits) ≈ mysignif(SpecificF_Ps(3.0u"MPa", 0.5u"kJ/kg/K"), TestDigits)

@test mysignif(SpecificF(80.0, 0.309_979_785E3), TestDigits) ≈ mysignif(SpecificF_Ps(80.0, 0.5), TestDigits)
@test mysignif(SpecificF(80.0u"MPa", 0.309_979_785E3u"K"), TestDigits) ≈ mysignif(SpecificF_Ps(80.0u"MPa", 0.5u"kJ/kg/K"), TestDigits)

@test mysignif(SpecificF(80.0, 0.565_899_909E3), TestDigits) ≈ mysignif(SpecificF_Ps(80.0, 3.0), TestDigits)
@test mysignif(SpecificF(80.0u"MPa", 0.565_899_909E3u"K"), TestDigits) ≈ mysignif(SpecificF_Ps(80.0u"MPa", 3.0u"kJ/kg/K"), TestDigits)

@test SpecificV(3.0, 0.307_842_258E3) ≈ SpecificV_Ps(3.0, 0.5)
@test SpecificV(3.0u"MPa", 0.307_842_258E3u"K") ≈ SpecificV_Ps(3.0u"MPa", 0.5u"kJ/kg/K")

@test SpecificV(80.0, 0.309_979_785E3) ≈ SpecificV_Ps(80.0, 0.5)
@test SpecificV(80.0u"MPa", 0.309_979_785E3u"K") ≈ SpecificV_Ps(80.0u"MPa", 0.5u"kJ/kg/K")

@test SpecificV(80.0, 0.565_899_909E3) ≈ SpecificV_Ps(80.0, 3.0)
@test SpecificV(80.0u"MPa", 0.565_899_909E3u"K") ≈ SpecificV_Ps(80.0u"MPa", 3.0u"kJ/kg/K")

@test SpecificU(3.0, 0.307_842_258E3) ≈ SpecificU_Ps(3.0, 0.5)
@test SpecificU(3.0u"MPa", 0.307_842_258E3u"K") ≈ SpecificU_Ps(3.0u"MPa", 0.5u"kJ/kg/K")

@test SpecificU(80.0, 0.309_979_785E3) ≈ SpecificU_Ps(80.0, 0.5)
@test SpecificU(80.0u"MPa", 0.309_979_785E3u"K") ≈ SpecificU_Ps(80.0u"MPa", 0.5u"kJ/kg/K")

@test SpecificU(80.0, 0.565_899_909E3) ≈ SpecificU_Ps(80.0, 3.0)
@test SpecificU(80.0u"MPa", 0.565_899_909E3u"K") ≈ SpecificU_Ps(80.0u"MPa", 3.0u"kJ/kg/K")

@test SpecificH(3.0, 0.307_842_258E3) ≈ SpecificH_Ps(3.0, 0.5)
@test SpecificH(3.0u"MPa", 0.307_842_258E3u"K") ≈ SpecificH_Ps(3.0u"MPa", 0.5u"kJ/kg/K")

@test SpecificH(80.0, 0.309_979_785E3) ≈ SpecificH_Ps(80.0, 0.5)
@test SpecificH(80.0u"MPa", 0.309_979_785E3u"K") ≈ SpecificH_Ps(80.0u"MPa", 0.5u"kJ/kg/K")

@test SpecificH(80.0, 0.565_899_909E3) ≈ SpecificH_Ps(80.0, 3.0)
@test SpecificH(80.0u"MPa", 0.565_899_909E3u"K") ≈ SpecificH_Ps(80.0u"MPa", 3.0u"kJ/kg/K")

@test SpecificCP(3.0, 0.307_842_258E3) ≈ SpecificCP_Ps(3.0, 0.5)
@test SpecificCP(3.0u"MPa", 0.307_842_258E3u"K") ≈ SpecificCP_Ps(3.0u"MPa", 0.5u"kJ/kg/K")

@test SpecificCP(80.0, 0.309_979_785E3) ≈ SpecificCP_Ps(80.0, 0.5)
@test SpecificCP(80.0u"MPa", 0.309_979_785E3u"K") ≈ SpecificCP_Ps(80.0u"MPa", 0.5u"kJ/kg/K")

@test SpecificCP(80.0, 0.565_899_909E3) ≈ SpecificCP_Ps(80.0, 3.0)
@test SpecificCP(80.0u"MPa", 0.565_899_909E3u"K") ≈ SpecificCP_Ps(80.0u"MPa", 3.0u"kJ/kg/K")

@test SpecificCV(3.0, 0.307_842_258E3) ≈ SpecificCV_Ps(3.0, 0.5)
@test SpecificCV(3.0u"MPa", 0.307_842_258E3u"K") ≈ SpecificCV_Ps(3.0u"MPa", 0.5u"kJ/kg/K")

@test SpecificCV(80.0, 0.309_979_785E3) ≈ SpecificCV_Ps(80.0, 0.5)
@test SpecificCV(80.0u"MPa", 0.309_979_785E3u"K") ≈ SpecificCV_Ps(80.0u"MPa", 0.5u"kJ/kg/K")

@test SpecificCV(80.0,0.565_899_909E3) ≈ SpecificCV_Ps(80.0, 3.0)
@test SpecificCV(80.0u"MPa", 0.565_899_909E3u"K") ≈ SpecificCV_Ps(80.0u"MPa", 3.0u"kJ/kg/K")

@test SpeedOfSound(3.0, 0.307_842_258E3) ≈ SpeedOfSound_Ps(3.0, 0.5)
@test SpeedOfSound(3.0u"MPa", 0.307_842_258E3u"K") ≈ SpeedOfSound_Ps(3.0u"MPa", 0.5u"kJ/kg/K")

@test SpeedOfSound(80.0, 0.309_979_785E3) ≈ SpeedOfSound_Ps(80.0, 0.5)
@test SpeedOfSound(80.0u"MPa", 0.309_979_785E3u"K") ≈ SpeedOfSound_Ps(80.0u"MPa", 0.5u"kJ/kg/K")

@test SpeedOfSound(80.0, 0.565_899_909E3) ≈ SpeedOfSound_Ps(80.0, 3.0)
@test SpeedOfSound(80.0u"MPa", 0.565_899_909E3u"K") ≈ SpeedOfSound_Ps(80.0u"MPa", 3.0u"kJ/kg/K")


text:
using ModelingToolkitStandardLibrary.Blocks
using ModelingToolkit, OrdinaryDiffEq

@parameters t

@testset "Gain" begin
    @named c = Constant(; k=1)
    @named gain = Gain(1;)
    @named int = Integrator(; k=1)
    @named model = ODESystem([connect(c.output, gain.input), connect(gain.output, int.input)], t, systems=[int, gain, c])
    sys = structural_simplify(model)

    prob = ODEProblem(sys, Pair[int.x=>1.0], (0.0, 1.0))

    sol = solve(prob, Rodas4())

    @test all(sol[c.output.u] .≈ 1)
    @test sol[int.output.u][end] ≈ 2 # expected solution after 1s
end

@testset "Feedback loop" begin
    @named c = Constant(; k=2)
    @named gain = Gain(1;)
    @named int = Integrator(; k=1)
    @named fb = Feedback(;)
    @named model = ODESystem(
        [
            connect(c.output, fb.input1),
            connect(fb.input2, int.output),
            connect(fb.output, gain.input),
            connect(gain.output, int.input),
        ],
        t,
        systems=[int, gain, c, fb]
    )
    sys = structural_simplify(model)

    prob = ODEProblem(sys, Pair[int.x=>0.0], (0.0, 100.0))

    sol = solve(prob, Rodas4())
    @test sol[int.output.u][end] ≈ 2 # expected solution after 1s
end

@testset "Add" begin
    @named c1 = Constant(; k=1)
    @named c2 = Sine(; frequency=1)
    @named add = Add(;)
    @named int = Integrator(; k=1)
    @named model = ODESystem(
        [
            connect(c1.output, add.input1),
            connect(c2.output, add.input2),
            connect(add.output, int.input),
        ],
        t,
        systems=[int, add, c1, c2]
    )
    sys = structural_simplify(model)
    prob = ODEProblem(sys, Pair[int.x=>0.0], (0.0, 1.0))
    sol = solve(prob, Rodas4())
    @test sol[add.output.u] ≈ 1 .+ sin.(2*pi*sol.t)

    @testset "weights" begin
        k1 = -1
        k2 = 2
        @named add = Add(;k1=k1, k2=k2)
        @named model = ODESystem(
            [
                connect(c1.output, add.input1),
                connect(c2.output, add.input2),
                connect(add.output, int.input),
            ],
            t,
            systems=[int, add, c1, c2]
        )
        sys = structural_simplify(model)
        prob = ODEProblem(sys, Pair[int.x=>0.0], (0.0, 1.0))
        sol = solve(prob, Rodas4())
        @test sol[add.output.u] ≈ k1 .* 1 .+ k2 .* sin.(2*pi*sol.t)
    end
end

@testset "Add3" begin
    @named c1 = Constant(; k=1)
    @named c2 = Sine(; frequency=1)
    @named c3 = Sine(; frequency=2)
    @named add = Add3(;)
    @named int = Integrator(; k=1)
    @named model = ODESystem(
        [
            connect(c1.output, add.input1),
            connect(c2.output, add.input2),
            connect(c3.output, add.input3),
            connect(add.output, int.input),
        ],
        t,
        systems=[int, add, c1, c2, c3]
    )
    sys = structural_simplify(model)
    prob = ODEProblem(sys, Pair[int.x=>0.0], (0.0, 1.0))
    sol = solve(prob, Rodas4())
    @test sol[add.output.u] ≈ 1 .+ sin.(2*pi*sol.t) .+ sin.(2*pi*2*sol.t)

    @testset "weights" begin
        k1 = -1
        k2 = 2
        k3 = -pi
        @named add = Add3(;k1=k1, k2=k2, k3=k3)
        @named model = ODESystem(
            [
                connect(c1.output, add.input1),
                connect(c2.output, add.input2),
                connect(c3.output, add.input3),
                connect(add.output, int.input),
            ],
            t,
            systems=[int, add, c1, c2, c3]
        )
        sys = structural_simplify(model)
        prob = ODEProblem(sys, Pair[int.x=>0.0], (0.0, 1.0))
        sol = solve(prob, Rodas4())
        @test sol[add.output.u] ≈ k1 .* 1 .+ k2 .* sin.(2*pi*sol.t) .+ k3 .* sin.(2*pi*2*sol.t)

    end
end

@testset "Product" begin
    @named c1 = Constant(; k=2)
    @named c2 = Sine(; frequency=1)
    @named prod = Product(;)
    @named int = Integrator(; k=1)
    @named model = ODESystem(
        [
            connect(c1.output, prod.input1),
            connect(c2.output, prod.input2),
            connect(prod.output, int.input),
        ],
        t,
        systems=[int, prod, c1, c2]
    )
    sys = structural_simplify(model)
    prob = ODEProblem(sys, Pair[int.x=>0.0], (0.0, 1.0))
    sol = solve(prob, Rodas4())
    @test sol[prod.output.u] ≈ 2 * sin.(2*pi*sol.t)
end

@testset "Division" begin
    @named c1 = Sine(; frequency=1)
    @named c2 = Constant(; k=2)
    @named div = Division(;)
    @named int = Integrator(; k=1)
    @named model = ODESystem(
        [
            connect(c1.output, div.input1),
            connect(c2.output, div.input2),
            connect(div.output, int.input),
        ],
        t,
        systems=[int, div, c1, c2]
    )
    sys = structural_simplify(model)
    prob = ODEProblem(sys, Pair[int.x=>0.0], (0.0, 1.0))
    sol = solve(prob, Rodas4())
    @test sol[div.output.u] ≈ sin.(2*pi*sol.t) ./ 2
end

@testset "Abs" begin
    @named c = Sine(; frequency=1)
    @named absb = Abs(;)
    @named int = Integrator(; k=1)
    @named model = ODESystem(
        [
            connect(c.output, absb.input),
            connect(absb.output, int.input),
        ],
        t,
        systems=[int, absb, c]
    )
    sys = structural_simplify(model)
    prob = ODEProblem(sys, Pair[int.x=>0.0], (0.0, 1.0))
    sol = solve(prob, Rodas4())
    @test sol[absb.output.u] ≈ abs.(sin.(2*pi*sol.t))
end

@testset "MatrixGain" begin
    K = [1 2; 3 4]
    @named gain = MatrixGain(K;)
    # TODO:
end

@testset "Sum" begin
    @named s = Sum(2;)
    # TODO:
end

@testset "Math" begin
    for (block, func) in [(Abs, abs), (Sign, sign), (Sin, sin), (Cos, cos), (Tan, tan), (Asin, asin), (Acos, acos), (Atan, atan), (Sinh, sinh), (Cosh, cosh), (Tanh, tanh), (Exp, exp)]
        @named source = Sine(frequency=1)
        @named b = block()
        @named int = Integrator()
        @named model = ODESystem([connect(source.output, b.input), connect(b.output, int.input)], t, systems=[int, b, source])
        sys = structural_simplify(model)
        prob = ODEProblem(sys, Pair[int.x=>0.0], (0.0, 1.0))
        sol = solve(prob, Rodas4())
        @test sol[b.output.u] ≈ func.(sol[source.output.u])
    end

    # input must be positive
    for (block, func) in [(Sqrt, sqrt), (Log, log), (Log10, log10)]
        @named source = Sine(; frequency=1, offset=2)
        @named b = block()
        @named int = Integrator()
        @named model = ODESystem([connect(source.output, b.input), connect(b.output, int.input)], t, systems=[int, b, source])
        sys = structural_simplify(model)
        prob = ODEProblem(sys, Pair[int.x=>0.0], (0.0, 1.0))
        sol = solve(prob, Rodas4())
        @test sol[b.output.u] ≈ func.(sol[source.output.u])
    end
end

@testset "Atan2" begin
    @named c1 = Constant(; k=1)
    @named c2 = Constant(; k=2)
    @named b = Atan2(;)
    @named int = Integrator(; k=1)
    @named model = ODESystem(
        [
            connect(c1.output, b.input1),
            connect(c2.output, b.input2),
            connect(b.output, int.input),
        ],
        t,
        systems=[int, b, c1, c2]
    )
    sys = structural_simplify(model)
    prob = ODEProblem(sys, Pair[int.x=>0.0], (0.0, 1.0))
    sol = solve(prob, Rodas4())
    @test sol[int.output.u][end] ≈ atan(1, 2)
end


text:
abstract type AbstractTreeArray{T, N} <: AbstractArray{T, N} end

Base.size(x::AbstractTreeArray) = x.dims
leaftype(x::AbstractTreeArray) = leaftype(x.tree)
leafeltype(x::AbstractTreeArray) = leafeltype(x.tree)
nleafblocks(x::AbstractTreeArray) = @. ($size(x) - 1) ÷ $leafblockunit(x) + 1
leafblockunit(x::AbstractTreeArray) = leafblockunit(x.tree)
rootnode(x::AbstractTreeArray) = rootnode(x.tree)

@inline function Base.getindex(x::AbstractTreeArray{<: Any, N}, i::Vararg{Int, N}) where {N}
    @boundscheck checkbounds(x, i...)
    @inbounds x.tree[i...]
end

@inline function Base.setindex!(x::AbstractTreeArray{<: Any, N}, v, i::Vararg{Int, N}) where {N}
    @boundscheck checkbounds(x, i...)
    @inbounds x.tree[i...] = v
    x
end

for f in (:isactive, :allocate!, :isallocated)
    @eval @inline function $f(x::AbstractTreeArray{<: Any, N}, I::Vararg{Int, N}) where {N}
        @boundscheck checkbounds(x, I...)
        @inbounds $f(x.tree, I...)
    end
end

function allocate!(x::AbstractTreeArray, mask::AbstractArray{Bool})
    promote_shape(x, mask)
    allocate!(rootnode(x.tree), mask)
end

cleanup!(x::AbstractTreeArray) = (cleanup!(x.tree); x)
Base.fill!(x::AbstractTreeArray, ::Nothing) = (fill!(x.tree, nothing); x)

struct TreeArray{T, N, Tnode} <: AbstractTreeArray{T, N}
    tree::TreeView{T, N, Tnode}
    dims::NTuple{N, Int}
    function TreeArray{T, N, Tnode}(tree::TreeView, dims::Tuple) where {T, N, Tnode}
        @assert (CartesianIndices(tree) ∩ CartesianIndices(dims)) == CartesianIndices(dims)
        @assert leaftype(Tnode) <: LeafNode
        new(tree, dims)
    end
end

function TreeArray{T}(dims::NTuple{N, Int}) where {T, N}
    TreeArray(DynamicHashNode{Node{LeafNode{T, N, 3}, N, 3}, N}, dims)
end

struct StructTreeArray{T, N, Tnode} <: AbstractTreeArray{T, N}
    tree::TreeView{T, N, Tnode}
    dims::NTuple{N, Int}
    function StructTreeArray{T, N, Tnode}(tree::TreeView, dims::Tuple) where {T, N, Tnode}
        @assert (CartesianIndices(tree) ∩ CartesianIndices(dims)) == CartesianIndices(dims)
        @assert leaftype(Tnode) <: StructLeafNode
        new(tree, dims)
    end
end

function StructTreeArray{T}(dims::NTuple{N, Int}) where {T, N}
    StructTreeArray(DynamicHashNode{Node{@StructLeafNode{T, N, 3}, N, 3}, N}, dims)
end

Base.propertynames(x::StructTreeArray{T}) where {T} = (:tree, :dims, fieldnames(T)...)
@inline function Base.getproperty(x::StructTreeArray{<: Any, N}, name::Symbol) where {N}
    name == :tree && return getfield(x, :tree)
    name == :dims && return getfield(x, :dims)
    T = fieldtype(leafeltype(x), name)
    PropertyArray{T, N, name}(x)
end

# construtors
for TreeArrayType in (:TreeArray, :StructTreeArray)
    @eval begin
        function $TreeArrayType(tree::TreeView{T, N, Tnode}, dims::Tuple) where {T, N, Tnode}
            $TreeArrayType{T, N, Tnode}(tree, dims)
        end

        function $TreeArrayType(::Type{Tnode}, dims::Tuple) where {Tnode <: Union{Node, HashNode}}
            $TreeArrayType(TreeView(Tnode()), dims)
        end

        function $TreeArrayType(::Type{Tnode}, dims::Tuple) where {Tnode <: Union{DynamicNode, DynamicHashNode}}
            t_child = childtype(Tnode)
            blockdims = block_index(totalsize(t_child), dims...)
            tree = TreeView(Tnode(blockdims...))
            $TreeArrayType(tree, dims)
        end
        $TreeArrayType(::Type{Tnode}, dims::Int...) where {Tnode} = $TreeArrayType(Tnode, dims)

        $TreeArrayType{T}(dims::Int...) where {T} = $TreeArrayType{T}(dims)
    end
end


text:
# ----------------------------------------------------------------------------------- #
# Copyright (c) 2020 Varnerlab
# Robert Frederick School of Chemical and Biomolecular Engineering
# Cornell University, Ithaca NY 14850

# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
# THE SOFTWARE.
# ----------------------------------------------------------------------------------- #

"""
    flux(S,[Lv,Uv],[Lx,Ux],c; min_flag=true)
    Computes the optimal metabolic flux distribution given the constraints.
    Inputs:
    `S` - stoichiometric_matrix (M x R)
    `[Lv,Uv]` - (R x 2) array of the flux lower (L) and upper (U) bounds
    `[Lx,Ux]` - (M x 2) array of species lower (L) and upper (U) bounds (if at steady-state, L = U = 0)
    `c` - R x 1 vector holding indexes for objective vector
    Outputs:
    `objective_value` - value of the objective function at the optimum
    `calculated_flux_array` - R x 1 flux array at the optimum
    `dual_value_array` - R x 1 dual values
    `uptake_array` - M x 1 array of S*v
    `exit_flag` = 0 if optimal
    `status_flag` = 5 if optimal
"""
function flux(stoichiometric_matrix::Array{Float64,2}, default_bounds_array::Array{Float64,2},
    species_bounds_array::Array{Float64,2}, objective_coefficient_array::Array{Float64,1}; min_flag::Bool = true)

    try

        # set some constants -
        TIME_RESTART_LIM = 60

        # Get the stoichiometric_matrix from data_dictionary -
        (number_of_species, number_of_fluxes) = size(stoichiometric_matrix)

        # # Setup the GLPK problem -
        lp_problem = GLPK.glp_create_prob()
        GLPK.glp_set_prob_name(lp_problem, "sample")
        GLPK.glp_set_obj_name(lp_problem, "objective")

        # Set solver parameters
        solver_parameters = GLPK.glp_smcp()
        GLPK.glp_init_smcp(solver_parameters)
        solver_parameters.msg_lev = GLPK.GLP_MSG_OFF

        # Are we doing min -or- max?
        if min_flag == true
            GLPK.glp_set_obj_dir(lp_problem, GLPK.GLP_MIN)
        else
            GLPK.glp_set_obj_dir(lp_problem, GLPK.GLP_MAX)
        end

        # Set the number of constraints and fluxes -
        GLPK.glp_add_rows(lp_problem, number_of_species)
        GLPK.glp_add_cols(lp_problem, number_of_fluxes)

        # Setup flux bounds, and objective function -
        (number_of_fluxes, number_of_bounds) = size(default_bounds_array)
        for flux_index = 1:number_of_fluxes

            flux_lower_bound = default_bounds_array[flux_index, 1]
            flux_upper_bound = default_bounds_array[flux_index, 2]

            # Check bounds type ... default is DB -
            if (flux_upper_bound == flux_lower_bound)
                flux_constraint_type = GLPK.GLP_FX
            else
                flux_constraint_type = GLPK.GLP_DB
            end

            # flux symbol? (later use name - for now, fake it)
            flux_symbol = "R_" * string(flux_index)

            # Set the bounds in GLPK -
            GLPK.glp_set_col_name(lp_problem, flux_index, flux_symbol)
            GLPK.glp_set_col_bnds(lp_problem, flux_index, flux_constraint_type, flux_lower_bound, flux_upper_bound)
        end

        # Setup objective function -
        for (flux_index, obj_coeff) in enumerate(objective_coefficient_array)

            # Set the objective function value in GLPK -
            GLPK.glp_set_obj_coef(lp_problem, flux_index, obj_coeff)
        end

        # Setup problem constraints for the metabolites -
        for species_index = 1:number_of_species

            species_lower_bound = species_bounds_array[species_index, 1]
            species_upper_bound = species_bounds_array[species_index, 2]

            # defualt
            species_constraint_type = GLPK.GLP_FX
            if (species_lower_bound != species_upper_bound)
                species_constraint_type = GLPK.GLP_DB
            end

            # set the symbol -
            species_symbol = "x_" * string(species_index)

            # Set the species bounds in GLPK -
            GLPK.glp_set_row_name(lp_problem, species_index, species_symbol)
            GLPK.glp_set_row_bnds(lp_problem, species_index, species_constraint_type, species_lower_bound, species_upper_bound)

        end

        # Setup the stoichiometric array -
        counter = 1
        row_index_array = zeros(Cint, number_of_species * number_of_fluxes)
        col_index_array = zeros(Cint, number_of_species * number_of_fluxes)
        species_index_vector = collect(1:number_of_species)
        flux_index_vector = collect(1:number_of_fluxes)
        flat_stoichiometric_array = zeros(Float64, number_of_species * number_of_fluxes)
        for species_index in species_index_vector
            for flux_index in flux_index_vector
                row_index_array[counter] = species_index
                col_index_array[counter] = flux_index
                flat_stoichiometric_array[counter] = stoichiometric_matrix[species_index, flux_index]
                counter += 1
            end
        end
        GLPK.glp_load_matrix(lp_problem, number_of_species * number_of_fluxes, GLPK.offset(row_index_array), GLPK.offset(col_index_array), GLPK.offset(flat_stoichiometric_array))

        # Call the solver -
        exit_flag = GLPK.glp_simplex(lp_problem, solver_parameters)

        # Get the objective function value -
        objective_value = GLPK.glp_get_obj_val(lp_problem)

        # Get the calculated flux values from GLPK -
        calculated_flux_array = zeros(Float64, number_of_fluxes)
        for flux_index in flux_index_vector
            calculated_flux_array[flux_index] = GLPK.glp_get_col_prim(lp_problem, flux_index)
        end

        # Get the dual values -
        dual_value_array = zeros(Float64, number_of_fluxes)
        for flux_index in flux_index_vector
            dual_value_array[flux_index] = GLPK.glp_get_col_dual(lp_problem, flux_index)
        end

        # is this solution optimal?
        status_flag = GLPK.glp_get_status(lp_problem)

        # Calculate the uptake array -
        uptake_array = stoichiometric_matrix * calculated_flux_array

        # return results_tuple -
        results_tuple = (objective_value = objective_value, calculated_flux_array = calculated_flux_array,
            dual_value_array = dual_value_array, uptake_array = uptake_array, exit_flag = exit_flag,
            status_flag = status_flag)

        # Formulate the return tuple -
        return results_tuple

    catch error

        # what is our error policy? => for now, just print the message
        error_message = sprint(showerror, error, catch_backtrace())
        println(error_message)
    end
end


text:
# Based on StackOverflow
function subscriptnumber(i::Int)
    if i < 0
        c = [Char(0x208B)]
    else
        c = []
    end
    for d in reverse(digits(abs(i)))
        push!(c, Char(0x2080+d))
    end
    return join(c)
end

function gather_plot_information(x::AbstractDataDrivenProblem{N, C, D}) where {N,C, D}

    X = getfield(x, :X)
    Y = getfield(x, :Y)
    DX = getfield(x, :DX)
    U = getfield(x, :U)

    t = x.t

    ylab = is_direct(x) ? "Sample ID" : "t"

    returns = []

    for (xi, si) in zip([X, DX, Y, U], ["x", "\U02202\U0209C"*"x", "y", "u"])
        isempty(xi) && continue
        push!(returns, (t, xi, si))
    end

    return returns, ylab
end

function gather_plot_information(x::AbstractDataDrivenSolution)
    p = get_problem(x)
    t = p.t

    t = is_discrete(p) ? t[2:end] : t
    ylab = is_direct(p) ? "Sample ID" : "t"

    Y = get_target(p)
    outsym = is_direct(p) ? "y" : (is_discrete(p) ? "x" : "\U02202\U0209C"*"x")
    est_sym = outsym * "\U00302"
    X, _, ts, u = get_oop_args(p)
    Ŷ = x.basis(X, parameters(x), ts, u)
    return [(t, Y, outsym), (t, Ŷ, est_sym)], ylab
end

## Problem

@recipe function probplot(p::AbstractDataDrivenProblem)
    plotins, ylab = gather_plot_information(p)

    layout := (length(plotins), 1)
    isdirec = is_direct(p)

    for (i, pins) in enumerate(plotins)
        yi, xi, lab = pins
        suff = isdirec ? "" : "(t)"
        @series begin
            if i > 1
                title := ""
            end
            label --> reduce(hcat, map(j->lab*subscriptnumber(j)*suff, 1:size(xi,1)))

            ylabel --> lab

            if length(plotins) == i
                xlabel --> ylab
            end

            subplot := i
            yi, permutedims(xi)
        end
    end

end

# Step 2...

## Solution !
# TODO Dispatch for SINDy to show a separate plot of the active
# terms in the candidates -> Need LatexStrings for that (most likely)
# Get the symbols of the states for plotting
# Add extreme values for error ?

@recipe function resplot(x::AbstractDataDrivenSolution; show_data::Bool = false)
    plotins, xlab = gather_plot_information(x)
    isdirec = is_direct(get_problem(x))
    layout := (2,1)

    suff = isdirec ? "" : "(t)"

    t, Y, lab = plotins[1]
    _ , Ŷ, lab_ = plotins[2]

    @series begin
        label --> reduce(hcat, map(j->lab_*subscriptnumber(j)*suff, 1:size(Y, 1)))
        ylabel --> lab_
        subplot := 1
        seriestype := :path
        t, permutedims(Ŷ)
    end

    if show_data
        @series begin
            label --> nothing
            ylabel --> lab
            subplot := 1
            seriestype --> :path
            primary := false
            linestyle --> :dash
            color --> :black
            t, permutedims(Y)
        end
    end

    elab = lab*suff*"-"*lab_*suff

    @series begin
        title := ""
        label --> reduce(hcat, map(j->"e"*subscriptnumber(j)*suff, 1:size(Y, 1)))
        ylabel --> elab
        subplot := 2
        xlabel --> xlab
        t, permutedims(Y-Ŷ)
    end

end


text:
using BinDeps
using Compat

@BinDeps.setup

version = "3370"
baseurl = "ftp://heasarc.gsfc.nasa.gov/software/fitsio/c/"

if is_unix()
    archivename = "cfitsio$(version).tar.gz"
elseif is_windows()
    archivename = "cfitsio_MSVC_$(WORD_SIZE)bit_DLL_$(version).zip"
end

libcfitsio = library_dependency("libcfitsio", aliases=["cfitsio"])
downloadsdir = BinDeps.downloadsdir(libcfitsio)
libdir = BinDeps.libdir(libcfitsio)
srcdir = BinDeps.srcdir(libcfitsio)

if is_apple()
    libfilename = "libcfitsio.dylib"
elseif is_unix()
    libfilename = "libcfitsio.so"
elseif is_windows()
    libfilename = "cfitsio.dll"
end

# Unix
prefix = joinpath(BinDeps.depsdir(libcfitsio), "usr")
provides(Sources, URI(baseurl*archivename), libcfitsio, unpacked_dir="cfitsio")
provides(BuildProcess,
         (@build_steps begin
             GetSources(libcfitsio)
             @build_steps begin
                 ChangeDirectory(joinpath(srcdir, "cfitsio"))
                 FileRule(joinpath(libdir, libfilename),
                          @build_steps begin
                              `./configure --prefix=$prefix --enable-reentrant`
                              `make shared install`
                          end)
             end
          end), libcfitsio, os = :Unix)

# Windows
provides(BuildProcess,
	 (@build_steps begin
             FileDownloader(baseurl*archivename,
                            joinpath(downloadsdir, archivename))
	     CreateDirectory(srcdir, true)
	     FileUnpacker(joinpath(downloadsdir, archivename), srcdir,
                          joinpath(srcdir, libfilename))
	     CreateDirectory(libdir, true)
	     @build_steps begin
		 ChangeDirectory(srcdir)
	         FileRule(joinpath(libdir, libfilename), @build_steps begin
		          `powershell -Command "cp $(libfilename) $(joinpath(libdir, libfilename))"`
			  end)
             end
	  end), libcfitsio, os = :Windows)

if is_windows()
    push!(BinDeps.defaults, BuildProcess)
end

@BinDeps.install @compat Dict(:libcfitsio => :libcfitsio)

if is_windows()
    pop!(BinDeps.defaults)
end


text:
const NCIterable = Union{BaseAttributes,AbstractDimensions,AbstractDataset,AbstractGroups}
Base.length(a::NCIterable) = length(keys(a))

############################################################
# Mappings
############################################################

# Mapping between NetCDF types and Julia types
const jlType = Dict(
    NC_BYTE   => Int8,
    NC_UBYTE  => UInt8,
    NC_SHORT  => Int16,
    NC_USHORT => UInt16,
    NC_INT    => Int32,
    NC_UINT   => UInt32,
    NC_INT64  => Int64,
    NC_UINT64 => UInt64,
    NC_FLOAT  => Float32,
    NC_DOUBLE => Float64,
    NC_CHAR   => Char,
    NC_STRING => String
)

# Inverse mapping
const ncType = Dict(value => key for (key, value) in jlType)

"""
    fillvalue(::Type{Int8})
    fillvalue(::Type{UInt8})
    fillvalue(::Type{Int16})
    fillvalue(::Type{UInt16})
    fillvalue(::Type{Int32})
    fillvalue(::Type{UInt32})
    fillvalue(::Type{Int64})
    fillvalue(::Type{UInt64})
    fillvalue(::Type{Float32})
    fillvalue(::Type{Float64})
    fillvalue(::Type{Char})
    fillvalue(::Type{String})

Default fill-value for the given type.
"""
@inline fillvalue(::Type{Int8})    = NC_FILL_BYTE
@inline fillvalue(::Type{UInt8})   = NC_FILL_UBYTE
@inline fillvalue(::Type{Int16})   = NC_FILL_SHORT
@inline fillvalue(::Type{UInt16})  = NC_FILL_USHORT
@inline fillvalue(::Type{Int32})   = NC_FILL_INT
@inline fillvalue(::Type{UInt32})  = NC_FILL_UINT
@inline fillvalue(::Type{Int64})   = NC_FILL_INT64
@inline fillvalue(::Type{UInt64})  = NC_FILL_UINT64
@inline fillvalue(::Type{Float32}) = NC_FILL_FLOAT
@inline fillvalue(::Type{Float64}) = NC_FILL_DOUBLE
@inline fillvalue(::Type{Char})    = NC_FILL_CHAR
@inline fillvalue(::Type{String})  = NC_FILL_STRING

iswritable(ds::NCDataset) = ds.iswritable

"Make sure that a dataset is in data mode"
function datamode(ds)
    if ds.isdefmode[]
        nc_enddef(ds.ncid)
        ds.isdefmode[] = false
    end
end

"Make sure that a dataset is in define mode"
function defmode(ds)
    if !ds.isdefmode[]
        nc_redef(ds.ncid)
        ds.isdefmode[] = true
    end
end

"Initialize the ds._boundsmap variable"
function initboundsmap!(ds)
    ds._boundsmap = Dict{String,String}()
    for vname in keys(ds)
        v = variable(ds,vname)
        bounds = get(v.attrib,"bounds",nothing)

        if bounds !== nothing
            ds._boundsmap[bounds] = vname
        end
    end
end

############################################################
# High-level
############################################################

"""
    NCDataset(filename::AbstractString, mode = "r";
              format::Symbol = :netcdf4,
              share::Bool = false,
              diskless::Bool = false,
              persist::Bool = false,
              attrib = [])

Load, create, or even overwrite a NetCDF file at `filename`, depending on `mode`

* `"r"` (default) : open an existing netCDF file or OPeNDAP URL
   in read-only mode.
* `"c"` : create a new NetCDF file at `filename` (an existing file with the same
  name will be overwritten).
* `"a"` : open `filename` into append mode (i.e. existing data in the netCDF
  file is not overwritten and a variable can be added).

If `share` is true, the `NC_SHARE` flag is set allowing to have multiple
processes to read the file and one writer process. Likewise setting `diskless`
or `persist` to `true` will enable the flags `NC_DISKLESS` or `NC_PERSIST` flag.
More information is available in the [NetCDF C-API](https://www.unidata.ucar.edu/software/netcdf/docs/).

Notice that this does not close the dataset, use `close` on the
result (or see below the `do`-block).

The optional parameter `attrib` is an iterable of attribute name and attribute
value pairs, for example a `Dict`, `DataStructures.OrderedDict` or simply a
vector of pairs (see example below).

# Supported `format` values:

* `:netcdf4` (default): HDF5-based NetCDF format.
* `:netcdf4_classic`: Only netCDF 3 compatible API features will be used.
* `:netcdf3_classic`: classic netCDF format supporting only files smaller than 2GB.
* `:netcdf3_64bit_offset`: improved netCDF format supporting files larger than 2GB.

Files can also be open and automatically closed with a `do` block.

```julia
NCDataset("file.nc") do ds
    data = ds["temperature"][:,:]
end
```

Here is an attribute example:
```julia
using DataStructures
NCDataset("file.nc", "c", attrib = OrderedDict("title" => "my first netCDF file")) do ds
   defVar(ds,"temp",[10.,20.,30.],("time",))
end;
```

`Dataset` is an alias of `NCDataset`.
"""
function NCDataset(filename::AbstractString,
                   mode::AbstractString = "r";
                   format::Symbol = :netcdf4,
                   share::Bool = false,
                   diskless::Bool = false,
                   persist::Bool = false,
                   attrib = [])

    ncid = -1
    isdefmode = Ref(false)

    ncmode =
        if mode == "r"
            NC_NOWRITE
        elseif mode == "a"
            NC_WRITE
        elseif mode == "c"
            ncmode  = NC_CLOBBER
        else
            throw(NetCDFError(-1, "Unsupported mode '$(mode)' for filename '$(filename)'"))
        end

    if diskless
        ncmode = ncmode | NC_DISKLESS

        if persist
            ncmode = ncmode | NC_PERSIST
        end
    end

    if share
        @show "share"
        ncmode = ncmode | NC_SHARE
    end

    @debug "ncmode: $ncmode"

    if (mode == "r") || (mode == "a")
        ncid = nc_open(filename,ncmode)
    elseif mode == "c"
        if format == :netcdf3_64bit_offset
            ncmode = ncmode | NC_64BIT_OFFSET
        elseif format == :netcdf4_classic
            ncmode = ncmode | NC_NETCDF4 | NC_CLASSIC_MODEL
        elseif format == :netcdf4
            ncmode = ncmode | NC_NETCDF4
        elseif format == :netcdf3_classic
            # do nothing
        else
            throw(NetCDFError(-1, "Unkown format '$(format)' for filename '$(filename)'"))
        end

        ncid = nc_create(filename,ncmode)
        isdefmode[] = true
    end

    iswritable = mode != "r"
    ds = NCDataset(ncid,iswritable,isdefmode)

    # set global attributes
    for (attname,attval) in attrib
        ds.attrib[attname] = attval
    end

    return ds
end

function NCDataset(f::Function,args...; kwargs...)
    ds = NCDataset(args...; kwargs...)
    try
        f(ds)
    finally
        @debug "closing netCDF NCDataset $(ds.ncid) $(NCDatasets.path(ds))"
        close(ds)
    end
end

export NCDataset, Dataset

############################################################
# High-level: user convenience
############################################################
"""
    keys(ds::NCDataset)

Return a list of all variables names in NCDataset `ds`.
"""
Base.keys(ds::NCDataset) = listVar(ds.ncid)

"""
    path(ds::NCDataset)
Return the file path (or the opendap URL) of the NCDataset `ds`
"""
path(ds::NCDataset) = nc_inq_path(ds.ncid)
export path

"""
    sync(ds::NCDataset)

Write all changes in NCDataset `ds` to the disk.
"""
function sync(ds::NCDataset)
    datamode(ds)
    nc_sync(ds.ncid)
end
export sync

"""
    close(ds::NCDataset)

Close the NCDataset `ds`. All pending changes will be written
to the disk.
"""
function Base.close(ds::NCDataset)
    @debug "closing netCDF NCDataset $(ds.ncid) $(NCDatasets.path(ds))"
    try
        nc_close(ds.ncid)
    catch err
        # like Base, allow close on closed file
        if isa(err,NetCDFError)
            if err.code == NC_EBADID
                return ds
            end
        end
        rethrow()
    end
    # prevent finalize to close file as ncid can reused for future files
    ds.ncid = -1
    return ds
end
export close

############################################################
# Common methods
############################################################
function Base.iterate(a::NCIterable, state = keys(a))
    if length(state) == 0
        return nothing
    end

    return (state[1] => a[popfirst!(state)], state)
end

"""
    varbyattrib(ds, attname = attval)

Returns a list of variable(s) which has the attribute `attname` matching the value `attval`
in the dataset `ds`.
The list is empty if the none of the variables has the match.
The output is a list of `CFVariable`s.

# Examples

Load all the data of the first variable with standard name "longitude" from the
NetCDF file `results.nc`.

```julia-repl
julia> ds = NCDataset("results.nc", "r");
julia> data = varbyattrib(ds, standard_name = "longitude")[1][:]
```

"""
function varbyattrib(ds::NCDataset; kwargs...)
    # Start with an empty list of variables
    varlist = []

    # Loop on the variables
    for v in keys(ds)
        var = ds[v]

        matchall = true

        for (attsym,attval) in kwargs
            attname = String(attsym)

            # Check if the variable has the desired attribute
            if haskey(var.attrib, attname)
                # Check if the attribute value is the selected one
                if var.attrib[attname] != attval
                    matchall = false
                    break
                end
            else
                matchall = false
                break
            end
        end

        if matchall
            push!(varlist, var)
        end
    end

    return varlist
end
export varbyattrib

"""
    haskey(ds::NCDataset,name)
    haskey(d::Dimensions,name)
    haskey(ds::Attributes,name)

Return true if the NCDataset `ds` (or dimension/attribute list) has a variable (dimension/attribute) with the name `name`.
For example:

```julia
ds = NCDataset("/tmp/test.nc","r")
if haskey(ds,"temperature")
    println("The file has a variable 'temperature'")
end

if haskey(ds.dim,"lon")
    println("The file has a dimension 'lon'")
end
```

This example checks if the file `/tmp/test.nc` has a variable with the
name `temperature` and a dimension with the name `lon`.
"""
Base.haskey(a::NCIterable,name::AbstractString) = name in keys(a)
Base.in(name::AbstractString,a::NCIterable) = name in keys(a)

function Base.show(io::IO,ds::AbstractDataset; indent="")
    try
        dspath = path(ds)
        printstyled(io, indent, "NCDataset: ",dspath,"\n", color=section_color())
    catch err
        if isa(err,NetCDFError)
            if err.code == NC_EBADID
                print(io,"closed NetCDF NCDataset")
                return
            end
        end
        rethrow()
    end

    print(io,indent,"Group: ",groupname(ds),"\n")
    print(io,"\n")

    dims = collect(ds.dim)

    if length(dims) > 0
        printstyled(io, indent, "Dimensions\n",color=section_color())

        for (dimname,dimlen) in dims
            print(io,indent,"   $(dimname) = $(dimlen)\n")
        end
        print(io,"\n")
    end

    varnames = keys(ds)

    if length(varnames) > 0

        printstyled(io, indent, "Variables\n",color=section_color())

        for name in varnames
            show(io,variable(ds,name); indent = "$(indent)  ")
            print(io,"\n")
        end
    end

    # global attribues
    if length(ds.attrib) > 0
        printstyled(io, indent, "Global attributes\n",color=section_color())
        show(io,ds.attrib; indent = "$(indent)  ");
    end

    # groups
    groupnames = keys(ds.group)

    if length(groupnames) > 0
        printstyled(io, indent, "Groups\n",color=section_color())
        for groupname in groupnames
            show(io,group(ds,groupname); indent = "  ")
        end
    end

end

"""
    write(dest_filename::AbstractString, src::AbstractDataset; include = keys(src), exclude = [], idimensions = Dict())
    write(dest::NCDataset, src::AbstractDataset; include = keys(src), exclude = [], idimensions = Dict())

Write the variables of `src` dataset into an empty `dest` dataset (which must be opened in mode `"a"` or `"c"`).
The keywords `include` and `exclude` configure which variable of `src` should be included
(by default all), or which should be `excluded` (by default none).

If the first argument is a file name, then the dataset is open in create mode (`"c"`).

This function is useful when you want to save the dataset from a multi-file dataset.

`idimensions` is a dictionary with dimension names mapping to a list of
indices if only a subset of the dataset should be saved.

## Example

```
NCDataset(fname_src) do ds
    write(fname_slice,ds,idimensions = Dict("lon" => 2:3))
end
```

All variables in the source file `fname_src` with a dimension `lon` will be sliced
along the indices `2:3` for the `lon` dimension. All attributes (and variables
without a dimension `lon`) will be copied over unmodified.

It is assumed that all the variable of the output file can be loaded in memory.
"""
function Base.write(dest::NCDataset, src::AbstractDataset;
                    include = keys(src),
                    exclude = String[],
                    idimensions = Dict())

    unlimited_dims = unlimited(src.dim)

    for (dimname,dimlength) in src.dim
        isunlimited = dimname in unlimited_dims

        # if haskey(dest.dim,dimname)
        #     # check length
        #     if (dest.dim[dimname] !== src.dim[dimname]) && !isunlimited
        #         throw(DimensionMismatch("length of the dimensions $dimname are inconstitent in files $(path(dest)) and $(path(src))"))
        #     end
        # else
            if isunlimited
                defDim(dest, dimname, Inf)
            else
                if haskey(idimensions,dimname)
                    dimlength = length(idimensions[dimname])
                    @debug "subset $dimname $(idimensions[dimname])"
                end
                defDim(dest, dimname, dimlength)
            end
        # end
    end

    # loop over variables
    for varname in include
        (varname ∈ exclude) && continue
        @debug "Writing variable $varname..."

        cfvar = src[varname]
        dimension_names = dimnames(cfvar)

        # indices for subset
        index = ntuple(i -> get(idimensions,dimension_names[i],:),length(dimension_names))

        destvar = defVar(dest, varname, eltype(cfvar.var), dimension_names; attrib = cfvar.attrib)
        # copy data
        destvar.var[:] = cfvar.var[index...]
    end

    # loop over all global attributes
    for (attribname,attribval) in src.attrib
        dest.attrib[attribname] = attribval
    end

    # loop over all groups
    for (groupname,groupsrc) in src.group
        groupdest = defGroup(dest,groupname)
        write(groupdest,groupsrc; idimensions = idimensions)
    end
    return dest
end

function Base.write(dest_filename::AbstractString, src::AbstractDataset; kwargs...)
    NCDataset(dest_filename,"c") do dest
        write(dest,src; kwargs...)
    end
    return nothing
end

get_chunk_cache() = nc_get_chunk_cache()

function set_chunk_cache(;size=nothing,nelems=nothing,preemption=nothing)
    size_orig,nelems_orig,preemption_orig = nc_get_chunk_cache()
    size = (isnothing(size) ? size_orig : size)
    nelems = (isnothing(nelems) ? nelems_orig : nelems)
    preemption = (isnothing(preemption) ? preemption_orig : preemption)
    nc_set_chunk_cache(size,nelems,preemption)
end


text:
module GeometryProcessingPlayground

# import FileIO
# using GLMakie
# include("visualize.jl")

end


text:
#----------------------------------------
# @__{slide()}__
# @__{theme("simple")}__
# @__{header("""
pgfplotsx()
""")}__
#
#
# The Plots pgfplotsx-backend
# =====
#----------------------------------------
# @__{slide()}__
#
# ## Why?
#
# - Missing features
# - (Hopefully) easier to maintain
#----
# @__{subslide()}__
#
# - PGFPlotsX has better documentation
# - it makes the right promises
#----------------------------------------
# @__{slide()}__
#
# ## How?
#----
# @__{subslide()}__
#
# ### Plots -> PGFPlotsX -> LuaLaTeX -> pdf -> convert
#
# - translates `plotattributes` to PGFPlotsX.Options + Types (e.g. `PGFPlotsX.Axis`)
# - creates a `.tex`-file that gets compiled to pdf
# - this gets then converted to the desired output format
#----
# @__{subslide()}__
#
# ### Issues
#
# - LaTeX dependency + binary dependencies for conversion
#----------------------------------------
# @__{slide()}__
#
# ## What's new?
#----
# @__{subslide()}__
#
# ### 3D contours
using Plots;
pgfplotsx();
f(x, y) = exp(-x^2 - y^2)
xs = range(-1, 1, length = 20)
ys = range(-1, 1, length = 20)
contour3d(xs, ys, [f(x, y) for x in xs, y in ys])
#----
# @__{subslide()}__
#
# ### Heatmaps
xs = [string("x", i) for i = 1:10]
ys = [string("y", i) for i = 1:4]
z = float((1:4) * reshape(1:10, 1, :))
pgfx_plot = heatmap(xs, ys, z, aspect_ratio = 1)
#----
# @__{subslide()}__
#
# ### 3D scatter
n = 100
ts = range(0, stop = 8π, length = n)
x = ts .* map(cos, ts)
y = (0.1ts) .* map(sin, ts)
z = 1:n
pl = plot(
        x,
        y,
        z,
        zcolor = reverse(z),
        m = (10, 0.8, :blues, Plots.stroke(0)),
        leg = false,
        cbar = true,
        w = 5,
)
pgfx_plot = plot!(pl, zeros(n), zeros(n), 1:n, w = 10)
#----
# @__{subslide()}__
#
# ### Surface plots
f(x, y) = exp(sqrt(x^2 + y^2))
xs = range(-1, 1, length = 20)
ys = range(-1, 1, length = 20)
surface(xs, ys, [f(x, y) for x in xs, y in ys], label = "Cosy surface")
#----
# @__{subslide()}__
#
# ### Quiver plots
x = (-2pi):0.2:(2 * pi)
y = sin.(x)
u = ones(length(x))
v = cos.(x)
plot(x, y, quiver = (u, v), arrow = true, xlims = (-7, 8), ylims = (-1.5, 1.5))
#----
# @__{subslide()}__
#
# ### Ribbons
aa = rand(10)
bb = rand(10)
cc = rand(10)
plot(collect(1:10), cc, ribbon = (bb, aa))
#----
# @__{subslide()}__
#
# ### Fill between functions
plot(sin, 0, 2pi)
plot!(cos, 0, 2pi, fillrange = (2, 1)) # numbers are the series indices
#----------------------------------------
# @__{slide()}__
#
# ## What's missing?
#
# - Layouting is a mess (especially setting margins)
# - [Meta-Issue](https://github.com/JuliaPlots/Plots.jl/issues/2290)
#----
# @__{subslide()}__
#
# ### Layouting
#
# - On the plots side, we have (right/left/top/bottom) margins, bounding box of the plotarea and boinding box of a subplot.
# - On the pgfplots side, we have (right-left/top-bottom) margins (inner and outer), height and width of the axis as well as x/y-shift.
# - ⇢ No 1:1 mapping
# - In addition to that pgfplots have a page (which is currently a4paper by default).
#----------------------------------------
# @__{slide()}__
#
# ## Lessons learned
#
# - translate first and "improve" afterwards
# - might make sense to directly target pgf/TikZ if used as backend
#----
# @__{subslide()}__
#
# Pgfplots is a convenience layer over pgf/TikZ. Targeting it for wrapper packages like PGFPlots/X.jl is a good thing.
# As a backend you don't really need the convenience, because that is already built in the frontend, you are more looking for fine grained control.
# Therefore, if I were to write a LaTeX Makie-backend it would be called `PGFMakie.jl`. ;)
#----------------------------------------


text:
wrap(content) = """
<?xml version="1.0" encoding="UTF-8"?>
<OME xmlns="http://www.openmicroscopy.org/Schemas/OME/2016-06"
     xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"
     xmlns:OME="http://www.openmicroscopy.org/Schemas/OME/2016-06"
     xsi:schemaLocation="http://www.openmicroscopy.org/Schemas/OME/2016-06
     http://www.openmicroscopy.org/Schemas/OME/2016-06/ome.xsd">
    $content
</OME>"""


text:
using DelayDiffEq, DiffEqProblemLibrary.DDEProblemLibrary, DiffEqDevTools
using LinearAlgebra
using Test

DDEProblemLibrary.importddeproblems()

const prob = DDEProblemLibrary.prob_dde_constant_2delays_long_ip
const prob_oop = DDEProblemLibrary.prob_dde_constant_2delays_long_oop
const prob_scalar = DDEProblemLibrary.prob_dde_constant_2delays_long_scalar

const testsol = TestSolution(solve(prob, MethodOfSteps(Vern9());
                             abstol = 1/10^14, reltol = 1/10^14))

@testset "NLFunctional" begin
  alg = MethodOfSteps(Tsit5(); fpsolve = NLFunctional(; max_iter = 10))

  ## in-place problem

  sol = solve(prob, alg)

  # check statistics
  @test sol.destats.nf > 3000
  @test sol.destats.nsolve == 0
  @test sol.destats.nfpiter > 300
  @test sol.destats.nfpconvfail > 50

  # compare it with the test solution
  sol2 = appxtrue(sol, testsol)
  @test sol2.errors[:L∞] < 2.7e-4

  ## out-of-place problem

  sol_oop = solve(prob_oop, alg)

  # compare it with the in-place solution
  @test sol_oop.destats.nf == sol.destats.nf
  @test sol_oop.destats.nsolve == sol.destats.nsolve
  @test sol_oop.destats.nfpiter == sol.destats.nfpiter
  @test sol_oop.destats.nfpconvfail == sol.destats.nfpconvfail
  @test sol_oop.t ≈ sol.t
  @test sol_oop.u ≈ sol.u
  @test isapprox(sol.u, sol_oop.u; atol = 1e-7)

  ## scalar problem

  sol_scalar = solve(prob_scalar, alg)

  # compare it with the in-place solution
  @test sol_scalar.destats.nf == sol.destats.nf
  @test sol_scalar.destats.nsolve == sol.destats.nsolve
  @test sol_scalar.destats.nfpiter == sol.destats.nfpiter
  @test sol_scalar.destats.nfpconvfail == sol.destats.nfpconvfail
  @test sol_scalar.t ≈ sol.t
  @test sol_scalar.u ≈ sol[1, :]
end

@testset "NLAnderson" begin
  alg = MethodOfSteps(Tsit5(); fpsolve = NLAnderson(; max_iter = 10))

  ## in-place problem

  sol = solve(prob, alg)

  # check statistics
  @test sol.destats.nf < 2500
  @test sol.destats.nsolve > 0
  @test sol.destats.nfpiter < 250
  @test sol.destats.nfpconvfail < 50

  # compare it with the test solution
  sol2 = appxtrue(sol, testsol)
  @test sol2.errors[:L∞] < 2.7e-4

  ## out-of-place problem

  sol_oop = solve(prob_oop, alg)

  # compare it with the in-place solution
  @test_broken sol_oop.destats.nf == sol.destats.nf
  @test_broken sol_oop.destats.nsolve == sol.destats.nsolve
  @test_broken sol_oop.destats.nfpiter == sol.destats.nfpiter
  @test_broken sol_oop.destats.nfpconvfail == sol.destats.nfpconvfail
  @test_broken sol_oop.t ≈ sol.t
  @test_broken sol_oop.u ≈ sol.u
  @test appxtrue(sol, sol_oop).errors[:L∞] < 3e-7

  ## scalar problem

  sol_scalar = solve(prob_scalar, alg)

  # compare it with the in-place solution
  @test_broken sol_scalar.destats.nf == sol.destats.nf
  @test_broken sol_scalar.destats.nsolve == sol.destats.nsolve
  @test_broken sol_scalar.destats.nfpiter == sol.destats.nfpiter
  @test_broken sol_scalar.destats.nfpconvfail == sol.destats.nfpconvfail
  @test_broken sol_scalar.t ≈ sol.t
  @test_broken sol_scalar.u ≈ sol[1, :]
end


text:
using Images
import ImageView

imgdy = grayim(UInt8[y for y=0x00:0xff, x=0:100])
imgdx = grayim(UInt8[x for y=1:100, x=0x00:0xff])
ImageView.view(imgdy)
ImageView.view(imgdx)

A = zeros(300, 200)
lbl = similar(A, Int)
fill!(lbl, 0)
A[80:120, 50:150] = 0.8
lbl[80:120, 50:150] = 1
cx, cy = 200, 100
for j = -30:30, i = -30:30
    if i^2+j^2 <= 900
        A[cx+i,cy+j] = rand()
        lbl[cx+i,cy+j] = 2
    end
end
ImageView.viewlabeled(A, lbl)