text:
@get_units @with_kw struct SimpleReservoir{T}
    Δt::T | "s"                                         # Model time step [s]
    maxvolume::Vector{T} | "m3"                         # maximum storage (above which water is spilled) [m³]
    area::Vector{T} | "m2"                              # reservoir area [m²]
    maxrelease::Vector{T} | "m3 s-1"                    # maximum amount that can be released if below spillway [m³ s⁻¹]
    demand::Vector{T} | "m3 s-1"                        # minimum (environmental) flow requirement downstream of the reservoir [m³ s⁻¹]
    targetminfrac::Vector{T} | "-"                      # target minimum full fraction (of max storage) [-]
    targetfullfrac::Vector{T} | "-"                     # target fraction full (of max storage) [-]
    volume::Vector{T} | "m3"                            # reservoir volume [m³]
    inflow::Vector{T} | "m3"                            # total inflow into reservoir [m³]
    outflow::Vector{T} | "m3 s-1"                       # outflow from reservoir [m³ s⁻¹]
    totaloutflow::Vector{T} | "m3"                      # total outflow from reservoir [m³]
    percfull::Vector{T} | "-"                           # fraction full (of max storage) [-]
    demandrelease::Vector{T} | "m3 s-1"                 # minimum (environmental) flow released from reservoir [m³ s⁻¹]
    precipitation::Vector{T}                            # average precipitation for reservoir area [mm Δt⁻¹]
    evaporation::Vector{T}                              # average evaporation for reservoir area [mm Δt⁻¹]

    function SimpleReservoir{T}(args...) where {T}
        equal_size_vectors(args)
        return new(args...)
    end
end

statevars(::SimpleReservoir) = (:volume,)

function initialize_simple_reservoir(config, nc, inds_riv, nriv, pits, Δt)
    # read only reservoir data if reservoirs true
    # allow reservoirs only in river cells
    # note that these locations are only the reservoir outlet pixels
    reslocs = ncread(
        nc,
        config.input,
        "lateral.river.reservoir.locs";
        optional = false,
        sel = inds_riv,
        type = Int,
        fill = 0,
    )

    # this holds the same ids as reslocs, but covers the entire reservoir
    rescoverage_2d = ncread(
        nc,
        config.input,
        "lateral.river.reservoir.areas";
        optional = false,
        allow_missing = true,
    )
    # for each reservoir, a list of 2D indices, needed for getting the mean precipitation
    inds_res_cov = Vector{CartesianIndex{2}}[]

    rev_inds_reservoir = zeros(Int, size(rescoverage_2d))

    # construct a map from the rivers to the reservoirs and
    # a map of the reservoirs to the 2D model grid
    resindex = fill(0, nriv)
    inds_res = CartesianIndex{2}[]
    rescounter = 0
    for (i, ind) in enumerate(inds_riv)
        res_id = reslocs[i]
        if res_id > 0
            push!(inds_res, ind)
            rescounter += 1
            resindex[i] = rescounter
            rev_inds_reservoir[ind] = rescounter

            # get all indices related to this reservoir outlet
            # done in this loop to ensure that the order is equal to the order in the
            # SimpleReservoir struct
            res_cov = findall(isequal(res_id), rescoverage_2d)
            push!(inds_res_cov, res_cov)
        end
    end

    resdemand = ncread(
        nc,
        config.input,
        "lateral.river.reservoir.demand";
        optional = false,
        sel = inds_res,
        type = Float,
        fill = 0,
    )
    resmaxrelease = ncread(
        nc,
        config.input,
        "lateral.river.reservoir.maxrelease";
        optional = false,
        sel = inds_res,
        type = Float,
        fill = 0,
    )
    resmaxvolume = ncread(
        nc,
        config.input,
        "lateral.river.reservoir.maxvolume";
        optional = false,
        sel = inds_res,
        type = Float,
        fill = 0,
    )
    resarea = ncread(
        nc,
        config.input,
        "lateral.river.reservoir.area";
        optional = false,
        sel = inds_res,
        type = Float,
        fill = 0,
    )
    res_targetfullfrac = ncread(
        nc,
        config.input,
        "lateral.river.reservoir.targetfullfrac";
        optional = false,
        sel = inds_res,
        type = Float,
        fill = 0,
    )
    res_targetminfrac = ncread(
        nc,
        config.input,
        "lateral.river.reservoir.targetminfrac";
        optional = false,
        sel = inds_res,
        type = Float,
        fill = 0,
    )

    # for surface water routing reservoir locations are considered pits in the flow network
    # all upstream flow goes to the river and flows into the reservoir
    pits[inds_res] .= true

    n = length(resarea)
    @info "Read `$n` reservoir locations."
    reservoirs = SimpleReservoir{Float}(
        Δt = Δt,
        demand = resdemand,
        maxrelease = resmaxrelease,
        maxvolume = resmaxvolume,
        area = resarea,
        targetfullfrac = res_targetfullfrac,
        targetminfrac = res_targetminfrac,
        volume = res_targetfullfrac .* resmaxvolume,
        inflow = fill(mv, n),
        outflow = fill(mv, n),
        totaloutflow = fill(mv, n),
        percfull = fill(mv, n),
        demandrelease = fill(mv, n),
        precipitation = fill(mv, n),
        evaporation = fill(mv, n),
    )

    return reservoirs,
    resindex,
    (
        indices_outlet = inds_res,
        indices_coverage = inds_res_cov,
        reverse_indices = rev_inds_reservoir,
    ),
    pits
end

"""
Update a single reservoir at position `i`.

This is called from within the kinematic wave loop, therefore updating only for a single
element rather than all at once.
"""
function update(res::SimpleReservoir, i, inflow, timestepsecs)

    vol = max(
        0.0,
        (
            res.volume[i] +
            (inflow * timestepsecs) +
            (res.precipitation[i] * (timestepsecs / res.Δt) / 1000.0) * res.area[i] -
            (res.evaporation[i] * (timestepsecs / res.Δt) / 1000.0) * res.area[i]
        ),
    )

    percfull = vol / res.maxvolume[i]
    # first determine minimum (environmental) flow using a simple sigmoid curve to scale for target level
    fac = scurve(percfull, res.targetminfrac[i], Float(1.0), Float(30.0))
    demandrelease = min(fac * res.demand[i] * timestepsecs, vol)
    vol = vol - demandrelease

    wantrel = max(0.0, vol - (res.maxvolume[i] * res.targetfullfrac[i]))
    # Assume extra maximum Q if spilling
    overflow_q = max((vol - res.maxvolume[i]), 0.0)
    torelease = min(wantrel, overflow_q + res.maxrelease[i] * timestepsecs - demandrelease)
    vol = vol - torelease
    outflow = torelease + demandrelease
    percfull = vol / res.maxvolume[i]

    # update values in place
    res.outflow[i] = outflow / timestepsecs
    res.inflow[i] += inflow * timestepsecs
    res.totaloutflow[i] += outflow
    res.demandrelease[i] = demandrelease / timestepsecs
    res.percfull[i] = percfull
    res.volume[i] = vol

    return res
end

@get_units @with_kw struct NaturalLake{T}
    Δt::T | "s"                             # Model time step [s]
    lowerlake_ind::Vector{Int} | "-"        # Index of lower lake (linked lakes)
    area::Vector{T} | "m2"                  # lake area [m²]
    threshold::Vector{T} | "m"              # water level threshold H₀ [m] below that level outflow is zero
    storfunc::Vector{Int} | "-"             # type of lake storage curve, 1: S = AH, 2: S = f(H) from lake data and interpolation
    outflowfunc::Vector{Int} | "-"          # type of lake rating curve, 1: Q = f(H) from lake data and interpolation, 2: General Q = b(H - H₀)ᵉ, 3: Case of Puls Approach Q = b(H - H₀)²
    b::Vector{T} | "m3/2 s-1 (if e=3/2)"    # rating curve coefficient
    e::Vector{T} | "-"                      # rating curve exponent
    sh::Vector{Union{SH,Missing}}           # data for storage curve
    hq::Vector{Union{HQ,Missing}}           # data for rating curve
    waterlevel::Vector{T} | "m"             # waterlevel H [m] of lake
    inflow::Vector{T} | "m3"                # inflow to the lake [m³]
    storage::Vector{T} | "m3"               # storage lake [m³]
    outflow::Vector{T} | "m3 s-1"           # outflow lake [m³ s⁻¹]
    totaloutflow::Vector{T} | "m3"          # total outflow lake [m³]
    precipitation::Vector{T}                # average precipitation for lake area [mm Δt⁻¹]
    evaporation::Vector{T}                  # average evaporation for lake area [mm Δt⁻¹]

    function NaturalLake{T}(args...) where {T}
        equal_size_vectors(args)
        return new(args...)
    end
end

function initialize_natural_lake(config, nc, inds_riv, nriv, pits, Δt)
    # read only lake data if lakes true
    # allow lakes only in river cells
    # note that these locations are only the lake outlet pixels
    lakelocs_2d = ncread(
        nc,
        config.input,
        "lateral.river.lake.locs";
        optional = false,
        type = Int,
        fill = 0,
    )
    lakelocs = lakelocs_2d[inds_riv]

    # this holds the same ids as lakelocs, but covers the entire lake
    lakecoverage_2d = ncread(
        nc,
        config.input,
        "lateral.river.lake.areas";
        optional = false,
        allow_missing = true,
    )
    # for each lake, a list of 2D indices, needed for getting the mean precipitation
    inds_lake_cov = Vector{CartesianIndex{2}}[]

    rev_inds_lake = zeros(Int, size(lakecoverage_2d))

    # construct a map from the rivers to the lakes and
    # a map of the lakes to the 2D model grid
    lakeindex = fill(0, nriv)
    inds_lake = CartesianIndex{2}[]
    lakecounter = 0
    for (i, ind) in enumerate(inds_riv)
        lake_id = lakelocs[i]
        if lake_id > 0
            push!(inds_lake, ind)
            lakecounter += 1
            lakeindex[i] = lakecounter
            rev_inds_lake[ind] = lakecounter

            # get all indices related to this lake outlet
            # done in this loop to ensure that the order is equal to the order in the
            # NaturalLake struct
            lake_cov = findall(isequal(lake_id), lakecoverage_2d)
            push!(inds_lake_cov, lake_cov)
        end
    end

    lakearea = ncread(
        nc,
        config.input,
        "lateral.river.lake.area";
        optional = false,
        sel = inds_lake,
        type = Float,
        fill = 0,
    )
    lake_b = ncread(
        nc,
        config.input,
        "lateral.river.lake.b";
        optional = false,
        sel = inds_lake,
        type = Float,
        fill = 0,
    )
    lake_e = ncread(
        nc,
        config.input,
        "lateral.river.lake.e";
        optional = false,
        sel = inds_lake,
        type = Float,
        fill = 0,
    )
    lake_threshold = ncread(
        nc,
        config.input,
        "lateral.river.lake.threshold";
        optional = false,
        sel = inds_lake,
        type = Float,
        fill = 0,
    )
    linked_lakelocs = ncread(
        nc,
        config.input,
        "lateral.river.lake.linkedlakelocs";
        sel = inds_lake,
        defaults = 0,
        type = Int,
        fill = 0,
    )
    lake_storfunc = ncread(
        nc,
        config.input,
        "lateral.river.lake.storfunc";
        optional = false,
        sel = inds_lake,
        type = Int,
        fill = 0,
    )
    lake_outflowfunc = ncread(
        nc,
        config.input,
        "lateral.river.lake.outflowfunc";
        optional = false,
        sel = inds_lake,
        type = Int,
        fill = 0,
    )
    lake_waterlevel = ncread(
        nc,
        config.input,
        "lateral.river.lake.waterlevel";
        optional = false,
        sel = inds_lake,
        type = Float,
        fill = 0,
    )

    # for surface water routing lake locations are considered pits in the flow network
    # all upstream flow goes to the river and flows into the lake
    pits[inds_lake] .= true

    # This is currently the same length as all river cells, but will be the
    # length of all lake cells. To do that we need to introduce a mapping.
    n_lakes = length(inds_lake)
    lakelocs = lakelocs_2d[inds_lake]

    @info "Read `$n_lakes` lake locations."

    sh = Vector{Union{SH,Missing}}(missing, n_lakes)
    hq = Vector{Union{HQ,Missing}}(missing, n_lakes)
    lowerlake_ind = fill(0, n_lakes)
    # lake CSV parameter files are expected in the same directory as path_static
    path = dirname(input_path(config, config.input.path_static))

    for i = 1:n_lakes
        lakeloc = lakelocs[i]
        if linked_lakelocs[i] > 0
            lowerlake_ind[i] = only(findall(x -> x == linked_lakelocs[i], lakelocs))
        end

        if lake_storfunc[i] == 2
            csv_path = joinpath(path, "lake_sh_$lakeloc.csv")
            sh[i] = read_sh_csv(csv_path)
        end

        if lake_outflowfunc[i] == 1
            csv_path = joinpath(path, "lake_hq_$lakeloc.csv")
            hq[i] = read_hq_csv(csv_path)
        end

        if lake_outflowfunc[i] == 3 && lake_storfunc[i] != 1
            @warn(
                "For the modified pulse approach (LakeOutflowFunc = 3) the LakeStorFunc should be 1"
            )
        end
    end
    n = length(lakearea)
    lakes = NaturalLake{Float}(
        Δt = Δt,
        lowerlake_ind = lowerlake_ind,
        area = lakearea,
        threshold = lake_threshold,
        storfunc = lake_storfunc,
        outflowfunc = lake_outflowfunc,
        b = lake_b,
        e = lake_e,
        waterlevel = lake_waterlevel,
        sh = sh,
        hq = hq,
        inflow = fill(mv, n),
        storage = initialize_storage(lake_storfunc, lakearea, lake_waterlevel, sh),
        outflow = fill(mv, n),
        totaloutflow = fill(mv, n),
        precipitation = fill(mv, n),
        evaporation = fill(mv, n),
    )

    return lakes,
    lakeindex,
    (
        indices_outlet = inds_lake,
        indices_coverage = inds_lake_cov,
        reverse_indices = rev_inds_lake,
    ),
    pits
end

"Determine the initial storage depending on the storage function"
function initialize_storage(storfunc, area, waterlevel, sh)
    storage = similar(area)
    for i in eachindex(storage)
        if storfunc[i] == 1
            storage[i] = area[i] * waterlevel[i]
        else
            storage[i] = interpolate_linear(waterlevel[i], sh[i].H, sh[i].S)
        end
    end
    return storage
end

statevars(::NaturalLake) = (:waterlevel,)

"""
    statevars(::SBM; snow=false)
    statevars(::LateralSSF)
    statevars(::SurfaceFlow)
    statevars(::NaturalLake)
    statevars(::SimpleReservoir)

Get a tuple of symbols representing the fields that are model states.
"""
function statevars end

function interpolate_linear(x, xp, fp)
    if x <= minimum(xp)
        return minimum(fp)
    elseif x >= maximum(xp)
        return maximum(fp)
    else
        idx = findall(i -> i <= x, xp)
        i1 = last(idx)
        i2 = i1 + 1
        return fp[i1] * (1.0 - (x - xp[i1]) / (xp[i2] - xp[i1])) +
               fp[i2] * (x - xp[i1]) / (xp[i2] - xp[i1])
    end
end

"""
Update a single lake at position `i`.

This is called from within the kinematic wave loop, therefore updating only for a single
element rather than all at once.
"""
function update(lake::NaturalLake, i, inflow, doy, timestepsecs)

    lo = lake.lowerlake_ind[i]
    has_lowerlake = lo != 0

    col = min(doy, 365)
    ### Modified Puls Approach (Burek et al., 2013, LISFLOOD) ###
    # outflowfunc = 3
    # Calculate lake factor and SI parameter
    if lake.outflowfunc[i] == 3
        lakefactor = lake.area[i] / (timestepsecs * pow(lake.b[i], 0.5))
        si_factor =
            (
                lake.storage[i] +
                (lake.precipitation[i] - lake.evaporation[i]) *
                (timestepsecs / lake.Δt) *
                lake.area[i] / 1000.0
            ) / timestepsecs + inflow
        #Adjust SIFactor for ResThreshold != 0
        si_factor_adj = si_factor - lake.area[i] * lake.threshold[i] / timestepsecs
        #Calculate the new lake outflow/waterlevel/storage

        if si_factor_adj > 0.0
            outflow = pow(
                -lakefactor + pow((pow(lakefactor, 2.0) + 2.0 * si_factor_adj), 0.5),
                2.0,
            )
        else
            outflow = 0.0
        end
        storage = (si_factor - outflow) * timestepsecs
        waterlevel = storage / lake.area[i]
    end

    ### Linearisation for specific storage/rating curves ###
    if lake.outflowfunc[i] == 1 || lake.outflowfunc[i] == 2

        diff_wl = has_lowerlake ? lake.waterlevel[i] - lake.waterlevel[lo] : 0.0

        if lake.outflowfunc[i] == 1
            outflow =
                interpolate_linear(lake.waterlevel[i], lake.hq[i].H, lake.hq[i].Q[:, col])
        else
            if diff_wl >= 0.0
                if lake.waterlevel[i] > lake.threshold[i]
                    outflow =
                        lake.b[i] * pow((lake.waterlevel[i] - lake.threshold[i]), lake.e[i])
                else
                    outflow = Float(0)
                end
            else
                if lake.waterlevel[lo] > lake.threshold[i]
                    outflow =
                        -1.0 *
                        lake.b[i] *
                        pow((lake.waterlevel[lo] - lake.threshold[i]), lake.e[i])
                else
                    outflow = Float(0)
                end
            end
        end

        storage =
            lake.storage[i] +
            inflow * timestepsecs +
            (lake.precipitation[i] / 1000.0) * (timestepsecs / lake.Δt) * lake.area[i] -
            (lake.evaporation[i] / 1000.0) * (timestepsecs / lake.Δt) * lake.area[i] -
            outflow * timestepsecs

        waterlevel = if lake.storfunc[i] == 1
            lake.waterlevel[i] + (storage - lake.storage[i]) / lake.area[i]
        else
            interpolate_linear(storage, lake.sh[i].S, lake.sh[i].H)
        end

        # update lower lake (linked lakes) in case flow from lower lake to upper lake occurs
        if diff_wl < 0.0
            lowerlake_storage = lake.storage[lo] + outflow * timestepsecs

            lowerlake_waterlevel = if lake.storfunc[lo] == 1
                lake.waterlevel[lo] + (lowerlake_storage - lake.storage[lo]) / lake.area[lo]
            else
                interpolate_linear(lowerlake_storage, lake.sh[lo].S, lake.sh[lo].H)
            end

            # update values for the lower lake in place
            lake.outflow[lo] = -outflow
            lake.storage[lo] = lowerlake_storage
            lake.waterlevel[lo] = lowerlake_waterlevel
        end

    end

    # update values in place
    lake.outflow[i] = max(outflow, 0.0) # for a linked lake flow can be negative
    lake.waterlevel[i] = waterlevel
    lake.inflow[i] += inflow * timestepsecs
    lake.totaloutflow[i] += outflow * timestepsecs
    lake.storage[i] = storage

    return lake
end


text:
@testset "Find Tokens" begin
    a = raw"""some markdown then `code` and @@dname block @@"""

    tokens = F.find_tokens(a, F.MD_TOKENS, F.MD_1C_TOKENS)

    @test tokens[1].name == :CODE_SINGLE
    @test tokens[2].name == :CODE_SINGLE
    @test tokens[3].name == :DIV_OPEN
    @test tokens[3].ss == "@@dname"
    @test tokens[4].ss == "@@"
    @test tokens[5].name == :EOS
end

@testset "Find blocks" begin
    st = raw"""
        some markdown then `code` and
        @@dname block @@
        then maybe an escape
        ~~~
        escape block
        ~~~
        and done {target} done.
        """

    steps = explore_md_steps(st)
    blocks, tokens = steps[:ocblocks]

    # inline code block
    β = blocks[1]
    @test β.name == :CODE_INLINE
    @test β.ss == "`code`"

    # div block
    β = blocks[2]
    @test β.name == :DIV
    @test β.ss == "@@dname block @@"

    # escape block
    β = blocks[3]
    @test β.name == :ESCAPE
    @test β.ss == "~~~\nescape block\n~~~"

    # escape block
    β = blocks[4]
    @test β.name == :LXB
    @test β.ss == "{target}"
end

@testset "Unicode lx" begin
    st = raw"""
    Call me “$x$”, not $🍕$.
    """

    steps = explore_md_steps(st)
    blocks, _ = steps[:ocblocks]

    # first math block
    β = blocks[1]
    @test β.name == :MATH_A
    @test β.ss == "\$x\$"

    # second math block
    β = blocks[2]
    @test β.name == :MATH_A
    @test β.ss == "\$🍕\$"
end

@testset "Lx defs+coms" begin
    st = raw"""
        \newcommand{\E}[1]{\mathbb E\left[#1\right]}blah de blah
        ~~~
        escape b1
        ~~~
        \newcommand{\eqa}[1]{\begin{eqnarray}#1\end{eqnarray}}
        \newcommand{\R}{\mathbb R}
        Then something like
        \eqa{ \E{f(X)} \in \R &\text{if}& f:\R\maptso\R }
        and we could try to show latex:
        ```latex
        \newcommand{\brol}{\mathbb B}
        ```
        """

    lxdefs, tokens, braces, blocks = explore_md_steps(st)[:latex]

    @test lxdefs[1].name == "\\E"
    @test lxdefs[1].narg == 1
    @test lxdefs[1].def  == "\\mathbb E\\left[#1\\right]"
    @test lxdefs[2].name == "\\eqa"
    @test lxdefs[2].narg == 1
    @test lxdefs[2].def  == "\\begin{eqnarray}#1\\end{eqnarray}"
    @test lxdefs[3].name == "\\R"
    @test lxdefs[3].narg == 0
    @test lxdefs[3].def  == "\\mathbb R"

    @test blocks[1].name == :ESCAPE
    @test blocks[2].name == :CODE_BLOCK_LANG

    lxcoms, tokens = F.find_lxcoms(tokens, lxdefs, braces)

    @test lxcoms[1].ss == "\\eqa{ \\E{f(X)} \\in \\R &\\text{if}& f:\\R\\maptso\\R }"
    lxd = getindex(lxcoms[1].lxdef)
    @test lxd.name == "\\eqa"
end

@testset "Lxdefs 2" begin
    st = raw"""
        \newcommand{\com}{blah}
        \newcommand{\comb}[ 2]{hello #1 #2}
        """

    lxdefs, tokens, braces, blocks, lxcoms = explore_md_steps(st)[:latex]

    @test lxdefs[1].name == "\\com"
    @test lxdefs[1].narg == 0
    @test lxdefs[1].def == "blah"
    @test lxdefs[2].name == "\\comb"
    @test lxdefs[2].narg == 2
    @test lxdefs[2].def == "hello #1 #2"

    #
    # Errors
    #

    # testing malformed newcommands
    st = raw"""abc \newcommand abc"""
    tokens = F.find_tokens(st, F.MD_TOKENS, F.MD_1C_TOKENS)
    blocks, tokens = F.find_all_ocblocks(tokens, F.MD_OCB_ALL)
    # Ill formed newcommand (needs two {...})
    @test_throws F.LxDefError F.find_lxdefs(tokens, blocks)
    st = raw"""abc \newcommand{abc} def"""
    tokens = F.find_tokens(st, F.MD_TOKENS, F.MD_1C_TOKENS)
    blocks, tokens = F.find_all_ocblocks(tokens, F.MD_OCB_ALL)
    # Ill formed newcommand (needs two {...})
    @test_throws F.LxDefError F.find_lxdefs(tokens, blocks)
end

@testset "Lxcoms 2" begin
    st = raw"""
        \newcommand{\com}{HH}
        \newcommand{\comb}[1]{HH#1HH}
        Blah \com and \comb{blah} etc
        ```julia
        f(x) = x^2
        ```
        etc \comb{blah} then maybe
        @@adiv inner part @@ final.
        """

    lxdefs, tokens, braces, blocks, lxcoms = explore_md_steps(st)[:latex]

    @test lxcoms[1].ss == "\\com"
    @test lxcoms[2].ss == "\\comb{blah}"

    @test blocks[1].name == :CODE_BLOCK_LANG
    @test blocks[1].ss == "```julia\nf(x) = x^2\n```"
    @test F.content(blocks[1]) == "\nf(x) = x^2\n"

    @test blocks[2].name == :DIV
    @test blocks[2].ss == "@@adiv inner part @@"
    @test F.content(blocks[2]) == " inner part "

    #
    # Errors
    #

    st = raw"""
        \newcommand{\comb}[1]{HH#1HH}
        etc \comb then.
        """

    tokens = F.find_tokens(st, F.MD_TOKENS, F.MD_1C_TOKENS)
    blocks, tokens = F.find_all_ocblocks(tokens, F.MD_OCB_ALL)
    lxdefs, tokens, braces, blocks = F.find_lxdefs(tokens, blocks)
    # Command comb expects 1 argument and there should be no spaces ...
    @test_throws F.LxComError F.find_lxcoms(tokens, lxdefs, braces)
end

@testset "lxcoms3" begin
    st = raw"""
        text A1 \newcommand{\com}{blah}text A2 \com and
        ~~~
        escape B1
        ~~~
        \newcommand{\comb}[ 1]{\mathrm{#1}} text C1 $\comb{b}$ text C2
        \newcommand{\comc}[ 2]{part1:#1 and part2:#2} then \comc{AA}{BB}.
        """

    lxdefs, tokens, braces, blocks, lxcoms = explore_md_steps(st)[:latex]

    @test lxdefs[1].name == "\\com" && lxdefs[1].narg == 0 &&  lxdefs[1].def == "blah"
    @test lxdefs[2].name == "\\comb" && lxdefs[2].narg == 1 && lxdefs[2].def == "\\mathrm{#1}"
    @test lxdefs[3].name == "\\comc" && lxdefs[3].narg == 2 && lxdefs[3].def == "part1:#1 and part2:#2"
    @test blocks[1].name == :ESCAPE
    @test blocks[1].ss == "~~~\nescape B1\n~~~"
end

@testset "Merge-blocks" begin
    st = raw"""
        @def title = "Convex Optimisation I"
        \newcommand{\com}[1]{⭒!#1⭒}
        \com{A}

        then some
        ## blah
        end \com{B}.
        """

    lxdefs, tokens, braces, blocks, lxcoms = explore_md_steps(st)[:latex]

    @test blocks[1].name == :MD_DEF
    @test F.content(blocks[1]) == " title = \"Convex Optimisation I\""
    @test blocks[2].name == :COMMENT
    @test F.content(blocks[2]) == " comment "
    @test blocks[3].name == :H2
    @test F.content(blocks[3]) == " blah "

    @test lxcoms[1].ss == "\\com{A}"
    @test lxcoms[2].ss == "\\com{B}"

    b2i = F.merge_blocks(lxcoms, blocks)

    @test b2i[1].ss == "@def title = \"Convex Optimisation I\"\n"
    @test b2i[2].ss == "\\com{A}"
    @test b2i[3].ss == ""
    @test b2i[4].ss == "## blah \n"
    @test b2i[5].ss == "\\com{B}"
end

@testset "Header blocks" begin
    st = raw"""
        # t1
        1
        ## t2
        2 ## trick
        ### t3
        3
        #### t4
        4
        ##### t5
        5
        ###### t6
        6
        """

    tokens, blocks = explore_md_steps(st)[:filter]

    @test blocks[1].name == :H1
    @test blocks[2].name == :H2
    @test blocks[3].name == :H3
    @test blocks[4].name == :H4
    @test blocks[5].name == :H5
    @test blocks[6].name == :H6

    set_curpath("index.md")

    h = """
        # t1
        1
        ## t2
        2
        ## t3 `blah` etc
        3
        ### t4
        4
        ### t2
        5
        ### t2
        6
        """ |> seval
    @test isapproxstr(h, """

t1

1

t2

2

t3 blah etc

3

t4

4

t2

5

t2

6

        """)

    # pathological issue 241
    h = raw"""
        ## example
        A
        ## example
        B
        ## example 2
        C
        """ |> seval
    @test  isapproxstr(h, """

example

A

example

B

example 2

C

        """)
end

@testset "Line skip" begin
    h = raw"""
        Hello \\ goodbye
        """ |> seval
    @test isapproxstr(h, """
Hello
 goodbye
""")
end

@testset "Header+lx" begin
    h = "# blah" |> fd2html_td
    @test h // """
blah
"""
    h = raw"""
        \newcommand{\foo}{foo}
        \newcommand{\header}{# hello}
        \foo
        \header
        """ |> fd2html_td
    @test h // """
foo
hello

"""
    h = raw"""
        \newcommand{\foo}{foo}
        \foo hello
        """ |> fd2html_td
    @test h // """
foo hello
"""
    h = raw"""
        \newcommand{\foo}{blah}
        # \foo hello
        """ |> fd2html_td
    @test h // """
blah hello
"""
    h = raw"""
        \newcommand{\foo}{foo}
        \newcommand{\header}[2]{!#1 \foo #2}
        \header{##}{hello}
        """ |> fd2html_td
    @test h // """
foo  hello
"""
end

ext: jl


text:
module Tomato

using JuMP
using XLSX
using DataFrames
using Distances
using Combinatorics
using Distributions

include("struct.jl")
include("utils.jl")
include("core.jl")

#struct
export vtx,veh
export delivery

#utils
export extract!
export V,K,T,d,dist
export set_optimizer!
export get_optimizer
export reset_optimizer

#core
export bnb_tolerance #method
export relaxed,directRounding,optimalRounding #method 2
export cluster,reduced #method 3
export optimizeVolume,iterOV #improvement method 1
export deliveryCost,inventoryCost,totalCost
export stockLevels,feasibility
export worst_removal,random_removal,repair,lns

end


text:
function hamming_distance(firefly_1::Firefly, firefly_2::Firefly)
    """Calculates the hamming distance between two fireflies.
    """
    dist = 0
    dist_info = []
    for (i, n) in enumerate(firefly_1.path)
        if n.id != firefly_2.path[i].id
            dist += 1
            push!(dist_info, true)
        else
            push!(dist_info, false)
        end
    end
    return dist, dist_info
end

max_stars_repo_name: zahidkizmaz/firefly-algorithm-tsp

max_issues_repo_name: zahidkizmaz/firefly-algorithm-tsp

max_forks_repo_name: zahidkizmaz/firefly-algorithm-tsp


text:
using Test, Genify, IRTools
using Random, TracedRandom, StatsBase
using Distributions, Gen

using Genify: unwrap, transform!
using MacroTools: isexpr
using IRTools: IR

include("loops.jl")
include("transform.jl")
include("primitives.jl")


text:
using DataStructures
import IterTools: subsets
import Base.Iterators: product

function rank!(oracle::Function, P, scores, lhs, rhs)

    # If one direction is impossible, they are implicitly ranked already:
    (LP.maybe_lt(P, lhs, rhs) && LP.maybe_gt(P, lhs, rhs)) || return

    # Otherwise, consult the oracle:
    choice = oracle(scores, LP.get_base(P), lhs, rhs)
    [LP.eq!, LP.gt!, LP.lt!][choice](P, lhs, rhs)

end

function lhs_iter(crit, bits)

    start = bits .+ 1
    stop = crit .+ bits .- 1

    product((:).(start, stop)...)

end

function rhs_iter(crit, bits, lhs)

    down = bits .== 1

    step = ones(Int, size(bits))
    step[down] .= -1

    start = lhs .+ step

    stop = copy(crit)
    stop[down] .= 1

    product((:).(start, step, stop)...)

end

function paprika(oracle::Function, scores; maxdegree=2, total=nothing)

    Y = shape(scores); P = LP.Model(Y); n = length(Y)

    for i = 1:n, j=1:Y[i]
        pts = getpts(scores, i, j)
        if pts ≢ nothing
            LP.fix(P, i, j, pts)
        end
    end

    if total ≢ nothing
        LP.fixtotal(P, total)
    end

    for degree = 2:min(maxdegree, n), base in subsets(1:n, degree)

        LP.set_base!(P, base)

        crit = Y[base]

        for i = 1 : 2^degree - 1

            bits = digits(i, base=2, pad=degree)

            for lhs in lhs_iter(crit, bits), rhs in rhs_iter(crit, bits, lhs)

                rank!(oracle, P, scores, lhs, rhs)

            end

        end

    end

    LP.optimize!(P)

    for i = 1:n, j = 1:Y[i]
        setpts(scores, i, j, Int(round(LP.value(P, i, j))))
    end

end


text:
using Flux.Losses: crossentropy, binarycrossentropy, logitbinarycrossentropy

@testset "Losses" begin

x = [1.,2.,3.]
cx = gpu(x)
@test crossentropy(x,x) ≈ crossentropy(cx,cx)
@test crossentropy(x,x, agg=identity) ≈ crossentropy(cx,cx, agg=identity) |> cpu
@test crossentropy(x,x, agg=x->mean([1.0;2.0;3.0].*x)) ≈ crossentropy(cx,cx, agg=x->mean(gpu([1.0;2.0;3.0]).*x))

x = [-1.1491, 0.8619, 0.3127]
y = [1, 1, 0.]
@test binarycrossentropy(σ.(x), y) ≈ binarycrossentropy(gpu(σ.(x)), gpu(y))
@test logitbinarycrossentropy(x, y) ≈ logitbinarycrossentropy(gpu(x), gpu(y))

@testset "GPU grad tests" begin
  x = rand(Float32, 3,3)
  y = rand(Float32, 3,3)

  for loss in ALL_LOSSES
    gpu_gradtest(loss, x, y)
  end
end

end #testset


text:
module Day2
export day2

function parseInput(inputfile="inputs/input02.txt")
    input = Any[]
    for line in eachline(inputfile)
        #if line == "" break end
        range_raw, letter_raw, password = split(line, ' ')
        range = let (min,max) = split(range_raw, '-')
            parse(Int, min):parse(Int, max)
        end
        letter = letter_raw[1]
        push!(input, (range, letter, password))
    end
    input
end

function day2(inputfile="inputs/input02.txt")
    parsedInput = parseInput(inputfile)
    (count(isValidPWD, parsedInput), count(isValidPWDPart2, parsedInput))
end

isValidPWD(data) = isValidPWD(data...)
function isValidPWD(range, letter, password)
    count(c->(c==letter), password) ∈ range
end

isValidPWDPart2(data) = isValidPWDPart2(data...)
function isValidPWDPart2(range, letter, password)
    xor(letter==password[range.start], letter==password[range.stop])
end

end #module


text:
__precompile__()

module ConicBenchmarkUtilities

using GZip

if VERSION < v"0.7.0-"
    import Compat: undef
    import Compat: @warn
end

if VERSION > v"0.7.0-"
    using SparseArrays
    using LinearAlgebra
    # this is required because findall return type changed in v0.7
    function SparseArrays.findnz(A::AbstractMatrix)
        I = findall(!iszero, A)
        return (getindex.(I, 1), getindex.(I, 2), A[I])
    end
end

export readcbfdata, cbftompb, mpbtocbf, writecbfdata
export remove_zero_varcones, socrotated_to_soc, remove_ints_in_nonlinear_cones, dualize

include("cbf_input.jl")
include("cbf_output.jl")
include("mpb.jl")
include("preprocess_mpb.jl")
include("convex_to_cbf.jl")
include("jump_to_cbf.jl")

end # module


text:
function testEgal(a::Array,b::Array)
	array_length(a)==array_length(b) || return false
	for i = 0:array_length(a) - 1 a[i]==b[i]||return false
	end
	return true
end

function testString(a::String) true end
function testString(a::Array) false end
function testString(a::Any) false end

assert (1==1) && !(2==1) && (1!=2) && !(1!=1)
assert !(2==true) && 1!=false && (1!=false) && 1==true

assert 3 < 4 && !(3 < 3) && !(4 < 3)
assert 3. < 4 && !(3. < 3) && !(4. < 3)
assert 3 < 4. && !(3 < 3.) && !(4 < 3.)
assert 3. < 4. && !(3. < 3.) && !(4. < 3.)
assert (0 < true) && (true < 2) && !(true < 1) && !(1 < true)
assert (-1 < false) && (false < 1) && !(false < 0) && !(0 < false)

assert 3 <= 4 && (3 <= 3) && !(4 <= 3)
assert 3. <= 4 && (3. <= 3) && !(4. <= 3)
assert 3 <= 4. && (3 <= 3.) && !(4 <= 3.)
assert 3. <= 4. && (3. <= 3.) && !(4. <= 3.)
assert (1 <= true) && !(true <= 0) && (true <= 1) && !(2 <= true)
assert !(1 <= false) && !(false <= -1) && (false <= 0) && (0 <= false)

assert !(3 > 4) && !(3 > 3) && (4 > 3)
assert !(3. > 4) && !(3. > 3) && (4. > 3)
assert !(3 > 4.) && !(3 > 3.) && (4 > 3.)
assert !(3. > 4.) && !(3. > 3.) && (4. > 3.)
assert (2 > true) && (true > 0) && !(true > 1) && !(1 > true)
assert (1 > false) && (false > -1) && !(false > 0) && !(0 > false)

assert !(3 >= 4) && (3 >= 3) && (4 >= 3)
assert !(3. >= 4) && (3. >= 3) && (4. >= 3)
assert !(3 >= 4.) && (3 >= 3.) && (4 >= 3.)
assert !(3. >= 4.) && (3. >= 3.) && (4. >= 3.)
assert (1 >= true) && (true >= 1) && !(true >= 2) && !(0 >= true)
assert !(-1 >= false) && !(false >= 1) && (false >= 0) && (0 >= false)

assert int(3)==3 && int(false)== 0 && int(true)==1 && int(1.1) == 1;
assert int(-3.1) == -3;

assert testEgal([1,2,3],[1,2,3])
assert !testEgal([1,2,3],[1,3,2])

a = "Je suis samuel"
b = "je suis Samuel"

assert testString(a)
assert !testString([1,2,3])
assert !testEgal(a,b)

a[8] = 'S'
b[0] = a[0]
assert testEgal(a,b)

a = "je suis"
b = "Samuel"
c = "je suis Samuel"

assert testEgal(a@" "@b,c)
assert !testEgal(a@b,c)
asssert testEgal([true,true]@[false,true],[true,true,false,true])
assert !testEgal([1,3,5]@[3,2],[1,3,3,5,2])

println("Success")


text:
# tracer components of the Redi rotation tensor

"""
    AbstractIsopycnalTensor

Abstract supertype for an isopycnal rotation model.
"""
abstract type AbstractIsopycnalTensor end

"""
    struct IsopycnalTensor{FT} <: AbstractIsopycnalTensor

A tensor that rotates a vector into the isopycnal plane using the local slopes
of the buoyancy field.

Slopes are computed via `slope_x = - ∂b/∂x / ∂b/∂z` and `slope_y = - ∂b/∂y / ∂b/∂z`,
with the negative sign to account for the stable stratification (`∂b/∂z < 0`).
Then, the components of the isopycnal rotation tensor are:

```
               ⎡     1 + slope_y²         - slope_x slope_y      slope_x ⎤
(1 + slope²)⁻¹ | - slope_x slope_y          1 + slope_x²         slope_y |
               ⎣       slope_x                 slope_y            slope² ⎦
```

where `slope² = slope_x² + slope_y²`.
"""
struct IsopycnalTensor{FT} <: AbstractIsopycnalTensor
    minimum_bz :: FT
end

"""
    struct SmallSlopeIsopycnalTensor{FT} <: AbstractIsopycnalTensor

A tensor that rotates a vector into the isopycnal plane using the local slopes
of the buoyancy field and employing the small-slope approximation, i.e., that
the horizontal isopycnal slopes, `slope_x` and `slope_y` are ``≪ 1``. Slopes are
computed via `slope_x = - ∂b/∂x / ∂b/∂z` and `slope_y = - ∂b/∂y / ∂b/∂z`, with
the negative sign to account for the stable stratification (`∂b/∂z < 0`). Then,
by utilizing the small-slope appoximation, the components of the isopycnal
rotation tensor are:

```
⎡   1            0         slope_x ⎤
|   0            1         slope_y |
⎣ slope_x      slope_y      slope² ⎦
```

where `slope² = slope_x² + slope_y²`.
"""
struct SmallSlopeIsopycnalTensor{FT} <: AbstractIsopycnalTensor
    minimum_bz :: FT
end

SmallSlopeIsopycnalTensor(; minimum_bz = 0) = SmallSlopeIsopycnalTensor(minimum_bz)

@inline function isopycnal_rotation_tensor_xz_fcc(i, j, k, grid::AbstractGrid{FT}, buoyancy, tracers, slope_model::SmallSlopeIsopycnalTensor) where FT
    bx = ∂x_b(i, j, k, grid, buoyancy, tracers)
    bz = ℑxzᶠᵃᶜ(i, j, k, grid, ∂zᶜᶜᶠ, buoyancy_perturbation, buoyancy.model, tracers)
    bz = max(bz, slope_model.minimum_bz)

    slope_x = - bx / bz

    return ifelse(bz == 0, zero(FT), slope_x)
end

@inline function isopycnal_rotation_tensor_xz_ccf(i, j, k, grid::AbstractGrid{FT}, buoyancy, tracers, slope_model::SmallSlopeIsopycnalTensor) where FT
    bx = ℑxzᶜᵃᶠ(i, j, k, grid, ∂xᶠᶜᶜ, buoyancy_perturbation, buoyancy.model, tracers)
    bz = ∂z_b(i, j, k, grid, buoyancy, tracers)
    bz = max(bz, slope_model.minimum_bz)

    slope_x = - bx / bz

    return ifelse(bz == 0, zero(FT), slope_x)
end

@inline function isopycnal_rotation_tensor_yz_cfc(i, j, k, grid::AbstractGrid{FT}, buoyancy, tracers, slope_model::SmallSlopeIsopycnalTensor) where FT
    by = ∂y_b(i, j, k, grid, buoyancy, tracers)
    bz = ℑyzᵃᶠᶜ(i, j, k, grid, ∂zᶜᶜᶠ, buoyancy_perturbation, buoyancy.model, tracers)
    bz = max(bz, slope_model.minimum_bz)

    slope_y = - by / bz

    return ifelse(bz == 0, zero(FT), slope_y)
end

@inline function isopycnal_rotation_tensor_yz_ccf(i, j, k, grid::AbstractGrid{FT}, buoyancy, tracers, slope_model::SmallSlopeIsopycnalTensor) where FT
    by = ℑyzᵃᶜᶠ(i, j, k, grid, ∂yᶜᶠᶜ, buoyancy_perturbation, buoyancy.model, tracers)
    bz = ∂z_b(i, j, k, grid, buoyancy, tracers)
    bz = max(bz, slope_model.minimum_bz)

    slope_y = - by / bz

    return ifelse(bz == 0, zero(FT), slope_y)
end

@inline function isopycnal_rotation_tensor_zz_ccf(i, j, k, grid::AbstractGrid{FT}, buoyancy, tracers, slope_model::SmallSlopeIsopycnalTensor) where FT
    bx = ℑxzᶜᵃᶠ(i, j, k, grid, ∂xᶠᶜᶜ, buoyancy_perturbation, buoyancy.model, tracers)
    by = ℑyzᵃᶜᶠ(i, j, k, grid, ∂yᶜᶠᶜ, buoyancy_perturbation, buoyancy.model, tracers)
    bz = ∂z_b(i, j, k, grid, buoyancy, tracers)
    bz = max(bz, slope_model.minimum_bz)

    slope_x = - bx / bz
    slope_y = - by / bz
    slope² = slope_x^2 + slope_y^2

    return ifelse(bz == 0, zero(FT), slope²)
end


text:
using Test
using LogicCircuits
using ProbabilisticCircuits
using DataFrames

@testset "BDD tests" begin
    # Set up a logic constraint ϕ as a BDD and scope size n. Sample m PSDDs.
    function case(ϕ::Bdd, n::Integer; atol::Real = 1e-3)
        # All possible valuations (including impossible ones).
        M = all_valuations(collect(1:n))
        # Get only possible worlds.
        W = M[findall(ϕ.(eachrow(M))),:]
        # Assign random probabilities for each world in W.
        R = rand(1:20, size(W, 1))
        # Construct a dataset that maps the distribution of R (world W[i] repeats R[i] times).
        D = DataFrame(vcat([repeat(W[i,:], 1, R[i])' for i ∈ 1:size(W, 1)]...), :auto)
        C = learn_bdd(ϕ, D; pseudocount = 0.0)
        T = DataFrame(M, :auto)
        # Test smoothness.
        @test issmooth(C)
        # Test decomposability.
        @test isdecomposable(C)
        # Test determinism.
        @test isdeterministic(C)
        # Tests if respects vtree.
        @test respects_vtree(C, C.vtree)
        # Test consistency.
        @test (EVI(C, T) .> -Inf) == ϕ.(eachrow(M))
    end
    case((1 ∧ 2) ∨ (3 ∧ ¬4) ∨ (¬1 ∧ 5), 5)
    case((1 → 3) ∧ (5 → ¬2), 5)
    case((1 ∧ 2 ∧ 3) ∨ (4 ∧ 5), 5)
    case(exactly(3, collect(1:5)), 5)
    case(atleast(3, collect(1:5)), 5)
    case(atmost(3, collect(1:5)), 5)
end


text:
using DiffEqBase, StaticArrays
using DiffEqBase: __init, ODEFunction, AbstractODEIntegrator

export CDS_KWARGS
#####################################################################################
#                                    Defaults                                       #
#####################################################################################
using SimpleDiffEq: SimpleATsit5
const DEFAULT_SOLVER = SimpleATsit5()
const DEFAULT_DIFFEQ_KWARGS = (abstol = 1e-6, reltol = 1e-6)
const CDS_KWARGS = (alg = DEFAULT_SOLVER, DEFAULT_DIFFEQ_KWARGS...)

_get_solver(a) = haskey(a, :alg) ? a[:alg] : DEFAULT_SOLVER

#####################################################################################
#                               Interface to DiffEq                                 #
#####################################################################################
function ContinuousDynamicalSystem(prob::ODEProblem, args...)
    return ContinuousDynamicalSystem(prob.f.f, prob.u0, prob.p, args...;
           t0 = prob.tspan[1])
end

"""
    ODEProblem(ds::ContinuousDynamicalSystem, tspan; u0 = ds.u0, callback=CallbackSet())
Transform a continuous dynamical system into an `ODEProblem`, optionally using a different
initial state and/or a callback.
"""
function DiffEqBase.ODEProblem(ds::CDS{IIP}, tspan; u0 = ds.u0, callback=CallbackSet()) where {IIP}
    return ODEProblem{IIP}(ODEFunction(ds.f; jac = ds.jacobian), u0, tspan, ds.p, callback)
end

#####################################################################################
#                                 Integrators                                       #
#####################################################################################
stateeltype(::AbstractODEIntegrator{A, IIP, S}) where {A, IIP, S} = eltype(S)
stateeltype(::AbstractODEIntegrator{A, IIP, S}) where {
    A, IIP, S<:Vector{<:AbstractArray{T}}} where {T} = T

function integrator(ds::CDS{iip}, u0 = ds.u0;
    tfinal = Inf, diffeq...) where {iip}

    u = safe_state_type(Val{iip}(), u0)
    prob = ODEProblem{iip}(ds.f, u, (ds.t0, typeof(ds.t0)(tfinal)), ds.p)

    solver = _get_solver(diffeq)
    integ = __init(prob, solver; DEFAULT_DIFFEQ_KWARGS...,
                   save_everystep = false, diffeq...)
    return integ
end

############################### Tangent #############################################
function tangent_integrator(ds::CDS, k::Int = dimension(ds); kwargs...)
    return tangent_integrator(ds, orthonormal(dimension(ds), k); kwargs...)
end
function tangent_integrator(ds::CDS{IIP}, Q0::AbstractMatrix;
    u0 = ds.u0, diffeq...) where {IIP}

    t0 = ds.t0
    Q = safe_matrix_type(Val{IIP}(), Q0)
    u = safe_state_type(Val{IIP}(), u0)

    k = size(Q)[2]
    k > dimension(ds) && throw(ArgumentError(
    "It is not possible to evolve more tangent vectors than the system's dimension!"
    ))

    tangentf = create_tangent(ds.f, ds.jacobian, ds.J, Val{IIP}(), Val{k}())
    tanprob = ODEProblem{IIP}(tangentf, hcat(u, Q), (t0, typeof(t0)(Inf)), ds.p)

    solver = _get_solver(diffeq)
    return __init(tanprob, solver; DEFAULT_DIFFEQ_KWARGS..., internalnorm = _tannorm,
                  save_everystep = false, diffeq...)
end

function _tannorm(u::AbstractMatrix, t)
    s = size(u)[1]
    x = zero(eltype(u))
    for i in 1:s
        @inbounds x += u[i, 1]^2
    end
    return sqrt(x/length(x))
end
_tannorm(u::Real, t) = abs(u)

# Auto-diffed in-place version
function tangent_integrator(ds::CDS{true, S, D, F, P, JAC, JM, true},
    Q0::AbstractMatrix;
    u0 = ds.u0, diffeq...) where {S, D, F, P, JAC, JM}

    t0 = ds.t0
    Q = safe_matrix_type(Val{true}(), Q0)
    u = safe_state_type(Val{true}(), u0)

    k = size(Q)[2]
    k > dimension(ds) && throw(ArgumentError(
    "It is not possible to evolve more tangent vectors than the system's dimension!"
    ))

    tangentf = create_tangent_iad(ds.f, ds.J, u, ds.p, t0, Val{k}())
    tanprob = ODEProblem{true}(tangentf, hcat(u, Q), (t0, typeof(t0)(Inf)), ds.p)

    solver = _get_solver(diffeq)
    return __init(tanprob, solver; DEFAULT_DIFFEQ_KWARGS..., save_everystep = false,
                  internalnorm = _tannorm, diffeq...)
end

############################### Parallel ############################################
# Vector-of-Vector does not work with DiffEq atm:
# This is a workaround currently, until DiffEq allows Vector[Vector]
function create_parallel(ds::CDS{true}, states)
    st = Matrix(hcat(states...))
    L = size(st)[2]
    paralleleom = (du, u, p, t) -> begin
        for i in 1:L
            ds.f(view(du, :, i), view(u, :, i), p, t)
        end
    end
    return paralleleom, st
end

# const STIFFSOLVERS = (ImplicitEuler, ImplicitMidpoint, Trapezoid, TRBDF2,
# GenericImplicitEuler,
# GenericTrapezoid, SDIRK2, Kvaerno3, KenCarp3, Cash4, Hairer4, Hairer42, Kvaerno4,
# KenCarp4, Kvaerno5, KenCarp5, Rosenbrock23,
# Rosenbrock32, ROS3P, Rodas3, RosShamp4, Veldd4, Velds4, GRK4T,
# GRK4A, Ros4LStab, Rodas4, Rodas42, Rodas4P)

function parallel_integrator(ds::CDS, states; diffeq...)
    peom, st = create_parallel(ds, states)
    pprob = ODEProblem(peom, st, (ds.t0, typeof(ds.t0)(Inf)), ds.p)
    solver = _get_solver(diffeq)
    # if typeof(solver) ∈ STIFFSOLVERS
    #     error("Stiff solvers can't support a parallel integrator.")
    # end
    if !(typeof(ds) <: CDS{true})
        return __init(pprob, solver; DEFAULT_DIFFEQ_KWARGS..., save_everystep = false,
                      internalnorm = _parallelnorm, diffeq...)
    else
        return __init(pprob, solver; DEFAULT_DIFFEQ_KWARGS..., save_everystep = false,
                      internalnorm = _tannorm, diffeq...)
    end
end

_parallelnorm(u::AbstractVector, t) = @inbounds DiffEqBase.ODE_DEFAULT_NORM(u[1], t)
_parallelnorm(u::Real, t) = abs(u)

#####################################################################################
#                                 Trajectory                                        #
#####################################################################################
function trajectory(ds::ContinuousDynamicalSystem, T, u = ds.u0;
    dt = 0.01, Ttr = 0.0, save_idxs = nothing, diffeq...)

    a = svector_access(save_idxs)
    integ = integrator(ds, u; dt = dt, diffeq...)
    trajectory(ds, integ, T, u, dt, Ttr, a; diffeq...)
end

function trajectory(ds::CDS{IIP, S, D}, integ, T, u, dt, Ttr, a; diffeq...) where {IIP, S, D}
    t0 = ds.t0
    tvec = (t0+Ttr):dt:(T+t0+Ttr)
    X = isnothing(a) ? D : length(a)
    sol = Vector{SVector{X, eltype(S)}}(undef, length(tvec))
    step!(integ, Ttr)
    for (i, t) in enumerate(tvec)
        while t > integ.t
            step!(integ)
        end
        sol[i] = SVector{X, eltype(S)}(obtain_access(integ(t), a))
    end
    return Dataset(sol)
end

#####################################################################################
#                                    Get States                                     #
#####################################################################################
get_state(integ::AbstractODEIntegrator{Alg, IIP, S}) where {Alg, IIP, S<:AbstractVector} =
    integ.u
get_state(integ::AbstractODEIntegrator{Alg, IIP, S}) where {Alg, IIP, S<:AbstractMatrix} =
    integ.u[:, 1]
get_state(integ::AbstractODEIntegrator{Alg, IIP, S}) where {Alg, IIP, S<:Vector{<:AbstractVector}} =
    integ.u[1]
get_state(integ::AbstractODEIntegrator{Alg, IIP, S}, k::Int) where {Alg, IIP, S<:Vector{<:AbstractVector}} =
    integ.u[k]
get_state(integ::AbstractODEIntegrator{Alg, IIP, S}, k::Int) where {Alg, IIP, S<:AbstractMatrix} =
    integ.u[:, k]

function set_state!(
    integ::AbstractODEIntegrator{Alg, IIP, S}, u::AbstractVector, k::Int = 1
    ) where {Alg, IIP, S<:Vector{<:AbstractVector}}
    integ.u[k] = u
    u_modified!(integ, true)
end
function set_state!(
    integ::AbstractODEIntegrator{Alg, IIP, S}, u::AbstractVector
    ) where {Alg, IIP, S<:Matrix}
    integ.u[:, 1] .= u
    u_modified!(integ, true)
end
function set_state!(
    integ::AbstractODEIntegrator{Alg, IIP, S}, u::AbstractVector
    ) where {Alg, IIP, S<:SMatrix{D, K}} where {D, K}
    integ.u = hcat(SVector{D}(u), integ.u[:, SVector{K-1}(2:K...)])
    u_modified!(integ, true)
end

get_deviations(integ::AbstractODEIntegrator{Alg, IIP, S}) where {Alg, IIP, S<:Matrix} =
    @view integ.u[:, 2:end]

@generated function get_deviations(
    integ::AbstractODEIntegrator{Alg, IIP, S}) where {Alg, IIP, S<:SMatrix{D,K}} where {D,K}
    gens = [:($k) for k=2:K]
    quote
        sind = SVector{$(K-1)}($(gens...))
        integ.u[:, sind]
    end
end

set_deviations!(integ::AbstractODEIntegrator{Alg, IIP, S}, Q) where {Alg, IIP, S<:Matrix} =
    (integ.u[:, 2:end] .= Q; u_modified!(integ, true))
set_deviations!(integ::AbstractODEIntegrator{Alg, IIP, S}, Q) where {Alg, IIP, S<:SMatrix} =
    (integ.u = hcat(integ.u[:,1], Q); u_modified!(integ, true))

function DiffEqBase.reinit!(integ::AbstractODEIntegrator, u0::AbstractVector,
    Q0::AbstractMatrix; kwargs...)

    set_state!(integ, u0)
    set_deviations!(integ, Q0)
    reinit!(integ, integ.u; kwargs...)
end


text:
# TODO rename testdense

function mm_spec_testdense_dimspec()::MPCCDimSpec
    dimspec = MPCCDimSpec(5, 2, 3, 2, 3, 2, 0)
    return dimspec
end

function mm_spec_testdense_defn(x, pr, ps)

    f = x[1]^2 * x[2]^3 * x[3]^4 * x[4]^5 * x[5]^6 * ( pr[1] + pr[2] )^2
            + pr[1]^2 * pr[2]^3 * exp(-2*x[1] - 3*x[2] - 5*x[3] - 7*x[4] - 11*x[5])

    ce = Vector{Num}(undef, 2)
    ce[1] = ( x[1] + x[2] + x[3] + x[4] + x[5] )^2 * sin(pr[1]) * cos(pr[2])
    ce[2] = pr[1] * ( x[1] + x[2] )^3 + pr[2] * ( x[3] + x[4] + x[5] )^2

    ci = Vector{Num}(undef, 3)
    ci[1] = x[1] + x[2] * ( pr[1] + pr[2] )^2
    ci[2] = x[3] + x[4] * x[5] * ( pr[1] - pr[2] )^2
    ci[3] = ( x[2] - x[1] ) * ( x[4] - x[3] ) * x[5]^2 * ( pr[2] - pr[1] )^2

    F = Matrix{Num}(undef, 3, 2)
    F[1,1] = ( x[1] + x[2] )^2
    F[1,2] = ( x[3] + x[4] + x[5] ) * ( pr[1] - pr[2] )^2
    F[2,1] = ( x[3] + x[4] ) * ( pr[1] - pr[2] )
    F[2,2] = ( x[4] - x[1] )
    F[3,1] = ( 1 + x[5] )^7
    F[3,2] = ( 1 + pr[1] + pr[2] ) * sin(x[1])

    return MPCCDefinition(f, ce, ci, F)
end

function mm_spec_testdense_nzmask()::MPCCModelNZMask
    nz_mask = MPCCModelNZMask(
        # ce
        [],
        # ci
        [],
        # F
        [
            # F1
            [ ],
            # F2
            [ ]
        ] )
    return nz_mask
end

function mm_spec_testdense_testvectors()

    x_val = Vector{Rational}(undef, 5)
    x_val = [4//3, 27//10, 7//2, 1//13, 5]

    pr_val = Vector{Rational}(undef, 2)
    pr_val = [1//7, 23//10]

    ps_val = Int64[]

    # f:
    f_chk_eval = (12167*exp(-16342/195))/49000 + 15667766415/11881376

    # ce:
    ce_chk_eval = Vector{Float64}(undef, 2)
    ce_chk_eval[1] = (6046681*cos(23/10)*sin(1/7))/38025
    ce_chk_eval[2] = 1425923221/7985250

    # ci:
    ci_chk_eval = Vector{Float64}(undef, 3)
    ci_chk_eval[1] = 2564521/147000
    ci_chk_eval[2] = 67391/12740
    ci_chk_eval[3] = -83200849/152880

    # F:
    F_chk_eval = Matrix{Float64}(undef, 3, 2)
    F_chk_eval[1, 1] = 14641/900
    F_chk_eval[1, 2] = 5084623/127400
    F_chk_eval[2, 1] = -14043/1820
    F_chk_eval[2, 2] = -49/39
    F_chk_eval[3, 1] = 279936
    F_chk_eval[3, 2] = (241*sin(4/3))/70

    # gradf:
    gradf_chk_eval = Vector{Float64}(undef, 5)
    gradf_chk_eval[1] = 47003299245/23762752 - (12167*exp(-16342/195))/24500
    gradf_chk_eval[2] = 8704314675/5940688 - (36501*exp(-16342/195))/49000
    gradf_chk_eval[3] = 2238252345/1485172 - (12167*exp(-16342/195))/9800
    gradf_chk_eval[4] = 78338832075/913952 - (12167*exp(-16342/195))/7000
    gradf_chk_eval[5] = 9400659849/5940688 - (133837*exp(-16342/195))/49000

    # ce:
    gradce_chk_eval = Matrix{Float64}(undef, 5, 2)
    gradce_chk_eval[1, 1] = (4918*cos(23/10)*sin(1/7))/195
    gradce_chk_eval[1, 2] = 14641/2100
    gradce_chk_eval[2, 1] = (4918*cos(23/10)*sin(1/7))/195
    gradce_chk_eval[2, 2] = 14641/2100
    gradce_chk_eval[3, 1] = (4918*cos(23/10)*sin(1/7))/195
    gradce_chk_eval[3, 2] = 5129/130
    gradce_chk_eval[4, 1] = (4918*cos(23/10)*sin(1/7))/195
    gradce_chk_eval[4, 2] = 5129/130
    gradce_chk_eval[5, 1] = (4918*cos(23/10)*sin(1/7))/195
    gradce_chk_eval[5, 2] = 5129/130

    # gradci:
    gradci_chk_eval = Matrix{Float64}(undef, 5, 3)
    gradci_chk_eval[1, 1] = 1
    gradci_chk_eval[1, 2] = 0
    gradci_chk_eval[1, 3] = 2029289/5096
    gradci_chk_eval[2, 1] = 29241/4900
    gradci_chk_eval[2, 2] = 0
    gradci_chk_eval[2, 3] = -2029289/5096
    gradci_chk_eval[3, 1] = 0
    gradci_chk_eval[3, 2] = 1
    gradci_chk_eval[3, 3] = -934841/5880
    gradci_chk_eval[4, 1] = 0
    gradci_chk_eval[4, 2] = 22801/980
    gradci_chk_eval[4, 3] = 934841/5880
    gradci_chk_eval[5, 1] = 0
    gradci_chk_eval[5, 2] = 22801/63700
    gradci_chk_eval[5, 3] = -83200849/382200

    # gradF:
    gradF_chk_eval = Matrix{Vector{Float64}}(undef, 3, 2)
    gradF_chk_eval[1, 1] = Vector{Float64}(undef, 5)
    gradF_chk_eval[1, 1][1] = 121/15
    gradF_chk_eval[1, 1][2] = 121/15
    gradF_chk_eval[1, 1][3] = 0
    gradF_chk_eval[1, 1][4] = 0
    gradF_chk_eval[1, 1][5] = 0
    gradF_chk_eval[1, 2] = Vector{Float64}(undef, 5)
    gradF_chk_eval[1, 2][1] = 0
    gradF_chk_eval[1, 2][2] = 0
    gradF_chk_eval[1, 2][3] = 22801/4900
    gradF_chk_eval[1, 2][4] = 22801/4900
    gradF_chk_eval[1, 2][5] = 22801/4900
    gradF_chk_eval[2, 1] = Vector{Float64}(undef, 5)
    gradF_chk_eval[2, 1][1] = 0
    gradF_chk_eval[2, 1][2] = 0
    gradF_chk_eval[2, 1][3] = -151/70
    gradF_chk_eval[2, 1][4] = -151/70
    gradF_chk_eval[2, 1][5] = 0
    gradF_chk_eval[2, 2] = Vector{Float64}(undef, 5)
    gradF_chk_eval[2, 2][1] = -1
    gradF_chk_eval[2, 2][2] = 0
    gradF_chk_eval[2, 2][3] = 0
    gradF_chk_eval[2, 2][4] = 1
    gradF_chk_eval[2, 2][5] = 0
    gradF_chk_eval[3, 1] = Vector{Float64}(undef, 5)
    gradF_chk_eval[3, 1][1] = 0
    gradF_chk_eval[3, 1][2] = 0
    gradF_chk_eval[3, 1][3] = 0
    gradF_chk_eval[3, 1][4] = 0
    gradF_chk_eval[3, 1][5] = 326592
    gradF_chk_eval[3, 2] = Vector{Float64}(undef, 5)
    gradF_chk_eval[3, 2][1] = (241*cos(4/3))/70
    gradF_chk_eval[3, 2][2] = 0
    gradF_chk_eval[3, 2][3] = 0
    gradF_chk_eval[3, 2][4] = 0
    gradF_chk_eval[3, 2][5] = 0

    # hessf:
    hessf_chk_eval = Matrix{Float64}(undef, 5, 5)
    hessf_chk_eval[1, 1] = (12167*exp(-16342/195))/12250 + 141009897735/95051008
    hessf_chk_eval[1, 2] = (36501*exp(-16342/195))/24500 + 26112944025/11881376
    hessf_chk_eval[1, 3] = (12167*exp(-16342/195))/4900 + 6714757035/2970344
    hessf_chk_eval[1, 4] = (12167*exp(-16342/195))/3500 + 235016496225/1827904
    hessf_chk_eval[1, 5] = (133837*exp(-16342/195))/24500 + 28201979547/11881376
    hessf_chk_eval[2, 1] = (36501*exp(-16342/195))/24500 + 26112944025/11881376
    hessf_chk_eval[2, 2] = (109503*exp(-16342/195))/49000 + 1611910125/1485172
    hessf_chk_eval[2, 3] = (36501*exp(-16342/195))/9800 + 1243473525/742586
    hessf_chk_eval[2, 4] = (36501*exp(-16342/195))/7000 + 43521573375/456976
    hessf_chk_eval[2, 5] = (401511*exp(-16342/195))/49000 + 5222588805/2970344
    hessf_chk_eval[3, 1] = (12167*exp(-16342/195))/4900 + 6714757035/2970344
    hessf_chk_eval[3, 2] = (36501*exp(-16342/195))/9800 + 1243473525/742586
    hessf_chk_eval[3, 3] = (12167*exp(-16342/195))/1960 + 959251005/742586
    hessf_chk_eval[3, 4] = (12167*exp(-16342/195))/1400 + 11191261725/114244
    hessf_chk_eval[3, 5] = (133837*exp(-16342/195))/9800 + 1342951407/742586
    hessf_chk_eval[4, 1] = (12167*exp(-16342/195))/3500 + 235016496225/1827904
    hessf_chk_eval[4, 2] = (36501*exp(-16342/195))/7000 + 43521573375/456976
    hessf_chk_eval[4, 3] = (12167*exp(-16342/195))/1400 + 11191261725/114244
    hessf_chk_eval[4, 4] = (12167*exp(-16342/195))/1000 + 78338832075/17576
    hessf_chk_eval[4, 5] = (133837*exp(-16342/195))/7000 + 47003299245/456976
    hessf_chk_eval[5, 1] = (133837*exp(-16342/195))/24500 + 28201979547/11881376
    hessf_chk_eval[5, 2] = (401511*exp(-16342/195))/49000 + 5222588805/2970344
    hessf_chk_eval[5, 3] = (133837*exp(-16342/195))/9800 + 1342951407/742586
    hessf_chk_eval[5, 4] = (133837*exp(-16342/195))/7000 + 47003299245/456976
    hessf_chk_eval[5, 5] = (1472207*exp(-16342/195))/49000 + 9400659849/5940688

    # hessce:
    hessce_chk_eval = Vector{Matrix{Float64}}(undef, 2)
    hessce_chk_eval[1] = Matrix{Float64}(undef, 5, 5)
    hessce_chk_eval[1][1, 1] = 2*cos(23/10)*sin(1/7)
    hessce_chk_eval[1][1, 2] = 2*cos(23/10)*sin(1/7)
    hessce_chk_eval[1][1, 3] = 2*cos(23/10)*sin(1/7)
    hessce_chk_eval[1][1, 4] = 2*cos(23/10)*sin(1/7)
    hessce_chk_eval[1][1, 5] = 2*cos(23/10)*sin(1/7)
    hessce_chk_eval[1][2, 1] = 2*cos(23/10)*sin(1/7)
    hessce_chk_eval[1][2, 2] = 2*cos(23/10)*sin(1/7)
    hessce_chk_eval[1][2, 3] = 2*cos(23/10)*sin(1/7)
    hessce_chk_eval[1][2, 4] = 2*cos(23/10)*sin(1/7)
    hessce_chk_eval[1][2, 5] = 2*cos(23/10)*sin(1/7)
    hessce_chk_eval[1][3, 1] = 2*cos(23/10)*sin(1/7)
    hessce_chk_eval[1][3, 2] = 2*cos(23/10)*sin(1/7)
    hessce_chk_eval[1][3, 3] = 2*cos(23/10)*sin(1/7)
    hessce_chk_eval[1][3, 4] = 2*cos(23/10)*sin(1/7)
    hessce_chk_eval[1][3, 5] = 2*cos(23/10)*sin(1/7)
    hessce_chk_eval[1][4, 1] = 2*cos(23/10)*sin(1/7)
    hessce_chk_eval[1][4, 2] = 2*cos(23/10)*sin(1/7)
    hessce_chk_eval[1][4, 3] = 2*cos(23/10)*sin(1/7)
    hessce_chk_eval[1][4, 4] = 2*cos(23/10)*sin(1/7)
    hessce_chk_eval[1][4, 5] = 2*cos(23/10)*sin(1/7)
    hessce_chk_eval[1][5, 1] = 2*cos(23/10)*sin(1/7)
    hessce_chk_eval[1][5, 2] = 2*cos(23/10)*sin(1/7)
    hessce_chk_eval[1][5, 3] = 2*cos(23/10)*sin(1/7)
    hessce_chk_eval[1][5, 4] = 2*cos(23/10)*sin(1/7)
    hessce_chk_eval[1][5, 5] = 2*cos(23/10)*sin(1/7)
    hessce_chk_eval[2] = Matrix{Float64}(undef, 5, 5)
    hessce_chk_eval[2][1, 1] = 121/35
    hessce_chk_eval[2][1, 2] = 121/35
    hessce_chk_eval[2][1, 3] = 0
    hessce_chk_eval[2][1, 4] = 0
    hessce_chk_eval[2][1, 5] = 0
    hessce_chk_eval[2][2, 1] = 121/35
    hessce_chk_eval[2][2, 2] = 121/35
    hessce_chk_eval[2][2, 3] = 0
    hessce_chk_eval[2][2, 4] = 0
    hessce_chk_eval[2][2, 5] = 0
    hessce_chk_eval[2][3, 1] = 0
    hessce_chk_eval[2][3, 2] = 0
    hessce_chk_eval[2][3, 3] = 23/5
    hessce_chk_eval[2][3, 4] = 23/5
    hessce_chk_eval[2][3, 5] = 23/5
    hessce_chk_eval[2][4, 1] = 0
    hessce_chk_eval[2][4, 2] = 0
    hessce_chk_eval[2][4, 3] = 23/5
    hessce_chk_eval[2][4, 4] = 23/5
    hessce_chk_eval[2][4, 5] = 23/5
    hessce_chk_eval[2][5, 1] = 0
    hessce_chk_eval[2][5, 2] = 0
    hessce_chk_eval[2][5, 3] = 23/5
    hessce_chk_eval[2][5, 4] = 23/5
    hessce_chk_eval[2][5, 5] = 23/5

    # hessci:
    hessci_chk_eval = Vector{Matrix{Float64}}(undef, 3)
    hessci_chk_eval[1] = Matrix{Float64}(undef, 5, 5)
    hessci_chk_eval[1][1, 1] = 0
    hessci_chk_eval[1][1, 2] = 0
    hessci_chk_eval[1][1, 3] = 0
    hessci_chk_eval[1][1, 4] = 0
    hessci_chk_eval[1][1, 5] = 0
    hessci_chk_eval[1][2, 1] = 0
    hessci_chk_eval[1][2, 2] = 0
    hessci_chk_eval[1][2, 3] = 0
    hessci_chk_eval[1][2, 4] = 0
    hessci_chk_eval[1][2, 5] = 0
    hessci_chk_eval[1][3, 1] = 0
    hessci_chk_eval[1][3, 2] = 0
    hessci_chk_eval[1][3, 3] = 0
    hessci_chk_eval[1][3, 4] = 0
    hessci_chk_eval[1][3, 5] = 0
    hessci_chk_eval[1][4, 1] = 0
    hessci_chk_eval[1][4, 2] = 0
    hessci_chk_eval[1][4, 3] = 0
    hessci_chk_eval[1][4, 4] = 0
    hessci_chk_eval[1][4, 5] = 0
    hessci_chk_eval[1][5, 1] = 0
    hessci_chk_eval[1][5, 2] = 0
    hessci_chk_eval[1][5, 3] = 0
    hessci_chk_eval[1][5, 4] = 0
    hessci_chk_eval[1][5, 5] = 0
    hessci_chk_eval[2] = Matrix{Float64}(undef, 5, 5)
    hessci_chk_eval[2][1, 1] = 0
    hessci_chk_eval[2][1, 2] = 0
    hessci_chk_eval[2][1, 3] = 0
    hessci_chk_eval[2][1, 4] = 0
    hessci_chk_eval[2][1, 5] = 0
    hessci_chk_eval[2][2, 1] = 0
    hessci_chk_eval[2][2, 2] = 0
    hessci_chk_eval[2][2, 3] = 0
    hessci_chk_eval[2][2, 4] = 0
    hessci_chk_eval[2][2, 5] = 0
    hessci_chk_eval[2][3, 1] = 0
    hessci_chk_eval[2][3, 2] = 0
    hessci_chk_eval[2][3, 3] = 0
    hessci_chk_eval[2][3, 4] = 0
    hessci_chk_eval[2][3, 5] = 0
    hessci_chk_eval[2][4, 1] = 0
    hessci_chk_eval[2][4, 2] = 0
    hessci_chk_eval[2][4, 3] = 0
    hessci_chk_eval[2][4, 4] = 0
    hessci_chk_eval[2][4, 5] = 22801/4900
    hessci_chk_eval[2][5, 1] = 0
    hessci_chk_eval[2][5, 2] = 0
    hessci_chk_eval[2][5, 3] = 0
    hessci_chk_eval[2][5, 4] = 22801/4900
    hessci_chk_eval[2][5, 5] = 0
    hessci_chk_eval[3] = Matrix{Float64}(undef, 5, 5)
    hessci_chk_eval[3][1, 1] = 0
    hessci_chk_eval[3][1, 2] = 0
    hessci_chk_eval[3][1, 3] = 22801/196
    hessci_chk_eval[3][1, 4] = -22801/196
    hessci_chk_eval[3][1, 5] = 2029289/12740
    hessci_chk_eval[3][2, 1] = 0
    hessci_chk_eval[3][2, 2] = 0
    hessci_chk_eval[3][2, 3] = -22801/196
    hessci_chk_eval[3][2, 4] = 22801/196
    hessci_chk_eval[3][2, 5] = -2029289/12740
    hessci_chk_eval[3][3, 1] = 22801/196
    hessci_chk_eval[3][3, 2] = -22801/196
    hessci_chk_eval[3][3, 3] = 0
    hessci_chk_eval[3][3, 4] = 0
    hessci_chk_eval[3][3, 5] = -934841/14700
    hessci_chk_eval[3][4, 1] = -22801/196
    hessci_chk_eval[3][4, 2] = 22801/196
    hessci_chk_eval[3][4, 3] = 0
    hessci_chk_eval[3][4, 4] = 0
    hessci_chk_eval[3][4, 5] = 934841/14700
    hessci_chk_eval[3][5, 1] = 2029289/12740
    hessci_chk_eval[3][5, 2] = -2029289/12740
    hessci_chk_eval[3][5, 3] = -934841/14700
    hessci_chk_eval[3][5, 4] = 934841/14700
    hessci_chk_eval[3][5, 5] = -83200849/1911000

    # hessF:
    hessF_chk_eval = Matrix{Matrix{Float64}}(undef, 3, 2)
    hessF_chk_eval[1, 1] = Matrix{Float64}(undef, 5, 5)
    hessF_chk_eval[1, 1][1, 1] = 2
    hessF_chk_eval[1, 1][1, 2] = 2
    hessF_chk_eval[1, 1][1, 3] = 0
    hessF_chk_eval[1, 1][1, 4] = 0
    hessF_chk_eval[1, 1][1, 5] = 0
    hessF_chk_eval[1, 1][2, 1] = 2
    hessF_chk_eval[1, 1][2, 2] = 2
    hessF_chk_eval[1, 1][2, 3] = 0
    hessF_chk_eval[1, 1][2, 4] = 0
    hessF_chk_eval[1, 1][2, 5] = 0
    hessF_chk_eval[1, 1][3, 1] = 0
    hessF_chk_eval[1, 1][3, 2] = 0
    hessF_chk_eval[1, 1][3, 3] = 0
    hessF_chk_eval[1, 1][3, 4] = 0
    hessF_chk_eval[1, 1][3, 5] = 0
    hessF_chk_eval[1, 1][4, 1] = 0
    hessF_chk_eval[1, 1][4, 2] = 0
    hessF_chk_eval[1, 1][4, 3] = 0
    hessF_chk_eval[1, 1][4, 4] = 0
    hessF_chk_eval[1, 1][4, 5] = 0
    hessF_chk_eval[1, 1][5, 1] = 0
    hessF_chk_eval[1, 1][5, 2] = 0
    hessF_chk_eval[1, 1][5, 3] = 0
    hessF_chk_eval[1, 1][5, 4] = 0
    hessF_chk_eval[1, 1][5, 5] = 0
    hessF_chk_eval[1, 2] = Matrix{Float64}(undef, 5, 5)
    hessF_chk_eval[1, 2][1, 1] = 0
    hessF_chk_eval[1, 2][1, 2] = 0
    hessF_chk_eval[1, 2][1, 3] = 0
    hessF_chk_eval[1, 2][1, 4] = 0
    hessF_chk_eval[1, 2][1, 5] = 0
    hessF_chk_eval[1, 2][2, 1] = 0
    hessF_chk_eval[1, 2][2, 2] = 0
    hessF_chk_eval[1, 2][2, 3] = 0
    hessF_chk_eval[1, 2][2, 4] = 0
    hessF_chk_eval[1, 2][2, 5] = 0
    hessF_chk_eval[1, 2][3, 1] = 0
    hessF_chk_eval[1, 2][3, 2] = 0
    hessF_chk_eval[1, 2][3, 3] = 0
    hessF_chk_eval[1, 2][3, 4] = 0
    hessF_chk_eval[1, 2][3, 5] = 0
    hessF_chk_eval[1, 2][4, 1] = 0
    hessF_chk_eval[1, 2][4, 2] = 0
    hessF_chk_eval[1, 2][4, 3] = 0
    hessF_chk_eval[1, 2][4, 4] = 0
    hessF_chk_eval[1, 2][4, 5] = 0
    hessF_chk_eval[1, 2][5, 1] = 0
    hessF_chk_eval[1, 2][5, 2] = 0
    hessF_chk_eval[1, 2][5, 3] = 0
    hessF_chk_eval[1, 2][5, 4] = 0
    hessF_chk_eval[1, 2][5, 5] = 0
    hessF_chk_eval[2, 1] = Matrix{Float64}(undef, 5, 5)
    hessF_chk_eval[2, 1][1, 1] = 0
    hessF_chk_eval[2, 1][1, 2] = 0
    hessF_chk_eval[2, 1][1, 3] = 0
    hessF_chk_eval[2, 1][1, 4] = 0
    hessF_chk_eval[2, 1][1, 5] = 0
    hessF_chk_eval[2, 1][2, 1] = 0
    hessF_chk_eval[2, 1][2, 2] = 0
    hessF_chk_eval[2, 1][2, 3] = 0
    hessF_chk_eval[2, 1][2, 4] = 0
    hessF_chk_eval[2, 1][2, 5] = 0
    hessF_chk_eval[2, 1][3, 1] = 0
    hessF_chk_eval[2, 1][3, 2] = 0
    hessF_chk_eval[2, 1][3, 3] = 0
    hessF_chk_eval[2, 1][3, 4] = 0
    hessF_chk_eval[2, 1][3, 5] = 0
    hessF_chk_eval[2, 1][4, 1] = 0
    hessF_chk_eval[2, 1][4, 2] = 0
    hessF_chk_eval[2, 1][4, 3] = 0
    hessF_chk_eval[2, 1][4, 4] = 0
    hessF_chk_eval[2, 1][4, 5] = 0
    hessF_chk_eval[2, 1][5, 1] = 0
    hessF_chk_eval[2, 1][5, 2] = 0
    hessF_chk_eval[2, 1][5, 3] = 0
    hessF_chk_eval[2, 1][5, 4] = 0
    hessF_chk_eval[2, 1][5, 5] = 0
    hessF_chk_eval[2, 2] = Matrix{Float64}(undef, 5, 5)
    hessF_chk_eval[2, 2][1, 1] = 0
    hessF_chk_eval[2, 2][1, 2] = 0
    hessF_chk_eval[2, 2][1, 3] = 0
    hessF_chk_eval[2, 2][1, 4] = 0
    hessF_chk_eval[2, 2][1, 5] = 0
    hessF_chk_eval[2, 2][2, 1] = 0
    hessF_chk_eval[2, 2][2, 2] = 0
    hessF_chk_eval[2, 2][2, 3] = 0
    hessF_chk_eval[2, 2][2, 4] = 0
    hessF_chk_eval[2, 2][2, 5] = 0
    hessF_chk_eval[2, 2][3, 1] = 0
    hessF_chk_eval[2, 2][3, 2] = 0
    hessF_chk_eval[2, 2][3, 3] = 0
    hessF_chk_eval[2, 2][3, 4] = 0
    hessF_chk_eval[2, 2][3, 5] = 0
    hessF_chk_eval[2, 2][4, 1] = 0
    hessF_chk_eval[2, 2][4, 2] = 0
    hessF_chk_eval[2, 2][4, 3] = 0
    hessF_chk_eval[2, 2][4, 4] = 0
    hessF_chk_eval[2, 2][4, 5] = 0
    hessF_chk_eval[2, 2][5, 1] = 0
    hessF_chk_eval[2, 2][5, 2] = 0
    hessF_chk_eval[2, 2][5, 3] = 0
    hessF_chk_eval[2, 2][5, 4] = 0
    hessF_chk_eval[2, 2][5, 5] = 0
    hessF_chk_eval[3, 1] = Matrix{Float64}(undef, 5, 5)
    hessF_chk_eval[3, 1][1, 1] = 0
    hessF_chk_eval[3, 1][1, 2] = 0
    hessF_chk_eval[3, 1][1, 3] = 0
    hessF_chk_eval[3, 1][1, 4] = 0
    hessF_chk_eval[3, 1][1, 5] = 0
    hessF_chk_eval[3, 1][2, 1] = 0
    hessF_chk_eval[3, 1][2, 2] = 0
    hessF_chk_eval[3, 1][2, 3] = 0
    hessF_chk_eval[3, 1][2, 4] = 0
    hessF_chk_eval[3, 1][2, 5] = 0
    hessF_chk_eval[3, 1][3, 1] = 0
    hessF_chk_eval[3, 1][3, 2] = 0
    hessF_chk_eval[3, 1][3, 3] = 0
    hessF_chk_eval[3, 1][3, 4] = 0
    hessF_chk_eval[3, 1][3, 5] = 0
    hessF_chk_eval[3, 1][4, 1] = 0
    hessF_chk_eval[3, 1][4, 2] = 0
    hessF_chk_eval[3, 1][4, 3] = 0
    hessF_chk_eval[3, 1][4, 4] = 0
    hessF_chk_eval[3, 1][4, 5] = 0
    hessF_chk_eval[3, 1][5, 1] = 0
    hessF_chk_eval[3, 1][5, 2] = 0
    hessF_chk_eval[3, 1][5, 3] = 0
    hessF_chk_eval[3, 1][5, 4] = 0
    hessF_chk_eval[3, 1][5, 5] = 326592
    hessF_chk_eval[3, 2] = Matrix{Float64}(undef, 5, 5)
    hessF_chk_eval[3, 2][1, 1] = -(241*sin(4/3))/70
    hessF_chk_eval[3, 2][1, 2] = 0
    hessF_chk_eval[3, 2][1, 3] = 0
    hessF_chk_eval[3, 2][1, 4] = 0
    hessF_chk_eval[3, 2][1, 5] = 0
    hessF_chk_eval[3, 2][2, 1] = 0
    hessF_chk_eval[3, 2][2, 2] = 0
    hessF_chk_eval[3, 2][2, 3] = 0
    hessF_chk_eval[3, 2][2, 4] = 0
    hessF_chk_eval[3, 2][2, 5] = 0
    hessF_chk_eval[3, 2][3, 1] = 0
    hessF_chk_eval[3, 2][3, 2] = 0
    hessF_chk_eval[3, 2][3, 3] = 0
    hessF_chk_eval[3, 2][3, 4] = 0
    hessF_chk_eval[3, 2][3, 5] = 0
    hessF_chk_eval[3, 2][4, 1] = 0
    hessF_chk_eval[3, 2][4, 2] = 0
    hessF_chk_eval[3, 2][4, 3] = 0
    hessF_chk_eval[3, 2][4, 4] = 0
    hessF_chk_eval[3, 2][4, 5] = 0
    hessF_chk_eval[3, 2][5, 1] = 0
    hessF_chk_eval[3, 2][5, 2] = 0
    hessF_chk_eval[3, 2][5, 3] = 0
    hessF_chk_eval[3, 2][5, 4] = 0
    hessF_chk_eval[3, 2][5, 5] = 0

    # fdp:
    fdp_chk_eval = Vector{Float64}(undef, 2)
    fdp_chk_eval[1] = (12167*exp(-16342/195))/3500 + 3206852775/2970344
    fdp_chk_eval[2] = (1587*exp(-16342/195))/4900 + 3206852775/2970344

    # cedp:
    cedp_chk_eval = Vector{Vector{Float64}}(undef, 2)
    cedp_chk_eval[1] = Vector{Float64}(undef, 2)
    cedp_chk_eval[1][1] = (6046681*cos(1/7)*cos(23/10))/38025
    cedp_chk_eval[1][2] = 1771561/27000
    cedp_chk_eval[2] = Vector{Float64}(undef, 2)
    cedp_chk_eval[2][1] = -(6046681*sin(1/7)*sin(23/10))/38025
    cedp_chk_eval[2][2] = 49729/676

    # cidp:
    cidp_chk_eval = Vector{Vector{Float64}}(undef, 2)
    cidp_chk_eval[1] = Vector{Float64}(undef, 3)
    cidp_chk_eval[1][1] = 4617/350
    cidp_chk_eval[1][2] = -151/91
    cidp_chk_eval[1][3] = 550999/1092
    cidp_chk_eval[2] = Vector{Float64}(undef, 3)
    cidp_chk_eval[2][1] = 4617/350
    cidp_chk_eval[2][2] = 151/91
    cidp_chk_eval[2][3] = -550999/1092

    # Fdp:
    Fdp_chk_eval = Vector{Matrix{Float64}}(undef, 2)
    Fdp_chk_eval[1] = Matrix{Float64}(undef, 3, 2)
    Fdp_chk_eval[1][1, 1] = 0
    Fdp_chk_eval[1][1, 2] = -33673/910
    Fdp_chk_eval[1][2, 1] = 93/26
    Fdp_chk_eval[1][2, 2] = 0
    Fdp_chk_eval[1][3, 1] = 0
    Fdp_chk_eval[1][3, 2] = sin(4/3)
    Fdp_chk_eval[2] = Matrix{Float64}(undef, 3, 2)
    Fdp_chk_eval[2][1, 1] = 0
    Fdp_chk_eval[2][1, 2] = 33673/910
    Fdp_chk_eval[2][2, 1] = -93/26
    Fdp_chk_eval[2][2, 2] = 0
    Fdp_chk_eval[2][3, 1] = 0
    Fdp_chk_eval[2][3, 2] = sin(4/3)

    # gradfdp:
    gradfdp_chk_eval = Vector{Vector{Float64}}(undef, 2)
    gradfdp_chk_eval[1] = Vector{Float64}(undef, 5)
    gradfdp_chk_eval[1][1] = 9620558325/5940688 - (12167*exp(-16342/195))/1750
    gradfdp_chk_eval[1][2] = 1781584875/1485172 - (36501*exp(-16342/195))/3500
    gradfdp_chk_eval[1][3] = 458121825/371293 - (12167*exp(-16342/195))/700
    gradfdp_chk_eval[1][4] = 16034263875/228488 - (12167*exp(-16342/195))/500
    gradfdp_chk_eval[1][5] = 1924111665/1485172 - (133837*exp(-16342/195))/3500
    gradfdp_chk_eval[2] = Vector{Float64}(undef, 5)
    gradfdp_chk_eval[2][1] = 9620558325/5940688 - (1587*exp(-16342/195))/2450
    gradfdp_chk_eval[2][2] = 1781584875/1485172 - (4761*exp(-16342/195))/4900
    gradfdp_chk_eval[2][3] = 458121825/371293 - (1587*exp(-16342/195))/980
    gradfdp_chk_eval[2][4] = 16034263875/228488 - (1587*exp(-16342/195))/700
    gradfdp_chk_eval[2][5] = 1924111665/1485172 - (17457*exp(-16342/195))/4900

    # gradcedp:
    gradcedp_chk_eval = Vector{Matrix{Float64}}(undef, 2)
    gradcedp_chk_eval[1] = Matrix{Float64}(undef, 5, 2)
    gradcedp_chk_eval[1][1, 1] = (4918*cos(1/7)*cos(23/10))/195
    gradcedp_chk_eval[1][1, 2] = 14641/300
    gradcedp_chk_eval[1][2, 1] = (4918*cos(1/7)*cos(23/10))/195
    gradcedp_chk_eval[1][2, 2] = 14641/300
    gradcedp_chk_eval[1][3, 1] = (4918*cos(1/7)*cos(23/10))/195
    gradcedp_chk_eval[1][3, 2] = 0
    gradcedp_chk_eval[1][4, 1] = (4918*cos(1/7)*cos(23/10))/195
    gradcedp_chk_eval[1][4, 2] = 0
    gradcedp_chk_eval[1][5, 1] = (4918*cos(1/7)*cos(23/10))/195
    gradcedp_chk_eval[1][5, 2] = 0
    gradcedp_chk_eval[2] = Matrix{Float64}(undef, 5, 2)
    gradcedp_chk_eval[2][1, 1] = -(4918*sin(1/7)*sin(23/10))/195
    gradcedp_chk_eval[2][1, 2] = 0
    gradcedp_chk_eval[2][2, 1] = -(4918*sin(1/7)*sin(23/10))/195
    gradcedp_chk_eval[2][2, 2] = 0
    gradcedp_chk_eval[2][3, 1] = -(4918*sin(1/7)*sin(23/10))/195
    gradcedp_chk_eval[2][3, 2] = 223/13
    gradcedp_chk_eval[2][4, 1] = -(4918*sin(1/7)*sin(23/10))/195
    gradcedp_chk_eval[2][4, 2] = 223/13
    gradcedp_chk_eval[2][5, 1] = -(4918*sin(1/7)*sin(23/10))/195
    gradcedp_chk_eval[2][5, 2] = 223/13

    # gradcidp:
    gradcidp_chk_eval = Vector{Matrix{Float64}}(undef, 2)
    gradcidp_chk_eval[1] = Matrix{Float64}(undef, 5, 3)
    gradcidp_chk_eval[1][1, 1] = 0
    gradcidp_chk_eval[1][1, 2] = 0
    gradcidp_chk_eval[1][1, 3] = -67195/182
    gradcidp_chk_eval[1][2, 1] = 171/35
    gradcidp_chk_eval[1][2, 2] = 0
    gradcidp_chk_eval[1][2, 3] = 67195/182
    gradcidp_chk_eval[1][3, 1] = 0
    gradcidp_chk_eval[1][3, 2] = 0
    gradcidp_chk_eval[1][3, 3] = 6191/42
    gradcidp_chk_eval[1][4, 1] = 0
    gradcidp_chk_eval[1][4, 2] = -151/7
    gradcidp_chk_eval[1][4, 3] = -6191/42
    gradcidp_chk_eval[1][5, 1] = 0
    gradcidp_chk_eval[1][5, 2] = -151/455
    gradcidp_chk_eval[1][5, 3] = 550999/2730
    gradcidp_chk_eval[2] = Matrix{Float64}(undef, 5, 3)
    gradcidp_chk_eval[2][1, 1] = 0
    gradcidp_chk_eval[2][1, 2] = 0
    gradcidp_chk_eval[2][1, 3] = 67195/182
    gradcidp_chk_eval[2][2, 1] = 171/35
    gradcidp_chk_eval[2][2, 2] = 0
    gradcidp_chk_eval[2][2, 3] = -67195/182
    gradcidp_chk_eval[2][3, 1] = 0
    gradcidp_chk_eval[2][3, 2] = 0
    gradcidp_chk_eval[2][3, 3] = -6191/42
    gradcidp_chk_eval[2][4, 1] = 0
    gradcidp_chk_eval[2][4, 2] = 151/7
    gradcidp_chk_eval[2][4, 3] = 6191/42
    gradcidp_chk_eval[2][5, 1] = 0
    gradcidp_chk_eval[2][5, 2] = 151/455
    gradcidp_chk_eval[2][5, 3] = -550999/2730

    # gradFdp:
    gradFdp_chk_eval = Vector{Matrix{Vector{Float64}}}(undef, 2)
    gradFdp_chk_eval[1] = Matrix{Vector{Float64}}(undef, 3, 2)
    gradFdp_chk_eval[1][1, 1] = Vector{Float64}(undef, 5)
    gradFdp_chk_eval[1][1, 1][1] = 0
    gradFdp_chk_eval[1][1, 1][2] = 0
    gradFdp_chk_eval[1][1, 1][3] = 0
    gradFdp_chk_eval[1][1, 1][4] = 0
    gradFdp_chk_eval[1][1, 1][5] = 0
    gradFdp_chk_eval[1][1, 2] = Vector{Float64}(undef, 5)
    gradFdp_chk_eval[1][1, 2][1] = 0
    gradFdp_chk_eval[1][1, 2][2] = 0
    gradFdp_chk_eval[1][1, 2][3] = -151/35
    gradFdp_chk_eval[1][1, 2][4] = -151/35
    gradFdp_chk_eval[1][1, 2][5] = -151/35
    gradFdp_chk_eval[1][2, 1] = Vector{Float64}(undef, 5)
    gradFdp_chk_eval[1][2, 1][1] = 0
    gradFdp_chk_eval[1][2, 1][2] = 0
    gradFdp_chk_eval[1][2, 1][3] = 1
    gradFdp_chk_eval[1][2, 1][4] = 1
    gradFdp_chk_eval[1][2, 1][5] = 0
    gradFdp_chk_eval[1][2, 2] = Vector{Float64}(undef, 5)
    gradFdp_chk_eval[1][2, 2][1] = 0
    gradFdp_chk_eval[1][2, 2][2] = 0
    gradFdp_chk_eval[1][2, 2][3] = 0
    gradFdp_chk_eval[1][2, 2][4] = 0
    gradFdp_chk_eval[1][2, 2][5] = 0
    gradFdp_chk_eval[1][3, 1] = Vector{Float64}(undef, 5)
    gradFdp_chk_eval[1][3, 1][1] = 0
    gradFdp_chk_eval[1][3, 1][2] = 0
    gradFdp_chk_eval[1][3, 1][3] = 0
    gradFdp_chk_eval[1][3, 1][4] = 0
    gradFdp_chk_eval[1][3, 1][5] = 0
    gradFdp_chk_eval[1][3, 2] = Vector{Float64}(undef, 5)
    gradFdp_chk_eval[1][3, 2][1] = cos(4/3)
    gradFdp_chk_eval[1][3, 2][2] = 0
    gradFdp_chk_eval[1][3, 2][3] = 0
    gradFdp_chk_eval[1][3, 2][4] = 0
    gradFdp_chk_eval[1][3, 2][5] = 0
    gradFdp_chk_eval[2] = Matrix{Vector{Float64}}(undef, 3, 2)
    gradFdp_chk_eval[2][1, 1] = Vector{Float64}(undef, 5)
    gradFdp_chk_eval[2][1, 1][1] = 0
    gradFdp_chk_eval[2][1, 1][2] = 0
    gradFdp_chk_eval[2][1, 1][3] = 0
    gradFdp_chk_eval[2][1, 1][4] = 0
    gradFdp_chk_eval[2][1, 1][5] = 0
    gradFdp_chk_eval[2][1, 2] = Vector{Float64}(undef, 5)
    gradFdp_chk_eval[2][1, 2][1] = 0
    gradFdp_chk_eval[2][1, 2][2] = 0
    gradFdp_chk_eval[2][1, 2][3] = 151/35
    gradFdp_chk_eval[2][1, 2][4] = 151/35
    gradFdp_chk_eval[2][1, 2][5] = 151/35
    gradFdp_chk_eval[2][2, 1] = Vector{Float64}(undef, 5)
    gradFdp_chk_eval[2][2, 1][1] = 0
    gradFdp_chk_eval[2][2, 1][2] = 0
    gradFdp_chk_eval[2][2, 1][3] = -1
    gradFdp_chk_eval[2][2, 1][4] = -1
    gradFdp_chk_eval[2][2, 1][5] = 0
    gradFdp_chk_eval[2][2, 2] = Vector{Float64}(undef, 5)
    gradFdp_chk_eval[2][2, 2][1] = 0
    gradFdp_chk_eval[2][2, 2][2] = 0
    gradFdp_chk_eval[2][2, 2][3] = 0
    gradFdp_chk_eval[2][2, 2][4] = 0
    gradFdp_chk_eval[2][2, 2][5] = 0
    gradFdp_chk_eval[2][3, 1] = Vector{Float64}(undef, 5)
    gradFdp_chk_eval[2][3, 1][1] = 0
    gradFdp_chk_eval[2][3, 1][2] = 0
    gradFdp_chk_eval[2][3, 1][3] = 0
    gradFdp_chk_eval[2][3, 1][4] = 0
    gradFdp_chk_eval[2][3, 1][5] = 0
    gradFdp_chk_eval[2][3, 2] = Vector{Float64}(undef, 5)
    gradFdp_chk_eval[2][3, 2][1] = cos(4/3)
    gradFdp_chk_eval[2][3, 2][2] = 0
    gradFdp_chk_eval[2][3, 2][3] = 0
    gradFdp_chk_eval[2][3, 2][4] = 0
    gradFdp_chk_eval[2][3, 2][5] = 0

    eval_test_pt = MPCCPointEval{Float64}(
        f_chk_eval,
        ce_chk_eval,
        ci_chk_eval,
        F_chk_eval,
        gradf_chk_eval,
        gradce_chk_eval,
        gradci_chk_eval,
        gradF_chk_eval,
        hessf_chk_eval,
        hessce_chk_eval,
        hessci_chk_eval,
        hessF_chk_eval,
        fdp_chk_eval,
        cedp_chk_eval,
        cidp_chk_eval,
        Fdp_chk_eval,
        gradfdp_chk_eval,
        gradcedp_chk_eval,
        gradcidp_chk_eval,
        gradFdp_chk_eval
    )

    return [ MPCCModelTestVector{Float64, Float64, Float64}(x_val, pr_val, ps_val, eval_test_pt) ]

end

function mm_spec_testdense_knownsols()
    return Vector{Function}(undef, 0)
end

function mm_spec_testdense_parameterisations()
    return Vector{MPCCModelParameterisation}(undef, 0)
end


text:
abstract type AbstractScalarTraceFunctionSpace{dim,T,FE,fdim} <: DiscreteFunctionSpace{dim,T,FE} end
# Scalar Trace Function Space (Scalar functions defined only on cell boundaries)
struct ScalarTraceFunctionSpace{dim,T<:Real,FE<:FiniteElement,fdim, FE2 <: FiniteElement,MM,N1,N2,N3} <: AbstractScalarTraceFunctionSpace{dim,T,FE,fdim}
    N::Matrix{T}
    dNdξ::Matrix{Vec{fdim,T}}
    fe::FE
    fs::ScalarFunctionSpace{dim,T, FE2, MM,N1,N2,N3,fdim}
end

#Constructor
function ScalarTraceFunctionSpace(psp::ScalarFunctionSpace{dim,T, FE, MM,N1,N2,N3,fdim},
    felem::FiniteElement{fdim,shape,order,gorder}) where {dim, T, FE, MM, N1,N2,N3,fdim,shape,order,gorder}
    qr_points = psp.qr_face_points
    qr_weights = psp.qr_face_weigths
    n_qpoints = length(qr_weights)

    # Function interpolation
    n_func_basefuncs = getnbasefunctions(felem)
    N    = fill(zero(T)          * T(NaN), n_func_basefuncs, n_qpoints)
    dNdξ = fill(zero(Vec{fdim,T}) * T(NaN), n_func_basefuncs, n_qpoints)

    for (qp, ξ) in enumerate(qr_points)
        for i in 1:n_func_basefuncs
            dNdξ[i, qp], N[i, qp]  = gradient(ξ -> value(felem, i, ξ), ξ, :all)
        end
    end
    ScalarTraceFunctionSpace{dim,T,typeof(felem),fdim, FE, MM,N1,N2,N3}(N, dNdξ,felem, psp)
end

#Data
@inline getnbasefunctions(fs::ScalarTraceFunctionSpace) = size(fs.N,1)
@inline getnquadpoints(fs::ScalarTraceFunctionSpace) = length(fs.fs.qr_face_weigths)
#@inline getdetJdS(fs::ScalarTraceFunctionSpace{dim,T,2}, cell::Int, face::Int, q_point::Int) where {dim,T} = fs.detJ[cell,face,q_point]*fs.qr_weights[q_point]
@inline getfacedetJdS(fs::ScalarTraceFunctionSpace, face::Int, q_point::Int) = getfacedetJdS(fs.fs, face, q_point)
@inline shape_value(fs::ScalarTraceFunctionSpace, q_point::Int, base_func::Int) = fs.N[base_func, q_point]
@inline getngeobasefunctions(fs::ScalarTraceFunctionSpace) = size(fs.fs.L,1)
@inline geometric_value(fs::ScalarTraceFunctionSpace, face::Int, q_point::Int, base_func::Int) = fs.fs.L[base_func, q_point, face]
@inline getnlocaldofs(fs::ScalarTraceFunctionSpace) = getnbasefunctions(fs)*size(fs.fs.L,3)
@inline getmesh(fs::ScalarTraceFunctionSpace) = getmesh(fs.fs)
@inline getfiniteelement(fs::ScalarTraceFunctionSpace) = fs.fe

"""
function spatial_coordinate(fs::ScalarTraceFunctionSpace{dim}, q_point::Int, x::AbstractVector{Vec{dim,T}}, orientation=true)
Map coordinates of quadrature point `q_point` of Scalar Trace Function Space `fs`
into domain with vertices `x`
"""
function spatial_coordinate(fs::ScalarTraceFunctionSpace{dim}, face::Int, q_point::Int, x::AbstractVector{Vec{dim,T}}, orientation=true) where {dim,T}
    n_base_funcs = getngeobasefunctions(fs)
    @assert length(x) == n_base_funcs
    vec = zero(Vec{dim,T})
    n = getnquadpoints(fs)
    @inbounds for i in 1:n_base_funcs
        or_q_point = orientation ? q_point : n - q_point + 1
        vec += geometric_value(fs, face, or_q_point, i) * x[i]
    end
    return vec
end


text:
#!/usr/bin/julia

# Trizen
# 15 November 2021
# https://github.com/trizen

# Compute the multiplicative order of `a` modulo `n`: znorder(a, n).
# This is the smallest positive integer k such that a^k == 1 (mod n).

using Primes

function divisors(n)

    d = Int64[1]

    for (p,e) in factor(n)
        t = Int64[]
        r = 1

        for i in 1:e
            r *= p
            for u in d
                push!(t, u*r)
            end
        end

        append!(d, t)
    end

    return sort(d)
end

function znorder(a, n)

    if isprime(n)
        for d in divisors(n-1)
            if (powermod(a, d, n) == 1)
                return d
            end
        end
    end

    f = factor(n)

    if (length(f) == 1)         # is prime power

        p = first(first(f))
        z = znorder(a, p)

        while (powermod(a, z, n) != 1)
            z *= p
        end

        return z
    end

    pp_orders = Int64[]

    for (p,e) in f
        push!(pp_orders, znorder(a, p^e))
    end

    return lcm(pp_orders)
end

isequal(znorder(97, factorial(14)), 25920)          || print("error")
isequal(znorder(53, factorial(15)), 2419200)        || print("error")
isequal(znorder(37, factorial(16)), 116121600)      || print("error")
isequal(znorder(31, factorial(17)), 6220800)        || print("error")
isequal(znorder(89, factorial(18)), 1045094400)     || print("error")
isequal(znorder(101,factorial(19)), 1254113280)     || print("error")
isequal(znorder(97, factorial(20)), 2239488000)     || print("error")

println(znorder(2, 341))        #> 10
println(znorder(97, 5040))      #> 12

max_stars_repo_name: trizen/julia-scripts

max_issues_repo_name: trizen/julia-scripts

max_forks_repo_name: trizen/julia-scripts


text:
#############################################################################
# Computes the bases that maxize a linear function over the base polyhedron
# using the greedy algorithm.
#############################################################################

export greedy, greedy_rand

TOL = 1e-2

# Computes one solution
function greedy(F::Function, x::AbstractArray,
                S = collect(1:length(x))::AbstractArray)
  @assert F([])[1] == 0 # f should be 0 at the empty set
  n = length(x)
  @assert length(S) == n
  i = sortperm(x, rev = true)
  V = sort(S)
  w = zeros(n)
  for ii = 1:n
    w[i[ii]] = F(V[i[1: ii]])[1] - F(V[i[1: ii-1]])[1]
  end
  return w
end

greedy(F::SubmodFunc, x::AbstractArray) = greedy(x -> evaluate(F, x), x)

greedy(f::LovaszExtAtom, x::AbstractArray) = greedy(f.func, x)

function greedy(f::LovaszExtAbsAtom, x::AbstractArray)
  w = greedy(f.func, abs.(x))
  for i = 1:length(w)
    if x[i] < 0
      w[i] = -w[i]
    end
  end
  return w
end

# Compute a solution with indexed on the same level set permutated
function greedy_rand(F::Function, x::AbstractArray, tol = 1e-3 ::Number,
                    S = collect(1:length(x))::AbstractArray)
  @assert F([])[1] == 0       # f should be 0 at the empty set
  n = length(x)
  @assert length(S) == n
  ordering = sortperm(x, rev = true)
  pernum = 1
  V = sort(S)
  ordering = V[ordering]      # the indexed reordered in a descending order
  # Mark level sets
  index = 1
  level_sets = zeros(2, n)    # the beginning and ending of each level set
  level_sets[1, 1] = 1
  for i = 2:n
    if x[ordering[i - 1]] - x[ordering[i]] >= tol
    level_sets[2, index] = i-1
    index += 1
    level_sets[1, index] = i
    end
  end
  level_sets[2, index] = n
  new_ordering = zeros(n)
  for i = 1:index
    new_ordering[Int(level_sets[1, i]) : Int(level_sets[2, i])] = shuffle(ordering[Int(level_sets[1, i]) : Int(level_sets[2, i])])
    pernum *= factorial(level_sets[2, i] - level_sets[1, i] + 1)
  end
  new_ordering = Int.(new_ordering)
  w = zeros(n)
  for ii = 1:n
    w[new_ordering[ii]] = F(V[new_ordering[1: ii]])[1] - F(V[new_ordering[1: ii-1]])[1]
  end
  return w, pernum, new_ordering
end

greedy_rand(F::SubmodFunc, x::AbstractArray, tol::Number) = greedy_rand(x -> evaluate(F, x), x, tol)


text:
# Tuples

struct ProbTuple <: Dist
    mgr
    left::Dist
    right::Dist
end

function prob_equals(x::ProbTuple, y::ProbTuple)
    if typeof(x.left) != typeof(y.left) || typeof(x.right) != typeof(y.right)
        # better to just define equality between different types to be false by default???
        false_constant(x.mgr)
    else
        # TODO order left and right checks by first fail principle
        left_eq = prob_equals(x.left, y.left)
        if !issat(left_eq)
            false_constant(x.mgr)
        else
            right_eq = prob_equals(x.right, y.right)
            left_eq & right_eq
        end
    end
end

function ite(cond::DistBool, then::T, elze::T) where {T <: ProbTuple}
    # TODO optimize of deterministic conditions
    left = ite(cond, then.left, elze.left)
    right = ite(cond, then.right, elze.right)
    ProbTuple(cond.mgr, left, right)
end

bools(t::ProbTuple) = [bools(t.left); bools(t.right)]


text:
using DelimitedFiles

content = readdlm("input.txt", '|')
nEntries = size(content)[1]

signals = split.(content[:, 1], " ", keepempty=false)
outputs = split.(content[:, 2], " ", keepempty=false)

function countP1()
    count = 0

    for i in 1:nEntries
        for output in outputs[i]
            len = length(output)

            if len == 2 || len ==  4 || len == 3 || len == 7
                count += 1
            end
        end
    end

    return count
end

countP1()

function countP2()
    value = 0

    numLengthMap = Dict(
        2 => ["1"],
        3 => ["7"],
        4 => ["4"],
        5 => ["2", "3", "5"],
        6 => ["0", "6", "9"],
        7 => ["8"]
    )

    for i in 1:nEntries
        numString = ""

        local overlapRef = Dict()

        for signal in signals[i]
            set = Set(signal)
            len = length(set)
            if len == 2 || len == 4
                overlapRef[len] = set
            end
        end

        for output in outputs[i]
            len = length(output)

            if len == 5
                digitSegs = Set(output)

                if length(digitSegs ∩ overlapRef[4]) == 2
                    numString *= numLengthMap[len][1] # 2
                elseif length(digitSegs ∩ overlapRef[2]) == 2
                    numString *= numLengthMap[len][2] # 3
                else
                    numString *= numLengthMap[len][3] # 5
                end

            elseif len == 6
                digitSegs = Set(output)

                if length(digitSegs ∩ overlapRef[2]) == 1
                    numString *= numLengthMap[len][2] # 6
                elseif length(digitSegs ∩ overlapRef[4]) == 4
                    numString *= numLengthMap[len][3] # 9
                else
                    numString *= numLengthMap[len][1] # 0
                end

            else
                numString *= numLengthMap[len][1]
            end

        end
        value += parse(Int, numString)
    end

    return value
end

countP2()

max_stars_repo_path: Day 8/main.jl

max_issues_repo_path: Day 8/main.jl

max_forks_repo_path: Day 8/main.jl