text:
using .PlotlyJS
import .PlotlyJS: plot
using LinearAlgebra: norm, dot
using StaticArrays

using .WignerSeitz: face_normal
# ---------------------------------------------------------------------------------------- #
# CONSTANTS

# default layout
const DEFAULT_PLOTLY_LAYOUT_3D  = Layout(
    showlegend=false,
    scene=attr(
        xaxis=attr(tickvals=[], zeroline=false,
                showgrid=false, showbackground=false,
                title=attr(text=""),
                ),
        yaxis=attr(tickvals=[], zeroline=false,
                showgrid=false, showbackground=false,
                title=attr(text=""),
                ),
        zaxis=attr(tickvals=[], zeroline=false,
                showgrid=false, showbackground=false,
                title=attr(text=""),
                ),
        aspectmode = "data",
        camera=attr(up = attr(x=0, z=1, y=0), center = attr(x=0, y=0, z=0)),
        dragmode = "turntable",
        ),
    margin=attr(l=0, r=0, b=0, t=0),
    autosize=false,
    #width=200, height=200,
    plot_bgcolor="rgba(255, 255, 255, 1)", paper_bgcolor="rgba(255, 255, 255, 1)",
    )

# colors from the "british" flatcolors color palette: https://flatuicolors.com/palette/gb
const BZ_COL          = Ref("rgb(47,54,64)")    # "electromagnetic"
const BASIS_COL       = Ref("rgb(39,60,117)")   # "pico void"
const BASIS_LIGHT_COL = Ref("rgb(212,216,227)") # 20% BASIS_COL, 80% white
const AXIS_COL        = Ref("rgb(194,54,22)")   # "harley davidson orange"
const AXIS_LIGHT_COL  = Ref("rgb(242,215,208)") # 20% AXIS_COL, 80% white

# ---------------------------------------------------------------------------------------- #
# 3D

function plot(c::Cell{3}, layout::Layout=DEFAULT_PLOTLY_LAYOUT_3D;
              config::PlotConfig = PlotConfig(responsive=true, displaylogo=false))

    setting(c) !== CARTESIAN && (c = cartesianize(c))
    scale = maximum(norm, basis(c))

    # BZ
    tbz = Vector{GenericTrace{Dict{Symbol,Any}}}(undef, length(c))
    for (i,poly) in enumerate(c)
        tbz[i] = PlotlyJS.scatter3d(
            x=push!(getindex.(poly, 1), poly[1][1]),
            y=push!(getindex.(poly, 2), poly[1][2]),
            z=push!(getindex.(poly, 3), poly[1][3]);
            mode="lines", hovertext="Cell", hoverinfo="text+x+y+z",
            line=attr(color=BZ_COL[], width=3)
            )
    end

    # lattice vectors
    tgs    = Vector{GenericTrace{Dict{Symbol,Any}}}(undef, 6)
    tgtips = Vector{GenericTrace{Dict{Symbol,Any}}}(undef, 3)
    intersects = axis_intersections(c, basis(c))
    for (i,V) in enumerate(basis(c))
        V′ = V./norm(V)
        name = "v$(i)"
        for j in (1,2) # in/outside BZ
            start = j == 1 ? 0.0           : intersects[i]
            stop  = j == 1 ? intersects[i] : 1.0
            V₀ = start*V .+ (j == 1 ? 0.0   : 0.025)*scale*V′
            V₁ = stop*V  .- (j == 1 ? 0.025 : 0.0)  *scale*V′
            tgs[i+(j-1)*3] = PlotlyJS.scatter3d(
                x=[V₀[1],V₁[1]], y=[V₀[2],V₁[2]], z=[V₀[3],V₁[3]];
                mode="lines", hovertext=name, hoverinfo="text",
                line=attr(color=ifelse(j==1, BASIS_LIGHT_COL[], BASIS_COL[]),
                          width=ifelse(j==1, 5, 6))
                )
        end
        tgtips[i] = PlotlyJS.cone(
            x=[V[1]], u=[V′[1]], y=[V[2]], v=[V′[2]], z=[V[3]], w=[V′[3]],
            sizeref=.1*scale, showscale=false, anchor="tail",
            colorscale=[[0, BASIS_COL[]], [1, BASIS_COL[]]],
            hovertext=name, hoverinfo="text+x+y+z")
    end

    # Cartesian axes
    cart_basis = cartesian_axes(Val(3))
    name_axs = ("x", "y", "z")
    taxs    = Vector{GenericTrace{Dict{Symbol,Any}}}(undef, 6)
    taxtips = Vector{GenericTrace{Dict{Symbol,Any}}}(undef, 3)
    intersects = axis_intersections(c, cart_basis)
    for (i, V) in enumerate(cart_basis)
        name = ""*string(name_axs[i])*""
        V′ = V./norm(V)
        for j in (1,2) # in/outside BZ
            start = j == 1 ? 0.0           : intersects[i]
            stop  = j == 1 ? intersects[i] : intersects[i]
            V₀ = start*V .+ (j == 1 ? 0.0   : 0.025)*scale*V′
            V₁ = stop*V  .- (j == 1 ? 0.025 : -0.2) *scale*V′
            taxs[i+(j-1)*3] = PlotlyJS.scatter3d(
                x=[V₀[1],V₁[1]], y=[V₀[2],V₁[2]], z=[V₀[3],V₁[3]];
                mode="lines", hovertext=name, hoverinfo="text",
                line=attr(color=ifelse(j==1, AXIS_LIGHT_COL[], AXIS_COL[]),
                          width=ifelse(j==1, 5, 6))
                )
            if j == 2
                taxtips[i] = PlotlyJS.cone(
                    x=[V₁[1]], u=[V[1]],
                    y=[V₁[2]], v=[V[2]],
                    z=[V₁[3]], w=[V[3]],
                    sizeref=.1*scale, showscale=false, anchor="tail",
                    colorscale=[[0, AXIS_COL[]], [1, AXIS_COL[]]],
                    hovertext=name, hoverinfo="text")
            end
        end
    end

    # combine traces and plot
    ts = vcat(tbz, tgs, tgtips, taxs, taxtips)
    return PlotlyJS.plot(ts, layout; config=config)
end

# ---------------------------------------------------------------------------------------- #
# 2D

# default layout
const DEFAULT_PLOTLY_LAYOUT_2D  = Layout(
    showlegend=false,
    xaxis=attr(tickvals=[], zeroline=false,
            showgrid=false, showbackground=false,
            title=attr(text=""),
            ),
    yaxis=attr(tickvals=[], zeroline=false,
            showgrid=false, showbackground=false,
            title=attr(text=""),
            scaleanchor="x", scaleratio=1
            ),
    aspectmode = "data",
    hovermode="closest",
    margin=attr(l=0, r=0, b=0, t=0),
    autosize=false,
    plot_bgcolor="rgba(255, 255, 255, 1)", paper_bgcolor="rgba(255, 255, 255, 1)",
    annotations=PlotlyBase.PlotlyAttribute[]
    )

function plot(c::Cell{2}, layout::Layout=DEFAULT_PLOTLY_LAYOUT_2D;
              config::PlotConfig = PlotConfig(responsive=true, displaylogo=false))
    layout = deepcopy(layout) # because we have to mutate to get arrows...

    setting(c) !== CARTESIAN && (c = cartesianize(c))

    scale = maximum(norm, basis(c))
    max_x, max_y = maximum(v->abs(v[1]), basis(c)), maximum(v->abs(v[2]), basis(c))
    get!(layout[:xaxis], :range, [-max_x-scale/15, max_x+scale/15])
    get!(layout[:yaxis], :range, [-max_y-scale/15, max_y+scale/15])

    # Cell boundaries
    tbz = Vector{GenericTrace{Dict{Symbol,Any}}}(undef, length(c))
    for (i,poly) in enumerate(c)
        tbz[i] = PlotlyJS.scatter(
            x=push!(getindex.(poly, 1), poly[1][1]),
            y=push!(getindex.(poly, 2), poly[1][2]);
            mode="lines", hovertext="Cell", hoverinfo="text+x+y",
            line=attr(color=BZ_COL[], width=3)
            )
    end

    # lattice vectors
    tgs    = Vector{GenericTrace{Dict{Symbol,Any}}}(undef, 4)
    tgtips = Vector{GenericTrace{Dict{Symbol,Any}}}(undef, 2)
    intersects = axis_intersections(c, basis(c))
    for (i,V) in enumerate(basis(c))
        V′ = V./norm(V)
        name = "v$(i)"
        for j in (1,2) # in/outside BZ
            start = j == 1 ? 0.0           : intersects[i]
            stop  = j == 1 ? intersects[i] : 1.0
            V₀ = start*V .+ (j == 1 ? 0.0   : 0.025)*scale*V′
            V₁ = stop*V  .- (j == 1 ? 0.025 : 0.0)  *scale*V′
            tgs[i+(j-1)*2] = PlotlyJS.scatter(
                x=[V₀[1],V₁[1]], y=[V₀[2],V₁[2]];
                mode="lines", hovertext=name, hoverinfo=ifelse(j==1,"text","text+x+y"),
                line=attr(color=ifelse(j==1, BASIS_LIGHT_COL[], BASIS_COL[]),
                          width=ifelse(j==1, 5, 6))
                )
        end
        # arrow heads have to be added as annotations to layout in 2D :/
        haskey(layout, :annotations) || (layout[:annotations] = PlotlyBase.PlotlyAttribute[])
        push!(layout[:annotations],
            attr(x=V[1]+.05V′[1]*scale, y=V[2]+.05V′[2]*scale,   # awful fidgeting; plotly's
                 ax=V[1]-.05V′[1]*scale, ay=V[2]-.05V′[2]*scale, # arrows are stupid
                 xref="ax", yref="ay", axref="x", ayref="y",
                 showarrow=true, arrowhead=2, arrowwidth=6, arrowsize=.5,
                 arrowcolor=BASIS_COL[]))
    end

    # Cartesian axes
    cart_basis = cartesian_axes(Val(2))
    name_axs = ("x", "y")
    taxs    = Vector{GenericTrace{Dict{Symbol,Any}}}(undef, 4)
    taxtips = Vector{GenericTrace{Dict{Symbol,Any}}}(undef, 2)
    intersects = axis_intersections(c, cart_basis)
    for (i, V) in enumerate(cart_basis)
        name = ""*string(name_axs[i])*""
        V′ = V./norm(V)
        for j in (1,2) # in/outside BZ
            start = j == 1 ? 0.0           : intersects[i]
            stop  = j == 1 ? intersects[i] : intersects[i]
            V₀ = start*V .+ (j == 1 ? 0.0   : 0.025)*scale*V′
            V₁ = stop*V  .- (j == 1 ? 0.025 : -0.2) *scale*V′

            taxs[i+2(j-1)] = PlotlyJS.scatter(
                x=[V₀[1],V₁[1]], y=[V₀[2],V₁[2]];
                mode="lines", hovertext=name, hoverinfo="text",
                line=attr(color=ifelse(j==1, AXIS_LIGHT_COL[], AXIS_COL[]),
                            width=ifelse(j==1, 5, 6))
                )
            if j == 2
                # arrow heads have to be added as annotations to layout in 2D :/
                haskey(layout, :annotations) || (layout[:annotations] = PlotlyBase.PlotlyAttribute[])
                push!(layout[:annotations],
                    attr(x=V₁[1]+.05V′[1]*scale, y=V₁[2]+.05V′[2]*scale,   # awful fidgeting; plotly's
                         ax=V₁[1]-.05V′[1]*scale, ay=V₁[2]-.05V′[2]*scale, # arrows are stupid
                         xref="ax", yref="ay", axref="x", ayref="y",
                         showarrow=true, arrowhead=2, arrowwidth=6, arrowsize=.5,
                         arrowcolor=AXIS_COL[]))
            end
        end
    end

    # combine traces and plot
    ts = vcat(tbz, tgs, taxs)
    return PlotlyJS.plot(ts, layout; config=config)
end

# ---------------------------------------------------------------------------------------- #
# UTILITIES

# the "outward" intersections of of lines with direction `Vs` though origo with a cell face
function axis_intersections(c::Cell{3},
            Vs::AbstractVector{<:AbstractVector}=cartesian_axes(Val(3)))

    intersects = MVector{3,Float64}(undef)
    fill!(intersects, Inf)
    # intersection between plane (r-rₒ) ⋅ n = 0 and line r(t) = l₀ + lt:
    #   t = (r₀ - l₀) ⋅ n / (l ⋅ n)     [https://wikipedia.org/wiki/Line–plane_intersection]
    for (i,poly) in enumerate(c)
        r₀ = sum(poly)/length(poly)
        n  = face_normal(c, i)
        for (j,V) in enumerate(Vs)
            t = dot(r₀, n)/dot(V, n)
            if t > 0.0 && t < intersects[j]
                intersects[j] = t
            end
        end
    end
    return intersects
end

# the "outward" intersections of of lines with direction `Vs` though origo with a cell segment
function axis_intersections(c::Cell{2},
    Vs::AbstractVector{<:AbstractVector}=cartesian_axes(Val(2)))

    intersects = MVector{2,Float64}(undef)
    fill!(intersects, Inf)
    # intersection between line rₐ(t) a + nₐt and line rᵦ(t) = β + nᵦt can be gotten by
    # linear algebra: [-nₐ|nᵦ][tₐ,tᵦ] = a-β. Here we, pick rᵦ for the axes `Vs` (β=0) and
    # let rₐ refer to cell segments. We're then only interested in tᵦ:
    fs = faces(c)
    vs = vertices(c)
    for i in eachindex(vs) # number of vertices and segments are equal in 2D
        if length(fs) == 1 # `merge = true`
            idxs = i ≠ length(vs) ? (i,i+1) : (i,1)
        else
            idxs = (fs[i][1], fs[i][2])
        end
        a  = vs[idxs[1]]
        nₐ = vs[idxs[2]] - a
        for (j,V) in enumerate(Vs)
            nᵦ = V
            tᵦ = (nₐ[2]*a[1]-nₐ[1]*a[2])/(-nₐ[1]*nᵦ[2] + nₐ[2]*nᵦ[1]) # ~ inv([-nₐ|nᵦ])
            if tᵦ > 0.0 && tᵦ < intersects[j]
                intersects[j] = tᵦ
            end
        end
    end
    return intersects
end

# a set of cartesian basis vectors in D dimensions
cartesian_axes(Dᵛ::Val{D}) where D = SVector(ntuple(i->SVector(ntuple(j->i==j ? 1.0 : 0.0, Dᵛ)), Dᵛ))


text:
@testset "test that synonyms work"  begin
    @test type_1_errors === false_positives
    @test type_2_errors === false_negatives
    @test misclassified === incorrectly_classified
end

# --------------------------------------------------------------------
# scalar

const fun_refs_scalar = [
    (true_positives,  (0,0,0,1)),
    (true_negatives,  (1,0,0,0)),
    (false_positives, (0,1,0,0)),
    (false_negatives, (0,0,1,0)),
    (condition_positive, (0,0,1,1)),
    (condition_negative, (1,1,0,0)),
    (predicted_positive, (0,1,0,1)),
    (predicted_negative, (1,0,1,0)),
    (correctly_classified,   (1,0,0,1)),
    (incorrectly_classified, (0,1,1,0)),
]

@testset "scalar fuzzy with/without encoding" begin
    # Test the simple non-array (i.e. scalar) methods that
    # expect scalar values (single observations).
    # The scalar methods require either an encoding or at least
    # one boolean parameter. The later results in FuzzyBinary.
    # Here we test all cases of scalar fuzzy matching.
    enc = LabelEnc.FuzzyBinary()
    for (fun, mask) in fun_refs_scalar
        @testset "$(fun): boolean parameters" begin
            @test @inferred(fun(false, false)) === mask[1]
            @test @inferred(fun(false, true )) === mask[2]
            @test @inferred(fun(true,  false)) === mask[3]
            @test @inferred(fun(true,  true )) === mask[4]
            @test @inferred(fun(false, false, enc)) === mask[1]
            @test @inferred(fun(false, true , enc)) === mask[2]
            @test @inferred(fun(true,  false, enc)) === mask[3]
            @test @inferred(fun(true,  true , enc)) === mask[4]
        end
        @testset "$(fun): numeric parameters" begin
            types = (Int32, Int64, Float32, Float64)
            for pos in (1, 2), neg in (0, -1, -2)
                for T1 in types, T2 in types
                    @testset "$T1 against $T2" begin
                        @test_throws ArgumentError fun(T1(neg), T2(pos))
                        @test @inferred(fun(T1(neg), T2(neg), enc)) === mask[1]
                        @test @inferred(fun(T1(neg), T2(pos), enc)) === mask[2]
                        @test @inferred(fun(T1(pos), T2(neg), enc)) === mask[3]
                        @test @inferred(fun(T1(pos), T2(pos), enc)) === mask[4]
                    end
                end
            end
        end
        @testset "$(fun): mixed parameters" begin
            for pos in (1, true, 0.5, 1.), neg in (0, -1, -2, false, 0.)
                @testset "pos = $pos, neg = $neg" begin
                    @test @inferred(fun(neg, neg, enc)) === mask[1]
                    @test @inferred(fun(neg, pos, enc)) === mask[2]
                    @test @inferred(fun(pos, neg, enc)) === mask[3]
                    @test @inferred(fun(pos, pos, enc)) === mask[4]
                    if typeof(pos) <: Bool || typeof(neg) <: Bool
                        @test @inferred(fun(neg, pos)) === mask[2]
                        @test @inferred(fun(pos, neg)) === mask[3]
                    end
                end
            end
        end
    end
end

@testset "scalar with encoding" begin
    for (fun, mask) in fun_refs_scalar
        @testset "$(fun): binary" begin
            enc = LabelEnc.ZeroOne()
            @test @inferred(fun(0,   0,   enc)) === mask[1]
            @test @inferred(fun(0,   0.2, enc)) === mask[1]
            @test @inferred(fun(0.0, 1.0, enc)) === mask[2]
            @test @inferred(fun(0.0, 0.6, enc)) === mask[2]
            @test @inferred(fun(1f0, 0,   enc)) === mask[3]
            @test @inferred(fun(0x1, 0x1, enc)) === mask[4]
            enc = LabelEnc.MarginBased()
            @test @inferred(fun(-1,    -1, enc)) === mask[1]
            @test @inferred(fun(-1,  -0.6, enc)) === mask[1]
            @test @inferred(fun(-1.0, 1.0, enc)) === mask[2]
            @test @inferred(fun(-1,   2.2, enc)) === mask[2]
            @test @inferred(fun(1f0,   -1, enc)) === mask[3]
            @test @inferred(fun(1,      5, enc)) === mask[4]
            enc = LabelEnc.TrueFalse()
            @test @inferred(fun(false, false, enc)) === mask[1]
            @test @inferred(fun(false, true , enc)) === mask[2]
            @test @inferred(fun(true,  false, enc)) === mask[3]
            @test @inferred(fun(true,  true , enc)) === mask[4]
            enc = LabelEnc.Indices(2)
            @test @inferred(fun(2, 2, enc)) === mask[1]
            @test @inferred(fun(2., 2, enc)) === mask[1]
            @test @inferred(fun(0x2, 2f0, enc)) === mask[1]
            @test @inferred(fun(2, 1, enc)) === mask[2]
            @test @inferred(fun(1, 2, enc)) === mask[3]
            @test @inferred(fun(1, 1, enc)) === mask[4]
            enc = LabelEnc.NativeLabels([:pos, :neg])
            @test_throws ArgumentError fun(:neg, :pos)
            @test @inferred(fun(:neg, :neg, enc)) === mask[1]
            @test @inferred(fun(:neg, :pos, enc)) === mask[2]
            @test @inferred(fun(:pos, :neg, enc)) === mask[3]
            @test @inferred(fun(:pos, :pos, enc)) === mask[4]
            @test fun(:neg, :neg, [:pos, :neg]) === mask[1]
            @test fun(:neg, :pos, [:pos, :neg]) === mask[2]
            @test fun(:pos, :neg, [:pos, :neg]) === mask[3]
            @test fun(:pos, :pos, [:pos, :neg]) === mask[4]
            enc = LabelEnc.OneVsRest("pos")
            @test_throws ArgumentError fun("neg", "pos")
            @test @inferred(fun("foo", "bar", enc)) === mask[1]
            @test @inferred(fun("baz", "pos", enc)) === mask[2]
            @test @inferred(fun("pos", "foo", enc)) === mask[3]
            @test @inferred(fun("pos", "pos", enc)) === mask[4]
            enc = LabelEnc.OneOfK(2)
            @test_throws ArgumentError fun([1,0], 1)
            @test_throws ArgumentError fun(1, [0,1])
            @test_throws ArgumentError fun(1, 1, enc)
            @test_throws ArgumentError fun(1, [0,1], enc)
            @test_throws ArgumentError fun([0,1], 1, enc)
            @test_throws ArgumentError fun([0,1], [0,1], enc)
            @test_throws ArgumentError fun([0 1], [0 1], enc)
        end
        @testset "$(fun): multiclass" begin
            enc = LabelEnc.Indices(3)
            @test @inferred(fun(1, 2, enc)) == Dict{Int,Int}(
                1 => mask[3],
                2 => mask[2],
                3 => mask[1],
            )
            @test @inferred(fun(2, 2, enc)) == Dict{Int,Int}(
                1 => mask[1],
                2 => mask[4],
                3 => mask[1],
            )
            @test @inferred(fun(3, 2, enc)) == Dict{Int,Int}(
                1 => mask[1],
                2 => mask[2],
                3 => mask[3],
            )
            enc = LabelEnc.NativeLabels([:a, :b, :c])
            @test @inferred(fun(:a, :b, enc)) == Dict{Symbol,Int}(
                :a => mask[3],
                :b => mask[2],
                :c => mask[1],
            )
            @test @inferred(fun(:b, :b, enc)) == Dict{Symbol,Int}(
                :a => mask[1],
                :b => mask[4],
                :c => mask[1],
            )
            @test @inferred(fun(:c, :b, enc)) == Dict{Symbol,Int}(
                :a => mask[1],
                :b => mask[2],
                :c => mask[3],
            )
            @test fun(:a, :b, [:a, :b, :c]) == Dict{Symbol,Int}(
                :a => mask[3],
                :b => mask[2],
                :c => mask[1],
            )
            @test fun(:b, :b, [:a, :b, :c]) == Dict{Symbol,Int}(
                :a => mask[1],
                :b => mask[4],
                :c => mask[1],
            )
            @test fun(:c, :b, [:a, :b, :c]) == Dict{Symbol,Int}(
                :a => mask[1],
                :b => mask[2],
                :c => mask[3],
            )
            enc = LabelEnc.OneOfK(3)
            @test_throws ArgumentError fun(1, 1, enc)
            @test_throws ArgumentError fun(1, [0,1], enc)
            @test_throws ArgumentError fun([0,1], 1, enc)
            @test_throws ArgumentError fun([0,1], [0,1], enc)
            @test_throws ArgumentError fun([0 1], [0 1], enc)
        end
    end
end

# --------------------------------------------------------------------
# arrays

y_true_p = [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0]
y_true_m = [1, 1, 1, 1, 1, 1, 1, 1, 1, 1,-1,-1,-1,-1,-1,-1,-1]
y_true_i = [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2]
y_true_b = Vector{Bool}(y_true_p)
y_true_pf = Vector{Float64}(y_true_p)
y_true_mf = Vector{Float64}(y_true_m)
y_true_sym = convertlabel([:pos,:neg], y_true_p, LabelEnc.ZeroOne())
y_true_sym2 = [:a, :a, :a, :a, :a, :a, :a, :a, :a, :a, :b, :b, :b, :b, :c, :c, :c]
y_true_c = [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3]

y_hat_p  = [1, 1, 0, 0, 1, 0, 1, 0, 0, 1, 0, 1, 0, 1, 0, 1, 0]
y_hat_m  = [1, 1,-1,-1, 1,-1, 1,-1,-1, 1,-1, 1,-1, 1,-1, 1,-1]
y_hat_i  = [1, 1, 2, 2, 1, 2, 1, 2, 2, 1, 2, 1, 2, 1, 2, 1, 2]
y_hat_b  = Vector{Bool}(y_hat_p)
y_hat_pf = Vector{Float64}(y_hat_p)
y_hat_mf = Vector{Float64}(y_hat_m)
y_hat_sym  = convertlabel([:pos,:neg], y_hat_p,  LabelEnc.ZeroOne())
y_hat_sym2  = [:a, :a, :b, :c, :a, :b, :a, :c, :c, :a, :c, :a, :c, :a, :b, :a, :c]
y_hat_c  = [1, 1, 2, 3, 1, 2, 1, 3, 3, 1, 3, 1, 3, 1, 2, 1, 3]

const fun_refs_arrays = [
    (true_positives,  5),
    (true_negatives,  4),
    (false_positives, 3),
    (false_negatives, 5),
    (condition_positive, 10),
    (condition_negative, 7),
    (predicted_positive, 8),
    (predicted_negative, 9),
    (correctly_classified, 9),
    (incorrectly_classified, 8),
]

@testset "arrays fuzzy with/without encoding" begin
    targets = [y_true_p, y_true_m, y_true_b, y_true_pf, y_true_mf]
    outputs = [y_hat_p, y_hat_m, y_hat_b, y_hat_pf, y_hat_mf, y_hat_b']
    for (fun, ref) in fun_refs_arrays
        @testset "$fun: binary" begin
            for target in targets, output in outputs
                @testset "$(typeof(target)) against $(typeof(output))" begin
                    @test @inferred(fun(target, output, LabelEnc.FuzzyBinary())) === ref
                    if eltype(target) <: Bool || eltype(output) <: Bool
                        @test @inferred(fun(target, output)) === ref
                    end
                end
            end
        end
    end
end

@testset "arrays with/without encoding" begin
    binary_data = [
        (y_true_p,    y_hat_p,    LabelEnc.ZeroOne()),
        (y_true_p,    y_hat_pf,   LabelEnc.ZeroOne()),
        (y_true_pf,   y_hat_p,    LabelEnc.ZeroOne()),
        (y_true_pf,   y_hat_pf,   LabelEnc.ZeroOne()),
        (y_true_m,    y_hat_m,    LabelEnc.MarginBased()),
        (y_true_m,    y_hat_mf,   LabelEnc.MarginBased()),
        (y_true_mf,   y_hat_m,    LabelEnc.MarginBased()),
        (y_true_mf,   y_hat_mf,   LabelEnc.MarginBased()),
        (y_true_i,    y_hat_i,    LabelEnc.Indices(2)),
        (y_true_b,    y_hat_b,    LabelEnc.TrueFalse()),
        (y_true_sym,  y_hat_sym,  LabelEnc.NativeLabels([:pos,:neg])),
        (y_true_sym,  y_hat_sym,  LabelEnc.OneVsRest(:pos)),
        (y_true_sym2, y_hat_sym2, LabelEnc.OneVsRest(:a)),
    ]
    multiclass_data = [
        (y_true_sym2,    y_hat_sym2,  LabelEnc.NativeLabels([:a,:b,:c])),
        (y_true_c,          y_hat_c,  LabelEnc.Indices(3)),
        (Int32.(y_true_c),  y_hat_c,  LabelEnc.Indices(3)),
        (y_true_c, Float32.(y_hat_c), LabelEnc.Indices(3)),
    ]
    for (fun, ref) in fun_refs_arrays
        @testset "$fun: binary" begin
            for (target, output, enc) in binary_data
                @testset "$(eltype(target)) against $(eltype(output)) (enc: $(typeof(enc)))" begin
                    @test @inferred(fun(target, output, enc)) === ref
                    if nlabel(target) == 2
                        @test fun(target, output) === ref
                    end
                end
            end
            # this is not type stable (number of labels not static)
            @test fun(y_true_sym, y_hat_sym, [:pos, :neg]) === ref
        end
        @testset "$fun: multiclass" begin
            for (target, output, enc) in multiclass_data
                @testset "$(eltype(target)) against $(eltype(output)) (enc: $(typeof(enc)))" begin
                    res = Dict(map(label(enc)) do l
                        l => fun(target, output, LabelEnc.OneVsRest(l))
                    end)
                    @test @inferred(fun(target, output, enc)) == res
                    @test fun(target, output) == res
                end
            end
        end
    end
end

ext: jl


text:
using DiffEqOperators, Test

@testset "Curl operation for a 3D Vector" begin

    s = x, y, z = (-5:1.25:5, -5:1.25:5, -5:1.25:5)
    dx = dy = dz = x[2] - x[1]

    # Vector u0 = (y^2 + z^2) ê₁ + (x^2 + z^2) ê₂  + (x^2 + y^2) ê₃

    u0 = zeros(Float64,length(x),length(y),length(z),3)

    for k in 1:length(z), j in 1:length(y), i in 1:length(x)
        u0[i,j,k,1] = y[j]^2 + z[k]^2
        u0[i,j,k,2] = x[i]^2 + z[k]^2
        u0[i,j,k,3] = x[i]^2 + y[j]^2
    end

    # Analytic Curl of the given vector given by u_analytic = 2(y-z) ê₁ + 2(z-x) ê₂  + 2(x-y) ê₃
    u_analytic = zeros(Float64,length(x)-2,length(y)-2,length(z)-2,3)
    for k in 1:length(z)-2, j in 1:length(y)-2, i in 1:length(x)-2
        u_analytic[i,j,k,1] = 2*(y[j+1] - z[k+1])
        u_analytic[i,j,k,2] = 2*(z[k+1] - x[i+1])
        u_analytic[i,j,k,3] = 2*(x[i+1] - y[j+1])
    end

    A = Curl(4,(dx,dy,dz),size(u0)[1:end-1].-2)

    u = A*u0

    for I in CartesianIndices(u)
        @test u[I] ≈ u_analytic[I] atol=1e-3
    end

    # check for multiplication with constant

    u1 = 2*A*u0

    for I in CartesianIndices(u)
        @test u1[I] ≈ 2*u[I] atol=1e-3
    end

    # check for non-uniform grid

    dx = dy = dz = 1.25*ones(length(s[1])-1)

    A = Curl(4,(dx,dy,dz),size(u0)[1:end-1].-2)

    u = A*u0

    for I in CartesianIndices(u)
        @test u[I] ≈ u_analytic[I] atol=1e-3
    end
end


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:3] <= 1)

                     @variable(m, -3.428579430434188 <= q <= 8.991868736935128)

                     add_NL_constraint(m, :(softplus(0.758903986598098 + 0.8106379052679875*softplus(0.47487042435595583 + 0.4488775132618512*softplus(0.1716391672055586 + 0.7301849089692007*softplus(0.33202390613731003 + -0.8931954910636541*$(x[1]) + 0.7142724585602056*$(x[2]) + -0.28852147983722487*$(x[3])) + 0.6276077641275091*softplus(-0.8935890604396688 + -0.29936605585485676*$(x[1]) + -0.5391102865466948*$(x[2]) + 0.8633503800224367*$(x[3]))) + 0.015398554931854491*softplus(0.8366427516313624 + 0.8654420505988725*softplus(0.33202390613731003 + -0.8931954910636541*$(x[1]) + 0.7142724585602056*$(x[2]) + -0.28852147983722487*$(x[3])) + -0.07448940361379819*softplus(-0.8935890604396688 + -0.29936605585485676*$(x[1]) + -0.5391102865466948*$(x[2]) + 0.8633503800224367*$(x[3])))) + -0.6668215058198093*softplus(0.8216722410540171 + 0.5667551555087149*softplus(0.1716391672055586 + 0.7301849089692007*softplus(0.33202390613731003 + -0.8931954910636541*$(x[1]) + 0.7142724585602056*$(x[2]) + -0.28852147983722487*$(x[3])) + 0.6276077641275091*softplus(-0.8935890604396688 + -0.29936605585485676*$(x[1]) + -0.5391102865466948*$(x[2]) + 0.8633503800224367*$(x[3]))) + 0.5912184082121308*softplus(0.8366427516313624 + 0.8654420505988725*softplus(0.33202390613731003 + -0.8931954910636541*$(x[1]) + 0.7142724585602056*$(x[2]) + -0.28852147983722487*$(x[3])) + -0.07448940361379819*softplus(-0.8935890604396688 + -0.29936605585485676*$(x[1]) + -0.5391102865466948*$(x[2]) + 0.8633503800224367*$(x[3]))))) + softplus(-0.09143860592972564 + 0.3407854070498746*softplus(0.47487042435595583 + 0.4488775132618512*softplus(0.1716391672055586 + 0.7301849089692007*softplus(0.33202390613731003 + -0.8931954910636541*$(x[1]) + 0.7142724585602056*$(x[2]) + -0.28852147983722487*$(x[3])) + 0.6276077641275091*softplus(-0.8935890604396688 + -0.29936605585485676*$(x[1]) + -0.5391102865466948*$(x[2]) + 0.8633503800224367*$(x[3]))) + 0.015398554931854491*softplus(0.8366427516313624 + 0.8654420505988725*softplus(0.33202390613731003 + -0.8931954910636541*$(x[1]) + 0.7142724585602056*$(x[2]) + -0.28852147983722487*$(x[3])) + -0.07448940361379819*softplus(-0.8935890604396688 + -0.29936605585485676*$(x[1]) + -0.5391102865466948*$(x[2]) + 0.8633503800224367*$(x[3])))) + 0.24612879529032217*softplus(0.8216722410540171 + 0.5667551555087149*softplus(0.1716391672055586 + 0.7301849089692007*softplus(0.33202390613731003 + -0.8931954910636541*$(x[1]) + 0.7142724585602056*$(x[2]) + -0.28852147983722487*$(x[3])) + 0.6276077641275091*softplus(-0.8935890604396688 + -0.29936605585485676*$(x[1]) + -0.5391102865466948*$(x[2]) + 0.8633503800224367*$(x[3]))) + 0.5912184082121308*softplus(0.8366427516313624 + 0.8654420505988725*softplus(0.33202390613731003 + -0.8931954910636541*$(x[1]) + 0.7142724585602056*$(x[2]) + -0.28852147983722487*$(x[3])) + -0.07448940361379819*softplus(-0.8935890604396688 + -0.29936605585485676*$(x[1]) + -0.5391102865466948*$(x[2]) + 0.8633503800224367*$(x[3]))))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
using ConvexTests, CSDP

using MathOptInterface
const MOI = MathOptInterface
using Dualization

@info "Starting CSDP tests"

exclusions =  [r"mip", r"exp", r"sdp_Complex_Semidefinite_constraint"]

do_tests("CSDP", () -> CSDP.Optimizer(printlevel=0); exclude = exclusions, first = true, last = false)

function dual_opt()
    DualOptimizer(CSDP.Optimizer(printlevel=0))
end

do_tests("CSDP", dual_opt; exclude = exclusions, first = false, last = true, variant="(dualized)")


text:
# automatically generated by the FlatBuffers compiler, do not modify

FlatBuffers.@with_kw mutable struct Geometry
    ends::Vector{UInt32} = []
    xy::Vector{Float64} = []
    z::Vector{Float64} = []
    m::Vector{Float64} = []
    t::Vector{Float64} = []
    tm::Vector{UInt64} = []
    type::GeometryType = 0
    parts::Vector{Geometry} = []
end
FlatBuffers.@ALIGN(Geometry, 1)
FlatBuffers.slot_offsets(::Type{T}) where {T <: Geometry} = [
    0x00000004, 0x00000006, 0x00000008, 0x0000000A,
    0x0000000C, 0x0000000E, 0x00000010, 0x00000012
]

Geometry(buf::AbstractVector{UInt8}) = FlatBuffers.read(Geometry, buf)
Geometry(io::IO) = FlatBuffers.deserialize(io, Geometry)

FlatBuffers.@with_kw mutable struct Feature
    geometry::Geometry = Geometry()
    properties::Vector{UInt8} = []
    columns::Vector{Column} = []
end
FlatBuffers.@ALIGN(Feature, 1)
FlatBuffers.slot_offsets(::Type{T}) where {T <: Feature} = [
    0x00000004, 0x00000006, 0x00000008
]
FlatBuffers.root_type(::Type{T}) where {T <: Feature} = true

Feature(buf::AbstractVector{UInt8}) = FlatBuffers.read(Feature, buf)
Feature(io::IO) = FlatBuffers.deserialize(io, Feature)


text:
global const cnst_c = 2.9979e10
global const slope_soft = 1.e-6

function readSingleCRShockDataFromOutputFile(file::String)

        # read file into memory
        f = open(file)
        lines = readlines(f)
        close(f)

        # filter only relevant lines. Output of every line starts with "CR DATA".
        g = occursin.("CR DATA", lines)
        lines = lines[g]

        # get number of output lines
        N = length(lines)

        # init datatype for analyze
        cr = CRShockData(N)

        for i ∈ 1:N
                line_content            = split(lines[i])
                cr.dt[i]                = parse(Float64,line_content[3])
                cr.Mach[i]              = parse(Float64,line_content[4])
                cr.Shock_Speed[i]       = parse(Float64,line_content[5])
                cr.Shock_Compress[i]    = parse(Float64,line_content[6])
                cr.Shock_Energy_In[i]   = parse(Float64,line_content[7])
                cr.Shock_Energy_Real[i] = parse(Float64,line_content[8])
                cr.Energy_P[i]          = parse(Float64,line_content[9])
                cr.Energy_e[i]          = parse(Float64,line_content[10])
        end

        return cr
end

function getCRMomentumDistributionFromPartID(snap_file::String, ID::Integer;
                                             pmin::Real=1.0, pmax::Real=1.0e6,
                                             Nbins::Integer=0, new::Bool=true)

    ymin = 1.e-20
    h = head_to_obj(snap_file)
    info = read_info(snap_file)

    if info == 1
        if Nbins == 0
            error("Can't read spectrum! No info block present!\nSupply number of momentum bins to proceed!")
        else
            info = Array{Info_Line,1}(undef,7)
            info[1] = Info_Line("ID", UInt32, Int32(1), [1, 0, 0, 0, 0, 0])
            info[2] = Info_Line("CRpN", Float32, Int32(24), [1, 0, 0, 0, 0, 0])
            info[3] = Info_Line("CRpS", Float32, Int32(24), [1, 0, 0, 0, 0, 0])
            info[4] = Info_Line("CRpC", Float32, Int32(1), [1, 0, 0, 0, 0, 0])
            info[5] = Info_Line("CReN", Float32, Int32(24), [1, 0, 0, 0, 0, 0])
            info[6] = Info_Line("CReS", Float32, Int32(24), [1, 0, 0, 0, 0, 0])
            info[7] = Info_Line("CReC", Float32, Int32(1), [1, 0, 0, 0, 0, 0])
        end
    end

    id = read_block_by_name(snap_file, "ID",
                            info=info[getfield.(info, :block_name) .== "ID"][1],
                            parttype=0)

    part = findfirst( id .== UInt32(ID) )[1]

    # protons
    CRpN = read_block_by_name(snap_file, "CRpN",
                            info=info[getfield.(info, :block_name) .== "CRpN"][1],
                            parttype=0)[part,:]
    CRpS = read_block_by_name(snap_file, "CRpS",
                            info=info[getfield.(info, :block_name) .== "CRpS"][1],
                            parttype=0)[part,:]
    CRpC = read_block_by_name(snap_file, "CRpC",
                            info=info[getfield.(info, :block_name) .== "CRpC"][1],
                            parttype=0)[part]

    # electrons
    CReN = read_block_by_name(snap_file, "CReN",
                            info=info[getfield.(info, :block_name) .== "CReN"][1],
                            parttype=0)[part,:]
    CReS = read_block_by_name(snap_file, "CReS",
                            info=info[getfield.(info, :block_name) .== "CReS"][1],
                            parttype=0)[part,:]
    CReC = read_block_by_name(snap_file, "CReC",
                            info=info[getfield.(info, :block_name) .== "CReC"][1],
                            parttype=0)[part]

    Nbins = length(CRpS)

    par = CRMomentumDistributionConfig(pmin, pmax, Nbins, new)

    cr = CRMomentumDistribution(par.Nbins)

    # compensate for io - needs to be converted to Float64!
    if new
        CRpN = 1.e-20 .* Float64.(CRpN)
        CReN = 1.e-20 .* Float64.(CReN)
    else
        CRpN = 1.e20 .* Float64.(CRpN)
        CReN = 1.e20 .* Float64.(CReN)
    end

    # get zeroth bin
    cr.CRp_bound[1] = pmin
    cr.CRp_dis[1] = CRpN[1]

    cr.CRe_bound[1] = pmin
    cr.CRe_dis[1] = CReN[1]

    # all other bins
    j = 2
    for i = 1:Nbins-1

        # upper boundary of bin
        cr.CRp_bound[j] = pmin * 10.0^(par.bin_width*i)
        cr.CRp_dis[j] = CRpN[i] * ( cr.CRp_bound[j]/cr.CRp_bound[j-1])^(-CRpS[i])
        if cr.CRp_bound[j] > CRpC/par.mc_p
            cr.CRp_bound[j] = CRpC/par.mc_p
        end

        cr.CRe_bound[j] = pmin * 10.0^(par.bin_width*i)
        cr.CRe_dis[j] = CReN[i] * ( cr.CRe_bound[j]/cr.CRe_bound[j-1])^(-CReS[i])
        if cr.CRe_bound[j] > CReC/par.mc_e
            cr.CRe_bound[j] = CReC/par.mc_e
        end

        # lower bound of next bin
        cr.CRp_bound[j+1] = cr.CRp_bound[j]
        cr.CRp_dis[j+1] = CRpN[i+1]
        if cr.CRp_bound[j] == CRpC/par.mc_p
            cr.CRp_bound[j+1] = CRpC/par.mc_p
        end

        cr.CRe_bound[j+1] = cr.CRe_bound[j]
        cr.CRe_dis[j+1] = CReN[i+1]
        if cr.CRe_bound[j] == CReC/par.mc_e
            cr.CRe_bound[j+1] = CReC/par.mc_e
        end

        j += 2

    end

    # last boundary
    cr.CRp_bound[j] = pmax
    if cr.CRp_bound[j-1] < CRpC/par.mc_p
        cr.CRp_bound[j] = CRpC/par.mc_p
    end
    cr.CRp_dis[j] = CRpN[Nbins] * ( cr.CRp_bound[j]/cr.CRp_bound[j-1])^(-CRpS[Nbins])
    cr.CRp_bound[j+1] = cr.CRp_bound[j]
    #cr.CRp_dis[j] = ymin

    cr.CRe_bound[j] = pmax
    if cr.CRe_bound[j-1] < CReC/par.mc_e
        cr.CRe_bound[j] = CReC/par.mc_e
    end
    cr.CRe_dis[j] = CReN[Nbins] * ( cr.CRe_bound[j]/cr.CRe_bound[j-1])^(-CReS[Nbins])
    cr.CRe_bound[j+1] = cr.CRe_bound[j]
    #cr.CRe_dis[j] = ymin * 1.e-2

    return cr
end

@inline function energy_integral(bound_low::Real, bound_up::Real,
                         norm::Real, slope::Real, ρ::Real)
    # energy integral (eq.21 M01)

    en = 4π * cnst_c * norm * bound_low^4 /  ρ
    energy = en * ( (bound_up/bound_low)^(4.0 - slope ) - 1.0 ) / ( 4.0 - slope )

    if ( 4.0 - slope_soft ) < slope < ( 4.0 + slope_soft )
        slope_var = (slope - 4.0)/slope_soft
        energy2 = en * log10(bound_up/bound_low)
        if slope_var != 0.0
            energy = energy * slope_var + energy2 * ( 1.0 - slope_var )
        else
            energy = energy2
        end
    end

    return energy
end

function calculateCREnergyInCGS(CR_N, CR_S, CR_Cut, ρ;
                                pmin::Real=10.0, pmax::Real=1.e7,
                                mc::Real=0.0, SelectBin::Int=-1,
                                units::GadgetPhysical=GadgetPhysical(),
                                verbose=true)
    # calculates engergy per bin in cgs units.

    if mc == 0.0
        error("Error! No mc specified! Please provide m_particle * c_light to continue!")
        return 1
    end

    CR_E = 0.0
    cnst_c = 2.9979e10
    ρ *= units.m_unit / units.l_unit^3  # convert rho in cgs units

    CR_N .*= 1.e20

    Nbins = length(CR_N)
    bin_width = log10(pmax/pmin)/Nbins
    bin = 10.0.^collect(log10(pmin):bin_width:log10(pmax)) .* mc
    above_cut = findall( bin .> CR_Cut/mc )
    if length(above_cut) > 0
        bin[above_cut] .= CR_Cut/mc
    end

    if SelectBin == -1
        CR_E = zeros(length(CR_N))
        for i = 1:length(CR_N)
            CR_E[i] = energy_integral(bin[i], bin[i+1], CR_N[i], CR_S[i], ρ)
        end
    elseif length(SelectBin) == 1
        CR_E = energy_integral(bin[SelectBin], bin[SelectBin+1], CR_N[SelectBin], CR_S[SelectBin], ρ)
    elseif length(SelectBin) > 1
        CR_E = zeros(length(SelectBin))
        if i ∈ SelectBin
            CR_E[i-SelectBin[1]+1] = energy_integral(bin[i], bin[i+1], CR_N[i], CR_S[i], ρ)
        end
    end
    return CR_E
end

@inline function density_integral(bound_low::Real, bound_up::Real,
                         norm::Real, slope::Real, ρ::Real)
    # density integral (eq. 9 M01)

    nb = 4π * norm * bound_low^3 / ρ
    density = nb * ( (bound_up/bound_low)^(3.0 - slope ) - 1.0 ) / ( 3.0 - slope )

    if ( 3.0 - slope_soft ) < slope < ( 3.0 + slope_soft )
        slope_var = (slope - 3.0)/slope_soft
        density2 = nb * log10(bound_up/bound_low)
        if slope_var != 0.0
            density = density * slope_var + density2 * ( 1.0 - slope_var )
        else
            density = density2
        end
    end

    return density
end

function calculateCRNumber(CR_N, CR_S, CR_Cut, ρ;pmin::Real=10.0, pmax::Real=1.e7,
                           mc::Real=0.0, SelectBin::Int=-1,
                           units::GadgetPhysical=GadgetPhysical(), verbose=true)
    # calculates engergy per bin in cgs units.

    CR_Num = 0.0
    cnst_c = 2.9979e10
    ρ *= units.m_unit / units.l_unit^3  # convert rho in cgs units

    CR_N .*= 1.e20  # compensate for io

    Nbins = length(CR_N)
    bin_width = log10(pmax/pmin)/Nbins
    bins = 10.0.^collect(log10(pmin):bin_width:log10(pmax)) .* mc
    above_cut = findall( bins .> CR_Cut/mc )

    if length(above_cut) > 0
        bins[above_cut] .= CR_Cut/mc
    end

    if SelectBin == -1
        CR_Num = zeros(length(CR_N))
        for i = 1:length(CR_N)
            CR_Num[i] = energy_integral(bin[i], bin[i+1], CR_N[i], CR_S[i], ρ)
        end
    elseif length(SelectBin) == 1
        CR_Num = energy_integral(bin[SelectBin], bin[SelectBin+1], CR_N[SelectBin], CR_S[SelectBin], ρ)
    elseif length(SelectBin) > 1
        CR_Num = zeros(length(SelectBin))
        if i ∈ SelectBin
            CR_Num[i-SelectBin[1]+1] = energy_integral(bin[i], bin[i+1], CR_N[i], CR_S[i], ρ)
        end
    end
    return CR_Num
end

max_stars_repo_name: LudwigBoess/GadJet

max_issues_repo_name: LudwigBoess/GadJet

max_forks_repo_name: LudwigBoess/GadJet


text:
#
# This file is part of the Simulate.jl Julia package, MIT license
#
# Paul Bayer, 2019
#
# This is a Julia package for discrete event simulation
#
# this implements the event handling
#

"""
    nextevent(c::Clock)

Return the next scheduled event.
"""
nextevent(c::Clock) = peek(c.sc.events)[1]

"""
    nextevtime(c::Clock)

Return the internal time (unitless) of next scheduled event.
"""
nextevtime(c::Clock) = peek(c.sc.events)[2]

"""
    evaluate(y,  m::Module)

Function barrier for arguments and keywords of a Fun. This allows Expr, Symbol
and Fun as arguments and keyword values.
"""
evaluate(y,  m::Module) = y    # catchall, gives back the argument
evaluate(y::Fun, m::Module) = sfExec(y, m)   # Fun: recursively call sfexec
function evaluate(y::Union{Symbol,Expr}, m::Module)  # Symbol,Expr: eval
    try
        ret = Core.eval(m, y)
        @warn "Evaluating expressions is slow, use `Fun` instead" maxlog=1
        return ret
    catch
        return y
    end
end

"""
    sfExec(x::Fun, m::Module)

Execute a Fun x.

- if x.f is an event!, its args and kws must be maintained for later evaluation,
    otherwise evaluate them now before passing them to x.f,
- branch to different invocation methods depending on args or kws,
- if x.f is not from Main and we are not on thread 1, call it with invokelatest
    to avoid a world age situation (it may be too new).
"""
function sfExec(x::Fun, m::Module)
    f_std = (parentmodule(x.f) != Main) || (threadid() == 1)
    if x.arg === nothing
        if x.kw === nothing
            f_std ? x.f() : Base.invokelatest(x.f)  # 1. no args and kws
        else
            kw = x.f === event! ? x.kw : (; zip(keys(x.kw), map(i->evaluate(i, m), values(x.kw)) )...)
            f_std ? x.f(; kw...) : Base.invokelatest(x.f; kw...)    # 2. only kws
        end
    else
        if x.kw === nothing
            arg = x.f === event! ? x.arg : map(i->evaluate(i, m), x.arg)
            f_std ? x.f(arg...) : Base.invokelatest(x.f, arg...)    # 3. only args
        else
            if x.f === event!
                arg = x.arg; kw = x.kw
            else
                arg = map(i->evaluate(i, m), x.arg)
                kw = (; zip(keys(x.kw), map(i->evaluate(i, m), values(x.kw)) )...)
            end
            f_std ? x.f(arg...; kw...) : Base.invokelatest(x.f, arg...; kw...) # 4. args and kws
        end
    end
end

# "Forward an expression to `evaluate`."
# sfExec(x::Expr, m::Module) = evaluate(x, m)

"""
    evExec(ex, m::Module=Main)

Function barrier for different ex: forward an event's `Fun`s or expressions
to further execution or evaluation.

# Return
the evaluated value or a tuple of evaluated values.
"""
evExec(ex::Fun, m::Module=Main) = sfExec(ex, m)
evExec(ex::Expr, m::Module=Main) = evaluate(ex, m)
evExec(ex::Tuple, m::Module=Main) = map(x->evExec(x, m), ex)

"""
    do_event!(c::Clock)

Execute or evaluate the next timed event on a clock c.
"""
function do_event!(c::Clock)
    c.time = c.tev
    ev = dequeue!(c.sc.events)
    evExec(ev.ex, ev.scope)
    c.evcount += 1
    if ev.Δt > 0.0  # schedule repeat event
        event!(c, ev.ex, c.time + ev.Δt, scope=ev.scope, cycle=ev.Δt)
    end
    c.tev = length(c.sc.events) ≥ 1 ? nextevtime(c) : c.time
end

"""
    do_tick!(s::Schedule)

First execute all sampling expressions in a schedule, then evaluate all
conditional events and if conditions are met, execute them.
"""
function do_tick!(c::Clock)
    c.time = c.tn
    foreach(x -> evExec(x.ex, x.scope), c.sc.samples)  # exec sampling
    # then lookup conditional events
    ix = findfirst(x->all(evExec(x.cond, x.scope)), c.sc.cevents)
    while ix !== nothing
        evExec(splice!(c.sc.cevents, ix).ex)
        if isempty(c.sc.cevents)
            if isempty(c.sc.samples) && isempty(c.ac) # no sampling and active clocks
                c.Δt = 0.0 # delete sample rate
            end
            break
        end
        ix = findfirst(x->all(evExec(x.cond, x.scope)), c.sc.cevents)
    end
    c.scount +=1
end

"""
    do_step!(c::Clock)

step forward to next tick or scheduled event.
"""
function do_step!(c::Clock)
    c.state = Busy()
    if (c.tev ≤ c.time) && (length(c.sc.events) > 0)
        c.tev = nextevtime(c)
    end

    if length(c.sc.events) > 0
        if c.Δt > 0.0
            if c.tn <= c.tev     # if t_next_tick  ≤ t_next_event
                do_tick!(c)
                if c.tn == c.tev # an event is scheduled at the same time
                    do_event!(c)
                end
                c.tn += c.Δt
            else
                do_event!(c)
            end
        else
            do_event!(c)
            c.tn = c.time
        end
    elseif c.Δt > 0.0
        do_tick!(c)
        c.tn += c.Δt
        c.tev = c.time
    else
        error("do_step!: nothing to evaluate")
    end
    length(c.processes) == 0 || yield() # let processes run
    (c.state == Busy()) && (c.state = Idle())
end


text:
function flush_suds!(S::GphysData,
  xc::Array{UnitRange{Int64},1},
  xn::Array{Int64,1},
  xt::Array{Int64,1},
  xz::Array{Int64,1},
  xj::Int64,
  v::Integer)

  (xj == 0) && return nothing

  # Post-read ================================================================
  resize!(xc, xj)
  resize!(xn, xj)
  resize!(xt, xj)
  resize!(xz, xj)

  if v > 1
    println("Done processing: ")
    println("Processed xj = ", xj, " segments")
    println("Segment ranges = ", xc)
    println("Segment starts = ", xn)
    println("Segment start times = ", xt)
    println("Channel segment lengths = ", xz)
  end

  # Determine length of each array in S.x
  Lx = zeros(Int64, S.n)
  for j = 1:xj
    nz = xz[j]
    for i in xc[j]
      Lx[i] += nz
    end
  end

  # Data assignment ==========================================================
  # Initialize arrays in S.x and counters sxi
  sxi = ones(Int64, S.n)
  for i = 1:S.n
    if isempty(S.x[i])
      S.x[i] = Array{Float32,1}(undef, Lx[i])
    else
      lxi = length(S.x[i])
      sxi[i] = lxi + 1
      resize!(S.x[i], lxi + Lx[i])
    end
  end

  # Loop again over xn; this time, copy each segment over
  for j = 1:xj
    xs = xn[j]
    si = xs
    for i in xc[j]
      t0 = xt[j]
      nz = xz[j]
      copyto!(S.x[i], sxi[i], BUF.x, si, nz)
      si += nz
      sxi[i] += nz

      # determine start time of channel
      if isempty(S.t[i])
        S.t[i] = Int64[1 t0; nz 0]
      else
        check_for_gap!(S, i, t0, nz, v)
      end
    end
  end
  return nothing
end


text:
abstract type AbstractBasis end

Base.size(basis::AbstractBasis) = size(basis.u)
getbeta(basis::AbstractBasis) = basis.β
getstatistics(basis::AbstractBasis) = basis.statistics

"""
    DimensionlessBasis <: AbstractBasis

Intermediate representation (IR) basis in reduced variables.

For a continuation kernel `K` from real frequencies, `ω ∈ [-ωmax, ωmax]`, to
imaginary time, `τ ∈ [0, β]`, this type stores the truncated singular
value expansion or IR basis:

    K(x, y) ≈ sum(u[l](x) * s[l] * v[l](y) for l in range(L))

The functions are given in reduced variables, `x = 2τ/β - 1` and
`y = ω/ωmax`, which scales both sides to the interval `[-1, 1]`.  The
kernel then only depends on a cutoff parameter `Λ = β * ωmax`.

# Examples
The following example code assumes the spectral function is a single
pole at `x = 0.2`. We first compute an IR basis suitable for fermions and `β*W ≤ 42`. Then we get G(iw) on the first few Matsubara frequencies:

```julia-repl
julia> using SparseIR

julia> basis = DimensionlessBasis(fermion, 42);

julia> gl = basis.s .* basis.v(0.2);

julia> giw = transpose(basis.uhat([1, 3, 5, 7])) * gl
```

# Fields
- `u::PiecewiseLegendrePolyVector`: Set of IR basis functions on the reduced imaginary time (`x`) axis. These functions are stored as piecewise Legendre polynomials.

  To obtain the value of all basis functions at a point or a array of
  points `x`, you can call the function `u(x)`.  To obtain a single
  basis function, a slice or a subset `l`, you can use `u[l]`.

- `uhat::PiecewiseLegendreFTVector`: Set of IR basis functions on the Matsubara frequency (`wn`) axis.
These objects are stored as a set of Bessel functions.

  To obtain the value of all basis functions at a Matsubara frequency
  or a array of points `wn`, you can call the function `uhat(wn)`.
  Note that we expect reduced frequencies, which are simply even/odd
  numbers for bosonic/fermionic objects. To obtain a single basis
  function, a slice or a subset `l`, you can use `uhat[l]`.

- `s`: Vector of singular values of the continuation kernel

- `v::PiecewiseLegendrePolyVector`: Set of IR basis functions on the reduced real frequency (`y`) axis.
These functions are stored as piecewise Legendre polynomials.

  To obtain the value of all basis functions at a point or a array of
  points `y`, you can call the function `v(y)`.  To obtain a single
  basis function, a slice or a subset `l`, you can use `v[l]`.

See also [`FiniteTempBasis`](@ref) for a basis directly in time/frequency.
"""
struct DimensionlessBasis{K<:AbstractKernel,T<:AbstractFloat} <: AbstractBasis
    kernel::K
    u::PiecewiseLegendrePolyVector{T}
    uhat::PiecewiseLegendreFTVector{T}
    s::Vector{T}
    v::PiecewiseLegendrePolyVector{T}
    statistics::Statistics
end

function Base.show(io::IO, a::DimensionlessBasis)
    return print(io, "DimensionlessBasis: statistics=$(getstatistics(a)), size=$(size(a))")
end

"""
    DimensionlessBasis(statistics, Λ, ε=nothing; kernel=LogisticKernel(Λ), sve_result=compute_sve(kernel; ε))

Construct an IR basis suitable for the given `statistics` and cutoff `Λ`.
"""
function DimensionlessBasis(
    statistics::Statistics, Λ, ε=nothing;
    kernel=LogisticKernel(Λ), sve_result=compute_sve(kernel; ε),
)
    u, s, v = sve_result
    size(u) == size(s) == size(v) || throw(DimensionMismatch("Mismatched shapes in SVE"))

    # The radius of convergence of the asymptotic expansion is Λ/2,
    # so for significantly larger frequencies we use the asymptotics,
    # since it has lower relative error.
    even_odd = Dict(fermion => :odd, boson => :even)[statistics]
    û = hat.(u, even_odd; n_asymp=conv_radius(kernel))
    return DimensionlessBasis(kernel, u, û, s, v, statistics)
end

"""
    Λ(basis)

Basis cutoff parameter `Λ = β * ωmax`.
"""
Λ(basis::DimensionlessBasis) = basis.kernel.Λ

function Base.getindex(basis::DimensionlessBasis, i)
    sve_result = basis.u[i], basis.s[i], basis.v[i]
    return DimensionlessBasis(basis.statistics, Λ(basis); kernel=basis.kernel, sve_result)
end

"""
    FiniteTempBasis <: AbstractBasis

Intermediate representation (IR) basis for given temperature.

For a continuation kernel `K` from real frequencies, `ω ∈ [-ωmax, ωmax]`, to
imaginary time, `τ ∈ [0, beta]`, this type stores the truncated singular
value expansion or IR basis:

    K(τ, ω) ≈ sum(u[l](τ) * s[l] * v[l](ω) for l in 1:L)

This basis is inferred from a reduced form by appropriate scaling of
the variables.

# Examples
The following example code assumes the spectral function is a single
pole at `ω = 2.5`. We first compute an IR basis suitable for fermions and `β = 10`, `W ≤ 4.2`. Then we get G(iw) on the first few Matsubara frequencies:

```julia-repl
julia> using SparseIR

julia> basis = FiniteTempBasis(fermion, 42, 4.2);

julia> gl = basis.s .* basis.v(2.5);

julia> giw = transpose(basis.uhat([1, 3, 5, 7])) * gl
```

# Fields
- `u::PiecewiseLegendrePolyVector`:
  Set of IR basis functions on the imaginary time (`tau`) axis.
  These functions are stored as piecewise Legendre polynomials.

  To obtain the value of all basis functions at a point or a array of
  points `x`, you can call the function `u(x)`.  To obtain a single
  basis function, a slice or a subset `l`, you can use `u[l]`.

- `uhat::PiecewiseLegendreFT`:
  Set of IR basis functions on the Matsubara frequency (`wn`) axis.
  These objects are stored as a set of Bessel functions.

  To obtain the value of all basis functions at a Matsubara frequency
  or a array of points `wn`, you can call the function `uhat(wn)`.
  Note that we expect reduced frequencies, which are simply even/odd
  numbers for bosonic/fermionic objects. To obtain a single basis
  function, a slice or a subset `l`, you can use `uhat[l]`.

- `s`: Vector of singular values of the continuation kernel

- `v::PiecewiseLegendrePoly`:
  Set of IR basis functions on the real frequency (`w`) axis.
  These functions are stored as piecewise Legendre polynomials.

  To obtain the value of all basis functions at a point or a array of
  points `w`, you can call the function `v(w)`.  To obtain a single
  basis function, a slice or a subset `l`, you can use `v[l]`.
"""
struct FiniteTempBasis{K,T} <: AbstractBasis
    kernel::K
    sve_result::Tuple{
        PiecewiseLegendrePolyVector{T},Vector{T},PiecewiseLegendrePolyVector{T}
    }
    statistics::Statistics
    β::T
    u::PiecewiseLegendrePolyVector{T}
    v::PiecewiseLegendrePolyVector{T}
    s::Vector{T}
    uhat::PiecewiseLegendreFTVector{T}
end

const _DEFAULT_FINITE_TEMP_BASIS = FiniteTempBasis{LogisticKernel,Float64}

function Base.show(io::IO, a::FiniteTempBasis{K,T}) where {K,T}
    return print(
        io, "FiniteTempBasis{$K, $T}($(getstatistics(a)), $(getbeta(a)), $(getwmax(a)))"
    )
end

"""
    FiniteTempBasis(statistics, β, wmax, ε=nothing; kernel=LogisticKernel(β * wmax), sve_result=compute_sve(kernel; ε))

Construct a finite temperature basis suitable for the given `statistics` and cutoffs `β` and `wmax`.
"""
function FiniteTempBasis(
    statistics::Statistics, β, wmax, ε=nothing;
    kernel=LogisticKernel(β * wmax), sve_result=compute_sve(kernel; ε),
)
    β > 0 || throw(DomainError(β, "Inverse temperature β must be positive"))
    wmax ≥ 0 || throw(DomainError(wmax, "Frequency cutoff wmax must be non-negative"))

    u, s, v = sve_result
    size(u) == size(s) == size(v) || throw(DimensionMismatch("Mismatched shapes in SVE"))

    # The polynomials are scaled to the new variables by transforming the
    # knots according to: tau = beta/2 * (x + 1), w = wmax * y.  Scaling
    # the data is not necessary as the normalization is inferred.
    wmax = kernel.Λ / β
    u_knots = β / 2 * (u.knots .+ 1)
    v_knots = wmax * v.knots
    u_ = PiecewiseLegendrePolyVector(u, u_knots; Δx=β / 2 * u.Δx, symm=u.symm)
    v_ = PiecewiseLegendrePolyVector(v, v_knots; Δx=wmax * v.Δx, symm=v.symm)

    # The singular values are scaled to match the change of variables, with
    # the additional complexity that the kernel may have an additional
    # power of w.
    s_ = √(β / 2 * wmax) * wmax^(-ypower(kernel)) * s

    # HACK: as we don't yet support Fourier transforms on anything but the
    # unit interval, we need to scale the underlying data.  This breaks
    # the correspondence between U.hat and Uhat though.
    û_base = scale.(u, √β)

    conv_radius = 40 * kernel.Λ
    even_odd = Dict(fermion => :odd, boson => :even)[statistics]
    û = hat.(û_base, even_odd; n_asymp=conv_radius)

    return FiniteTempBasis(kernel, sve_result, statistics, float(β), u_, v_, s_, û)
end

Base.firstindex(::AbstractBasis) = 1
Base.length(basis::AbstractBasis) = length(basis.s)

"""
    iswellconditioned(basis)

Return `true` if the sampling is expected to be well-conditioned.
"""
iswellconditioned(::DimensionlessBasis) = true
iswellconditioned(::FiniteTempBasis) = true

function Base.getindex(basis::FiniteTempBasis, i)
    u, s, v = basis.sve_result
    sve_result = u[i], s[i], v[i]
    return FiniteTempBasis(
        basis.statistics, getbeta(basis), getwmax(basis); kernel=basis.kernel, sve_result
    )
end

"""
    getwmax(basis::FiniteTempBasis)

Real frequency cutoff.
"""
getwmax(basis::FiniteTempBasis) = basis.kernel.Λ / getbeta(basis)

"""
    finite_temp_bases(β, wmax, ε, sve_result=compute_sve(LogisticKernel(β * wmax); ε))

Construct FiniteTempBasis objects for fermion and bosons using the same LogisticKernel instance.
"""
function finite_temp_bases(
    β::AbstractFloat, wmax::AbstractFloat, ε,
    sve_result=compute_sve(LogisticKernel(β * wmax); ε),
)
    basis_f = FiniteTempBasis(fermion, β, wmax, ε; sve_result)
    basis_b = FiniteTempBasis(boson, β, wmax, ε; sve_result)
    return basis_f, basis_b
end

"""
    default_tau_sampling_points(basis)

Default sampling points on the imaginary time/`x` axis.
"""
default_tau_sampling_points(basis::AbstractBasis) = _default_sampling_points(basis.u)

"""
    _default_matsubara_sampling_points(basis; mitigate=true)

Default sampling points on the imaginary frequency axis.
"""
function default_matsubara_sampling_points(basis::AbstractBasis; mitigate=true)
    return _default_matsubara_sampling_points(basis.uhat, mitigate)
end

"""
    default_omega_sampling_points(basis)

Default sampling points on the real-frequency axis.
"""
default_omega_sampling_points(basis::AbstractBasis) = _default_sampling_points(basis.v)

function _default_sampling_points(u)
    poly = last(u)
    maxima = roots(deriv(poly))
    left = (first(maxima) + poly.xmin) / 2
    right = (last(maxima) + poly.xmax) / 2
    return [left; maxima; right]
end

function _default_matsubara_sampling_points(uhat, mitigate=true)
    # Use the (discrete) extrema of the corresponding highest-order basis
    # function in Matsubara.  This turns out to be close to optimal with
    # respect to conditioning for this size (within a few percent).
    polyhat = last(uhat)
    wn = findextrema(polyhat)

    # While the condition number for sparse sampling in tau saturates at a
    # modest level, the conditioning in Matsubara steadily deteriorates due
    # to the fact that we are not free to set sampling points continuously.
    # At double precision, tau sampling is better conditioned than iwn
    # by a factor of ~4 (still OK). To battle this, we fence the largest
    # frequency with two carefully chosen oversampling points, which brings
    # the two sampling problems within a factor of 2.
    if mitigate
        wn_outer = [first(wn), last(wn)]
        wn_diff = 2 * round.(Int, 0.025 * wn_outer)
        length(wn) ≥ 20 && append!(wn, wn_outer - wn_diff)
        length(wn) ≥ 42 && append!(wn, wn_outer + wn_diff)
        unique!(wn)
    end

    if iseven(first(wn))
        pushfirst!(wn, 0)
        unique!(wn)
    end

    sort!(wn)

    return wn
end

function _get_kernel(Λ, kernel)
    if isnothing(kernel)
        kernel = LogisticKernel(Λ)
    elseif kernel.Λ ≉ Λ
        error("kernel.Λ ≉ Λ")
    end
    return kernel
end


text:
module NedelecRefFEsTest

using Test
using Gridap.Polynomials
using Gridap.Fields
using Gridap.TensorValues
using Gridap.Fields: MockField
using Gridap.ReferenceFEs

p = QUAD
D = num_dims(QUAD)
et = Float64
order = 0

reffe = NedelecRefFE(et,p,order)
test_reference_fe(reffe)
@test num_terms(get_prebasis(reffe)) == 4
@test get_order(get_prebasis(reffe)) == 0
@test num_dofs(reffe) == 4

@test Conformity(reffe) == CurlConformity()

p = QUAD
D = num_dims(QUAD)
et = Float64
order = 1

reffe = NedelecRefFE(et,p,order)
test_reference_fe(reffe)
@test num_terms(get_prebasis(reffe)) == 12
@test num_dofs(reffe) == 12
@test get_order(get_prebasis(reffe)) == 1

prebasis = get_prebasis(reffe)
dof_basis = get_dof_basis(reffe)

v = VectorValue(3.0,0.0)
field = GenericField(x->v*x[1])

cache = return_cache(dof_basis,field)
r = evaluate!(cache, dof_basis, field)
test_dof_array(dof_basis,field,r)

cache = return_cache(dof_basis,prebasis)
r = evaluate!(cache, dof_basis, prebasis)
test_dof_array(dof_basis,prebasis,r)

p = TET
D = num_dims(p)
et = Float64
order = 0

reffe = NedelecRefFE(et,p,order)
test_reference_fe(reffe)
@test num_terms(get_prebasis(reffe)) == 6
@test get_order(get_prebasis(reffe)) == 0
@test num_dofs(reffe) == 6
@test Conformity(reffe) == CurlConformity()

p = TRI
D = num_dims(p)
et = Float64
order = 0

reffe = NedelecRefFE(et,p,order)
test_reference_fe(reffe)
@test num_terms(get_prebasis(reffe)) == 3
@test get_order(get_prebasis(reffe)) == 0
@test num_dofs(reffe) == 3
@test Conformity(reffe) == CurlConformity()
dof_basis = get_dof_basis(reffe)

#using Gridap.Geometry
#using Gridap.Visualization
#grid = compute_reference_grid(p,10)
#x = get_node_coordinates(grid)
#shapes = get_shapefuns(reffe)
#ux = evaluate(shapes,x)
#writevtk(grid,"nede",nodaldata=["s$i"=>ux[:,i] for i in 1:num_dofs(reffe)])

# Factory function
reffe = ReferenceFE(QUAD,nedelec,0)
@test num_terms(get_prebasis(reffe)) == 4
@test get_order(get_prebasis(reffe)) == 0
@test num_dofs(reffe) == 4
@test Conformity(reffe) == CurlConformity()

reffe = ReferenceFE(QUAD,nedelec,Float64,0)
@test num_terms(get_prebasis(reffe)) == 4
@test get_order(get_prebasis(reffe)) == 0
@test num_dofs(reffe) == 4
@test Conformity(reffe) == CurlConformity()

@test Conformity(reffe,:L2) == L2Conformity()
@test Conformity(reffe,:Hcurl) == CurlConformity()
@test Conformity(reffe,:HCurl) == CurlConformity()

@test Nedelec() == nedelec

end # module


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct KeyVaultProperties <: SwaggerModel
    keyname::Any # spec type: Union{ Nothing, String } # spec name: keyname
    keyversion::Any # spec type: Union{ Nothing, String } # spec name: keyversion
    keyvaulturi::Any # spec type: Union{ Nothing, String } # spec name: keyvaulturi

    function KeyVaultProperties(;keyname=nothing, keyversion=nothing, keyvaulturi=nothing)
        o = new()
        validate_property(KeyVaultProperties, Symbol("keyname"), keyname)
        setfield!(o, Symbol("keyname"), keyname)
        validate_property(KeyVaultProperties, Symbol("keyversion"), keyversion)
        setfield!(o, Symbol("keyversion"), keyversion)
        validate_property(KeyVaultProperties, Symbol("keyvaulturi"), keyvaulturi)
        setfield!(o, Symbol("keyvaulturi"), keyvaulturi)
        o
    end
end # type KeyVaultProperties

const _property_map_KeyVaultProperties = Dict{Symbol,Symbol}(Symbol("keyname")=>Symbol("keyname"), Symbol("keyversion")=>Symbol("keyversion"), Symbol("keyvaulturi")=>Symbol("keyvaulturi"))
const _property_types_KeyVaultProperties = Dict{Symbol,String}(Symbol("keyname")=>"String", Symbol("keyversion")=>"String", Symbol("keyvaulturi")=>"String")
Base.propertynames(::Type{ KeyVaultProperties }) = collect(keys(_property_map_KeyVaultProperties))
Swagger.property_type(::Type{ KeyVaultProperties }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_KeyVaultProperties[name]))}
Swagger.field_name(::Type{ KeyVaultProperties }, property_name::Symbol) =  _property_map_KeyVaultProperties[property_name]

function check_required(o::KeyVaultProperties)
    true
end

function validate_property(::Type{ KeyVaultProperties }, name::Symbol, val)
end


text:
# #  Fast Lagrangian top.

# Source code: [`fast_top_tut.jl`](fast_top_tut.jl)

# ## Description

# Fast-spinning Lagrange  top simulated by a beam model. The reference
# solution is available in Krysl, P., Endres, L.: Explicit Newmark/Verlet
# algorithm for time integration of the rotational dynamics of rigid bodies.
# Int. J. Numer. Meth. Eng. 62, 2154–2177 (2005).

# ## Goals

# - Illustrate integration of the nonlinear equations
#   of motion with the Newmark algorithm.

using LinearAlgebra
using PlotlyJS

# ## Definition of the basic inputs

# The finite element code realize on the basic functionality implemented in this
# package.

using FinEtools

# The material parameters may be defined with the specification of the units.
# The elastic properties and the mass density are:
E = 71240.0 * phun("MPa")
nu = 0.31
rho = 2.7e3 * phun("kg/m^3")

# The top is a block of square cross-section with dimensions
Width = 60 * phun("mm"); Height = 60 * phun("mm");
Length = 4*Width

Mass = Width*Height*Length*rho;
IL =  (1/12)*Mass*(Width^2+Height^2);
Ib =  (1/12)*Mass*(Length^2+Height^2)+Mass*(Length/2)^2;
Ih =  (1/12)*Mass*(Width^2+Length^2)+Mass*(Length/2)^2;
Omega0 = 313*pi;
R0 = rotmat3([0.05 0 0]);
g = 9.81 * phun("m/sec^2");
q = [0,0,-g*Width*Height*rho];
utol = 1e-3 * phun("m");
maxit = 12
dt = min(2*pi/norm(Omega0)/10, 0.005);
tend = 0.8 * phun("sec");
ng = 1/2; nb = 1/4*(1/2+ng)^2;
# Choose the mass formulation:
mass_type=1;

# ## Cross-section

# Cross-sectional properties are incorporated in the cross-section property. The
# three arguments supplied are functions. All are returning "constants". In
# particular the first two functions each return the dimension of the
# cross-section as a constant(the beam has a uniform cross-section); the third
# function defines the orientation of the cross-section in the global Cartesian
# coordinates. `[1.0, 0.0, 0.0]` is the vector that together with the tangent
# to the midline curve of the beam spans the $x_1x_2$ plane of the local
# coordinates for the beam.
using FinEtoolsFlexStructures.CrossSectionModule: CrossSectionRectangle
cs = CrossSectionRectangle(s -> Width, s -> Width, s -> [1.0, 0.0, 0.0])

# Select the number of elements per leg.
spin_vector = R0*[0, 0, Omega0];
X=[0 0 0;    reshape(R0*[0,0,Length], 1, 3)];
n=2;
tolerance=Length/n/100;

using FinEtoolsFlexStructures.MeshFrameMemberModule: frame_member
members = []
push!(members, frame_member(X, n, cs))

using FinEtoolsFlexStructures.MeshFrameMemberModule: merge_members
fens, fes = merge_members(members; tolerance = tolerance);

# Material properties
using FinEtoolsDeforLinear
material = MatDeforElastIso(DeforModelRed3D, rho, E, nu, 0.0)

# ## Fields

# Now we start constructing the discrete finite element model.
# We begin by constructing the requisite fields, geometry and displacement.
# These are the so-called "configuration variables", all initialized to 0.
# This is that geometry field.
geom0 = NodalField(fens.xyz)
# This is the displacement field, three unknown displacements per node.
u0 = NodalField(zeros(size(fens.xyz, 1), 3))
# This is the rotation field, three unknown rotations per node are represented
# with a rotation matrix, in total nine numbers. The utility function
# `initial_Rfield`
using FinEtoolsFlexStructures.RotUtilModule: initial_Rfield
Rfield0 = initial_Rfield(fens)

# Finally, this is the displacement and rotation field for incremental changes,
# incremental displacements and incremental rotations. In total, 6 unknowns per
# node.
dchi = NodalField(zeros(size(fens.xyz,1), 6))

# ## Support conditions

# The "bottom" of the top is pinned.
box = fill(0.0, 6)
initbox!(box, X[1,:])
supportn = selectnode(fens; box = box, inflate = tolerance)
initbox!(box, X[2,:])
tipn = selectnode(fens; box = box, inflate = tolerance)
for i in [1, 2, 3]
    setebc!(dchi, supportn, true, i)
end
applyebc!(dchi)
numberdofs!(dchi);

# Initial conditions
v0 = deepcopy(dchi) # need the numbering of the degrees of freedom
for i in 1:size(v0.values, 1)
    v0.values[i, 4:6] = spin_vector
end

# For disambiguation we will refer to the stiffness and mass functions by qualifying them with the corotational-beam module, `FEMMCorotBeamModule`.
using FinEtoolsFlexStructures.FEMMCorotBeamModule
CB = FEMMCorotBeamModule
femm = CB.FEMMCorotBeam(IntegDomain(fes, GaussRule(1, 2)), material)
fi = ForceIntensity(q);

using FinEtoolsFlexStructures.RotUtilModule: initial_Rfield, update_rotation_field!
using DelimitedFiles

function integrate(tend, CB, geom0, u0, Rfield0, dchi, v0, report)
    # Make sure we don't clobber any of these variable fields
    geom0, u0, Rfield0, dchi, v0 = deepcopy((geom0, u0, Rfield0, dchi, v0))
    # Additional fields
    stepdchi = deepcopy(dchi)
    u1 = deepcopy(u0)
    v1 = deepcopy(dchi)
    a1 = deepcopy(dchi) # zero out the Acceleration
    a0 = deepcopy(dchi) # zero out the Acceleration
    Rfield1 = deepcopy(Rfield0)
    v0v = gathersysvec(v0)
    a0v = gathersysvec(a0)
    vpv = gathersysvec(v0);
    dchipv = gathersysvec(dchi);
    stepdchiv = gathersysvec(dchi);
    rhs = gathersysvec(dchi);
    TMPv = deepcopy(rhs)
    utol = 1e-13*dchi.nfreedofs;

    t = 0.0;
    step = 0;
    while (t <= tend)
        t = t + dt;
        (mod(step, 50)==0) && println("Time $(t)"); # pause
        # Initialization
        applyebc!(dchi) # Apply boundary conditions
        u1.values[:] = u0.values[:]; # guess
        Rfield1.values[:] = Rfield0.values[:]; # guess
        stepdchi.values[:] .= 0.0;# Total increment in current step
        a1.values[:] = -(1/nb/dt)*v0.values[:] -(1/2-nb)/nb*a0.values[:];
        v1.values[:] = v0.values[:] + dt*((1-ng)*a0.values[:] + ng*a1.values[:]);
        gathersysvec!(v0, v0v)
        gathersysvec!(a0, a0v)
        dchipv = dt*v0v + (dt^2/2*(1-2*nb))*a0v
        vpv = v0v +(dt*(1-ng))*a0v;

        iter = 1;
        while true
            F = CB.distribloads_global(femm, geom0, u1, Rfield1, dchi, fi)
            Fr = CB.restoringforce(femm, geom0, u1, Rfield1, dchi);       # Internal forces
            @. rhs = F + Fr;
            K = CB.stiffness(femm, geom0, u1, Rfield1, dchi);
            M = CB.mass(femm, geom0, u1, Rfield1, dchi);
            G = CB.gyroscopic(femm, geom0, u1, Rfield1, v1, dchi);
            gathersysvec!(stepdchi, stepdchiv)
            @. TMPv = ((-1/(nb*dt^2))*stepdchiv+(1/(nb*dt^2))*dchipv)
            rhs .+= M*TMPv
            @. TMPv = ((-ng/nb/dt)*stepdchiv+(ng/nb/dt)*dchipv - vpv)
            rhs .+= G*TMPv;
            dchi = scattersysvec!(dchi, (K+(ng/nb/dt)*G+(1/(nb*dt^2))*M)\rhs); # Disp. incr
            u1.values[:] += (dchi.values[:,1:3])[:];   # increment displacement
            stepdchi.values[:] += dchi.values[:]
            v1.values[:] += (ng/nb/dt)*dchi.values[:];
            a1.values[:] += (1/nb/dt^2)*dchi.values[:];
            update_rotation_field!(Rfield1, dchi)
            if maximum(abs.(dchi.values[:])) < utol# convergence check
                break;
            end
            if (iter > maxit)# bailout for failed convergence
                error("Possible failed convergence");
            end
            iter += 1;
        end
        u0.values[:] = u1.values[:];       # update the displacement
        Rfield0.values[:] = Rfield1.values[:]; # update the rotations
        v0.values[:] = v1.values[:];       # update the velocities
        a0.values[:] = a1.values[:];       # update the accelerations

        report(step, u1, Rfield1)

        step=step+1;
    end
end

# The visualization utilities take advantage of the PlotlyJS library.
using PlotlyJS
using VisualStructures: plot_space_box, plot_midline, plot_solid, render, react!, default_layout_3d, save_to_json

# Display the graph of the motion of the tip of the top.

function updategraph(step, u1, Rfield1)
    if (mod(step,20)==0)
        push!(tipx, X[2,1]+u1.values[tipn[1], 1])
        push!(tipy, X[2,2]+u1.values[tipn[1], 2])
        plots = cat(tbox, tref, scatter(;x=tipx./Length, y=tipy./Length, mode="markers", name = "Sol", line_color = "rgb(155, 15, 15)"); dims = 1)
        react!(pl, plots, pl.plot.layout)
        sleep(0.01)
    end
end

# integrate(tend, CB, geom0, u0, Rfield0, dchi, v0, updategraph)

tipx = Float64[]
tipy = Float64[]
tipz = Float64[]
layout = Layout(; scene=attr(
    xaxis = attr(title="X"),
    yaxis = attr(title="Y"),
    zaxis = attr(title="Z"),
    camera = attr(
        up=attr(x=-0.8, y=-0.6, z=0.07),
        center=attr(x=0.0, y=0.0, z=0.0),
        eye=attr(x=0.1, y=0.12, z=2.17),
        projection = attr(type = "orthographic")
        )), showlegend = false)
tbox = plot_space_box([[-1.1*Width -1.1*Width 0]; [1.1*Width 1.1*Width 1.1*Length]])
tshape0s = plot_solid(fens, fes; x = geom0.values, u = 0.0.*dchi.values[:, 1:3], R = Rfield0.values, facecolor = "rgb(125, 155, 125)", opacity = 0.3);
tshape0m = plot_midline(fens, fes; x = geom0.values, u = 0.0.*dchi.values[:, 1:3], color = "rgb(125, 105, 175)", lwidth = 4)
plots = cat(tbox,  tshape0s, tshape0m; dims = 1)
pl2 = render(plots; layout = layout)
sleep(0.5)

function updateplot(step, u1, Rfield1)
    push!(tipx, X[2,1]+u1.values[tipn[1], 1])
    push!(tipy, X[2,2]+u1.values[tipn[1], 2])
    push!(tipz, X[2,3]+u1.values[tipn[1], 3])
    if (mod(step,13)==0)
        curv = scatter3d(;x=tipx, y=tipy, z=tipz, mode="lines", name = "Sol", line_color = "rgb(155, 15, 15)")
        tshape1s = plot_solid(fens, fes; x = geom0.values, u = u1.values, R = Rfield1.values, facecolor = "rgb(125, 15, 15)");
        tshape1m = plot_midline(fens, fes; x = geom0.values, u = u1.values, color = "rgb(125, 15, 15)", lwidth = 4)
        plots = cat(tbox,  tshape0s, tshape0m,  curv, tshape1s, tshape1m; dims = 1)
        react!(pl2, plots, pl2.plot.layout)
        sleep(0.12)
    end
end

integrate(4.5*tend, CB, geom0, u0, Rfield0, dchi, v0, updateplot)


text:
using PyPlot

function plot_network(; flow=nothing)
    for i in 1:n
        i0, i1 = orig[i], dest[i]
        plot([absn[i0], absn[i1]], [ordn[i0], ordn[i1]],
             lw= isa(flow, Array{Float64}) ? 20*flow[i] : 5,
             zorder=1, color="darkblue")
        # If the flow is null, plot a red edge.
        if(isa(flow, Array{Float64}) && flow[i] == 0)
            plot([absn[i0], absn[i1]], [ordn[i0], ordn[i1]],
             lw= 5, zorder=1, color="red")
        end
    end

    # Plot reservoirs.
    scatter(absn[1:mr], ordn[1:mr], c="steelblue", s=200, zorder=2, label="Tank")
    # Plot demands.
    scatter(absn[mr+1:m], ordn[mr+1:m], c="k", s=200, zorder=2, label="Demand")

    legend()
end

max_stars_repo_name: frapac/Knitro-Zoo

max_issues_repo_name: frapac/Knitro-Zoo

max_forks_repo_name: frapac/Knitro-Zoo


text:
### A Pluto.jl notebook ###
# v0.18.2

using Markdown
using InteractiveUtils

# ╔═╡ 14f719cd-ffdc-4d9b-83f4-efd7e7b54c82
using Pkg; Pkg.activate(".")

# ╔═╡ c27879d8-9702-11eb-3050-f1c4477273fe
using PlotlyLight

# ╔═╡ 301cbfcf-88f6-4eed-9b75-6cc193242d12
PlotlyLight.Defaults.parent_style[] = "height: 400px;"

# ╔═╡ e2dedee2-9702-11eb-1364-118bcaff2607
Plot(Config(x=1:50,y=randn(50), type="bar"))

# ╔═╡ Cell order:
# ╠═14f719cd-ffdc-4d9b-83f4-efd7e7b54c82
# ╠═c27879d8-9702-11eb-3050-f1c4477273fe
# ╠═301cbfcf-88f6-4eed-9b75-6cc193242d12
# ╠═e2dedee2-9702-11eb-1364-118bcaff2607


text:
module WienerHopf
using ApproxFun, SingularIntegralEquations, RiemannHilbert, ForwardDiff

import ApproxFunBase: SegmentDomain, fromcanonical, tocanonical, angle,
                        reverseorientation, isambiguous, arclength,
                        leftendpoint, rightendpoint, complexlength,
                        AnyDomain, tocanonicalD, fromcanonicalD,
                        setcanonicaldomain, mappoint
import SingularIntegralEquations: stieltjes, stieltjesmoment!, undirected
import Base: convert, isempty, +, -, *, ==
import ForwardDiff: derivative
import RiemannHilbert: RiemannDual, finitepart,
            fpleftstieltjesmoment!, fprightstieltjesmoment!,
            collocationpoints, fpcauchymatrix, evaluationmatrix

export SqrtLine, collocationpoints, fpcauchymatrix, evaluationmatrix

"""
    SqrtLine{a}(c)

represents the line at angle `a` in the complex plane, centred at `c`.
This is encoded by a map from -1..1 that allows rapid expansion of
functions with sqrt-decay.
"""
struct SqrtLine{angle,T<:Number} <: SegmentDomain{T}
    center::T

    #TODO get this inner constructor working again.
    SqrtLine{angle,T}(c) where {angle,T} = new{angle,T}(c)
    SqrtLine{angle,T}() where {angle,T} = new{angle,T}(zero(T))
end

SqrtLine{a}(c::Number) where {a} = SqrtLine{a,typeof(c)}(c)
SqrtLine{a}() where {a} = SqrtLine{a,Float64}()
SqrtLine() = SqrtLine{false}()

angle(d::SqrtLine{a}) where {a} = a*π

reverseorientation(d::SqrtLine{true}) = SqrtLine{false}(d.center)
reverseorientation(d::SqrtLine{false}) = SqrtLine{true}(d.center)
reverseorientation(d::SqrtLine{a}) where {a} = SqrtLine{a-1}(d.center)

isambiguous(d::SqrtLine) = isnan(d.center)
convert(::Type{Domain{T}},d::SqrtLine{a}) where {a,T<:Number} = SqrtLine{a,T}(d.center)
convert(::Type{SqrtLine{a,T}},::AnyDomain) where {a,T<:Number} = SqrtLine{a,T}(NaN)
convert(::Type{IT},::AnyDomain) where {IT<:SqrtLine}=SqrtLine(NaN,NaN)

## Map interval

##TODO non-1 alpha,beta

isempty(::SqrtLine) = false

fourtoone(x) = (x+x^3)/(1-x^2)^2
ifourtoone(α) = if abs(α) ≤ 0.001
        α - 3α^3 + 22α^5 # Taylor series
    else
        (1 + sqrt(1 + 16α^2) - sqrt(2)sqrt(1 + sqrt(1 + 16α^2)))/(4α)
    end

function ifourtoone(k, α)
    x = ifourtoone(α)
    k == 1 && return x
    k == 2 && return inv(x)
    c = x + inv(x) - inv(α)
    k == 3 && return (-c+sqrt(c^2-4))/2
    k == 4 && return (-c-sqrt(c^2-4))/2
    throw(ArgumentError("$k"))
end

sqrtline_tocanonical(x) = ifourtoone(x)
sqrtline_fromcanonical(x) = fourtoone(x)

sqrtline_tocanonicalD(x) = derivative(ifourtoone,x)
sqrtline_fromcanonicalD(x) = derivative(fourtoone,x)

tocanonical(d::SqrtLine,x) = sqrtline_tocanonical(cis(-angle(d)).*(x-d.center))
tocanonical(d::SqrtLine{false},x) = sqrtline_tocanonical(x-d.center)
tocanonical(d::SqrtLine{true},x) = sqrtline_tocanonical(d.center-x)

tocanonicalD(d::SqrtLine,x) = cis(-angle(d)).*sqrtline_tocanonicalD(cis(-angle(d)).*(x-d.center))
tocanonicalD(d::SqrtLine{false},x) = sqrtline_tocanonicalD(x-d.center)
tocanonicalD(d::SqrtLine{true},x) = -sqrtline_tocanonicalD(d.center-x)

fromcanonical(d::SqrtLine,x) = cis(angle(d))*sqrtline_fromcanonical(x)+d.center
fromcanonical(d::SqrtLine{false},x) = sqrtline_fromcanonical(x)+d.center
fromcanonical(d::SqrtLine{true},x) = -sqrtline_fromcanonical(x)+d.center

fromcanonicalD(d::SqrtLine,x) = cis(angle(d))*sqrtline_fromcanonicalD(x)
fromcanonicalD(d::SqrtLine{false},x) = sqrtline_fromcanonicalD(x)
fromcanonicalD(d::SqrtLine{true},x) = -sqrtline_fromcanonicalD(x)

invfromcanonicalD(d::SqrtLine,x) = cis(-angle(d))*sqrtline_invfromcanonicalD(x)
invfromcanonicalD(d::SqrtLine{false},x) = sqrtline_invfromcanonicalD(x)
invfromcanonicalD(d::SqrtLine{true},x) = -sqrtline_invfromcanonicalD(x)

==(d::SqrtLine{a},m::SqrtLine{a}) where {a} = d.center == m.center

# algebra
*(c::Real,d::SqrtLine{false}) = SqrtLine{sign(c)>0 ? false : true}(isapprox(d.center,0) ? d.center : c*d.center,d.α,d.β)
*(c::Real,d::SqrtLine{true}) = SqrtLine{sign(c)>0 ? true : false}(isapprox(d.center,0) ? d.center : c*d.center,d.α,d.β)
*(c::Number,d::SqrtLine) = SqrtLine(isapprox(d.center,0) ? d.center : c*d.center,angle(d)+angle(c),d.α,d.β)
*(d::SqrtLine,c::Number) = c*d
for OP in (:+,:-)
    @eval begin
        $OP(c::Number,d::SqrtLine{a}) where {a} = SqrtLine{a}($OP(c,d.center),d.α,d.β)
        $OP(d::SqrtLine{a},c::Number) where {a} = SqrtLine{a}($OP(d.center,c),d.α,d.β)
    end
end

# algebra
arclength(d::SqrtLine) = Inf
leftendpoint(d::SqrtLine) = -Inf
rightendpoint(d::SqrtLine) = Inf
complexlength(d::SqrtLine) =Inf

# Sum over all inverses of fromcanonical, see [Olver,2014]

fpstieltjes(s,f,z) = finitepart(stieltjes(s,f,RiemannDual(z,1)))

function stieltjes(S::Space{<:SqrtLine},f,z)
    if domain(S) == SqrtLine()
        # TODO: rename tocanonical
        s = setcanonicaldomain(S)
        u = undirected(z)
        stieltjes(s,f,ifourtoone(1,z))+stieltjes(s,f,ifourtoone(2,u))+
            stieltjes(s,f,ifourtoone(3,u))+stieltjes(s,f,ifourtoone(4,u))-
            2fpstieltjes(s,f,1)-2fpstieltjes(s,f,-1)
    else
        stieltjes(setdomain(S,SqrtLine()),f,mappoint(domain(S),SqrtLine(),z))
    end
end

function stieltjesmoment!(ret,S::PolynomialSpace{<:SqrtLine},z,filter=identity)
    if domain(S) == SqrtLine()
        s = setcanonicaldomain(S)
        u = complex(undirected(z))
        tmp = similar(ret)
        stieltjesmoment!(ret, s, ifourtoone(1,z),filter)
        stieltjesmoment!(tmp, s, ifourtoone(2,u),filter); ret .+= tmp
        stieltjesmoment!(tmp, s, ifourtoone(3,u),filter); ret .+= tmp
        stieltjesmoment!(tmp, s, ifourtoone(4,u),filter); ret .+= tmp
        fpleftstieltjesmoment!(tmp, s); ret .-= 2 .* tmp
        fprightstieltjesmoment!(tmp, s); ret .-= 2 .* tmp
        ret
    else
        stieltjesmoment!(ret,setdomain(S,SqrtLine()),mappoint(domain(S),SqrtLine(),z),filter)
    end
end

fprightstieltjesmoment!(V, sp::PolynomialSpace{<:SqrtLine}) = fill!(V,0)
fpleftstieltjesmoment!(V, sp::PolynomialSpace{<:SqrtLine}) = fill!(V,0)

end # module


text:
"""
`COM.prettyprint(x)` prints nicely the object `x` for which `COM.prettify(x)` method is defined.
"""
module COM
prettify(x) = sprint(io -> show(io, "text/plain", x))
prettyprint(io::IO, x) = print(io, prettify(x))
prettyprint(x) = prettyprint(stdout, x)
end

max_stars_repo_name: genkuroki/public

max_issues_repo_name: genkuroki/public

max_forks_repo_name: genkuroki/public


text:
"""
This is where module will be documented

#TODO
# plot l(tau) given AmP parameters as LongParams
# convert to plot of L(t)
# getE0() [solves eb = f, given  uHb]
# runDEBstd(f, p::LongParams, timespan)
  ## runs getE0(), then then plots l(t), sets its own timespan
# add reporting tau_b, l_b, maybe uHb also

"""
module DEBStandardModel

  using DifferentialEquations, Plots, DelimitedFiles
  gr()             #specify GR backend

  struct DEBstd  # specifies std model
    ## matching of function to u0 and parameters rests on user
    fname::Function
    growthparams::Vector{Float64}
    maturityparams::Vector{Float64}
    u0::Vector{Float64}
  end

  function stdDEBscaled(du, u, p, t)  # eqns 2.26 to 2.28, DEB book 3rd ed
    # requires the parameter uE0 for initial condition
    # and juvreached better reset via callback at tau_b
    # and adultreached via callback on uH
    # and perhaps u[4] for reproductive output?
    g, kap, k, lT, f, juvreached, adultreached  = p
    uE, l, uH = u
    du[1] = juvreached*f*l^2 - uE*l^2 * (g+l)/(uE+l^3)
    du[2] = (1.0/3.0)*(g*uE - l^4 - lT)/(uE + l^3)
    du[3] = (1-adultreached)*((1-kap) * uE*l^2 * (g+l)/(uE+l^3) - k*uH )
  end

function makeparsdict(csvfilename)
  lines = readdlm(csvfilename,',')
  parsdict = Dict(lines[i,1] => lines[i,2] for i = 1:size(lines)[1])
  return parsdict
  end

  function packDEBscaled_std(parsdict, fname = stdDEBscaled)
    # IN:
    ## parsdict     Dict, keys are AmP symbols (as strings)
    ##
    #
    # OUT:
    ## p::DEBstd

    #unpack parsdict to parameter names
    T_ref = parsdict["T_ref"]
    E_Hb = parsdict[ "E_Hb"]
    s_G = parsdict[ "s_G"]
    p_Am = parsdict["p_Am"]
    T_body = parsdict["T_body"]
    E_Hp = parsdict["E_Hp"]
    p_M = parsdict["p_M"]
    T_A = parsdict["T_A"]
    z = parsdict["z"]
    kap_P = parsdict["kap_P"]
    h_a = parsdict["h_a"]
    v = parsdict["v"]
    E_G = parsdict["E_G"]
    kap_R = parsdict["kap_R"]
    kap_X = parsdict["kap_X"]
    k_J = parsdict["k_J"]
    F_m = parsdict["F_m"]
    p_T = parsdict["p_T"]
    kap = parsdict["kap"]
    f   = parsdict["f"]
    E0  = parsdict["E0"]

    # compute derived parameters
    k_M = p_M/E_G
    E_m = p_Am/v
    g   = E_G/(kap*E_m)
    L_m = v/(g*k_M)
    L_T = p_T/p_M
    k   = k_J/k_M
    lT  = L_T/L_m
    feedingyes, adultno = 0.0, 0.0
    uHb = E_Hb/(g*E_m*L_m^3)
    uHp = E_Hp/(g*E_m*L_m^3)
    uE0 = E0/(g*E_m*L_m^3)

    #pack the relevant vectors
    growthpars   = [f, g, kap, k, lT, feedingyes, adultno]
    maturitypars = [uHb, uHp]
    u0           = [uE0, 0.0, 0.0]

    return DEBstd(fname, growthpars, maturitypars, u0)

  end
  ## Hons course version: stdAmPparams = [f, g, kap, k, lT, feedingyes, adultno]
  ## while uE0 is used to set initial condition
  ## and uHb, uHb are used to set callbacks for solve
  ## 11 params in all
  ## maybe make LongParams named tuple (stdAmParams, maturitylevels, uE0)
  ## and have a dummy/output function that sets up problem, solves it, plots l(t)

function output(p::DEBstd, tspan)
  f, growthpars, maturitypars, u0 = p.fname, p.growthparams, p.maturityparams, p.u0
  prob = ODEProblem(f, u0, tspan, growthpars)

  #create stdDEBcallback
  uhb, uhp = maturitypars
  function reached_uHb(u, t, integrator)     u[3] - uhb         end
  function startfeeding!(integrator)
    integrator.p[6] = 1.0
    @show integrator.t    #uncomment this to see age at birth
    end
  function reached_uHp(u, t, integrator) u[3] - uhp end
  function nowadult!(integrator) integrator.p[7] = 0.0 end
  cbonuHb = ContinuousCallback(reached_uHb, startfeeding!)
  cbonuHp = ContinuousCallback(reached_uHp, nowadult!)
  stdDEBcallback = CallbackSet(cbonuHb, cbonuHp)

  #solve and plot
  sol = solve(prob, cb=stdDEBcallback)
  plot(sol, vars=2 )
  # plot(title = "labelslater" ) #maybe add timestamp?

  end

function rundemo(datafile, tspan) # and another method for f as function?
  parsdict = makeparsdict(datafile)
  runparams = packDEBscaled_std(parsdict)
  output(runparams, tspan)
  end

# test: rundemo("Msplendens_AmPfittedparams.csv", [0.0, 20.0]) should produce a specific figure

end # module

#=
NOTES

First target: plot entire l(t) given the AmP parameter values, also
return age at birth and length at birth.
Ten LongParams values: seven for growth model, two for maturity switches, one for initial state.
Environment via f is constant.

Second target: determine E0, given E_H^b and assuming e0 = f (or, if f is variable, then
equal to average of f)

Long term target: fit all parameters,

=#

max_stars_repo_name: HenriLaurie/DEBStandardModel

max_issues_repo_name: HenriLaurie/DEBStandardModel

max_forks_repo_name: HenriLaurie/DEBStandardModel