text:
_scip = JuMP.optimizer_with_attributes(SCIP.Optimizer, "display/verblevel" => 0)
_glpk = GLPK.Optimizer

@testset "solve master problem: ground structure, discrete diameters" begin
    #       7    9      even arc numbers for
    #   () - d2 - ()    reversed arcs
    #   /1   /3   /5
    #  s1 - () - d1
    #    11   13
    inst = Instance([Node(0,0), Node(40,0), Node(20, 20)],
                    [-50, 20, 30],
                    fill(Bounds(60,80), 3),
                    [Diameter(t...) for t in [(0.8, 1.0),(1.0, 1.2)]])
    topo = squaregrid(2, 3, 20.0, antiparallel=true)
    model, y, z, q = GndStr.make_master(inst, topo, _scip)

    JuMP.optimize!(model)
	status = JuMP.termination_status(model)
    @test status == MOI.OPTIMAL

    ysol = JuMP.value.(y)
    zsol = JuMP.value.(z)
    qsol = JuMP.value.(q)

    # shortest tree is obvious:
    @test ysol[11] ≈ 1.0
    @test ysol[13] ≈ 1.0
    @test ysol[4] ≈ 1.0
    @test sum(ysol) ≈ 3.0    # all others 0

    # only the smallest diameter is chosen
    @test zsol[11,1] ≈ 1.0
    @test zsol[13,1] ≈ 1.0
    @test zsol[4,1] ≈ 1.0
    @test sum(zsol) ≈ 3.0    # all others 0

    # uniq flow solution
    @test qsol[11] ≈ 50.0
    @test qsol[13] ≈ 20.0
    @test qsol[4] ≈ 30.0
    @test sum(qsol) ≈ 100.0  # all others 0
end

@testset "solve subproblem: ground structure, discrete diameters" begin
    #       7    9      even arc numbers for
    #   () - d2 - ()    reversed arcs
    #   /1   /3   /5
    #  s1 - () - d1
    #    11   13
    inst = Instance([Node(0,0), Node(40,0), Node(20, 20)],
                    [-50, 20, 30],
                    fill(Bounds(60,80), 3),
                    [Diameter(t...) for t in [(0.8, 1.0),(1.0, 1.2)]])
    topo = squaregrid(2, 3, 20.0, antiparallel=true)

    nnodes = 2*3
    narcs = 2*7
    ndiam = length(inst.diameters)

    # recreate solution from master problem
    zsol = fill(false, (narcs, ndiam))
    zsol[11,1] = true
    zsol[13,1] = true
    zsol[4,1] = true

    qsol = fill(0.0, narcs)
    qsol[11] = 50.0
    qsol[13] = 20.0
    qsol[4] = 30.0

    @testset "feasible subproblem" begin
        cand = GndStr.CandSol(zsol, qsol, fill(0.0, narcs))
        model, π, Δl, Δu, ploss, plb, pub = GndStr.make_sub(
            inst, topo, cand, _glpk)
        JuMP.optimize!(model)
	    status = JuMP.termination_status(model)
        @test status == MOI.OPTIMAL

        # primal solution
        πsol = JuMP.value.(π)
        Δlsol = JuMP.value.(Δl)
        Δusol = JuMP.value.(Δu)

        @test Δlsol ≈ zeros(nnodes)
        @test Δusol ≈ zeros(nnodes)

        # dual solution
        μsol = JuMP.dual.(ploss)
        λlsol = JuMP.dual.(plb)
        λusol = JuMP.dual.(pub)

        @test μsol ≈ zeros(length(ploss))
        @test λlsol ≈ zeros(nnodes)
        @test λusol ≈ zeros(nnodes)
    end

    @testset "infeasible subproblem" begin
        cand = GndStr.CandSol(zsol, 10 * qsol, fill(0.0, narcs)) # scaled
        model, π, Δl, Δu, ploss, plb, pub = GndStr.make_sub(
            inst, topo, cand, _glpk)
        JuMP.optimize!(model)
	    status = JuMP.termination_status(model)
        @test status == MOI.OPTIMAL

        # primal solution
        πsol = JuMP.value.(π)
        Δlsol = JuMP.value.(Δl)
        Δusol = JuMP.value.(Δu)

        # complementarity
        @test Δlsol .* Δusol ≈ zeros(nnodes)

        # innodes should not have slack
        terms = [2, 3, 6]
        innodes = filter(v -> !(v in terms),  1:nnodes)
        @test Δlsol[innodes] ≈ zeros(length(innodes))
        @test Δusol[innodes] ≈ zeros(length(innodes))

        # at least two vertices have slack
        slack = Δlsol + Δusol
        @test sum(slack[terms] .> 0) >= 2

        # dual solution, TODO: fix sign of dual multipliers
        μsol = abs.(JuMP.dual.(ploss))
        λlsol = abs.(JuMP.dual.(plb))
        λusol = abs.(JuMP.dual.(pub))

        # at least two bounds active
        @test sum(λlsol[terms] .> 0) >= 1
        @test sum(λusol[terms] .> 0) >= 1

        # at least one path active
        @test sum(μsol .> 0) >= 2

        @test λlsol[innodes] ≈ zeros(length(innodes))
        @test λusol[innodes] ≈ zeros(length(innodes))
    end
end

@testset "compare relaxation and exact for subproblem" begin
    #     __s1__  __s2
    #   t3      t4
    nodes = [Node(100,0), Node(300,0), Node(0,0), Node(200,0)]
    arcs = [Arc(1,3), Arc(1,4), Arc(2,4)]
    topo = Topology(nodes, arcs)

    demand = [-600, -400, 400, 600]
    bounds = fill(Bounds(40, 80), 4)
    diams = [Diameter(1.0, 1.0), Diameter(2.0, 2.0)]
    inst = Instance(nodes, demand, bounds, diams, ploss_coeff_nice)

    zsol = [true false; true false; true false]
    qsol = [400, 200, 400]
    cand = GndStr.CandSol(zsol, qsol, qsol.^2)

    @testset "solving the exact subproblem" begin
        model, π, Δl, Δu, ploss, plb, pub =
            GndStr.make_sub(inst, topo, cand, _glpk, relaxed=false)
        JuMP.optimize!(model)
	    status = JuMP.termination_status(model)
        @test status == MOI.OPTIMAL
        @test JuMP.objective_value(model) ≈ 3600
    end

    @testset "solving the relaxation" begin
        model, π, Δl, Δu, ploss, plb, pub =
            GndStr.make_sub(inst, topo, cand, _glpk, relaxed=true)
        JuMP.optimize!(model)
	    status = JuMP.termination_status(model)
        @test status == MOI.OPTIMAL
        @test JuMP.objective_value(model) ≈ 0
    end
end

@testset "run GBD iterations" begin
    #       7    9      even arc numbers for
    #   () - d2 - ()    reversed arcs
    #   /1   /3   /5
    #  s1 - () - d1
    #    11   13

    # slightly irregular grid
    nodes = [Node(0,0), Node(45,0), Node(25, 22)]
    demand = [-50, 20, 30]
    bounds = fill(Bounds(60,80), 3)
    diams = [Diameter(t...) for t in [(0.8, 1.0),(1.0, 1.2)]]

    topo = Topology([Node(t...) for t in [(0,22), (0,0), (25,22), (25,0), (45,22), (45,0)]],
                    [Arc(t...) for t in [(1,2), (2,1), (3,4), (4,3), (5,6), (6,5),
                                         (1,3), (3,1), (3,5), (5,3),
                                         (2,4), (4,2), (4,6), (6,4)]])

    @testset "low flow: very easy instance" begin
        inst = Instance(nodes, 1*demand, bounds, diams)

        result = optimize(inst, topo, GndStr.IterGBD(_scip, _glpk))
        @test result.status == MOI.OPTIMAL

        zsol = result.solution.zsol
        @test zsol[4,1] == true
        @test zsol[11,1] == true
        @test zsol[13,1] == true
        @test sum(zsol) == 3

        @test result.dualbound ≈ 67.0
        @test result.niter == 1
    end

    @testset "medium flow: difficult instance" begin
        inst = Instance(nodes, 5*demand, bounds, diams)

        result = optimize(inst, topo, GndStr.IterGBD(_scip, _glpk))
        @test result.status == MOI.OPTIMAL

        zsol = result.solution.zsol
        @test zsol[4,1] == true
        @test zsol[11,2] == true
        @test zsol[13,1] == true
        @test sum(zsol) == 3

        @test result.dualbound ≈ 72.0
        @test result.niter == 3
    end

    @testset "high flow: iteration limit instance" begin
        inst = Instance(nodes, 30*demand, bounds, diams)

        result = optimize(inst, topo, GndStr.IterGBD(_scip, _glpk, maxiter=3))
        @test result.status == MOI.ITERATION_LIMIT
        @test result.solution == nothing
        @test result.dualbound ≈ 156.0
        @test result.niter == 3
    end

    @testset "high flow: time limit instance" begin
        inst = Instance(nodes, 30*demand, bounds, diams)

        result = optimize(inst, topo, GndStr.IterGBD(_scip, _glpk, timelimit=5.0))
        @test result.status == MOI.TIME_LIMIT
        @test result.solution == nothing
    end

    @testset "high flow on triangle: infeasible" begin
        inst3 = Instance([Node(0,0), Node(50,0)],
                         20*[-50, 50],
                         [Bounds(60,80), Bounds(60,80)],
                         [Diameter(t...) for t in [(0.8, 1.0),(1.0, 1.2)]])
        topo3 = Topology([Node(0,0), Node(50,0), Node(30, 40)],
                         [Arc(1,3), Arc(1,2), Arc(2,3)])

        result = optimize(inst3, topo3, GndStr.IterGBD(_scip, _glpk))
        @test result.status == MOI.INFEASIBLE
        @test result.solution == nothing
        @test result.dualbound == Inf
        @test result.niter == 2
    end

    @testset "on the example with subproblem/relaxation gap" begin
        #     __s1__  __s2
        #   t3      t4
        nodes = [Node(100,0), Node(300,0), Node(0,0), Node(200,0)]
        arcs = [Arc(1,3), Arc(1,4), Arc(2,4)]
        topo = Topology(nodes, arcs)

        demand = [-600, -400, 400, 600]
        bounds = fill(Bounds(40, 80), 4)
        diams = [Diameter(1.0, 1.0), Diameter(2.0, 2.0)]
        inst = Instance(nodes, demand, bounds, diams, ploss_coeff_nice)

        result = optimize(inst, topo, GndStr.IterGBD(_scip, _glpk))
        @test result.status == MOI.OPTIMAL
        zsol = result.solution.zsol
        @test sum(zsol[:,2]) == 1
        @test sum(zsol[:,1]) == 2 # solution not quite uniqe
        @test result.dualbound ≈ 400.0
        @test result.niter == 2
    end

    @testset "difficult instance with disconnected candidates" begin
        inst = Instance([Node(0,0), Node(100,0), Node(200,100)],
                        [800, -900, 100],
                        fill(Bounds(40,80), 3),
                        [Diameter(t...) for t in [(1.0, 1.0), (2.0, 3.2)]],
                        ploss_coeff_nice)
        topo = squaregrid(2, 3, 100.0, antiparallel=true)

        # trigger the cuts for disconnected candidate
        solver = GndStr.IterGBD(_scip, _glpk, addnogoods=true, addcritpath=false)
        result = optimize(inst, topo, solver)
        @test result.status == MOI.OPTIMAL

        zsol = result.solution.zsol
        @test sum(zsol) == 3

        @test result.dualbound ≈ 520.0
    end

end

@testset "Linear overestimation of supremum terms" begin
    values = Float64[i*j - 1 for i=1:4, j=1:4]
    values = 2*tril(values) - triu(values)

    cand_i, cand_j = 3, 2
    a, b, c = GndStr.linear_overest(values, cand_i, cand_j, _glpk)
    @test size(a) == (4,)
    @test size(b) == (4,)
    @test size(c) == ()

    @test a[cand_i] + b[cand_j] + c == values[cand_i, cand_j]
end

@testset "Solve semimaster without z vars" begin
    #       7    9      even arc numbers for
    #   () - d2 - ()    reversed arcs
    #   /1   /3   /5
    #  s1 - () - d1
    #    11   13
    inst = Instance([Node(0,0), Node(40,0), Node(20, 20)],
                    [-50, 20, 30],
                    fill(Bounds(60,80), 3),
                    [Diameter(t...) for t in [(0.8, 1.0),(1.0, 1.2)]])
    topo = squaregrid(2, 3, 20.0, antiparallel=true)
    model, y, q = GndStr.make_semimaster(inst, topo, _scip)

    JuMP.optimize!(model)
	status = JuMP.termination_status(model)
    @test status == MOI.OPTIMAL

    ysol = JuMP.value.(y)
    qsol = JuMP.value.(q)

    # shortest tree is obvious:
    @test ysol[11] ≈ 1.0
    @test ysol[13] ≈ 1.0
    @test ysol[4] ≈ 1.0
    @test sum(ysol) ≈ 3.0 # all others 0

    # uniq flow solution
    @test qsol[11] ≈ 50.0
    @test qsol[13] ≈ 20.0
    @test qsol[4] ≈ 30.0
    @test sum(qsol) ≈ 100.0 # all others 0
end

@testset "Solve semisubproblem with free z vars" begin
    #       7    9      even arc numbers for
    #   () - d2 - ()    reversed arcs
    #   /1   /3   /5
    #  s1 - () - d1
    #    11   13

    # slightly irregular grid
    nodes = [Node(0,0), Node(45,0), Node(25, 22)]
    demand = [-50, 20, 30]
    bounds = fill(Bounds(60,80), 3)
    diams = [Diameter(t...) for t in [(0.8, 1.0),(1.0, 1.2)]]

    topo = Topology([Node(t...) for t in [(0,22), (0,0), (25,22), (25,0), (45,22), (45,0)]],
                    [Arc(t...) for t in [(1,2), (2,1), (3,4), (4,3), (5,6), (6,5),
                                         (1,3), (3,1), (3,5), (5,3),
                                         (2,4), (4,2), (4,6), (6,4)]])

    zsol = fill(false, length(topo.arcs), length(diams))
    zsol[4,1]  = true
    zsol[11,1] = true
    zsol[13,1] = true
    qsol = fill(0.0, length(topo.arcs))
    qsol[[4, 11, 13]] = [30.0, 50.0, 20.0]

    @testset "low flow: very easy instance" begin
        factor = 1.0
        inst = Instance(nodes, factor*demand, bounds, diams)

        sol = GndStr.CandSol(zsol, factor*qsol, qsol.^2)
        model, candarcs, z = GndStr.make_semisub(
            inst, topo, sol, _scip)
        @test length(candarcs) == 3
        @test size(z) == (3, 2)

        JuMP.optimize!(model)
	    status = JuMP.termination_status(model)
        @test status == MOI.OPTIMAL

        znew = fill(false, length(topo.arcs), length(diams))
        znew[candarcs,:] = (JuMP.value.(z) .> 0.5)
        @test znew[4,1] == true
        @test znew[11,1] == true
        @test znew[13,1] == true
        @test sum(znew) == 3
    end

    @testset "medium flow: difficult instance" begin
        factor = 5.0
        inst = Instance(nodes, factor*demand, bounds, diams)

        sol = GndStr.CandSol(zsol, factor*qsol, qsol.^2)
        model, candarcs, z = GndStr.make_semisub(
            inst, topo, sol, _scip)
        @test length(candarcs) == 3
        @test size(z) == (3, 2)

        JuMP.optimize!(model)
	    status = JuMP.termination_status(model)
        @test status == MOI.OPTIMAL

        znew = fill(false, length(topo.arcs), length(diams))
        znew[candarcs,:] = (JuMP.value.(z) .> 0.5)
        @test znew[4,1] == true
        @test znew[11,2] == true
        @test znew[13,1] == true
        @test sum(znew) == 3

    end

    @testset "high flow on triangle: infeasible" begin
        inst = Instance([Node(0,0), Node(50,0)],
                        [-1000, 1000],
                        [Bounds(60,80), Bounds(60,80)],
                        [Diameter(t...) for t in [(0.8, 1.0),(1.0, 1.2)]])
        topo = Topology([Node(0,0), Node(50,0), Node(30, 40)],
                        [Arc(1,3), Arc(1,2), Arc(2,3)])

        zsol = fill(false, 3, 2)
        zsol[1,1] = true
        qsol = 1000*[1.0, 0.0, 0.0]
        sol = GndStr.CandSol(zsol, qsol, qsol.^2)
        model, candarcs, z = GndStr.make_semisub(
            inst, topo, sol, _scip)
        @test length(candarcs) == 1
        @test size(z) == (1, 2)

        JuMP.optimize!(model)
	    status = JuMP.termination_status(model)
        @test status == MOI.INFEASIBLE
    end

end

@testset "Solve semi decomposition with nogoods on y" begin
    #       7    9      even arc numbers for
    #   () - d2 - ()    reversed arcs
    #   /1   /3   /5
    #  s1 - () - d1
    #    11   13

    # slightly irregular grid
    nodes = [Node(0,0), Node(45,0), Node(25, 22)]
    demand = [-50, 20, 30]
    bounds = fill(Bounds(60,80), 3)
    diams = [Diameter(t...) for t in [(0.8, 1.0),(1.0, 1.2)]]

    topo = Topology([Node(t...) for t in [(0,22), (0,0), (25,22), (25,0), (45,22), (45,0)]],
                    [Arc(t...) for t in [(1,2), (2,1), (3,4), (4,3), (5,6), (6,5),
                                         (1,3), (3,1), (3,5), (5,3),
                                         (2,4), (4,2), (4,6), (6,4)]])

    @testset "low flow: very easy instance" begin
        inst = Instance(nodes, 1 * demand, bounds, diams)
        result = optimize(inst, topo, GndStr.IterTopo(_scip, _scip))
        @test result.status == MOI.OPTIMAL

        zsol = result.solution.zsol
        qsol = result.solution.qsol

        # shortest tree is obvious, smallest diameter enough:
        @test zsol[11, 1] ≈ 1.0
        @test zsol[13, 1] ≈ 1.0
        @test zsol[4, 1] ≈ 1.0
        @test sum(zsol) ≈ 3.0 # all others 0

        # uniq flow solution
        @test qsol[11] ≈ 50.0
        @test qsol[13] ≈ 20.0
        @test qsol[4] ≈ 30.0
        @test sum(qsol) ≈ 100.0 # all others 0
    end

    @testset "medium flow: difficult instance" begin
        inst = Instance(nodes, 10 * demand, bounds, diams)
        result = optimize(inst, topo, GndStr.IterTopo(_scip, _scip))
        @test result.status == MOI.OPTIMAL

        zsol = result.solution.zsol
        qsol = result.solution.qsol

        # shortest tree is obvious, need one large diameter
        @test sum(zsol[2, :]) == 1
        @test sum(zsol[7, :]) == 1
        @test sum(zsol[11, :]) == 1
        @test sum(zsol[13, :]) == 1
        @test sum(zsol[:, 1]) == 1
        @test sum(zsol[:, 2]) == 3
        @test sum(zsol) ≈ 4.0 # all others 0
    end

    @testset "high flow on triangle: infeasible instance" begin
        inst = Instance([Node(0,0), Node(50,0)],
                        [-1000, 1000],
                        [Bounds(60,80), Bounds(60,80)],
                        [Diameter(t...) for t in [(0.8, 1.0),(1.0, 1.2)]])
        topo = Topology([Node(0,0), Node(50,0), Node(30, 40)],
                        [Arc(1,3), Arc(1,2), Arc(2,3),
                         Arc(3,1), Arc(2,1), Arc(3,2)])

        result = optimize(inst, topo, GndStr.IterTopo(_scip, _scip))
        @test result.status == MOI.INFEASIBLE
    end
end


text:
"""
  Kinetic energy
"""
function Eₖ(uh,h,dΩ)
  0.5*sum(∫(uh⋅uh*h)dΩ)
end

"""
  Potential energy
"""
function Eₚ(a,b,c,dΩ)
  c*sum(∫(a*b)dΩ)
end

"""
  Total energy
"""
function Eₜ(uh,H,hh,g,dΩ)
  Eₖ(uh,H,dΩ)+Eₚ(hh,hh,0.5*g,dΩ)
end

"""
  Kinetic to potential
"""
function compute_kin_to_pot!(w,unv,divvh,hnv)
  mul!(w,divvh,hnv)
  unv⋅w
end

"""
  Potential to kinetic
"""
function compute_pot_to_kin!(w,hnv,qdivu,unv)
   mul!(w,qdivu,unv)
   hnv⋅w
end

"""
  Total Mass
"""
function compute_total_mass!(w,L2MM,hh)
  mul!(w,L2MM,hh)
  sum(w)
end

"""Write scalar diagnostics to csv.
Diagnostics should be added as kwargs: field=values"""
function write_to_csv(csv_file_path; kwargs...)
  df = DataFrame(kwargs... )
  header = names(df)
  CSV.write(csv_file_path, df, header=header)
end

"""Append line to existing csv file.
row should be given as kwargs: field=value"""
function append_to_csv(csv_file_path; kwargs...)
  df = DataFrame(kwargs )
  CSV.write(csv_file_path, df, delim=",", append=true)
end

"""Create a .cvs file header. header should be specified as args"""
function initialize_csv(csv_file_path, args...)
  header = [String(_) for _ in args]
  CSV.write(csv_file_path,[], writeheader=true, header=header)
end

"""Wrapper to get a vector from a csv field. The fieldname argument should be passed as a symbol,
i.e. :<fieldname>"""
function get_scalar_field_from_csv(csv_file_path, fieldname)
  t = CSV.read(csv_file_path, DataFrame)
  getproperty(t, fieldname)
end

"""
  Full diagnostics for the shallow water equations (mass, vorticity, kinetic energy, potential energy, power)
"""

function compute_diagnostics_shallow_water!(h_tmp, w_tmp,
                                         model, dΩ, dω, S, L2MM, H1MM,
					 h, u, w, ϕ, F, h2, c)
  mass_i = compute_total_mass!(h_tmp, L2MM, get_free_dof_values(h))
  vort_i = compute_total_mass!(w_tmp, H1MM, get_free_dof_values(w))
  kin_i  = Eₖ(u,h,dΩ)
  pot_i  = Eₚ(h,h2,c,dΩ)
  pow_i  = sum(∫(ϕ*DIV(F))dω)

  mass_i, vort_i, kin_i, pot_i, pow_i
end

function dump_diagnostics_shallow_water!(h_tmp, w_tmp,
                                         model, dΩ, dω, S, L2MM, H1MM,
                                         h, u, w, ϕ, F, g, step, dt,
                                         output_file, dump_on_screen)

  mass_i, vort_i, kin_i, pot_i, pow_i = compute_diagnostics_shallow_water!(
                                          h_tmp, w_tmp,
                                          model, dΩ, dω, S, L2MM, H1MM,
                                          h, u, w, ϕ, F, h, 0.5*g)

  append_to_csv(output_file;
                time       = step*dt/24/60/60,
                mass       = mass_i,
                vorticity  = vort_i,
                kinetic    = kin_i,
                potential  = pot_i,
                power      = pow_i)

  if dump_on_screen
    @printf("%5d %14.9e %14.9e %14.9e %14.9e %14.9e %14.9e\n",
             step, mass_i, vort_i, kin_i, pot_i, kin_i+pot_i, pow_i)
  end
end

function dump_diagnostics_thermal_shallow_water!(h_tmp, w_tmp,
                                         model, dΩ, dω, S, L2MM, H1MM,
                                         h, u, E, w, ϕ, F, eF, step, dt,
                                         output_file, dump_on_screen)

  mass_i, vort_i, kin_i, pot_i, pow_k2p_i = compute_diagnostics_shallow_water!(
                                              h_tmp, w_tmp,
                                              model, dΩ, dω, S, L2MM, H1MM,
                                              h, u, w, ϕ, F, E, 0.5)

  buoy_i    = compute_total_mass!(h_tmp, L2MM, get_free_dof_values(E))
  pow_k2i_i = 0.5*sum(∫(h*DIV(eF))dω)

  append_to_csv(output_file;
                time       = step*dt/24/60/60,
                mass       = mass_i,
                vorticity  = vort_i,
                buoyancy   = buoy_i,
                kinetic    = kin_i,
                internal   = int_i,
                power_k2p  = pow_k2p_i,
                power_k2i  = pow_k2i_i)

  if dump_on_screen
    @printf("%5d %14.9e %14.9e %14.9e %14.9e %14.9e %14.9e %14.9e %14.9e\n",
             step, mass_i, vort_i, buoy_i, kin_i, int_i, kin_i+int_i, pow_k2p_i, pow_k2i_i)
  end
end

ext: jl


text:
import RvSpectMLBase
using Test

@testset "Spectra1DBasic" begin
    using RvSpectMLBase.TheoreticalInstrument
    n1 = TheoreticalInstrument1D()
    lambdas = TheoreticalInstrument.calc_λs(n1)
    num_pixels = length(lambdas)
    spec = Spectra1DBasic(lambdas, ones(num_pixels), ones(num_pixels), n1)
end

@testset "Spectra2DBasic" begin
    num_pixels = 1000
    num_orders = 10
    lambda2d = collect(reshape(range(5000.0,stop=6000.0,length=num_pixels*num_orders),num_pixels,num_orders))
    n2 = TheoreticalInstrument2D(lambda2d)
    spec = TheoreticalInstrument.demo_generate_spectrum_line(n2)
    @test get_λ_range(spec).min ≈ 5000
    @test get_λ_range(spec).max ≈ 6000

    df = DataFrame(:lambda_lo=>[5501], :lambda_hi=>[5505])
    cl = make_chunk_list_around_lines(spec, df)
    @test_nowarn calc_normalization(cl)
    norm = calc_normalization(cl)
    @test_nowarn normalize_spectrum!(spec,norm )
    @test calc_normalization(cl) ≈ 1
    @test_nowarn calc_normalization_var_weighted(cl)
    norm = calc_normalization(cl)
    @test_nowarn normalize_spectrum!(spec,norm )
    @test calc_normalization_var_weighted(cl) ≈ 1
    @test_nowarn apply_doppler_boost!(spec,1.001)
    @test_nowarn RvSpectMLBase.discard_large_metadata(spec)
end

max_stars_repo_name: RvSpectML/RvSpectMLTemplate

max_issues_repo_name: RvSpectML/RvSpectMLTemplate

max_forks_repo_name: RvSpectML/RvSpectMLTemplate


text:
"""
    struct Shortcut{T}
        chain::T
        f::Function
        σ::Function
    end
Shortcut connection for ResNet-type blocks
"""
struct Shortcut{T}
    chain::T
    f::Function
    σ::Function
end

Shortcut(chain::T) where {T} = Shortcut{typeof(chain)}(chain, +, tanh)

Flux.@functor Shortcut

(nn::Shortcut)(x) = nn.σ.(nn.f(nn.chain(x), x))

function Base.show(io::IO, model::Shortcut{T}) where {T}
    print(
        io,
        "Shortcut{$T}\n",
        "chain: $(model.chain)\n",
        "connection: $(model.f)\n",
        "activation: $(model.σ)\n",
    )
end


text:
"""
An unofficial [Julia](https://julialang.org/) library to communicate
with [Itead Nextion display](https://nextion.itead.cc/).

The Nextion.jl package is licensed under the MIT "Expat" License:

Copyright (c) 2018: [Sébastien Celles](https://github.com/scls19fr/).

```
Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in all
copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE.
```

Source: [https://github.com/scls19fr/Nextion.jl](https://github.com/scls19fr/Nextion.jl)

Contributors: [https://github.com/scls19fr/Nextion.jl/graphs/contributors](https://github.com/scls19fr/Nextion.jl/graphs/contributors)

Reference documents for implementation:
- [https://nextion.itead.cc/resources/documents/instruction-set/](https://nextion.itead.cc/resources/documents/instruction-set/)
- [https://www.itead.cc/wiki/Nextion_Instruction_Set](https://www.itead.cc/wiki/Nextion_Instruction_Set)
"""
module Nextion
    import Base: show, read, write, close, reset,
                    setproperty!, getproperty,
                    getindex, push!
                    Int32, UInt16, length, first
    import Sockets: send

    # hardware
    export NexSerial
    #export PageID, ComponentID, Name
    export init
    #  Base: close, reset
    export send  # to avoid WARNING: Base.send is deprecated

    export Name, PageID, ComponentID

    # drawing
    export Draw

    # contants
    export Return
    export Alignment, Background,
        Colour, Scroll, Baudrate, GPIO

    # events
    export Event

    # resources
    export Font, Picture

    # controls
    export NexButton
    export NexCheckbox
    export NexCrop
    export NexDualStateButton
    export NexGauge
    export NexHotspot
    export NexNumber
    export NexPage, isshown
    #  Base: show
    export NexPicture
    export NexProgressBar
    export NexQRcode
    export NexRadio
    export NexScrollText
    export NexSlider
    export NexText
    export NexWaveform

    # invisibles
    export NexTimer
    export NexUpload
    export NexVariable

    # enhanced
    export NexEEPROM
    export NexGPIO
    export NexRTC

    # more
    export Backlight

    # package code goes here
    include("constants.jl")
    include("resources.jl")
    include("hardware.jl")

    include("draw/draw.jl")

    include("object.jl")
    include("range_number.jl")
    include("events.jl")
    include("nexproperties.jl")
    include("interfaces/interfaces.jl")

    include("controls/button.jl")
    include("controls/checkbox.jl")
    include("controls/crop.jl")
    include("controls/dualstatebutton.jl")
    include("controls/gauge.jl")
    include("controls/hotspot.jl")
    include("controls/number.jl")
    include("controls/page.jl")
    include("controls/picture.jl")
    include("controls/progressbar.jl")
    include("controls/qrcode.jl")
    include("controls/radio.jl")
    include("controls/scrolltext.jl")
    include("controls/slider.jl")
    include("controls/text.jl")
    include("controls/waveform.jl")

    include("invisibles/variable.jl")
    include("invisibles/timer.jl")
    include("invisibles/upload.jl")

    include("enhanced/eeprom.jl")
    include("enhanced/gpio.jl")
    include("enhanced/rtc.jl")

end # module


text:
"""
# SeisModels

Give various physical properties of the interior of the Earth for a range of
models.

### Units

Throughout, the module uses the following unit conventions:

- Distances are in km
- Velocities are in km/s
- Densities are in g/cm³
- Accelerations are in m/s²
- Pressures (from `pressure`) are in Pa
- Moduli (e.g., from `youngs_modulus`) are in GPa
- Masses in are in kg

### Physical properties

All models should have:

- P-wave velocity (accessed with the `vp` function) in km/s
- S-wave velocity (`vs`) in km/s
- (Average) Earth surface radius (`surface_radius`) in km

Models may also provide:

- Density (`density`) in g/cm^3
- Radial anisotropy:
  - VPV, VPH, VSV and VSH for radially-anistropic models (`vpv`, etc.) in km/s
  - Radial anisotropy parameter η (`eta`)
- Attenuation:
  - Shear wave quality factor (`Qμ`)
  - Bulk sound quality factor (`Qκ`)

(For anisotropic models, the values of VP and VS rather than VPV etc. should
be some average—perhaps the Voigt average—of the anisotropic velocities.)

Derived properties for all models:

- Radius from centre of Earth to a given depth (`depth`) in km
- Bulk and shear moduli (`bulk_modulus`, `shear_modulus`) in GPa
- Young's modulus (`youngs_modulus`) in GPa
- Poisson's ratio

Derived properties can be computed for all models with density:

- Mass up to radius r (`mass`) in kg
- Mass from radius r to the surface (`surface_mass`) in kg
- Pressure (`pressure`) in Pa
- Acceleration due to gravity g (`gravity`) in m/s^2

### Inbuilt models

- Earth
  - `PREM`
  - `IASP91`
  - `AK135`
- Moon
  - `MOON_WEBER_2011`
"""
module SeisModels

using DelimitedFiles: readdlm, writedlm
using Printf: @printf

using QuadGK: quadgk

export
    # Model types
    SeisModel,
    SeisModel1D,
    LinearLayeredModel,
    PREMPolyModel,
    SteppedLayeredModel,

    # Model properties
    depth,
    hasattenuation,
    hasdensity,
    hasreffrequency,
    isanisotropic,
    radius,
    reffrequency,
    surface_radius,

    # Evaluation functions
    evaluate,
    vp,
    vs,
    density,
    vph,
    vpv,
    vsh,
    vsv,
    eta,
    Qμ, Qmu,
    Qκ, Qkappa,

    # Derived properties
    bulk_modulus,
    gravity,
    mass,
    moment_of_inertia,
    poissons_ratio,
    pressure,
    shear_modulus,
    surface_mass,
    youngs_modulus,

    # IO
    read_mineos,
    write_mineos,
    read_tvel,
    write_tvel

@static if VERSION < v"1.1"
    isnothing(x) = x === nothing
end

"""
Abstract supertype of all models of the Earth in the `SeisModels` module.
All models should be a subtype of this or one of `SeisModel`'s abstract
subtypes.
"""
abstract type SeisModel end

# Allow models to be broadcasted as scalars
Base.broadcastable(x::SeisModel) = Ref(x)

# Default equality and approximate equality comparisons
Base.:(==)(m1::T, m2::T) where {T<:SeisModel} =
    all(isequal(getfield(m1, f), getfield(m2, f)) for f in fieldnames(T))
Base.:(≈)(m1::T, m2::T; kwargs...) where {T<:SeisModel} =
    all(≈(getfield(m1, f), getfield(m2, f); kwargs...) for f in fieldnames(T))

## 1D models
include("models_1d.jl")

## Basic properties of models
include("basic_properties.jl")

## Conversion between models
include("conversion.jl")

## IO
include("io.jl")

## Predefined models
# For each `planet`, models are defined in a module in "$(planet)/$(planet).jl".
# Each module must include a list of model symbols in ALL_MODELS.
for planet in (:Earth, :Moon)
    include(joinpath(@__DIR__, string(planet), string(planet, ".jl")))
    for model in getfield(SeisModels, planet).ALL_MODELS
        @eval begin
            using .$(planet): $(model)
            export $(model)
        end
    end
end

end # module


text:
using RvSpectMLBase
import EchelleInstruments
import EchelleInstruments: EXPRES
using EchelleCCFs
using Experimental
#using DataFrames, Query

all_spectra = include(joinpath(pkgdir(EchelleInstruments),"examples/read_expres_data_101501.jl"))

order_list_timeseries = extract_orders(all_spectra,pipeline_plan, recalc=true )

line_list_df = prepare_line_list_pass1(linelist_for_ccf_filename, all_spectra, pipeline_plan,  v_center_to_avoid_tellurics=ccf_mid_velocity, Δv_to_avoid_tellurics = 30e3, recalc=true)

(ccfs, v_grid) = ccf_total(order_list_timeseries, line_list_df, pipeline_plan,  mask_scale_factor=10.0, ccf_mid_velocity=ccf_mid_velocity, recalc=true)

line_width = RvSpectMLBase.calc_line_width(v_grid,view(ccfs,:,1),frac_depth=0.05)

line_list_df = prepare_line_list_pass1(linelist_for_ccf_filename, all_spectra, pipeline_plan,  v_center_to_avoid_tellurics=ccf_mid_velocity, Δv_to_avoid_tellurics = line_width)

(ccfs, v_grid) = ccf_total(order_list_timeseries, line_list_df, pipeline_plan,  mask_type=:gaussian, mask_scale_factor=10.0, ccf_mid_velocity=ccf_mid_velocity, recalc=true)

(order_ccfs, v_grid_order_ccfs) = ccf_orders(order_list_timeseries, line_list_df, pipeline_plan)

#rvs_ccf = calc_rvs_from_ccf_total(ccfs, pipeline_plan, v_grid=v_grid, times = order_list_timeseries.times, recalc=true)

#=
if need_to(pipeline_plan,:scalpels)
   rvs_scalpels = map(n->Scalpels.clean_rvs_scalpels(rvs_ccf, ccfs, num_basis=n), 1:5)
   println("RMS RVs cleaned by Scalpels: ",std.(rvs_scalpels) )
   dont_need_to!(pipeline_plan,:scalpels)
end
=#


text:
using LinearAlgebra: norm
using Printf: @printf

"""
    QuadES(T::Type{<:AbstractFloat}; maxlevel::Integer=10, h0::Real=one(T)/8)

A callable object to integrate a function over the range [0, ∞) using the
*exp-sinh quadrature*. It utilizes the change of variables to transform the
integrand into a form well-suited to the trapezoidal rule.

`QuadES` tries to calculate integral values `maxlevel` times at a maximum;
the step size of a trapezoid is started from `h0` and is halved in each
following repetition for finer accuracy. The repetition is terminated when the
difference from the previous estimation gets smaller than a certain threshold.
The threshold is determined by the runtime parameters, see below.

The type `T` represents the accuracy of interval. The integrand should accept
values `x<:T` as its parameter.

---

    I, E = (q::QuadES)(f::Function;
                       atol::Real=zero(T),
                       rtol::Real=atol>0 ? zero(T) : sqrt(eps(T)))
                       where {T<:AbstractFloat}

Numerically integrate `f(x)` over the interval [0, ∞) and return the integral
value `I` and an estimated error `E`. The `E` is not exactly equal to the
difference from the true value. However, one can expect that the integral value
`I` is converged if `E <= max(atol, rtol*norm(I))` is true. Otherwise, the
obtained `I` would be unreliable; the number of repetitions exceeds the
`maxlevel` before converged.

The integrand `f` can also return any value other than a scalar, as far as
`+`, `-`, multiplication by real values, and `LinearAlgebra.norm`, are
implemented. For example, `Vector` or `Array` of numbers are acceptable
although, unfortunately, it may not be very performant.

# Examples
```jldoctest
julia> using DoubleExponentialFormulas

julia> using LinearAlgebra: norm

julia> qes = QuadES(Float64);

julia> f(x) = 2/(1 + x^2);

julia> I, E = qes(f);

julia> I ≈ π
true

julia> E ≤ sqrt(eps(Float64))*norm(I)
true

julia> I, E = qes(x -> [1/(1 + x^2), 2/(1 + x^2)]);

julia> I ≈ [π/2, π]
true

julia> E ≤ sqrt(eps(Float64))*norm(I)
true
```
"""
struct QuadES{T<:AbstractFloat,N}
    h0::T
    origin::Tuple{T,T}
    table0⁺::Vector{Tuple{T,T}}
    table0⁻::Vector{Tuple{T,T}}
    tables⁺::NTuple{N,Vector{Tuple{T,T}}}
    tables⁻::NTuple{N,Vector{Tuple{T,T}}}
end
function QuadES(T::Type{<:AbstractFloat}; maxlevel::Integer=12, h0::Real=one(T))
    @assert maxlevel > 0
    @assert h0 > 0
    origin  = samplepoint(QuadES, zero(T))
    table0⁺ = generate_table⁺(QuadES, h0,  1)
    table0⁻ = generate_table⁻(QuadES, h0, -1)
    tables⁺ = Vector{Tuple{T,T}}[]
    tables⁻ = Vector{Tuple{T,T}}[]
    for level in 1:maxlevel
        h = h0/2^level
        push!(tables⁺, generate_table⁺(QuadES, h,  2))
        push!(tables⁻, generate_table⁻(QuadES, h, -2))
    end
    return QuadES{T,maxlevel}(T(h0), origin, table0⁺, table0⁻,
                              Tuple(tables⁺), Tuple(tables⁻))
end

function (q::QuadES{T,N})(f::Function; atol::Real=zero(T),
                          rtol::Real=atol>0 ? zero(T) : sqrt(eps(T))) where {T<:AbstractFloat,N}
    eval(t) = f(t[1])*t[2]
    x0, w0 = q.origin
    Σ = f(x0)*w0
    istart⁺ = startindex(f, q.table0⁺, 1)
    Σ += mapsum(eval, q.table0⁺, istart⁺)
    Σ += mapsum(eval, q.table0⁻)
    h0 = q.h0
    I = h0*Σ
    E = zero(eltype(I))
    for level in 1:N
        table⁺ = q.tables⁺[level]
        table⁻ = q.tables⁻[level]
        istart⁺ = startindex(f, table⁺, 2*istart⁺ - 1)
        Σ += mapsum(eval, table⁺, istart⁺)
        Σ += mapsum(eval, table⁻)
        h = h0/2^level
        prevI = I
        I = h*Σ
        E = estimate_error(T, prevI, I)
        tol = max(norm(I)*rtol, atol)
        !(E > tol) && break
    end
    return I, E
end

function Base.show(io::IO, ::MIME"text/plain", q::QuadES{T,N}) where {T<:AbstractFloat,N}
    @printf(io, "DoubleExponentialFormulas.QuadES{%s}: maxlevel=%d, h0=%.3e",
            string(T), N, q.h0)
end

function generate_table⁺(::Type{QuadES}, h::T, step::Int) where {T<:AbstractFloat}
    table = Tuple{T,T}[]
    k = 1
    while true
        t = k*h
        xk, wk = samplepoint(QuadES, t)
        xk ≥ floatmax(T) && break
        wk ≥ floatmax(T) && break
        push!(table, (xk, wk))
        k += step
    end
    reverse!(table)
    return table
end

function generate_table⁻(::Type{QuadES}, h::T, step::Int) where {T<:AbstractFloat}
    table = Tuple{T,T}[]
    k = -1
    while true
        t = k*h
        xk, wk = samplepoint(QuadES, t)
        # xk ≤ eps(zero(T)) may be better for accuracy but it may hurts
        # singular point durability for some kind of integrands.
        xk ≤ eps(T) && break
        wk ≤ floatmin(T) && break
        push!(table, (xk, wk))
        k += step
    end
    reverse!(table)
    return table
end

function samplepoint(::Type{QuadES}, t)
    expsinht = exp(sinh(t)*π/2)
    ϕ = expsinht
    ϕ′ = (cosh(t)*π/2)*expsinht
    return ϕ, ϕ′
end


text:
export macau

using CSV
using DelimitedFiles
using Distributed: workers
using JLD
using LinearAlgebra.BLAS
using Compat: eachrow
#### updated
using DataFrames

function macau(data::RelationData;
              num_latent::Int = 10,
              lambda_beta     = NaN,
              burnin          = 500,
              psamples        = 200,
              verbose::Bool   = true,
              full_lambda_u   = true,
              reset_model     = true,
              compute_ff_size = 6500,
              latent_pids     = workers(),
              latent_blas_threads = 1,
              cg_pids         = workers(),
              full_prediction = false,
              rmse_train      = false,
              tol             = NaN,
              output          = "",
              output_beta     = false,
              output_type     = "csv",
              clamp::Vector{Float64}  = Float64[],
              f               = false)
  correct = Float64[]
  local yhat_full

  if output_beta && isempty(output)
    error("To output samples of beta ('output_beta = true') you have to set also output prefix, e.g., output = \"my_model\".")
  end

  output_type in ["csv", "binary"] || error("output_type must be either \"csv\" or \"binary\".")

  verbose && println("Model setup")
  if reset_model
    reset!(data, num_latent, lambda_beta = lambda_beta, compute_ff_size = compute_ff_size, cg_pids=cg_pids)
  end

  if full_prediction
    yhat_full = zeros(size(data.relations[1]))
  end

  latent_multi_threading = false
  local latent_data_refs

  if length(latent_pids) >= 1
    ## initializing multi-threaded latent sampling:
    if length(data.relations) == 1 &&
       ! hasFeatures(data.relations[1])

       latent_multi_threading = true
       verbose && println("Setting up multi-threaded sampling of latent vectors. Using $(length(latent_pids)) threads.")
       fastidf = FastIDF(data.relations[1].data)
       latent_data_refs = map( i -> @spawnat( latent_pids[i], fetch(fastidf)), 1:length(latent_pids) )
       ## setting blas threads
       if latent_blas_threads >= 1
         for p in latent_pids
           if VERSION >= v"0.5.0"
             remotecall_wait(BLAS.set_num_threads, p, latent_blas_threads)
           else
             remotecall_wait(blas_set_num_threads, p, latent_blas_threads)
           end
         end
       end
    else
      verbose && println("Cannot use multi-threaded sampling of latent vectors, currently only works for 1 relation.")
    end
  end

  verbose && println("Sampling")
  err_avg  = 0.0
  roc_avg  = 0.0
  rmse_avg = 0.0
  counter_prob = 1

  local probe_rat_all, clamped_rat_all
  local probe_stdev::Vector{Float64}
  local train_rat_all
  f_output = Any[]

  ## Gibbs sampling loop
  for i in 1 : burnin + psamples
    time0 = time()

    # sample relation model (alpha)
    for r in data.relations
      if r.model.alpha_sample
        err = pred(r) - getValues(r.data)
        r.model.alpha = sample_alpha(r.model.alpha_lambda0, r.model.alpha_nu0, err)
      end
      if hasFeatures(r)
        r.model.beta = sample_beta_rel(r)
        r.temp.linear_values = r.model.mean_value .+ r.F * r.model.beta
      end
    end

    # Sampling latent vectors
    for j in 1:length(data.entities)
      en = data.entities[j]
      mj = en.model

      if latent_multi_threading
        sample_v = data.entities[j == 1 ? 2 : 1].model.sample
        if hasFeatures(en)
          mj.uhat = F_mul_beta(en)'
          mu_matrix = mj.mu .+ mj.uhat
          sample_latent_all2!(data.relations[1], latent_data_refs, latent_pids, en.modes[1], mu_matrix, mj.Lambda)
        else
          sample_latent_all2!(data.relations[1], latent_data_refs, latent_pids, en.modes[1], mj.mu, mj.Lambda)
        end
      else
        ## single thread
        if hasFeatures(en)
          mj.uhat = F_mul_beta(en)'
          mu_matrix = mj.mu .+ mj.uhat
          sample_user2_all!(data.entities[j], mu_matrix, en.modes, en.modes_other)
        else
          sample_user2_all!(data.entities[j], en.modes, en.modes_other)
        end
      end

      # update prior
      local U::Matrix{Float64}
      nu = mj.nu0
      Tinv = mj.WI
      if hasFeatures(en)
        U = mj.sample - mj.uhat
        if full_lambda_u
          nu   += size(mj.beta, 1)
          Tinv += mj.beta' * mj.beta * en.lambda_beta
        end
      else
        U = mj.sample
      end

      mj.mu, mj.Lambda = rand( ConditionalNormalWishart(U, mj.mu0, mj.b0, Tinv, nu) )

    end

    for en in data.entities
      update_beta!(en, tol)
    end

    rel = data.relations[1]
    probe_rat = pred(rel, rel.test_vec, rel.test_F)

    if full_prediction && i > burnin
      yhat_full += pred_all( data.relations[1] )
    end

    if i > burnin
      if length(output) > 0
        ## saving latent vectors to disk
        for en in data.entities
          ndigits = convert( Int, floor(log10(psamples)) ) + 1
          nstr    = lpad(string(i-burnin), ndigits, "0")
          fn = "$(output)-$(en.name)-$nstr"
          if output_type == "binary"
            save(string(fn, ".jld"), "sample", en.model.sample)
          else output_type == "csv"
            CSV.write(string(fn, ".csv"), DataFrame(en.model.sample), header=false)
          end
          if output_beta && hasFeatures(en)
            if output_type == "binary"
              save(string(fn, ".beta.jld"), "beta", en.model.beta)
            else output_type == "csv"
              CSV.write(string(fn, ".beta.csv"), DataFrame(en.model.beta), header=false)
            end
          end
        end
      end
      if rmse_train
        train_rat = pred(rel)
      end
      if i == burnin + 1
        verbose && println("--------- Burn-in complete, averaging posterior samples ----------")
        counter_prob  = 1
        probe_rat_all = probe_rat
        probe_stdev   = probe_rat .^ 2
        if rmse_train
          train_rat_all = train_rat
        end
      else
        probe_rat_all  = (counter_prob*probe_rat_all + probe_rat)/(counter_prob+1)
        probe_stdev   += probe_rat .^ 2
        if rmse_train
          train_rat_all = (counter_prob*train_rat_all + train_rat)/(counter_prob+1)
        end
        counter_prob  = counter_prob + 1
      end
    else
      probe_rat_all = probe_rat
    end

    if typeof(f) <: Function && i > burnin
      f_out = f( data )
      push!(f_output, f_out)
    end

    time1    = time()
    haveTest = numTest(rel) > 0

    correct  = (rel.test_label .== (probe_rat_all .< rel.class_cut) )
    err_avg  = mean(correct)
    err      = mean(rel.test_label .== (probe_rat .< rel.class_cut))

    clamped_rat     = isempty(clamp) ? probe_rat     : makeClamped(probe_rat, clamp)
    clamped_rat_all = isempty(clamp) ? probe_rat_all : makeClamped(probe_rat_all, clamp)

    rmse_avg = haveTest ? sqrt(mean( (rel.test_vec[:,end] - clamped_rat_all) .^ 2 )) : NaN
    rmse     = haveTest ? sqrt(mean( (rel.test_vec[:,end] - clamped_rat) .^ 2 ))     : NaN
    roc_avg  = haveTest ? AUC_ROC(rel.test_label, -vec(probe_rat_all))             : NaN
    if verbose
      estr = join(map(en -> toStr(en), data.entities), " ")
      rstr = join(map(r  -> toStr(r),  data.relations), " ")
      if i <= burnin
        verbose && @printf("%3d: ROC=%6.4f RMSE=%6.4f | %s | %s [%1.1fs]\n", i, roc_avg, rmse_avg, estr, rstr, time1 - time0)
      else
        verbose && @printf("%3d: ROC=%6.4f RMSE=%6.4f | %s | %s [%1.1fs]\n", i, roc_avg, rmse_avg, estr, rstr, time1 - time0)
      end
    end
  end

  result = Dict{AbstractString,Any}()
  result["num_latent"]  = num_latent
  result["burnin"]      = burnin
  result["psamples"]    = psamples
  result["lambda_beta"] = data.entities[1].lambda_beta
  result["RMSE"]        = rmse_avg
  result["accuracy"]    = err_avg
  result["ROC"]         = roc_avg
  if rmse_train
    ## calculating prediction on training set
    train_cl = isempty(clamp) ? train_rat_all : makeClamped(train_rat_all, clamp)
    result["RMSE_train"]  = sqrt(mean( (getValues(data.relations[1].data) - train_cl) .^ 2 ))
  end
  if full_prediction
    result["predictions_full"] = yhat_full / psamples
  end
  if numTest(data.relations[1]) > 0
    rel = data.relations[1]
    result["predictions"] = copy(rel.test_vec)
    result["predictions"][!, :pred]= vec(clamped_rat_all)
    if psamples >= 3
      tmp = vec(probe_stdev - probe_rat_all.^2 * psamples) / (psamples - 1)
      tmp[ tmp .< 0 ] .= 0
      result["predictions"][!, :stdev] = sqrt.(tmp)
    else
      result["predictions"][!, :stdev] = repeat([NaN], inner=[length(probe_rat_all)])
    end
    train_count = zeros(Int, numTest(rel), length(size(rel)) )
    for i in 1:numTest(rel)
      for mode in 1:length(size(rel))
        train_count[i,mode] = getCount(rel.data, mode, rel.test_vec[i,mode])
      end
    end
    # result["train_counts"] = convert( DataFrame, train_count )
    ### updated 3
    train_counts = DataFrame(train_counts_1 = train_count[:,1], train_counts_2 = train_count[:,2])
    result["train_counts"] = train_counts
  end
  if typeof(f) <: Function
    result["f_output"] = f_output
  end
  result["latent_multi_threading"] = latent_multi_threading
  return result
end


text:
# give nicer names to entities that had name clashes
const PARType = _Type
export PARType

# server impleneation intentionally empty as this is just the client part


text:
function intrinsic_value(contract::PSCallOptionContract, currentPriceValue::Float64)::PSResult

    # initialize -
    payoffValue = 0.0
    profitLossValue = 0.0

    # get data from the contract object -
    sense = contract.sense
    strikePrice = contract.strikePrice
    premiumValue = contract.premimumValue
    number_of_contracts = contract.numberOfContracts
    contract_multiplier = contract.contractMultiplier

    # compute the iv -
    if (sense == :buy)

        # compute the P/L -
        payoffValue = max(0.0, (currentPriceValue - strikePrice))
        profitLossValue = (payoffValue - premiumValue)

    elseif (sense == :sell)

        # compute the P/L -
        payoffValue = min(0.0,-1.0*(currentPriceValue - strikePrice))
        profitLossValue = (payoffValue + premiumValue)
    end

    # compute the intrinsic value -
    iv = (contract_multiplier)*(number_of_contracts)*payoffValue
    pl = (contract_multiplier)*(number_of_contracts)*profitLossValue

    # make a named tuple -
    named_tuple = (intrinsic_value=iv, pl_value=pl)

    # return -
    return PSResult(named_tuple)
end

function intrinsic_value(contract::PSPutOptionContract, currentPriceValue::Float64)::PSResult

    # initialize -
    payoffValue = 0.0
    profitLossValue = 0.0

    # get data from the contract object -
    sense = contract.sense
    strikePrice = contract.strikePrice
    premiumValue = contract.premimumValue
    number_of_contracts = contract.numberOfContracts
    contract_multiplier = contract.contractMultiplier

    # compute the iv -
    if (sense == :buy)

        # compute the P/L -
        payoffValue = max(0.0, (strikePrice - currentPriceValue))
        profitLossValue = (payoffValue - premiumValue)

    elseif (sense == :sell)

        # compute the P/L -
        payoffValue = min(0.0,-1.0*(strikePrice - currentPriceValue))
        profitLossValue = (payoffValue + premiumValue)
    end

    # compute the intrinsic value -
    pl = (contract_multiplier)*(number_of_contracts)*profitLossValue
    iv = (contract_multiplier)*(number_of_contracts)*payoffValue

    # make a named tuple -
    named_tuple = (intrinsic_value=iv, pl_value=pl)

    # return -
    return PSResult(named_tuple)
end

function intrinsic_value(equityObject::PSEquityAsset, currentPriceValue::Float64)::PSResult

    # initialize -
    iv = 0.0

    # get data from equityObject -
    purchasePricePerShare = equityObject.purchasePricePerShare
    numberOfShares = equityObject.numberOfShares

    # compute the intrinsic value -
    iv = numberOfShares*(currentPriceValue - purchasePricePerShare)

    # make a named tuple -
    named_tuple = (intrinsic_value=iv, pl_value=iv)

    # return -
    return PSResult(named_tuple)
end

function intrinsic_value(contractSet::Set{PSAbstractAsset}, underlyingPriceValue::Float64)::PSResult

    # initialize -
    tmp_iv_array = Float64[]

    # go through each contract, compute the iv, store -
    for (index,contract) in enumerate(contractSet)

        # compute -
        result = intrinsic_value(contract,underlyingPriceValue)
        if (isa(result.value, Exception) == true)
            return result
        end
        iv_value = result.value.intrinsic_value

        # store -
        push!(tmp_iv_array,iv_value)
    end

    # sum -
    total_iv_value = sum(tmp_iv_array)

    # return -
    return PSResult{Float64}(total_iv_value)
end

function compute_profit_loss_value(contractSet::Set{PSAbstractAsset}, underlyingPriceValue::Float64)::PSResult

    # initialize -
    tmp_pl_array = Float64[]

    # go through each contract, compute the iv, store -
    for (index,contract) in enumerate(contractSet)

        # compute -
        result = intrinsic_value(contract,underlyingPriceValue)
        if (isa(result.value, Exception) == true)
            return result
        end
        pl_value = result.value.pl_value

        # store -
        push!(tmp_pl_array,pl_value)
    end

    # sum -
    total_pl_value = sum(tmp_pl_array)

    # return -
    return PSResult{Float64}(total_pl_value)
end


text:
module LogExpFunctions

using DocStringExtensions: SIGNATURES
using Base: Math.@horner, @irrational

export loghalf, logtwo, logπ, log2π, log4π
export xlogx, xlogy, logistic, logit, log1psq, log1pexp, log1mexp, log2mexp, logexpm1,
    softplus, invsoftplus, log1pmx, logmxp1, logaddexp, logsubexp, logsumexp, softmax,
    softmax!

include("constants.jl")
include("basicfuns.jl")
include("logsumexp.jl")

end # module


text:
# Autogenerated wrapper script for openbabel_jll for i686-w64-mingw32-cxx03
export libinchi, libopenbabel, obabel, obconformer, obdistgen, obenergy, obfit, obfitall, obgen, obgrep, obminimize, obmm, obprobe, obprop, obrms, obrotamer, obrotate, obspectrophore, obsym, obtautomer, obthermo, roundtrip

using XML2_jll
using Zlib_jll
JLLWrappers.@generate_wrapper_header("openbabel")
JLLWrappers.@declare_library_product(libinchi, "libinchi.dll")
JLLWrappers.@declare_library_product(libopenbabel, "libopenbabel.dll")
JLLWrappers.@declare_executable_product(obabel)
JLLWrappers.@declare_executable_product(obconformer)
JLLWrappers.@declare_executable_product(obdistgen)
JLLWrappers.@declare_executable_product(obenergy)
JLLWrappers.@declare_executable_product(obfit)
JLLWrappers.@declare_executable_product(obfitall)
JLLWrappers.@declare_executable_product(obgen)
JLLWrappers.@declare_executable_product(obgrep)
JLLWrappers.@declare_executable_product(obminimize)
JLLWrappers.@declare_executable_product(obmm)
JLLWrappers.@declare_executable_product(obprobe)
JLLWrappers.@declare_executable_product(obprop)
JLLWrappers.@declare_executable_product(obrms)
JLLWrappers.@declare_executable_product(obrotamer)
JLLWrappers.@declare_executable_product(obrotate)
JLLWrappers.@declare_executable_product(obspectrophore)
JLLWrappers.@declare_executable_product(obsym)
JLLWrappers.@declare_executable_product(obtautomer)
JLLWrappers.@declare_executable_product(obthermo)
JLLWrappers.@declare_executable_product(roundtrip)
function __init__()
    JLLWrappers.@generate_init_header(XML2_jll, Zlib_jll)
    JLLWrappers.@init_library_product(
        libinchi,
        "bin\\libinchi.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libopenbabel,
        "bin\\libopenbabel.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_executable_product(
        obabel,
        "bin\\obabel.exe",
    )

    JLLWrappers.@init_executable_product(
        obconformer,
        "bin\\obconformer.exe",
    )

    JLLWrappers.@init_executable_product(
        obdistgen,
        "bin\\obdistgen.exe",
    )

    JLLWrappers.@init_executable_product(
        obenergy,
        "bin\\obenergy.exe",
    )

    JLLWrappers.@init_executable_product(
        obfit,
        "bin\\obfit.exe",
    )

    JLLWrappers.@init_executable_product(
        obfitall,
        "bin\\obfitall.exe",
    )

    JLLWrappers.@init_executable_product(
        obgen,
        "bin\\obgen.exe",
    )

    JLLWrappers.@init_executable_product(
        obgrep,
        "bin\\obgrep.exe",
    )

    JLLWrappers.@init_executable_product(
        obminimize,
        "bin\\obminimize.exe",
    )

    JLLWrappers.@init_executable_product(
        obmm,
        "bin\\obmm.exe",
    )

    JLLWrappers.@init_executable_product(
        obprobe,
        "bin\\obprobe.exe",
    )

    JLLWrappers.@init_executable_product(
        obprop,
        "bin\\obprop.exe",
    )

    JLLWrappers.@init_executable_product(
        obrms,
        "bin\\obrms.exe",
    )

    JLLWrappers.@init_executable_product(
        obrotamer,
        "bin\\obrotamer.exe",
    )

    JLLWrappers.@init_executable_product(
        obrotate,
        "bin\\obrotate.exe",
    )

    JLLWrappers.@init_executable_product(
        obspectrophore,
        "bin\\obspectrophore.exe",
    )

    JLLWrappers.@init_executable_product(
        obsym,
        "bin\\obsym.exe",
    )

    JLLWrappers.@init_executable_product(
        obtautomer,
        "bin\\obtautomer.exe",
    )

    JLLWrappers.@init_executable_product(
        obthermo,
        "bin\\obthermo.exe",
    )

    JLLWrappers.@init_executable_product(
        roundtrip,
        "bin\\roundtrip.exe",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
"""
    integral(f, t = echotime(f))

Integral of time profile `f` between `0` and `t`. Unless specified, the echotime
is used as the upper integral limit.

For the `PGSE`, `SinOGSE`, `CosOGSE` and `DoublePGSE` sequences, analytical
expressions are available. Otherwise a numerical integral is computed.
"""
function integral(f::TimeProfile, t = echotime(f))
    quadgk(f, zero(t), t)
end

function integral(f::PGSE, t = echotime(f))
    δ, Δ = f.δ, f.Δ
    (zero(t) ≤ t < δ) * t + (δ ≤ t < Δ + δ) * δ - (Δ ≤ t < Δ + δ) * (t - Δ)
end

function integral(f::DoublePGSE, t = echotime(f))
    δ, Δ, p = f.δ, f.Δ, f.tpause
    tmid = p + Δ + δ
    if zero(t) ≤ t < δ
        t
    elseif δ ≤ t < Δ
        δ
    elseif Δ ≤ t < Δ + δ
        δ - (t - Δ)
    elseif tmid ≤ t < tmid + δ
        t - tmid
    elseif tmid + δ ≤ t < tmid + Δ
        δ
    elseif tmid + Δ ≤ t < tmid + Δ + δ
        δ - (t - Δ - tmid)
    else
        zero(t)
    end
end

function integral(f::CosOGSE, t = echotime(f))
    δ, Δ, n = f.δ, f.Δ, f.nperiod
    ((t < δ) * sin(2π * n * t / δ) - (Δ ≤ t) * sin(2π * n * (t - Δ) / δ)) * δ / (2π * n)
end

function integral(f::SinOGSE, t = echotime(f))
    δ, Δ, n = f.δ, f.Δ, f.nperiod
    ((t < δ) * (1 - cos(2π * n * t / δ)) - (Δ ≤ t) * (1 - cos(2π * n * (t - Δ) / δ))) * δ /
    (2π * n)
end


text:
module CommunalHelperCassetteSwapBlock
using ..Ahorn, Maple

function swapFinalizer(entity)
    x, y = Ahorn.position(entity)

    width = Int(get(entity.data, "width", 8))
    height = Int(get(entity.data, "height", 8))

    entity.data["nodes"] = [(x + width, y)]
end

@mapdef Entity "CommunalHelper/CassetteSwapBlock" CassetteSwapBlock(x::Integer,
                                                                    y::Integer,
                                                                    width::Integer=Maple.defaultBlockWidth,
                                                                    height::Integer=Maple.defaultBlockHeight,
                                                                    index::Integer=0,
                                                                    tempo::Number=1.0,
                                                                    noReturn::Bool=false)

const colorNames = Dict{String, Int}(
    "Blue" => 0,
    "Rose" => 1,
    "Bright Sun" => 2,
    "Malachite" => 3
)

const colors = Dict{Int, Ahorn.colorTupleType}(
    1 => (240, 73, 190, 255) ./ 255,
	2 => (252, 220, 58, 255) ./ 255,
	3 => (56, 224, 78, 255) ./ 255
)

const ropeColors = Dict{Int, Ahorn.colorTupleType}(
    1 => (194, 116, 171, 255) ./ 255,
	2 => (227, 214, 148, 255) ./ 255,
	3 => (128, 224, 141, 255) ./ 255
)

const defaultRopeColor = (110, 189, 245, 255) ./ 255
const defaultColor = (73, 170, 240, 255) ./ 255

const placements = Ahorn.PlacementDict(
    "Cassette Swap Block ($index - $color) (Communal Helper)" => Ahorn.EntityPlacement(
        CassetteSwapBlock,
        "rectangle",
        Dict{String, Any}(
            "index" => index,
        ),
        swapFinalizer
    ) for (color, index) in colorNames
)

Ahorn.editingOptions(entity::CassetteSwapBlock) = Dict{String, Any}(
    "index" => colorNames
)

Ahorn.nodeLimits(entity::Maple.SwapBlock) = 1, 1

Ahorn.minimumSize(entity::CassetteSwapBlock) = 16, 16
Ahorn.resizable(entity::CassetteSwapBlock) = true, true

function Ahorn.selection(entity::CassetteSwapBlock)
    x, y = Ahorn.position(entity)
    stopX, stopY = Int.(entity.data["nodes"][1])

    width = Int(get(entity.data, "width", 8))
    height = Int(get(entity.data, "height", 8))

    return [Ahorn.Rectangle(x, y, width, height), Ahorn.Rectangle(stopX, stopY, width, height)]
end

const block = "objects/cassetteblock/solid"
const crossSprite = "objects/CommunalHelper/cassetteMoveBlock/x"

function renderTrail(ctx, x::Number, y::Number, width::Number, height::Number, trail::String, index::Int)
    tilesWidth = div(width, 8)
    tilesHeight = div(height, 8)

    ropeColor = get(ropeColors, index, defaultRopeColor)

    for i in 2:tilesWidth - 1
        Ahorn.drawImage(ctx, trail, x + (i - 1) * 8, y + 2, 6, 0, 8, 6, tint=ropeColor)
        Ahorn.drawImage(ctx, trail, x + (i - 1) * 8, y + height - 8, 6, 14, 8, 6, tint=ropeColor)
    end

    for i in 2:tilesHeight - 1
        Ahorn.drawImage(ctx, trail, x + 2, y + (i - 1) * 8, 0, 6, 6, 8, tint=ropeColor)
        Ahorn.drawImage(ctx, trail, x + width - 8, y + (i - 1) * 8, 14, 6, 6, 8, tint=ropeColor)
    end

    for i in 2:tilesWidth - 1, j in 2:tilesHeight - 1
        Ahorn.drawImage(ctx, trail, x + (i - 1) * 8, y + (j - 1) * 8, 6, 6, 8, 8, tint=ropeColor)
    end

    Ahorn.drawImage(ctx, trail, x + width - 8, y + 2, 14, 0, 6, 6, tint=ropeColor)
    Ahorn.drawImage(ctx, trail, x + width - 8, y + height - 8, 14, 14, 6, 6, tint=ropeColor)
    Ahorn.drawImage(ctx, trail, x + 2, y + 2, 0, 0, 6, 6, tint=ropeColor)
    Ahorn.drawImage(ctx, trail, x + 2, y + height - 8, 0, 14, 6, 6, tint=ropeColor)
end

function Ahorn.renderSelectedAbs(ctx::Ahorn.Cairo.CairoContext, entity::CassetteSwapBlock, room::Maple.Room)
    startX, startY = Int(entity.data["x"]), Int(entity.data["y"])
    stopX, stopY = Int.(entity.data["nodes"][1])

    width = Int(get(entity.data, "width", 32))
    height = Int(get(entity.data, "height", 32))

    tileWidth = ceil(Int, width / 8)
    tileHeight = ceil(Int, height / 8)

    index = Int(get(entity.data, "index", 0))
    color = get(colors, index, defaultColor)

    for i in 1:tileWidth, j in 1:tileHeight
        tx = (i == 1) ? 0 : ((i == tileWidth) ? 16 : 8)
        ty = (j == 1) ? 0 : ((j == tileHeight) ? 16 : 8)

        Ahorn.drawImage(ctx, block, stopX + (i - 1) * 8, stopY + (j - 1) * 8, tx, ty, 8, 8, tint=color)
    end

    Ahorn.drawArrow(ctx, startX + width / 2, startY + height / 2, stopX + width / 2, stopY + height / 2, Ahorn.colors.selection_selected_fc, headLength=6)
end

function Ahorn.renderAbs(ctx::Ahorn.Cairo.CairoContext, entity::CassetteSwapBlock, room::Maple.Room)
    x, y = Ahorn.position(entity)
    stopX, stopY = Int.(entity.data["nodes"][1])

    width = Int(get(entity.data, "width", 32))
    height = Int(get(entity.data, "height", 32))

    tileWidth = ceil(Int, width / 8)
    tileHeight = ceil(Int, height / 8)

    index = Int(get(entity.data, "index", 0))
    color = get(colors, index, defaultColor)

    renderTrail(ctx, min(x, stopX), min(y, stopY), abs(x - stopX) + width, abs(y - stopY) + height, "objects/swapblock/target", index)

    for i in 1:tileWidth, j in 1:tileHeight
        tx = (i == 1) ? 0 : ((i == tileWidth) ? 16 : 8)
        ty = (j == 1) ? 0 : ((j == tileHeight) ? 16 : 8)

        Ahorn.drawImage(ctx, block, x + (i - 1) * 8, y + (j - 1) * 8, tx, ty, 8, 8, tint=color)
    end

    if Bool(get(entity.data, "noReturn", false))
        noReturnSprite = Ahorn.getSprite(crossSprite, "Gameplay")
        Ahorn.drawImage(ctx, noReturnSprite, x + div(width - noReturnSprite.width, 2), y + div(height - noReturnSprite.height, 2), tint=color)
    end
end

end

max_stars_repo_name: FlynnSC/CommunalHelper

max_issues_repo_name: FlynnSC/CommunalHelper

max_forks_repo_name: FlynnSC/CommunalHelper


text:
"""
    GeneralizedGaussian( μ, α, β )

The "Generalized Normal" or "Generalized Gaussian" with shape parameter 'β',
scale 'α', and location 'μ' has probability density function

```math
f(x; \\mu \\alpha \\beta ) = \\frac{\\beta}{2\\alpha\\Gamma(1/\\beta)}
\\text{exp}\\left( -|x - \\mu|/\\alpha \\right),

\\alpha = \\sigma \\sqrt{ \\frac{\\Gamma(1/\\beta)}{\\Gamma(3/\\beta)} }

```
where `β = 1` incorporates the Laplacian distribution, `β = 2` is the Normal distribution,
and as `β → ∞`, the distribution approaches Uniform on `[μ-α, μ+α]`.

```julia
GeneralizedGaussian(m, a, b)    # Generalized Gaussian with mu, alpha, and beta.
params(d)                       # Get the parameters (m,a,b)
shape(d)                        # returns β shape
scale(d)                        # returns the scale parameter, α
```

External Links
* [Generalized Normal Distribution on Wikipedia](https://en.wikipedia.org/wiki/Generalized_normal_distribution)
* [Reference implementation paper](https://www.researchgate.net/publication/254282790_Simulation_of_the_p-generalized_Gaussian_distribution)

References
[1] Nadarajah, Saralees (September 2005). "A generalized normal distribution".
    Journal of Applied Statistics. 32 (7): 685&ndash, 694.
[2]  Gonzalez-Farias, G., Molina, J. A. D., & Rodríguez-Dagnino, R. M. (2009).
    Efficiency of the approximated shape parameter estimator in the generalized
    Gaussian distribution. IEEE Transactions on Vehicular Technology, 58(8),
    4214-4223.
"""
struct GeneralizedGaussian{T<:Real}
    μ::T
    α::T
    β::T
    function GeneralizedGaussian{T}(μ::T, α::T, β::T) where T
        α > zero(α) || error("standard deviation must be positive.")
        β > zero(β) || error("β must be greater than zero.")
        new{T}(μ, α, β)
    end
end

GeneralizedGaussian(μ::T, α::T, β::T) where {T<:Real} = GeneralizedGaussian{T}(μ, α, β)
GeneralizedGaussian(μ::Real, α::Real, β::Real) = GeneralizedGaussian(promote(μ, α, β)...)
GeneralizedGaussian(μ::Integer, σ::Integer, β::Integer) = GeneralizedGaussian(Float64(μ), Float64(σ), Float64(β))

"""
    GeneralizedGaussian(β)

Builds a default distribution with shape `β` and `μ=0, α=1`.
"""
GeneralizedGaussian(β::Real) = GeneralizedGaussian(zero(β), one(β), β)

"""
    GeneralizedGaussian()

Builds the Normal distribution case, where `μ=0, α=√2, β=2`, or N(0, 1).
"""
GeneralizedGaussian() = GeneralizedGaussian(0.0, √2, 2)

### Conversion
function convert(::Type{GeneralizedGaussian{T}}, μ::Real, σ::Real, β::Real) where T <: Real
    GeneralizedGaussian(T(μ), T(α), T(β))
end
function convert(::Type{GeneralizedGaussian{T}}, d::GeneralizedGaussian{S}) where {T <: Real, S <: Real}
    GeneralizedGaussian(T(d.μ), T(d.α), T(d.β))
end

### Parameters
params(d::GeneralizedGaussian) = (d.μ, d.α, d.β)
location(d::GeneralizedGaussian) = d.μ
scale(d::GeneralizedGaussian) = d.α
shape(d::GeneralizedGaussian) = d.β

var(d::GeneralizedGaussian) = (d.α^2) * (gamma(3.0 * inv(d.β)) / gamma(inv(d.β)))
std(d::GeneralizedGaussian) = (d.α) * sqrt(gamma(3.0 * inv(d.β)) / gamma(inv(d.β)))

mean(d::GeneralizedGaussian) = d.μ
median(d::GeneralizedGaussian) = d.μ
mode(d::GeneralizedGaussian) = d.μ

### Statistics
skewness(d::GeneralizedGaussian{T}) where T = zero(T)
kurtosis(d::GeneralizedGaussian) = gamma(5.0 * inv(d.β)) * gamma(inv(d.β)) / (gamma(3.0 * inv(d.β))^2) - 3.0
entropy(d::GeneralizedGaussian) = inv(d.β) - log( d.β / (2.0 * d.α * gamma(inv(d.β))))

### Evaluation
function pdf(d::GeneralizedGaussian{T}, x::Real) where T<:Real
    if x == -Inf || x == Inf
        return zero(T)
    else
        A = d.β / (2.0 * d.α * gamma(1. / d.β) )
        return A * exp( -(abs(x - d.μ) / d.α)^d.β )
    end
end
logpdf(d::GeneralizedGaussian, x::Real) = log(pdf(d, x))

"""
    cdf(d, x)
Calculates the CDF of the distribution. To determine the CDF, the incomplete
gamma function is required. The CDF  of the Gamma distribution provides this,
with the necessary 1/Γ(a) normalization.
"""
function cdf(d::GeneralizedGaussian{T}, x::Real) where T<:Real
    if x == -Inf || x == Inf
        return zero(T)
    else
        v = cdf(Gamma(inv(d.β), 1), (abs(x - μ) / d.α)^d.β) * inv(2)
        return typeof(v)(1/2) + sign(x - μ) * v
    end
end

"""
    rand(rng, d)
Extract a sample from the Generalized Gaussian distribution `d`. The sampling
procedure is implemented from from [2].
"""
function rand(rng::AbstractRNG, d::GeneralizedGaussian)
    # utilizing the sampler from the Gamma distribution.
    g = Gamma(inv(d.β), 1)
    # random variable with value -1 or 1 with probability (1/2).
    b = 2.0 * rand(Bernoulli()) - 1
    return d.μ + inv(sqrt(d.α)) * rand(rng, g)^inv(d.β) * b
end

# multi-sample case
function rand(rng::AbstractRNG, d::GeneralizedGaussian, dims::Dims)
    out = Array{eltype(params(d))}(undef, dims)
    @inbounds for i in eachindex(out)
        out[i] = rand(rng, d)
    end
    out
end

# various function signatures for sampling the distribution.
rand(d::GeneralizedGaussian) = rand(GLOBAL_RNG, d)
rand(d::GeneralizedGaussian, dims::Dims) = rand(GLOBAL_RNG, d, dims)
rand(d::GeneralizedGaussian, dims::Int) = rand(GLOBAL_RNG, d, (dims))
rand(rng::AbstractRNG, d::GeneralizedGaussian, dim::Int, dims::Int...) = rand(rng, d, (dim, dims...))
rand(d::GeneralizedGaussian, dims::Int...) = rand(GLOBAL_RNG, d, dims...)

function Zn(X::AbstractArray, β::Real)
    n = size(X,1)
    S1 = 1/n * sum(abs.(X).^(2β))
    S2 = 1/n * sum(abs.(X).^β)
    return (S1/(S2^2)) - (β + 1)
end

function Zn′(X::AbstractArray, β::Real)
    n = size(X,1)
    S1 = sum(abs.(X).^β) # sum
    S2 = sum(abs.(X).^(2β))
    L1 = sum(abs.(X).^β .* log.(abs.(X)))    # log
    L2 = sum(abs.(X).^(2β) .* log.(abs.(X)))
    num1 = (2/n * L2) * (1/n * S1)^2
    denom = (1/n * S1)^4
    num2 = (1/n * L1) * (1/n * S2) * (2/n * S1)
    return (num1/denom) - (num2/denom) - 1
end

"""
    β = gcmsearch(X::AbstractArray, βi::Real)

Globally convergent method for β-parameter estimation using Newton-Raphson
iterative search from [1].

[1] Song, Kai-Sheng. "A globally convergent and consistent method
    for estimating the shape parameter of a generalized Gaussian distribution."
    IEEE Transactions on Information Theory 52.2 (2006): 510-527.
"""
function gcmsearch(X::AbstractArray{T}, βi::Real) where T<:AbstractFloat
    ϵ = √eps(T)
    N_MAX = 1_000
    βp = copy(βi)
    num, denom = Zn(X, βp), Zn′(X, βp)
    βn = βp - (num/denom)
    βp = copy(βn)
    for _ = 2:N_MAX
        num, denom = Zn(X, βp), Zn′(X, βp)
        βn = βp - (num/denom)
        if (abs(βn - βp) < ϵ) # check if converged?
            break
        else
            βp = copy(βn) # prep for next iteration.
        end
    end
    βn
end


text:
using CSV, DataFrames
using ScikitLearn
using Random, Statistics
using JLD2

# Import the necessary modules and libraries from ScikitLearn
@sk_import ensemble:RandomForestClassifier

# Import our training function which trains a random forest
include("rf.jl")

# Load the data for training and validation
train = CSV.read("data/listings_clean_train.csv", header=false)
val = CSV.read("data/listings_clean_val.csv", header=false)

# Load the parameters dataframe
@load "params.jld2" params

# Get the experiment ID from the arguments that were passed to julia
experiment_id = 1
if length(ARGS) > 0
	experiment_id = parse(Int, ARGS[1])
end

# Set a random seed - always a good idea when using functions you didn't write
Random.seed!(1776)
for i = 1:nrow(params)
	# skip all experiment ids except for the one we care about
	if params[i, :experiment_id] != experiment_id
		continue
	end
	R2tr, R2va = trainrandomforest(train, val, params[i, :max_depth],
	                               params[i, :n_estimators], params[i, :max_features])
	@save "results/results-$experiment_id.jld2" R2tr R2va
end


text:
using Flux, Mill, BenchmarkTools, Random

X1 = Flux.onehotbatch([1,2,3,4,5], 1:10)
X2 = Mill.maybehotbatch([1,2,3,4,5], 1:10)
X3 = Mill.maybehotbatch(fill(missing, 5), 1:10)
X4 = Mill.maybehotbatch([1,missing,3,missing,5], 1:10)

copyinto(xs...; n=20÷length(xs)) = vcat([collect(deepcopy.(xs)) for _ in 1:n]...)

# should be equally fast
@btime hcat(X1, X1);
@btime hcat(X2, X2);
@btime hcat(X3, X3);
@btime hcat(X4, X4);

# slow, but possible
@btime hcat(X2, X4);
@btime hcat(X3, X4);

# should be equally fast
@btime reduce(hcat, $(copyinto(X1)));
@btime reduce(hcat, $(copyinto(X2)));
@btime reduce(hcat, $(copyinto(X3)));
@btime reduce(hcat, $(copyinto(X4)));

# slow, but possible
@btime reduce(hcat, $(copyinto(X2, X4)));
@btime reduce(hcat, $(copyinto(X3, X4)));

X1 = NGramMatrix([randstring(10) for _ in 1:10])
X2 = NGramMatrix(fill(missing, 10))
X3 = NGramMatrix([isodd(i) ? missing : randstring(10) for i in 1:10])

# should be equally fast
@btime hcat(X1, X1);
@btime hcat(X2, X2);
@btime hcat(X3, X3);

# slow, but possible
@btime hcat(X2, X3);

# should be equally fast
@btime reduce(hcat, $(copyinto(X1)));
@btime reduce(hcat, $(copyinto(X2)));
@btime reduce(hcat, $(copyinto(X3)));

# slow, but possible
@btime reduce(hcat, $(copyinto(X2, X3)));


text:
"""
    struct Output

The composite type for output.

$(TYPEDFIELDS)
"""
struct Output
    "the relevant problem for the output"
      prob :: Problem
    "the path for the output file"
      path :: String
    "the fields to be saved; the relevant problem for this diagnostic"
    fields :: Dict{Symbol, Function}

  @doc """
      Output(prob, filename, fieldtuples...)

  Define output for `prob` with fields and functions that calculate
  the output in the list of tuples `fieldtuples = (fldname, func)...`.
  """
  function Output(prob, path, fieldtuples::Dict{Symbol,Function})
    truepath = uniquepath(withoutjld2(path) * ".jld2") # ensure path ends in ".jld2"
    return new(prob, truepath, fieldtuples)
  end
end

Output(prob, path, fields...) = Output(prob, path, Dict{Symbol,Function}([fields...]))
Output(prob, path, field::Tuple{Symbol,T}) where T = Output(prob, path, Dict{Symbol,Function}([field]))

withoutjld2(path) = (length(path)>4 && path[end-4:end] == ".jld2") ? path[1:end-5] : path

"""
    uniquepath(path)

Return `path` with a number appended if `isfile(path) == true`. The number is incremented
until `path` does not exist.
"""
function uniquepath(path)
  n = 1

  if isfile(path)
    path = withoutjld2(path) * "_$n.jld2"
  end

  while isfile(path)
    n += 1
    path = withoutjld2(path)[1:end-length("_$(n-1)")] * "_$n.jld2"
  end

  return path
end

getindex(out::Output, key) = out.fields[key](out.prob)

"""
    saveoutput(out)

Save the fields in `out.fields` to `out.path`.
"""
function saveoutput(out)
  groupname = "snapshots"
  jldopen(out.path, "a+") do path
    path["$groupname/t/$(out.prob.clock.step)"] = out.prob.clock.t
    for fieldname in keys(out.fields)
      path["$groupname/$fieldname/$(out.prob.clock.step)"] = out[fieldname]
    end
  end

  return nothing
end

"""
    savefield(file, location, data)

Saves a particular field's `data` to `file`.
"""
savefield(file, location, data) = file[location] = data
savefield(file, location, data::AbstractArray) = file[location] = Array(data)

"""
    savefields(file, grid)

Saves some parameters of problem's `grid` to `file`.
"""
function savefields(file::JLD2.JLDFile{JLD2.MmapIO}, grid::TwoDGrid)
  for field in [:nx, :ny, :Lx, :Ly, :x, :y]
    savefield(file, "grid/$field", getfield(grid, field))
  end

  return nothing
end

function savefields(file::JLD2.JLDFile{JLD2.MmapIO}, grid::OneDGrid)
  for field in [:nx, :Lx, :x]
    savefield(file, "grid/$field", getfield(grid, field))
  end

  return nothing
end

forbidden_output_types = [Function, FFTW.FFTWPlan, CUDA.CUFFT.CuFFTPlan]

function savefields(file::JLD2.JLDFile{JLD2.MmapIO}, params::AbstractParams)
  for name in fieldnames(typeof(params))
    field = getfield(params, name)
    if !(any([typeof(field) <: forbidden for forbidden in forbidden_output_types]))
      savefield(file, "params/$name", field)
    end
  end

  return nothing
end

function savefields(file::JLD2.JLDFile{JLD2.MmapIO}, clock::Clock)
  file["clock/dt"] = clock.dt

  return nothing
end

function savefields(file::JLD2.JLDFile{JLD2.MmapIO}, eqn::Equation)
  savefield(file, "eqn/L", eqn.L)
  file["eqn/dims"] = eqn.dims
  file["eqn/T"] = eqn.T

  return nothing
end

"""
    saveproblem(prob, filename)

Save certain aspects of a problem `prob` to `filename`.
"""
function saveproblem(prob, filename)
  file = jldopen(filename, "a+")

  for field in [:eqn, :clock, :grid, :params]
    savefields(file, getfield(prob, field))
  end
  close(file)

  return nothing
end

saveproblem(out::Output) = saveproblem(out.prob, out.path)

"""
    savediagnostic(diag, diagname, filename)

Save diagnostics in `diag` to `filename`, labeled by `diagname`.
"""
function savediagnostic(diag, diagname, filename)
  jldopen(filename, "a+") do file
    file["diags/$diagname/steps"] = diag.steps
    file["diags/$diagname/t"] = diag.t
    file["diags/$diagname/data"] = diag.data
  end

  return nothing
end

show(io::IO, out::Output) =
     print(io, "Output\n",
               "  ├──── prob: ", summary(out.prob), '\n',
               "  ├──── path: ", out.path, '\n',
               "  └── fields: ", out.fields)


text:
module KnapsackPlates

module Args
	import ...Utilities
	import ...Utilities.Args: Arg
	#export accepted_arg_list, throw_if_incompatible_options
	# No extra flags for this model yet.
	Utilities.Args.accepted_arg_list(::Val{:KnapsackPlates}) = Vector{Arg}()
	function Utilities.Args.throw_if_incompatible_options(
		::Val{:KnapsackPlates}, p_args
	)
		nothing
	end
end

using JuMP

import ..build_model
"""
    build_model(::Val{:KnapsackPlates}, model, d, p, l, w, L, W, [options])

The KnapsackPlates models each plate as a 1D knapsack of pieces and smaller
subplates. Basically, a plate is allowed to be cut in just one orientation
(except the original plate, that has a hack to allow the first stage to
be any of the two orientations) and, this way, the plate is treated as a
1D knapsack problem using the size of pieces and plates in the considered
dimension.

NOTE: THIS MODEL IS BROKEN. The source has the skeleton of the implementation
but it fails compilation as it is.

```
TODO: add high-level description of the mathematical model here.
```

"""
function build_model(
	::Val{:KnapsackPlates}, model, d :: Vector{D}, p :: Vector{P},
	l :: Vector{S}, w :: Vector{S}, L :: S, W :: S,
	options :: Dict{String, Any} = Dict{String, Any}()
) where {D, S, P}
	@assert length(d) == length(l) && length(l) == length(w)
	num_fplate_types = convert(D, length(d))

	@variable(model, 0 <= fplate[i = 1:num_fplate_types] <= d[i], Int)
	@variable(model, 0 <= iplate[i = 1:num_iplate_types] <= qt_fit_op[i], Int)

	@objective(model, Max,
		sum(p[fpi] * fplate[fpi] for fpi = 1:num_fplate_types)
	)

	# The original plate issue (no way to cut allow cutting both horizonatlly and
	# vertically) may be solved by having a default cut orientation for the
	# original plate, and allowing to cut a same-size-but-opposite-orientation
	# plate. Such plate has length L and is contained in the list of plates to
	# be cut vertically.
	@constraint(model,
		sum(lfp .* fp) + sum(lip2cv .* ip2cv) <= L
	)
	# iip2ch: Index of Intermediary Plates To Cut Horizontally
	for iip2ch = 1:num_ip2ch_types
		fitting_fp = fp_fitting_iipch[iip2ch]
		fitting_ip2cv = iip2cv_fitting_iip2ch[iip2ch]
		isempty(fitting_fp) || isempty(fitting_ip2cv) && continue
		@constraint(model,
			sum(lfp[fitting_fp] .* fp[fitting_fp]) +
			sum(lip2cv[fitting_ip2cv] .* ip2cv[fitting_ip2cv]) <=
			lip2ch[iip2ch] * ip2ch[iip2ch]
		)
	end

	for iip2cv = 1:num_ip2cv_types
		fitting_fp = fp_fitting_iipcv[iip2cv]
		fitting_ip2ch = iip2ch_fitting_iip2cv[iip2cv]
		isempty(fitting_fp) || isempty(fitting_ip2ch) && continue
		@constraint(model,
			sum(wfp[fitting_fp] .* fp[fitting_fp]) +
			sum(wip2ch[fitting_ip2ch] .* ip2ch[fitting_ip2ch]) <=
			wip2cv[iip2cv] * ip2cv[iip2cv]
		)
	end

	model
end

end # module


text:
export Study, Experiment, Reconstruction, Visualization, DatasetStore,
       studydir, BrukerDatasetStore, BrukerStore, getStudy, getStudies, getExperiment,
       getExperiments, MDFDatasetStore, MDFStore, addReco, getReco, getRecons, findReco,
       findBrukerFiles, id, getVisus, getVisuPath, remove, addStudy, getNewExperimentNum,
       exportToMDFStore, generateSFDatabase, loadSFDatabase, addVisu, readonly, getNewCalibNum,
       calibdir, try_chmod, getMDFStudyFolderName

########################################

# The following are base types describing
# the dataset store at a certain level
struct Study
  path::String
  name::String
  subject::String
  date::DateTime
end

id(s::Study) = s.name

# might not be so clever to use explicity type fields here
# maybe better a dict
struct Experiment
  path::String
  num::Int64
  name::String
  numFrames::Int64
  df::Vector{Float64}
  sfGradient::Float64
  numAverages::Int64
  operator::String
  time::String
  # more ...
end

mutable struct Reconstruction
  path::String
  num::Int64
  params::Dict
end

mutable struct Visualization
  path::String
  num::Int64
  params::Dict
end

abstract type DatasetStore end

struct BrukerDatasetStore <: DatasetStore
   path::String
end

const BrukerStore = BrukerDatasetStore("/opt/mpidata")

function try_chmod(path, mode; recursive=true)
  try
    chmod(path,mode,recursive=recursive)
  catch
  end
  return
end

struct MDFDatasetStore <: DatasetStore
  path::String

  function MDFDatasetStore(path::String)
    if ispath(path)
      mkpath(joinpath(path,"measurements"))
      try_chmod(joinpath(path,"measurements"), 0o777, recursive=true)
      mkpath(joinpath(path,"reconstructions"))
      try_chmod(joinpath(path,"reconstructions"), 0o777, recursive=true)
      mkpath(joinpath(path,"calibrations"))
      try_chmod(joinpath(path,"calibrations"), 0o777, recursive=true)
    end
    return new(path)
  end
end

const MDFStore = MDFDatasetStore("/opt/data/Bruker")

### generic functions ###
function ishidden(filename::AbstractString)
  @static if Sys.isunix()
    s = basename(filename)
    return (!isempty(s) && s[1] == '.')
  else
    attr = ccall((:GetFileAttributesA), stdcall, Cint, (Ptr{UInt8},),bytestring(filename))
    return attr & 0x2 > 0
  end
end

function getStudies(d::DatasetStore)
  s = Study[]

  files = readdir( studydir(d) )
  for file in files
    fullpath = joinpath(studydir(d),file)
    if isdir(fullpath) && !ishidden(fullpath)
      try
      study = getStudy(d, file)
      if study != nothing
        push!(s, study)
      end
      catch e
        @warn e
      end
    end
  end
  return s
end

function remove(study::Study)
  if isdir(study.path)
    rm(study.path, recursive=true)

    #TODO remove recos!
  end
end

function getExperiment(path::String)

  prefix, ext = splitext(path)
  if isdir(path) #Ugly
    p = path
    b = MPIFiles.BrukerFileFast(path) #use fast path for BrukerFiles
  else
    p = string(prefix,".mdf")
    if isfile(p)
      b = MDFFile(p)
    else
      return nothing
    end
  end

  exp = Experiment(p, parse(Int64,last(splitdir(prefix))),
                      string(experimentName(b)), acqNumFrames(b),
                      round.(1000 .* vec(dfStrength(b)[1,:,1]),digits=2), maximum(abs.(acqGradient(b))),
                      acqNumAverages(b), scannerOperator(b), string(acqStartTime(b)))

  return exp
end

function getExperiment(s::Study, numExp::Integer)
   getExperiment(joinpath(s.path,string(numExp)))
end

function remove(exp::Experiment)
  if isfile(exp.path)
    rm(exp.path)
  end
end

function exportToMDFStore(path::String, mdf::MDFDatasetStore)
  # pretend to be a measurement to enforce loading data from time domain in case post processed data is not availible
  b = BrukerFile(path,isCalib=false)
  exportpath = ""

  if _iscalib(path)
    exportpath = getNewCalibPath(mdf)
    saveasMDF(exportpath,b,applyCalibPostprocessing=true)
    @info "Calibration data from $path sucessfully exported to $exportpath."
  else
    name = studyName(b)
    subject = experimentSubject(b)
    date = studyTime(b)
    mdfPath = joinpath(studydir(mdf),name)
    s = Study(mdfPath,name,subject,date)
    exportpath = getNewExperimentPath(mdf,s)
    saveasMDF(exportpath, b)
    @info "Measurment data from $path sucessfully exported to $exportpath."
  end
  # Store export path in Bruker directory
  open(joinpath(path,"mdf"),write=true) do io
    write(io, exportpath)
  end
  # log action
  open("/opt/DataArchiveOptmpidata/convertBrukerToMDF/log.csv", append=true) do io
    write(io,"$path, $exportpath\n")
  end

  return exportpath
end

###  Implementations of abstract interfaces ###

readonly(::BrukerDatasetStore) = true
readonly(::MDFDatasetStore) = false

studydir(d::BrukerDatasetStore) = d.path
studydir(d::MDFDatasetStore) = joinpath(d.path,"measurements")
calibdir(d::MDFDatasetStore) = joinpath(d.path,"calibrations")

function getStudy(d::BrukerDatasetStore, studyfolder::String)
  study = nothing
  studypath = joinpath(d.path,studyfolder)
  if !ishidden(studypath) && isdir(studypath)
    w = split(studyfolder,'_')
    if length(w) >= 5 && length(w[1])==8 # only these can be study folders
      # w_ = w[1:end-2]
      # date = w[1]
      # date = string(date[1:4],"/",date[5:6],"/",date[7:8])

      w = split(studyfolder,'_')
      dateStr = w[1]
      timeStr = w[2]
      date = DateTime(string(dateStr[1:4],"-",dateStr[5:6],"-",dateStr[7:8],"T",
			   timeStr[1:2],":",timeStr[3:4],":",timeStr[5:6]))

      j = JcampdxFile()
      subjfile = string(studypath,"/subject")
      if isfile(subjfile)
        read(j,string(studypath,"/subject"),maxEntries=14) #magic number...
        name = latin1toutf8(j["SUBJECT_study_name"])
        # name = string(latin1toutf8(j["SUBJECT_name_string"]),
        #              "_",latin1toutf8(j["SUBJECT_study_name"]),
        #              "_",latin1toutf8(j["SUBJECT_study_nr"]))
        s1 = latin1toutf8(j["SUBJECT_id"])
	s2 = latin1toutf8(j["SUBJECT_name_string"])
        subject = (s1 == s2) ? s1 : s1*s2
      else
        # Workaround if no subject file is present => use first dataset
        # and derive the study from the Brukerfile
        r = readdir(studypath)

        found = false
        for file in r

          if tryparse(Int64,file) != nothing
            b = BrukerFileFast(joinpath(studypath, file ))
            name = studyName(b)
            subject = experimentSubject(b)
            found = true
            break
          end
        end
        if !found
          return nothing
        end
      end
      study = Study(studypath, name, subject, date )
    end
  end
  return study
end

function getStudy(d::MDFDatasetStore, studyfolder::String)
  study = nothing
  studypath = joinpath( studydir(d), studyfolder)
  if length(studyfolder) >= 15 &&
     isascii(studyfolder[1:15]) &&
     all([tryparse(Int,studyfolder[l:l])!=nothing for l=union(1:8,10:15)])

    w = split(studyfolder,'_')
    dateStr = w[1]
    timeStr = w[2]
    date = DateTime(string(dateStr[1:4],"-",dateStr[5:6],"-",dateStr[7:8],"T",
			   timeStr[1:2],":",timeStr[3:4],":",timeStr[5:6]))
    name = join(w[3:end],"_")
  else
    date = Dates.unix2datetime(stat(studypath).mtime)
    name = studyfolder
  end

  subject = ""
  study = Study(studypath, name, subject, date )
  return study
end

getMDFStudyFolderName(study::Study) = getMDFStudyFolderName(study.name, study.date)

function getMDFStudyFolderName(name::String, date::DateTime)
  return string(split(string(date),"T")[1][union(1:4,6:7,9:10)],"_",
                split(string(date),"T")[2][union(1:2,4:5,7:8)],"_",name)
end

function addStudy(d::MDFDatasetStore, study::Study)
  studypath = joinpath( studydir(d), getMDFStudyFolderName(study))
  mkpath(studypath)
  try_chmod(studypath, 0o770, recursive=true)

  nothing
end

@static if Sys.isunix()
  function findBrukerFiles(path::AbstractString, mindepth::Int=1, maxdepth::Int=2)
    candidatePaths = split(read(`find $path -maxdepth $maxdepth -mindepth $mindepth -type d`,String),"\n")[1:end-1]
    mask = zeros(Bool,length(candidatePaths))
    for (i,candidatePath) in enumerate(candidatePaths)
      if isfile(joinpath(candidatePath,"acqp")) &&
         isfile(joinpath(candidatePath,"method")) &&
         isfile(joinpath(candidatePath,"visu_pars"))
        mask[i] = true
      end
    end
    return String.(candidatePaths[mask])
  end
else
  function findBrukerFiles(path::AbstractString)
    files = readdir(path)
    bfiles = String[]
    for file in files
      if isdir(joinpath(path,file))
       try
        if isfile(joinpath(path,file,"acqp")) &&
           isfile(joinpath(candidatePath,"method")) &&
           isfile(joinpath(candidatePath,"visu_pars"))
          push!(bfiles, joinpath(path,file))
        else
          rfiles = findBrukerFiles(joinpath(path,file))
          if rfiles != nothing && length(rfiles) > 0
            push!(bfiles, rfiles...)
          end
        end
       catch
        continue
       end
      end
    end
  return bfiles
  end
end

function findSFFiles(d::BrukerDatasetStore)
  studies = readdir(d.path)
  bfiles = String[]

  for study in studies
    studypath = joinpath(d.path,study)
    if isdir(studypath)
      experiments = readdir(studypath)
      for exp in experiments
        path = joinpath(d.path,study,exp)
        if _iscalib(path)
          push!(bfiles, path)
        end
      end
    end
  end
  BrukerMDFSFs = readdir("/opt/data/MDF_SFs/")
  for BrukerMDFSF in BrukerMDFSFs
    push!(bfiles,joinpath("/opt/data/MDF_SFs/",BrukerMDFSF))
  end
  return bfiles
end

function findSFFiles(d::MDFDatasetStore)
  bfiles = String[]

  path = joinpath(s.path,"calibrations/")

  files = readdir(path)

  for file in files
    prefix, ext = splitext(file)
    if !isdir(file) && tryparse(Int64,prefix) != nothing &&
       (ext == ".mdf" || ext == ".hdf" || ext == ".h5") && !occursin("td.mdf",file)
      try
        push!(bfiles, joinpath(path,file))
      catch e
        @debug "" e
      end
    end
  end

  bfiles
end

####

function generateSFDatabase(d::DatasetStore, filename::AbstractString)
  fileList = findSFFiles(d)
  A = generateSFDatabase(fileList)
  writedlm(filename, A, ',')
end

function generateSFDatabase(fileList::Vector)

  A = Array{Any}(undef,length(fileList)+1,16)

  # Headerrow
  A[1,1] = "Name"
  A[1,2] = "Gradient"
  A[1,3] = "DFx"
  A[1,4] = "DFy"
  A[1,5] = "DFz"
  A[1,6] = "Size x"
  A[1,7] = "Size y"
  A[1,8] = "Size z"
  A[1,9] = "Bandwidth"
  A[1,10] = "Tracer"
  A[1,11] = "TracerBatch"
  A[1,12] = "DeltaSampleConcentration"
  A[1,13] = "DeltaSampleVolume"
  A[1,14] = "Path"
  A[1,15] = "StartDate"
  A[1,16] = "MeasurementTime"

  for (k,sf) in enumerate(fileList)
    i=k+1
    b = MPIFile(sf)
    _innerGenerateSFDatabase(A,i,sf,b)
  end
  return A
end

function _innerGenerateSFDatabase(A,i,sf,b)
  A[i,1] = experimentName(b)
  A[i,2] = maximum(acqGradient(b))
  df = vec(dfStrength(b)).*1e3
  A[i,3:5] .= 0.0
  for l=1:min(length(df),3)
    A[i,l+2] = df[l]
  end
  N = calibSize(b)
  A[i,6] = N[1]
  A[i,7] = N[2]
  A[i,8] = N[3]
  A[i,9] = rxBandwidth(b) / 1e6
  A[i,10] = tracerName(b)[1]
  A[i,11] = tracerBatch(b)[1]
  A[i,12] = 0.0#deltaSampleConcentration(b)
  A[i,13] = 0.0#deltaSampleVolume(b)
  A[i,14] = sf #filepath(b)
  A[i,15] = string(acqStartTime(b))
  A[i,16] = 0.0#b["PVM_ScanTimeStr"]
end

function generateSFDatabase(d::MDFDatasetStore)
  oldfile = joinpath(d.path,"SF_DatabaseOld.csv")
  newfile = joinpath(d.path,"SF_Database.csv")
  generateSFDatabase_(d, oldfile, newfile)
end

# HAAACKKK
function generateSFDatabase(d::BrukerDatasetStore)
  oldfile = "/opt/data/SF_DatabaseOld.csv"
  newfile = "/opt/data/SF_Database.csv"
  generateSFDatabase_(d, oldfile, newfile)
end

function generateSFDatabase_(d::DatasetStore, oldfile, newfile)

  if isfile(newfile)
    if isfile(oldfile)
      cp(newfile, oldfile, force=true)
    else
      cp(newfile, oldfile, force=false)
    end
  end

  generateSFDatabase(d, newfile)
end

function loadSFDatabase(d::BrukerDatasetStore)
  if isfile("/opt/data/SF_Database.csv")
    A = readdlm("/opt/data/SF_Database.csv",',')
    if size(A,2) != 16
      A = readdlm("/opt/data/SF_Database.csv",'\t')
    end
    return A
  else
    return nothing
  end
end

function loadSFDatabase(d::MDFDatasetStore)
  files = readdir(calibdir(d))
  @debug "system function database" files
  mdffiles = files[endswith.(files,".mdf") .& .!endswith.(files,"td.mdf")]
  fileList = calibdir(d).*"/".*mdffiles
  A = generateSFDatabase(fileList)
  return A
end

####

function getExperiments(d::BrukerDatasetStore, s::Study)

  files = findBrukerFiles(s.path) # make me fast

  experiments = Experiment[]

  for file in files
    try
      exp = getExperiment(file)

      push!(experiments, exp)
    catch e
      @debug "" e
    end
  end
  return experiments
end

function getExperiments(d::MDFDatasetStore, s::Study)

  files = readdir(s.path)

  experiments = Experiment[]

  @debug "Time for get Experiments"
  @time for file in files
    prefix, ext = splitext(file)
    if !isdir(file) && tryparse(Int64,prefix) != nothing &&
       (ext == ".mdf" || ext == ".hdf" || ext == ".h5") &&
       isfile(joinpath(s.path,file))
       exp = getExperiment(joinpath(s.path,file))
       push!(experiments, exp)
    end
  end
  sort!(experiments,lt=(a,b)->(a.num < b.num))
  return experiments
end

function getNewNumInFolder(d::MDFDatasetStore, path)
  if !isdir(path)
    mkpath(path)
    try_chmod(path, 0o777, recursive=true)
    return 1
  end

  files = readdir(path)
  num = 1
  if length(files) > 0
    for i=1:length(files)
      pref, ext = splitext(files[i])
      num_ = tryparse(Int64, pref)
      if num_ != nothing && num_+1>num
        num = num_+1
      end
    end
  end

  return num
end

function getNewExperimentNum(d::MDFDatasetStore, s::Study)
  return getNewNumInFolder(d, s.path)
end

function getNewExperimentPath(d::MDFDatasetStore, s::Study)
  addStudy(d,s)
  expNum = getNewExperimentNum(d,s)
  path = joinpath(studydir(d),s.name,string(expNum)*".mdf")
  # touch new mdf file
  touch(path)
  try_chmod(path, 0o660)
  return path
end

function getNewCalibNum(d::MDFDatasetStore)
  return getNewNumInFolder(d, calibdir(d))
end

function getNewCalibPath(d::MDFDatasetStore)
    calibNum = getNewCalibNum(d)
    path = joinpath(calibdir(d),string(calibNum)*".mdf")
    # touch new mdf file
    touch(path)
    try_chmod(path, 0o660)
    return path
end

####### Reconstruction Store MDF ###################

function getReco(d::MDFDatasetStore, study::Study, exp::Experiment, recoNum::Int64)
  path = joinpath(d.path, "reconstructions", getMDFStudyFolderName(study), string(exp.num), string(recoNum))
  filename = path*".mdf"
  if !isfile(filename)
    filename = path*".hdf"
  end
  r = Reconstruction(filename, recoNum, Dict())
  loadParams(r)
  return r
end

# This functions searches for recoparams and returns the corresponding recoNumber
# 0 indicates that the set of parameters was not found
function findReco(d::MDFDatasetStore, study::Study, exp::Experiment, recoParams::Dict)
  recoNum = 0

  recoParams_ = deepcopy(recoParams)
  # We do not care if the reconstruction has been done by a different
  # user. Therefore, we remove the :reconstructor field
  if haskey(recoParams_, :reconstructor)
    delete!(recoParams_, :reconstructor)
  end
  recons = getRecons(d, study, exp)
  for reco in recons
    if haskey(reco.params, :reconstructor)
      delete!(reco.params, :reconstructor)
    end
    if recoParams_ == reco.params
      recoNum = reco.num
    end
  end

  return recoNum
end

function getReco(d::MDFDatasetStore, study::Study, exp::Experiment, recoParams::Dict)
  getReco(d, study, exp, findReco(d,study,exp,recoParams) )
end

# The following function is certainly not ideal when considering a "getReco" scenario
function addReco(d::MDFDatasetStore, study::Study, exp::Experiment, image)

  outputpath = joinpath(d.path, "reconstructions", getMDFStudyFolderName(study), string(exp.num))
  # create data directory
  mkpath(outputpath)
  try_chmod(outputpath, 0o777, recursive=true)

  recoNum = getNewNumInFolder(d, outputpath)

  filepath = joinpath(outputpath, string(recoNum))

  saveRecoData(filepath*".mdf", image)
  #save(filepath*".jld","recoParams",recoParams)
end

function remove(reco::Reconstruction)
  if isfile(reco.path)
    rm(reco.path)
  end
  visufile = getVisuPath(reco)
  if isfile(visufile)
    rm(visufile)
  end
end

function save(reco::Reconstruction)
  h5open(reco.path, "r+") do file
    if exists(file, "/reconstruction/_parameters")
      o_delete(file, "/reconstruction/_parameters")
    end
    saveParams(file, "/reconstruction/_parameters", reco.params)
  end
end

function loadParams(reco::Reconstruction)

  if isfile(reco.path)
   h5open(reco.path, "r") do file
    g = file["/reconstruction"]
    if exists(g, "_parameters") #new world order
      reco.params = loadParams(reco.path, "/reconstruction/_parameters")
    else #this needs to go
      @debug "opening legacy file"
      prefix, ext = splitext(reco.path)
      reco.params = load(prefix*".jld","recoParams")
    end
   end
  end
  nothing
end

function getVisuPath(reco::Reconstruction)
  prefix, ext = splitext(reco.path)
  return prefix*".visu"
end

function getRecons(d::MDFDatasetStore, study::Study, exp::Experiment)

  recons = Reconstruction[]

  datadir = joinpath(d.path, "reconstructions", getMDFStudyFolderName(study), string(exp.num))

  if isdir(datadir)
    files = readdir(datadir)
    for file in files
        prefix, ext = splitext(file)
        fullfile = joinpath(datadir,file)

        if ext == ".hdf"  || ext == ".mdf"

        filename = joinpath(datadir,prefix*".mdf")
        if !isfile(filename)
            filename = joinpath(datadir,prefix*".hdf")
        end

        num = parse(Int64,prefix)
        r = Reconstruction(filename, num,  Dict())
        loadParams(r)

        push!(recons, r )
        end
    end
  end

  return recons
end

####### Visualization Store #######

function getVisu(d::MDFDatasetStore, study::Study, exp::Experiment, reco::Reconstruction, numVisu)

  filename = joinpath(d.path, "reconstructions", getMDFStudyFolderName(study),
			       string(exp.num), string(reco.num)*".visu")

  if isfile(filename)

  #  file = h5open(filename, "r")
 #   g = file[string(numVisu)]
    params = loadParams(filename, string(numVisu))

#    close(file)
  else
    params = Dict{Symbol,Any}()
  end
  return Visualization(filename, numVisu, params)
end

function getVisus(d::MDFDatasetStore, study::Study, exp::Experiment, reco::Reconstruction)

  visus = Visualization[]

  filename = joinpath(d.path, "reconstructions", getMDFStudyFolderName(study), string(exp.num),
			      string(reco.num)*".visu")

  if isfile(filename)

    file = h5open(filename, "r")
    g = file["/"]
    for obj in g
      key = HDF5.name(obj)
      numVisu = parse(Int64, last(splitdir(key)) )
      params = loadParams(file, key) #this is the dictionarry

      push!(visus, Visualization(filename, numVisu, params))
    end
    close(file)
  end
  return visus
end

function remove(visu::Visualization)
  h5open(visu.path, "r+") do file
    if exists(file, string(visu.num))
      o_delete(file, string(visu.num))
    end
  end
end

function save(visu::Visualization)
  if isfile(visu.path)
    file = h5open(visu.path, "r+")
  else
    file = h5open(visu.path, "w")
  end
  if exists(file, string(visu.num))
    o_delete(file, string(visu.num))
  end
  saveParams(file, string(visu.num), visu.params)
  close(file)
end

function addVisu(d::MDFDatasetStore, study::Study, exp::Experiment, reco::Reconstruction, visuParams)

  filename = joinpath(d.path, "reconstructions", getMDFStudyFolderName(study), string(exp.num),
			      string(reco.num)*".visu" )

  visus = getVisus(d, study, exp, reco)

  if isempty(visus)
    num = 1
  else
    num = last(visus).num + 1
  end

  visu =  Visualization(filename, num, visuParams)

  save(visu)
end