text:
# Copyright (c) 2013: Cbc.jl contributors
#
# Use of this source code is governed by an MIT-style license that can be found
# in the LICENSE.md file or at https://opensource.org/licenses/MIT.

module TestMOIWrapper

using Test
using MathOptInterface
import Cbc

const MOI = MathOptInterface

function runtests()
    for name in names(@__MODULE__; all = true)
        if startswith("$(name)", "test_")
            @testset "$(name)" begin
                getfield(@__MODULE__, name)()
            end
        end
    end
    return
end

function test_SolverName()
    @test MOI.get(Cbc.Optimizer(), MOI.SolverName()) ==
          "COIN Branch-and-Cut (Cbc)"
end

function test_supports_incremental_interface()
    @test !MOI.supports_incremental_interface(Cbc.Optimizer())
    return
end

function test_runtests()
    model = MOI.Utilities.CachingOptimizer(
        MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()),
        MOI.instantiate(Cbc.Optimizer; with_bridge_type = Float64),
    )
    MOI.set(model, MOI.Silent(), true)
    MOI.Test.runtests(
        model,
        MOI.Test.Config(
            exclude = Any[
                MOI.ConstraintDual,
                MOI.DualObjectiveValue,
                MOI.ConstraintBasisStatus,
                MOI.VariableBasisStatus,
            ],
        ),
        exclude = [
            # TODO(odow): upstream bug in Cbc
            "test_linear_Indicator_",
            "test_linear_SOS1_integration",
            "test_linear_SOS2_integration",
            "test_solve_SOS2_add_and_delete",
            # Can't prove infeasible.
            "test_conic_NormInfinityCone_INFEASIBLE",
            "test_conic_NormOneCone_INFEASIBLE",
            "test_solve_TerminationStatus_DUAL_INFEASIBLE",
        ],
    )
    return
end

function test_params()
    # Note: we generate a non-trivial problem to ensure that Cbc struggles to
    # find a solution at the root node.
    knapsack_model = MOI.Utilities.Model{Float64}()
    N = 100
    x = MOI.add_variables(knapsack_model, N)
    MOI.add_constraint.(knapsack_model, x, MOI.ZeroOne())
    MOI.add_constraint(
        knapsack_model,
        MOI.ScalarAffineFunction(
            MOI.ScalarAffineTerm.([1 + sin(i) for i in 1:N], x),
            0.0,
        ),
        MOI.LessThan(10.0),
    )
    MOI.set(
        knapsack_model,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
        MOI.ScalarAffineFunction(
            MOI.ScalarAffineTerm.([cos(i) for i in 1:N], x),
            0.0,
        ),
    )
    model = Cbc.Optimizer()
    MOI.set(model, MOI.RawOptimizerAttribute("maxSol"), 1)
    MOI.set(model, MOI.RawOptimizerAttribute("presolve"), "off")
    MOI.set(model, MOI.RawOptimizerAttribute("cuts"), "off")
    MOI.set(model, MOI.RawOptimizerAttribute("heur"), "off")
    MOI.set(model, MOI.RawOptimizerAttribute("logLevel"), 0)
    MOI.copy_to(model, knapsack_model)
    MOI.optimize!(model)
    @test MOI.get(model, MOI.TerminationStatus()) == MOI.SOLUTION_LIMIT
    @test MOI.get(model, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT
    MOI.empty!(model)
    MOI.copy_to(model, knapsack_model)
    MOI.optimize!(model)
    @test MOI.get(model, MOI.TerminationStatus()) == MOI.SOLUTION_LIMIT
    @test MOI.get(model, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT
    return
end

"""
    test_threads()

Test solving a model with the threads parameter set.

See issues #112 and #186.
"""
function test_threads()
    model = MOI.Utilities.CachingOptimizer(
        MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()),
        MOI.instantiate(Cbc.Optimizer; with_bridge_type = Float64),
    )
    MOI.set(model, MOI.RawOptimizerAttribute("presolve"), "off")
    MOI.set(model, MOI.RawOptimizerAttribute("cuts"), "off")
    MOI.set(model, MOI.RawOptimizerAttribute("heur"), "off")
    MOI.set(model, MOI.RawOptimizerAttribute("threads"), 4)
    MOI.set(model, MOI.RawOptimizerAttribute("logLevel"), 3)
    N = 100
    x = MOI.add_variables(model, N)
    MOI.add_constraint.(model, x, MOI.ZeroOne())
    w = [1 + sin(i) for i in 1:N]
    c = [1 + cos(i) for i in 1:N]
    MOI.add_constraint(
        model,
        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(w, x), 0.0),
        MOI.LessThan(10.0),
    )
    obj = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(c, x), 0.0)
    MOI.set(model, MOI.ObjectiveFunction{typeof(obj)}(), obj)
    MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE)
    MOI.optimize!(model)
    @test MOI.get(model, MOI.TerminationStatus()) == MOI.OPTIMAL
    return
end

function test_PrimalStatus()
    model = MOI.Utilities.Model{Float64}()
    x = MOI.add_variable(model)
    MOI.add_constraint(model, x, MOI.GreaterThan(1.0))
    MOI.add_constraint(model, x, MOI.LessThan(0.0))
    cbc = Cbc.Optimizer()
    MOI.copy_to(cbc, model)
    MOI.optimize!(cbc)
    MOI.get(cbc, MOI.PrimalStatus()) == MOI.NO_SOLUTION
    return
end

function test_issue_187()
    if Sys.iswindows() || Sys.islinux()
        # This test segfaults on Windows and Linux
        @test_broken 1 == 2
        return
    end
    model = MOI.Utilities.CachingOptimizer(
        MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()),
        MOI.instantiate(Cbc.Optimizer; with_bridge_type = Float64),
    )
    MOI.set(model, MOI.Silent(), true)
    x = MOI.add_variables(model, 2)
    MOI.add_constraint.(model, x, MOI.ZeroOne())
    MOI.set.(model, MOI.VariablePrimalStart(), x, 0.0)
    y = MOI.add_variables(model, 2)
    MOI.add_constraint.(model, y, MOI.ZeroOne())

    MOI.add_constraint(
        model,
        MOI.Utilities.operate(vcat, Float64, x[1], 1.0 * y[1]),
        MOI.Indicator{MOI.ACTIVATE_ON_ONE}(MOI.EqualTo(1.0)),
    )
    MOI.add_constraint(
        model,
        MOI.Utilities.operate(vcat, Float64, x[2], 1.0 * y[1]),
        MOI.Indicator{MOI.ACTIVATE_ON_ONE}(MOI.EqualTo(0.0)),
    )
    MOI.set.(model, MOI.VariablePrimalStart(), y, 0.0)
    MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE)
    f = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(1.0, x), 0.0)
    MOI.set(model, MOI.ObjectiveFunction{typeof(f)}(), f)
    MOI.optimize!(model)
    @test MOI.get(model, MOI.TerminationStatus()) == MOI.OPTIMAL
    @test ≈(sum(MOI.get(model, MOI.VariablePrimal(), x)), 1.0, atol = 1e-4)
    return
end

# The test_linear_SOS1_integration test with the additional requirement that all
# variables are integer.
function test_SOS1()
    model = MOI.Bridges.full_bridge_optimizer(
        MOI.Utilities.CachingOptimizer(
            MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()),
            Cbc.Optimizer(),
        ),
        Float64,
    )
    config = MOI.Test.Config()
    MOI.set(model, MOI.Silent(), true)
    @test MOI.supports_constraint(
        model,
        MOI.VectorOfVariables,
        MOI.SOS1{Float64},
    )
    @test MOI.supports_constraint(
        model,
        MOI.VariableIndex,
        MOI.LessThan{Float64},
    )
    v = MOI.add_variables(model, 3)
    MOI.add_constraint.(model, v, MOI.Integer())
    @test MOI.get(model, MOI.NumberOfVariables()) == 3
    vc1 = MOI.add_constraint(model, v[1], MOI.LessThan(1.0))
    @test vc1.value == v[1].value
    vc2 = MOI.add_constraint(model, v[2], MOI.LessThan(1.0))
    @test vc2.value == v[2].value
    vc3 = MOI.add_constraint(model, v[3], MOI.LessThan(2.0))
    @test vc3.value == v[3].value
    c1 = MOI.add_constraint(
        model,
        MOI.VectorOfVariables([v[1], v[2]]),
        MOI.SOS1([1.0, 2.0]),
    )
    c2 = MOI.add_constraint(
        model,
        MOI.VectorOfVariables([v[1], v[3]]),
        MOI.SOS1([1.0, 2.0]),
    )
    @test MOI.get(
        model,
        MOI.NumberOfConstraints{MOI.VectorOfVariables,MOI.SOS1{Float64}}(),
    ) == 2
    #=
        To allow for permutations in the sets and variable vectors
        we're going to sort according to the weights
    =#
    cs_sos = MOI.get(model, MOI.ConstraintSet(), c2)
    cf_sos = MOI.get(model, MOI.ConstraintFunction(), c2)
    p = sortperm(cs_sos.weights)
    @test isapprox(cs_sos.weights[p], [1.0, 2.0], config)
    @test cf_sos.variables[p] == v[[1, 3]]
    objf =
        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([2.0, 1.0, 1.0], v), 0.0)
    MOI.set(
        model,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
        objf,
    )
    MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE)
    @test MOI.get(model, MOI.ObjectiveSense()) == MOI.MAX_SENSE
    @test MOI.get(model, MOI.TerminationStatus()) == MOI.OPTIMIZE_NOT_CALLED
    MOI.optimize!(model)
    @test MOI.get(model, MOI.TerminationStatus()) == config.optimal_status
    @test MOI.get(model, MOI.ResultCount()) >= 1
    @test MOI.get(model, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT
    @test isapprox(MOI.get(model, MOI.ObjectiveValue()), 3, config)
    @test isapprox(MOI.get(model, MOI.VariablePrimal(), v), [0, 1, 2], config)
    MOI.delete(model, c1)
    MOI.delete(model, c2)
    MOI.optimize!(model)
    @test MOI.get(model, MOI.TerminationStatus()) == config.optimal_status
    @test MOI.get(model, MOI.ResultCount()) >= 1
    @test MOI.get(model, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT
    @test isapprox(MOI.get(model, MOI.ObjectiveValue()), 5, config)
    @test isapprox(MOI.get(model, MOI.VariablePrimal(), v), [1, 1, 2], config)
    return
end

# The test_linear_SOS2_integration test with the additional requirement that all
# variables are integer.
function test_SOS2()
    model = MOI.Bridges.full_bridge_optimizer(
        MOI.Utilities.CachingOptimizer(
            MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()),
            Cbc.Optimizer(),
        ),
        Float64,
    )
    config = MOI.Test.Config()
    MOI.set(model, MOI.Silent(), true)
    @test MOI.supports_constraint(
        model,
        MOI.VectorOfVariables,
        MOI.SOS1{Float64},
    )
    @test MOI.supports_constraint(
        model,
        MOI.VectorOfVariables,
        MOI.SOS2{Float64},
    )
    v = MOI.add_variables(model, 10)
    @test MOI.get(model, MOI.NumberOfVariables()) == 10
    bin_constraints = []
    for i in 1:8
        vc = MOI.add_constraint(model, v[i], MOI.Interval(0.0, 2.0))
        @test vc.value == v[i].value
        push!(bin_constraints, MOI.add_constraint(model, v[i], MOI.ZeroOne()))
        @test bin_constraints[i].value == v[i].value
    end
    MOI.add_constraint(
        model,
        MOI.ScalarAffineFunction(
            MOI.ScalarAffineTerm.([1.0, 2.0, 3.0, -1.0], v[[1, 2, 3, 9]]),
            0.0,
        ),
        MOI.EqualTo(0.0),
    )
    MOI.add_constraint(
        model,
        MOI.ScalarAffineFunction(
            MOI.ScalarAffineTerm.(
                [5.0, 4.0, 7.0, 2.0, 1.0, -1.0],
                v[[4, 5, 6, 7, 8, 10]],
            ),
            0.0,
        ),
        MOI.EqualTo(0.0),
    )
    MOI.add_constraint(
        model,
        MOI.VectorOfVariables(v[[1, 2, 3]]),
        MOI.SOS1([1.0, 2.0, 3.0]),
    )
    vv = MOI.VectorOfVariables(v[[4, 5, 6, 7, 8]])
    sos2 = MOI.SOS2([5.0, 4.0, 7.0, 2.0, 1.0])
    c = MOI.add_constraint(model, vv, sos2)
    #=
        To allow for permutations in the sets and variable vectors
        we're going to sort according to the weights
    =#
    cs_sos = MOI.get(model, MOI.ConstraintSet(), c)
    cf_sos = MOI.get(model, MOI.ConstraintFunction(), c)
    p = sortperm(cs_sos.weights)
    @test isapprox(cs_sos.weights[p], [1.0, 2.0, 4.0, 5.0, 7.0], config)
    @test cf_sos.variables[p] == v[[8, 7, 5, 4, 6]]
    objf = MOI.ScalarAffineFunction(
        MOI.ScalarAffineTerm.([1.0, 1.0], [v[9], v[10]]),
        0.0,
    )
    MOI.set(
        model,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
        objf,
    )
    MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE)
    @test MOI.get(model, MOI.ObjectiveSense()) == MOI.MAX_SENSE
    @test MOI.get(model, MOI.TerminationStatus()) == MOI.OPTIMIZE_NOT_CALLED
    MOI.optimize!(model)
    @test MOI.get(model, MOI.TerminationStatus()) == config.optimal_status
    @test MOI.get(model, MOI.ResultCount()) >= 1
    @test MOI.get(model, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT
    @test isapprox(MOI.get(model, MOI.ObjectiveValue()), 15.0, config)
    @test isapprox(
        MOI.get(model, MOI.VariablePrimal(), v),
        [0.0, 0.0, 1.0, 1.0, 0.0, 1.0, 0.0, 0.0, 3.0, 12.0],
        config,
    )
    for cref in bin_constraints
        MOI.delete(model, cref)
    end
    MOI.add_constraint.(model, v, MOI.Integer())
    MOI.optimize!(model)
    @test MOI.get(model, MOI.TerminationStatus()) == config.optimal_status
    @test MOI.get(model, MOI.ResultCount()) >= 1
    @test MOI.get(model, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT
    @test isapprox(MOI.get(model, MOI.ObjectiveValue()), 30.0, config)
    @test isapprox(
        MOI.get(model, MOI.VariablePrimal(), v),
        [0.0, 0.0, 2.0, 2.0, 0.0, 2.0, 0.0, 0.0, 6.0, 24.0],
        config,
    )
    return
end

"""
    test_VariablePrimalStart()

Testing that VariablePrimalStart is actually applied is a little convoluted.

We formulate a MIP with various setttings turned off to avoid a trivial solve in
presolve.

Then we solve and return the optimal primal solution and the number of nodes
visited.

For the second pass, we rebuild the same MIP, but this time we pass the optimal
solution as the VariablePrimalStart, and we set maxSol=1 to force Cbc to exit
after finding a single solution. Because we passed a primal feasible point, it
should return the optimal solution after exploring 0 nodes.
"""
function test_VariablePrimalStart()
    function formulate_and_solve(start)
        model = MOI.Utilities.CachingOptimizer(
            MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()),
            MOI.instantiate(Cbc.Optimizer; with_bridge_type = Float64),
        )
        MOI.set(model, MOI.RawOptimizerAttribute("presolve"), "off")
        MOI.set(model, MOI.RawOptimizerAttribute("cuts"), "off")
        MOI.set(model, MOI.RawOptimizerAttribute("heur"), "off")
        MOI.set(model, MOI.RawOptimizerAttribute("logLevel"), 0)
        N = 100
        x = MOI.add_variables(model, N)
        MOI.add_constraint.(model, x, MOI.ZeroOne())
        w = [1 + sin(i) for i in 1:N]
        c = [1 + cos(i) for i in 1:N]
        if start !== nothing
            MOI.set(model, MOI.RawOptimizerAttribute("maxSol"), 1)
            MOI.set.(model, MOI.VariablePrimalStart(), x, start)
        end
        MOI.add_constraint(
            model,
            MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(w, x), 0.0),
            MOI.LessThan(10.0),
        )
        obj = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(c, x), 0.0)
        MOI.set(model, MOI.ObjectiveFunction{typeof(obj)}(), obj)
        MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE)
        MOI.optimize!(model)
        sol = MOI.get.(model, MOI.VariablePrimal(), x)
        return sol, MOI.get(model, MOI.NodeCount())
    end
    x_sol, nodes = formulate_and_solve(nothing)
    y_sol, nodes_start = formulate_and_solve(x_sol)
    @test x_sol == y_sol
    @test nodes > 0
    @test nodes_start == 0
    return
end

end

TestMOIWrapper.runtests()


text:
"""
`mer_caltarget(color::Symbol)`

Returns the spectrum of a given MER Caltarget color patch.

Allowed color patch names:

`:white, :grey, :black, :yellow, :red, :green, :blue`

Result: `RSpec`
"""
function mer_caltarget(color::Symbol)
   color == :white  ? reflectance_spec(MER_caltarget_white[:, 1], MER_caltarget_white[:, 2]) :
   color == :grey   ? reflectance_spec(MER_caltarget_grey[:, 1], MER_caltarget_grey[:, 2]) :
   color == :black  ? reflectance_spec(MER_caltarget_black[:, 1], MER_caltarget_black[:, 2]) :
   color == :yellow ? reflectance_spec(MER_caltarget_yellow[:, 1], MER_caltarget_yellow[:, 2]) :
   color == :red    ? reflectance_spec(MER_caltarget_red[:, 1], MER_caltarget_red[:, 2]) :
   color == :green  ? reflectance_spec(MER_caltarget_green[:, 1], MER_caltarget_green[:, 2]) :
   color == :blue   ? reflectance_spec(MER_caltarget_blue[:, 1], MER_caltarget_blue[:, 2]) :
   throw(DomainError(color, "Nonexistent MER calibration target color."))
end

"""
`mer_caltarget(color::Int)`

Returns the spectrum of a given MER Caltarget color patch.

Allowed color patch indices: `1 ≤ color ≤ 7`.

Result: `RSpec`
"""
function mer_caltarget(num::Int)
   if num < 1 || num > 7
      error("Color does not exist.")
   else
      i=[:white, :grey, :black, :yellow, :red, :green, :blue]
      mer_caltarget(i[num])
   end
end

const MER_caltarget_yellow=
   [400.0	      0.052336;
    430.841121	   0.059813;
    478.504673	   0.069782;
    491.121495	   0.079751;
    502.336449	   0.08972;
    514.953271	   0.114642;
    524.766355	   0.142056;
    530.373832	   0.161994;
    535.981308	   0.186916;
    538.006231	   0.195777;
    543.925234	   0.231499;
    550.0	      0.265282;
    556.074766	   0.301004;
    561.993769	   0.332572;
    568.068536	   0.358048;
    573.987539	   0.373555;
    580.062305	   0.381585;
    585.981308	   0.387124;
    592.056075	   0.386847;
    597.975078	   0.383801;
    604.049844	   0.377155;
    609.968847	   0.374663;
    616.043614	   0.370786;
    621.962617	   0.368294;
    628.037383	   0.365247;
    633.956386	   0.364694;
    640.031153	   0.361648;
    645.950156	   0.360263;
    652.024922	   0.363863;
    657.943925	   0.364694;
    664.018692	   0.364971;
    669.937695	   0.371616;
    676.012461	   0.377432;
    693.925234	   0.3982;
    700.0	      0.404292;
    706.074766	   0.415092;
    711.993769	   0.425337;
    718.068536	   0.433091;
    723.987539	   0.440568;
    730.062305	   0.449706;
    735.981308	   0.459121;
    742.056075	   0.462167;
    747.975078	   0.46549;
    754.049844	   0.467982;
    759.968847	   0.468259;
    766.043614	   0.462721;
    771.962617	   0.459121;
    778.037383	   0.452198;
    783.956386	   0.448044;
    807.943925	   0.414815;
    814.018692	   0.40983;
    819.937695	   0.400969;
    826.012461	   0.394323;
    831.931464	   0.388785;
    838.006231	   0.381862;
    843.925234	   0.376324;
    850.0	      0.37134;
    856.074766	   0.366078;
    861.993769	   0.362201;
    868.068536	   0.358602;
    873.987539	   0.358602;
    880.062305	   0.355832;
    885.981308	   0.352233;
    892.056075	   0.351956;
    897.975078	   0.351679;
    904.049844	   0.349463;
    909.968847	   0.350294;
    916.043614	   0.351679;
    921.962617	   0.352786;
    928.037383	   0.352786;
    933.956386	   0.356386;
    940.031153	   0.356663;
    945.950156	   0.359709;
    952.024922	   0.361648;
    957.943925	   0.363863;
    964.018692	   0.366632;
    976.012461	   0.372447;
    981.931464	   0.376047;
    988.006231	   0.377985;
    993.925234	   0.381862;
   1000.0	      0.38657;
   1012.616822	   0.393769;
   1021.028037	   0.398754;
   1047.663551	   0.42866]

const MER_caltarget_red=
   [400.0         0.047352;
    433.64486	   0.049844;
    444.859813	   0.049844;
    450.46729	   0.052336;
    457.476636	   0.049844;
    463.084112	   0.054829;
    465.88785	   0.054829;
    468.691589	   0.052336;
    530.373832	   0.054829;
    544.392523	   0.057321;
    551.401869	   0.062305;
    557.009346	   0.06729;
    561.214953	   0.069782;
    566.82243	   0.079751;
    573.831776	   0.097196;
    585.046729	   0.129595;
    597.663551	   0.161994;
    604.672897	   0.179439;
    610.280374	   0.181931;
    615.88785	   0.184424;
    618.691589	   0.186916;
    625.700935	   0.189408;
    629.906542	   0.1919;
    634.11215	   0.189408;
    642.523364	   0.194393;
    646.728972	   0.199377;
    657.943925	   0.201869;
    667.757009	   0.206854;
    678.971963	   0.219315;
    683.17757	   0.224299;
    685.981308	   0.224299;
    691.588785	   0.231776;
    704.205607	   0.249221;
    708.411215	   0.256698;
    712.616822	   0.261682;
    719.626168	   0.266667;
    725.233645	   0.276636;
    737.850467	   0.279128;
    746.261682	   0.276636;
    757.476636	   0.264174;
    760.280374	   0.264174;
    764.485981	   0.261682;
    772.897196	   0.249221;
    778.504673	   0.241745;
    781.308411	   0.241745;
    792.523364	   0.226791;
    803.738318	   0.216822;
    814.953271	   0.206854;
    820.560748	   0.201869;
    828.971963	   0.199377;
    841.588785	   0.194393;
    862.616822	   0.1919;
    885.046729	   0.194393;
    890.654206	   0.199377;
    894.859813	   0.199377;
    899.065421	   0.196885;
    908.878505	   0.206854;
    913.084112	   0.206854;
    920.093458	   0.216822;
    924.299065	   0.216822;
    938.317757	   0.234268;
    943.925234	   0.239252;
    946.728972	   0.249221;
    955.140187	   0.256698;
    959.345794	   0.269159;
    970.560748	   0.289097;
    977.570093	   0.30405;
    985.981308	   0.321495;
    992.990654	   0.343925;
    998.598131	   0.353894;
   1001.401869	   0.358879;
   1005.607477	   0.37134;
   1014.018692	   0.383801;
   1018.224299	   0.391277;
   1028.037383	   0.413707;
   1033.64486	   0.42866;
   1037.850467	   0.436137;
   1047.663551	   0.45109]

const MER_caltarget_green=
   [400.0	      0.099688;
    406.074766	   0.102735;
    411.993769	   0.096642;
    418.068536	   0.089997;
    423.987539	   0.079751;
    430.062305	   0.072274;
    442.056075	   0.06369;
    459.968847	   0.063967;
    466.043614	   0.064244;
    471.962617	   0.065351;
    483.956386	   0.072828;
    490.031153	   0.076982;
    495.950156	   0.081689;
    502.024922	   0.088335;
    507.943925	   0.099135;
    514.018692	   0.118795;
    519.937695	   0.147594;
    526.012461	   0.166424;
    531.931464	   0.17224;
    538.006231	   0.168086;
    543.925234	   0.156179;
    550.0	      0.140948;
    556.074766	   0.12378;
    561.993769	   0.10855;
    568.068536	   0.094704;
    573.987539	   0.088889;
    580.062305	   0.083904;
    585.981308	   0.081689;
    592.056075	   0.081412;
    597.975078	   0.080582;
    604.049844	   0.081689;
    609.968847	   0.081966;
    616.043614	   0.084181;
    621.962617	   0.088889;
    628.037383	   0.093319;
    633.956386	   0.098858;
    640.031153	   0.105781;
    645.950156	   0.111596;
    652.024922	   0.119349;
    669.937695	   0.142333;
    676.012461	   0.15064;
    681.931464	   0.157286;
    688.006231	   0.16504;
    693.925234	   0.174178;
    700.0	      0.181101;
    706.074766	   0.190239;
    711.993769	   0.204638;
    718.068536	   0.220422;
    723.987539	   0.239806;
    730.062305	   0.255867;
    735.981308	   0.28965;
    742.056075	   0.306819;
    747.975078	   0.332572;
    754.049844	   0.358048;
    759.968847	   0.379093;
    766.043614	   0.396816;
    771.962617	   0.409553;
    778.037383	   0.422291;
    783.956386	   0.426445;
    790.031153	   0.432537;
    795.950156	   0.435583;
    802.024922	   0.436414;
    807.943925	   0.43946;
    814.018692	   0.441398;
    819.937695	   0.442229;
    826.012461	   0.442783;
    831.931464	   0.44306;
    838.006231	   0.444444;
    843.925234	   0.445552;
    850.0	      0.445829;
    856.074766	   0.445829;
    861.993769	   0.445829;
    868.068536	   0.446106;
    873.987539	   0.448321;
    880.062305	   0.448875;
    885.981308	   0.448875;
    892.056075	   0.448321;
    897.975078	   0.446383;
    904.049844	   0.446106;
    909.968847	   0.448321;
    916.043614	   0.449152;
    921.962617	   0.45109;
    928.037383	   0.454136;
    933.956386	   0.453306;
    940.031153	   0.452198;
    945.950156	   0.454136;
    952.024922	   0.454413;
    957.943925	   0.454967;
    964.018692	   0.454967;
    969.937695	   0.451367;
    976.012461	   0.45469;
    981.931464	   0.453859;
    988.006231	   0.453859;
    993.925234	   0.452752;
   1000.0	      0.45469;
   1006.074766	   0.450813;
   1018.068536	   0.45109;
   1023.987539	   0.450537;
   1030.062305	   0.451644;
   1035.981308	   0.452475;
   1042.056075	   0.452752;
   1047.663551	   0.453583]

const MER_caltarget_blue=
   [400.0	      0.173624;
    406.074766	   0.18027;
    411.993769	   0.206854;
    418.068536	   0.222084;
    423.987539	   0.234545;
    430.062305	   0.246452;
    435.981308	   0.253375;
    442.056075	   0.253375;
    447.975078	   0.253375;
    454.049844	   0.252821;
    459.968847	   0.25116;
    466.043614	   0.241468;
    471.962617	   0.223745;
    478.037383	   0.209069;
    483.956386	   0.227345;
    490.031153	   0.233714;
    495.950156	   0.218761;
    502.024922	   0.178055;
    507.943925	   0.138456;
    519.937695	   0.080028;
    526.012461	   0.067013;
    531.931464	   0.062028;
    550.0	      0.050675;
    556.074766	   0.050398;
    561.993769	   0.049844;
    568.068536	   0.049844;
    573.987539	   0.049844;
    580.062305	   0.049844;
    585.981308	   0.049844;
    592.056075	   0.049567;
    597.975078	   0.049567;
    604.049844	   0.049567;
    609.968847	   0.049567;
    616.043614	   0.049567;
    621.962617	   0.050952;
    628.037383	   0.052613;
    633.956386	   0.055659;
    640.031153	   0.061198;
    645.950156	   0.070059;
    652.024922	   0.085566;
    657.943925	   0.10855;
    676.012461	   0.271097;
    681.931464	   0.342818;
    688.006231	   0.449983;
    693.925234	   0.539979;
    700.0	      0.594808;
    706.074766	   0.63856;
    711.993769	   0.667636;
    718.068536	   0.684804;
    723.987539	   0.695327;
    730.062305	   0.698927;
    735.981308	   0.706127;
    742.056075	   0.706957;
    747.975078	   0.711665;
    754.049844	   0.716096;
    759.968847	   0.716649;
    766.043614	   0.716649;
    771.962617	   0.716096;
    778.037383	   0.71748;
    783.956386	   0.718865;
    790.031153	   0.719142;
    795.950156	   0.717757;
    802.024922	   0.716926;
    807.943925	   0.716649;
    814.018692	   0.716096;
    819.937695	   0.717757;
    826.012461	   0.716372;
    831.931464	   0.716372;
    838.006231	   0.714988;
    843.925234	   0.707511;
    850.0	      0.706957;
    856.074766	   0.704465;
    861.993769	   0.703911;
    868.068536	   0.699481;
    873.987539	   0.695604;
    880.062305	   0.694496;
    885.981308	   0.689235;
    892.056075	   0.675666;
    897.975078	   0.652129;
    904.049844	   0.639668;
    909.968847	   0.656282;
    916.043614	   0.671236;
    921.962617	   0.674005;
    928.037383	   0.672343;
    933.956386	   0.665697;
    940.031153	   0.659882;
    945.950156	   0.653236;
    952.024922	   0.642991;
    957.943925	   0.63136;
    964.018692	   0.620838;
    969.937695	   0.606715;
    976.012461	   0.58567;
    981.931464	   0.570717;
    988.006231	   0.554379;
    993.925234	   0.530564;
   1000.0	      0.510903;
   1006.074766	   0.484874;
   1018.068536	   0.433645;
   1035.981308	   0.36857;
   1042.056075	   0.34974;
   1047.663551	   0.331464]

const MER_caltarget_white=
   [400.0	      0.272222;
    402.790698	   0.311111;
    408.372093	   0.35;
    412.55814	   0.386728;
    416.744186	   0.414815;
    419.534884	   0.425617;
    425.116279	   0.43642;
    433.488372	   0.458025;
    448.837209	   0.483951;
    461.395349	   0.507716;
    479.534884	   0.531481;
    497.674419	   0.550926;
    514.418605	   0.557407;
    521.395349	   0.563889;
    531.162791	   0.566049;
    549.302326	   0.574691;
    564.651163	   0.574691;
    578.604651	   0.576852;
    589.767442	   0.574691;
    593.953488	   0.572531;
    598.139535	   0.576852;
    619.069767	   0.572531;
    631.627907	   0.57037;
    640.0	      0.57037;
    648.372093	   0.566049;
    667.906977	   0.563889;
    681.860465	   0.559568;
    697.209302	   0.559568;
    709.767442	   0.555247;
    720.930233	   0.553086;
    730.697674	   0.550926;
    734.883721	   0.550926;
    743.255814	   0.548765;
    762.790698	   0.544444;
    771.162791	   0.542284;
    776.744186	   0.542284;
    785.116279	   0.537963;
    790.697674	   0.540123;
    794.883721	   0.540123;
    811.627907	   0.533642;
    828.372093	   0.531481;
    840.930233	   0.529321;
    850.697674	   0.52284;
    884.186047	   0.518519;
    891.162791	   0.514198;
    899.534884	   0.512037;
    910.697674	   0.512037;
    919.069767	   0.514198;
    935.813953	   0.512037;
    949.767442	   0.509877;
    955.348837	   0.505556;
    967.906977	   0.505556;
    984.651163	   0.501235;
   1000.0	      0.496914;
   1012.55814	   0.494753;
   1020.930233	   0.490432;
   1030.697674	   0.492593;
   1047.44186	   0.490432]

const MER_caltarget_grey=
   [400.0         0.239815;
    404.186047	   0.259259;
    411.162791	   0.295988;
    416.744186	   0.315432;
    433.488372	   0.337037;
    453.023256	   0.356481;
    469.767442	   0.375926;
    483.72093	   0.375926;
    506.046512	   0.384568;
    539.534884	   0.384568;
    563.255814	   0.382407;
    584.186047	   0.378086;
    598.139535	   0.375926;
    617.674419	   0.367284;
    630.232558	   0.365123;
    644.186047	   0.365123;
    651.162791	   0.360802;
    672.093023	   0.354321;
    686.046512	   0.35;
    702.790698	   0.345679;
    726.511628	   0.337037;
    740.465116	   0.337037;
    747.44186	   0.332716;
    760.0         0.328395;
    780.930233	   0.326235;
    800.465116	   0.319753;
    814.418605	   0.319753;
    836.744186	   0.313272;
    856.27907	   0.30679;
    868.837209	   0.30463;
    892.55814	   0.300309;
    926.046512	   0.289506;
    946.976744	   0.285185;
    969.302326	   0.280864;
    979.069767	   0.280864;
   1000.0         0.276543;
   1012.55814	   0.270062;
   1022.325581	   0.270062;
   1033.488372	   0.267901;
   1048.837209	   0.267901]

const MER_caltarget_black=
   [400.0     	   0.168519;
    408.372093	   0.17284;
    430.697674	   0.183642;
    439.069767	   0.187963;
    461.395349	   0.190123;
    469.767442	   0.190123;
    493.488372	   0.190123;
    501.860465	   0.187963;
    522.790698	   0.185802;
    533.953488	   0.183642;
    553.488372	   0.181481;
    564.651163	   0.181481;
    584.186047	   0.17716;
    595.348837	   0.17716;
    613.488372	   0.17284;
    626.046512	   0.170679;
    642.790698	   0.166358;
    672.093023	   0.162037;
    701.395349	   0.155556;
    730.697674	   0.151235;
    760.0         0.151235;
    790.697674	   0.144753;
    821.395349	   0.142593;
    829.767442	   0.138272;
    849.302326	   0.138272;
    877.209302	   0.138272;
    888.372093	   0.133951;
    905.116279	   0.133951;
    935.813953	   0.13179;
    966.511628	   0.13179;
    995.813953	   0.127469;
   1025.116279	   0.120988;
   1034.883721	   0.127469;
   1048.837209	   0.12963]

ext: jl


text:
using Serialization, DataStructures, Random

export computePredictionMetricsSparseVectorData

# Returns a priority queue of ordered item counts, for a dataset in sparse vector format
function getOrderedItemCountsSparseVectorData(trainingBasketsDictFileName,
                                              trainingBasketsDictObjectName)

  # Load training data
  trainingUsersBasketsDict = open(deserialize, trainingBasketsDictFileName);
  println("Loaded $trainingBasketsDictFileName")

  # Build set of training instances
  numTrainingInstances = length(collect(keys(trainingUsersBasketsDict)))
  trainingInstances = fill(Vector{Int}(), numTrainingInstances)
  trainingInstanceIndex = 1
  numItems = 0
  for trainingInstanceUserId in collect(keys(trainingUsersBasketsDict))
    trainingInstanceItems = deepcopy(trainingUsersBasketsDict[trainingInstanceUserId].basketItems)
    heldOutItem = trainingUsersBasketsDict[trainingInstanceUserId].heldOutItem
    push!(trainingInstanceItems, heldOutItem)

    trainingInstances[trainingInstanceIndex] = trainingInstanceItems

    trainingInstanceIndex += 1

    numItems = trainingUsersBasketsDict[trainingInstanceUserId].numItemsInCatalog
  end

  # Compute ordered item counts, by descending popularity
  orderedItemCounts = PriorityQueue(Dict{Int, Int}(), Base.Order.Reverse)
  for i = 1:numTrainingInstances
    trainingInstance = trainingInstances[i]

    for itemId in trainingInstance
      if haskey(orderedItemCounts, itemId)
        orderedItemCounts[itemId] += 1
      else
        # Create a new entry for this item
        orderedItemCounts[itemId] = 1
      end
    end
  end

  return orderedItemCounts
end

# Compute MPR and precision metrics from DPP prediction results for a dataset in sparse vector format
function computePredictionMetricsSparseVectorData(trainingBasketsDictFileName, trainingBasketsDictObjectName,
                                                  testBasketsDictFileName, testBasketsDictObjectName,
                                                  resultsForTestInstancesDictFileName)
  resultsForTestInstancesDictObjectName = "resultsForTestInstancesDict"
  Random.seed!(1234)

  # Load test data
  testUsersBasketsDict = open(deserialize, testBasketsDictFileName);
  println("Loaded $testBasketsDictFileName")

  # Build set of test instances
  numTestInstances = length(collect(keys(testUsersBasketsDict)))
  testInstances = fill(Vector{Int}(), numTestInstances)
  testInstanceIndex = 1
  numItems = 0
  for testInstanceUserId in collect(keys(testUsersBasketsDict))
    testInstanceItems = deepcopy(testUsersBasketsDict[testInstanceUserId].basketItems)
    heldOutItem = testUsersBasketsDict[testInstanceUserId].heldOutItem
    push!(testInstanceItems, heldOutItem)

    testInstances[testInstanceIndex] = testInstanceItems

    testInstanceIndex += 1

    numItems = testUsersBasketsDict[testInstanceUserId].numItemsInCatalog
  end

  # Compute popularity-stratification weights for each item, where the weight is
  # proportional to 1 / (itemCount ^ beta), for beta = 0.5
  orderedItemCounts = getOrderedItemCountsSparseVectorData(trainingBasketsDictFileName,
                                                           trainingBasketsDictObjectName)
  sumItemCounts = sum(values(orderedItemCounts))
  itemPopWeights = PriorityQueue(Dict{Int, Float64}(), Base.Order.Reverse)
  beta = 0.5
  startingItemId = 1
  for i = startingItemId:numItems
    if haskey(orderedItemCounts, i)
      itemPopWeights[i] = 1 / (orderedItemCounts[i] ^ beta)
    else
      # If an item id is not available in orderedItemCounts, then its popularity is 0.
      itemPopWeights[i] = 0
      println("itemPopWeights[$i]: 0")
    end
  end
  # Normalize itemPopWeights, so they sum to 1
  sumItemPopWeights = sum(collect(values(itemPopWeights)))
  for i = startingItemId:numItems
    itemPopWeights[i] = itemPopWeights[i] / sumItemPopWeights
  end

  # Load DPP prediction results from disk
  resultsForTestInstancesDict = open(deserialize, resultsForTestInstancesDictFileName);
  println("Loaded $resultsForTestInstancesDictFileName")

  # Compute mean percentile rank and precision across predictions for all test instances
  shuffle!(testInstances)
  numTestInstances = length(testInstances)
  sumPredictedRanks = 0
  basketSizeResultsDict = Dict{Int, BasketSizeResult}()
  predictedRanks = Int[]
  for i = 1:numTestInstances
    testInstance = testInstances[i]

    # Get predicted rank for test instance
    testResult = resultsForTestInstancesDict[testInstance]
    sumPredictedRanks += testResult.predictedRankActualNextItem
    push!(predictedRanks, testResult.predictedRankActualNextItem)

    # Last item in the test instance is not part of the basket
    basketSize = length(testInstance) - 1
    # Update metrics for the size of the current basket
    if basketSize >= 2
      if haskey(basketSizeResultsDict, basketSize)
        basketSizeResult = basketSizeResultsDict[basketSize]
        basketSizeResult.numBasketInstances += 1
        basketSizeResult.sumPredictedRanks += testResult.predictedRankActualNextItem
        push!(basketSizeResult.predictedRanks, testResult.predictedRankActualNextItem)
        basketSizeResultsDict[basketSize] = basketSizeResult
      else
        # Create a new entry for this basket size
        basketSizeResult = BasketSizeResult(basketSize, 1, testResult.predictedRankActualNextItem,
                                            [testResult.predictedRankActualNextItem])
        basketSizeResultsDict[basketSize] = basketSizeResult
      end

      if i % 100 == 0
        println("Processed $i test instances")
      end
    end
  end

  # Compute overall mean percentile rank
  meanRank = sumPredictedRanks / numTestInstances
  meanRankPercent = meanRank / numItems
  meanPercentileRank = 100.0 - (100.0 * meanRankPercent)
  println("meanPercentileRank for all baskets: $meanPercentileRank")

  # Compute mean percentile rank for each basket size
  for basketSize in sort(collect(keys(basketSizeResultsDict)))
    basketSizeResult = basketSizeResultsDict[basketSize]
    meanRank = basketSizeResult.sumPredictedRanks / basketSizeResult.numBasketInstances
    meanRankPercent = meanRank / numItems
    meanPercentileRank = 100.0 - (100.0 * meanRankPercent)
    println("For basket size $basketSize, meanPercentileRank: $meanPercentileRank")
  end

  # Compute overall precision@1 (accuracy)
  hits = length(findall(x -> x < 1, predictedRanks))
  precisionAtOnePercent = hits / numTestInstances * 100
  println("")
  println("precision@1 for all baskets: $precisionAtOnePercent")

  # Compute overall popularity-weighted precision@1
  popWeightedPrecisionAtOnePercent =
    computePopularityWeightedPrecisionAtK(itemPopWeights, resultsForTestInstancesDict,
                                          testInstances, predictedRanks, 1)
  println("")
  println("Popularity-weighted precision@1 for all baskets: $popWeightedPrecisionAtOnePercent")

  # Compute precision@1 (accuracy) for each basket size
  for basketSize in sort(collect(keys(basketSizeResultsDict)))
    basketSizeResult = basketSizeResultsDict[basketSize]
    hits = length(findall(x -> x < 1, basketSizeResult.predictedRanks))
    precisionAtOnePercent = hits / basketSizeResult.numBasketInstances * 100
    println("For basket size $basketSize, precision@1: $precisionAtOnePercent")
  end

  # Compute overall precision@5
  hits = length(findall(x -> x < 5, predictedRanks))
  precisionAtFivePercent = hits / numTestInstances * 100
  println("")
  println("precision@5 for all baskets: $precisionAtFivePercent")

  # Compute overall popularity-weighted precision@5
  popWeightedPrecisionAtFivePercent =
    computePopularityWeightedPrecisionAtK(itemPopWeights, resultsForTestInstancesDict,
                                          testInstances, predictedRanks, 5)
  println("")
  println("Popularity-weighted precision@5 for all baskets: $popWeightedPrecisionAtFivePercent")

  # Compute precision@5 for each basket size
  for basketSize in sort(collect(keys(basketSizeResultsDict)))
    basketSizeResult = basketSizeResultsDict[basketSize]
    hits = length(findall(x -> x < 5, basketSizeResult.predictedRanks))
    precisionAtFivePercent = hits / basketSizeResult.numBasketInstances * 100
    println("For basket size $basketSize, precision@5: $precisionAtFivePercent")
  end

  # Compute overall precision@10
  hits = length(findall(x -> x < 10, predictedRanks))
  precisionAtTenPercent = hits / numTestInstances * 100
  println("")
  println("precision@10 for all baskets: $precisionAtTenPercent")

  # Compute overall popularity-weighted precision@10
  popWeightedPrecisionAtTenPercent =
    computePopularityWeightedPrecisionAtK(itemPopWeights, resultsForTestInstancesDict,
                                          testInstances, predictedRanks, 10)
  println("")
  println("Popularity-weighted precision@10 for all baskets: $popWeightedPrecisionAtTenPercent")

  # Compute precision@10 for each basket size
  for basketSize in sort(collect(keys(basketSizeResultsDict)))
    basketSizeResult = basketSizeResultsDict[basketSize]
    hits = length(findall(x -> x < 10, basketSizeResult.predictedRanks))
    precisionAtTenPercent = hits / basketSizeResult.numBasketInstances * 100
    println("For basket size $basketSize, precision@10: $precisionAtTenPercent")
  end

  # Compute overall precision@20
  hits = length(findall(x -> x < 20, predictedRanks))
  precisionAtTwentyPercent = hits / numTestInstances * 100
  println("")
  println("precision@20 for all baskets: $precisionAtTwentyPercent")

  # Compute overall popularity-weighted precision@20
  popWeightedPrecisionAtTwentyPercent =
    computePopularityWeightedPrecisionAtK(itemPopWeights, resultsForTestInstancesDict,
                                          testInstances, predictedRanks, 20)
  println("")
  println("Popularity-weighted precision@20 for all baskets: $popWeightedPrecisionAtTwentyPercent")

  # Compute precision@20 for each basket size
  for basketSize in sort(collect(keys(basketSizeResultsDict)))
    basketSizeResult = basketSizeResultsDict[basketSize]
    hits = length(findall(x -> x < 20, basketSizeResult.predictedRanks))
    precisionAtTwentyPercent = hits / basketSizeResult.numBasketInstances * 100
    println("For basket size $basketSize, precision@20: $precisionAtTwentyPercent")
  end
end

# Computes the popularity weighted (stratified) precision@k
function computePopularityWeightedPrecisionAtK(itemPopWeights, resultsForTestInstancesDict,
                                               testInstances, predictedRanks, k)
  popWeightedHits = 0
  numTestInstances = length(testInstances)
  sumPopWeightsForTestRatings = 0

  startingHeldOutItemId = 0

  for i = 1:numTestInstances
    testInstance = testInstances[i]
    heldOutTestItemId = resultsForTestInstancesDict[testInstance].actualNextItemId - startingHeldOutItemId

    sumPopWeightsForTestRatings += itemPopWeights[heldOutTestItemId]

    if predictedRanks[i] < k
      popWeightedHits += itemPopWeights[heldOutTestItemId]
    end
  end

  popWeightedPrecisionAtKPercent = popWeightedHits / sumPopWeightsForTestRatings * 100
  return popWeightedPrecisionAtKPercent
end


text:
@testset "Viewing" begin
  @testset "Domain" begin
    g = CartesianGrid{T}(10,10)
    b = Box(P2(1,1), P2(5,5))
    v = view(g, b)
    @test v == CartesianGrid(P2(1,1), P2(5,5), dims=(4,4))

    p = PointSet(collect(vertices(g)))
    v = view(p, b)
    @test centroid(v, 1) == P2(1,1)
    @test centroid(v, nelements(v)) == P2(5,5)
  end

  @testset "Data" begin
    dummydata(domain, table) = DummyData(domain, Dict(paramdim(domain) => table))
    dummymeta(domain, table) = meshdata(domain, Dict(paramdim(domain) => table))

    for dummy in [dummydata, dummymeta]
      g = CartesianGrid{T}(10,10)
      t = (a=1:100, b=1:100)
      d = dummy(g, t)
      b = Box(P2(1,1), P2(5,5))
      v = view(d, b)
      @test domain(v) == CartesianGrid(P2(1,1), P2(5,5), dims=(4,4))
      @test Tables.columntable(values(v)) == (a=[12,13,14,15,22,23,24,25,32,33,34,35,42,43,44,45],
                                              b=[12,13,14,15,22,23,24,25,32,33,34,35,42,43,44,45])

      p = PointSet(collect(vertices(g)))
      d = dummy(p, t)
      v = view(d, b)
      dd = domain(v)
      @test centroid(dd, 1) == P2(1,1)
      @test centroid(dd, nelements(dd)) == P2(5,5)
      tt = Tables.columntable(values(v))
      @test tt == (a=[13,14,15,16,17,24,25,26,27,28,35,36,37,38,39,46,47,48,49,50,57,58,59,60,61],
                  b=[13,14,15,16,17,24,25,26,27,28,35,36,37,38,39,46,47,48,49,50,57,58,59,60,61])

      g = CartesianGrid{T}(250,250)
      t = (a=rand(250*250), b=rand(250*250))
      d = dummy(g, t)
      s1 = slice(d, T(50.5):T(100.2), T(41.7):T(81.3))
      pts1 = [centroid(s1, i) for i in 1:nelements(s1)]
      X1 = reduce(hcat, coordinates.(pts1))
      @test all(T[50.5,41.7] .≤ minimum(X1, dims=2))
      @test all(maximum(X1, dims=2) .≤ T[100.2,81.3])

      p = sample(d, 100)
      s2 = slice(p, T(50.5):T(150.7), T(175.2):T(250.3))
      pts2 = [centroid(s2, i) for i in 1:nelements(s2)]
      X2 = reduce(hcat, coordinates.(pts2))
      @test all(T[50.5,175.2] .≤ minimum(X2, dims=2))
      @test all(maximum(X2, dims=2) .≤ T[150.7,250.3])
    end
  end
end


text:
using Test
using EquationEditor
const EE = EquationEditor
const server, task = serveEE()
@test !isready(server.out)
@test !istaskdone(task)
close(server)
#= For debugging
server
task
put!(server.in, "die!")
if isready(server.out)
    println(String(take!(server.out)))
end
if isready(server.out)
    println(String(take!(server.out)))
end
=#

max_stars_repo_name: hustf/EquationEditor

max_issues_repo_name: hustf/EquationEditor

max_forks_repo_name: hustf/EquationEditor


text:
using StatsBase, Plots; pyplot()

names = ["Mary","Mel","David","John","Kayley","Anderson"]
randomName() = rand(names)
X = 3:8
N = 10^6
sampleLengths = [length(randomName()) for _ in 1:N]

bar(X,counts(sampleLengths)/N, ylims=(0,0.35),
    xlabel="Name length", ylabel="Estimated p(x)", legend=:none)

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook


text:
module FilteredMatrices

using LinearAlgebra, LinearMaps, UnPack

export delta, order_estimate, Kernels

"""
    delta(H::AbstractMatrix{<:Number}, ε::Number, range; order = 50, kernel = missing)

Construct a `LinearMap` that acts like matrix `δ_K(ε-H)`, where `δ_K` is a Chebyshev
expansion of a Delta function to `order = K`. See `order_estimate` for estimators based on
the desired energy resolution. A `range = (εmin, εmax)` of the spectrum of `H` is needed to
ensure convergence. The `kernel` is a function `k(n, order)`. Built-in kernels are:
    - missing                 : no kernel
    - Kernels.Jackson         : Jackson kernel
    - Kernels.Fejer           : Fejer kernel
    - Kernels.Lanczos(M = 3)  : Lanczos kernel `sinc(n/order)^M` with integer `M`
    - Kernels.Lorentz(λ = 4)  : Lorentz kernel `sinh(λ(1-n/order))/sinh(λ)` with real `λ`
"""
function delta(H::AbstractMatrix{T}, ε::Number, range; order = 50, kernel = (n, order) -> 1) where {T<:Number}
    M, N = size(H)
    M == N || throw(ArgumentError("Only square matrices are supported"))
    bracket = bandbracket(range)
    v = Vector{T}(undef, N)
    v´ = Vector{T}(undef, N)
    args = (H = H, ε = ε, bracket = bracket, order = order, v = v, v´ = v´, kernel = kernel)
    LinearMap{T}((vdst, v0) -> delta_mul!(vdst, v0, args), N, N; ismutating = true, ishermitian = true)
end

bandbracket((εmin, εmax)) = (εmax + εmin)/2, abs(εmax - εmin)

"""
    order_estimate(ε, Δε, range)

Estimate the Chebyshev expansion order to achieve a resolution `Δε` at energy `ε` within a
`range = (εmin, εmax)`. In terms of normalized energies `σ` and `Δσ` in a `(-1,1)` window,
the required order must be greater than 4π*sqrt(1-σ^2)/Δσ. The prefactor `4` is approximate
and can depend on the kernel used and the definition of resolution.
"""
function order_estimate(ε, Δε, range)
    (center, halfwidth) = bandbracket(range)
    σ = (ε - center)/halfwidth
    abs(σ) < 1 || throw(ArgumentError("Energy outside band range"))
    Δσ = Δε/halfwidth
    return ceil(Int, 4π*sqrt(1-σ^2)/Δσ)
end

function delta_mul!(vdst::AbstractVector{T}, v0, args) where {T}
    @unpack H, ε, bracket, order, v, v´, kernel = args

    # We use adjoint of CSC matrix to accelerate multiplication
    H´ = H'
    center, halfwidth = bracket
    α = 2/halfwidth
    β = 2center/halfwidth
    σ = (ε - center)/halfwidth
    ρ = 1/(pi*sqrt(1-σ^2))

    # t and t´ are used to build tn = ρ Tn(σ) recurseively
    t  = ρ
    t´ = σ * t

    # v and v´ are used to build vn = Tn(h)⋅v recursively, where h = (H-center)/halfwidth
    @. v = v0
    fill!(v´, zero(T))
    nextCheby!(v´, H´, v, α, β)

    # vdst accumulates delta(σ-H) * v ≈ ρ (v + 2∑ Tn(σ)Tn(h)⋅v) = t0 v0 + 2 ∑ tn vn
    @. vdst = kernel(0, order) * t * v + kernel(1, order) * 2t´ * v´

    for n in 2:order
        nextCheby!(v, H´, v´, α, β)
        t = nextCheby(t, σ, t´)
        @. vdst += kernel(n, order) * 2t * v
        v´, v = v, v´
        t´, t = t, t´
    end
    return vdst
end

function nextCheby!(v´, H´, v, α, β)
    mul!(v´, H´, v, α, -1)
    @. v´ -= β * v
    return v´
end

nextCheby(t´, σ, t) = 2σ * t - t´

module Kernels

Jackson(n, order) =
    ((order - n + 1) * cos(π * n / (order + 1)) + sin(π * n / (order + 1)) * cot(π / (order + 1))) / (order + 1)

Fejer(n, order) = 1 - n/order

Lorentz(λ::Real = 4) = (n, order) -> sinh(λ*(1-n/order))/sinh(λ)

Lanczos(M::Integer = 3) = (n, order) -> sinc(n/order)^M

end # module

end # module


text:
"""
    LRFinderResult(lrs, losses)

Result of the learning rate finder [`lrfind`](#). Use `plot`
to visualize.
"""
struct LRFinderResult
    losses
    lrs
    estimators
    estimates
end

"""
    lrfind(learner[, dataiter; kwargs...]) -> LRFinderResult

Run the learning rate finder. Exponentially increases the learning rate from a very
low value to a very high value and uses the losses to estimate an optimal learning
rate. Return a [`LRFinderResult`](#).

## Keyword arguments

- `nsteps = 100`: maximum number of steps to run the learning rate finder for
- `startlr = 1e-7`: minimum learning rate
- `endlr = 10`: maximum learning rate
- `divergefactor`: stop finder early if loss goes higher than lowest loss times
    this factor
- `estimators = [Steepest(), MinDivByTen()]`: list of [`LREstimator`](#)s
"""
function lrfind(
        learner,
        dataiter = learner.data.training;
        nsteps = 100,
        startlr = 1e-7,
        endlr = 10,
        divergefactor = 4,
        estimators = [Steepest(), MinDivByTen()])
    losses = Float64[]
    lrs = Float64[]
    bestloss = Inf
    scheduler = FluxTraining.removecallback!(learner, Scheduler)  # remove current `Scheduler` so it does not interfere
    modelcheckpoint = deepcopy(cpu(learner.model))

    withfields(
        learner,
        model = modelcheckpoint,
        params = Flux.params(modelcheckpoint),
        optimizer = deepcopy(learner.optimizer)
        ) do

        FluxTraining.runepoch(learner, TrainingPhase()) do _
            for (i, batch) in zip(1:nsteps, dataiter)
                lr = startlr * (endlr / startlr) ^ (i / nsteps)
                learner.optimizer.eta = lr

                state = step!(learner, TrainingPhase(), batch)

                push!(losses, state.loss)
                push!(lrs, learner.optimizer.eta)
                bestloss = min(state.loss, bestloss)
                if state.loss > bestloss * divergefactor
                    throw(FluxTraining.CancelEpochException("Learning rate finder diverged."))
                end
            end
        end
    end

    # Restore previous `Learner` state
    isnothing(scheduler) || FluxTraining.addcallback!(learner, scheduler)
    FluxTraining.model!(learner, modelcheckpoint)
    return LRFinderResult(losses, lrs, estimators)
end

# `LREstimator`s allow giving back different suggestions for a learning rate.

"""
    abstract type LREstimator

Estimator for an optimal learning rate. Needs to implement [`estimatelr`](#).

See [`Steepest`](#) and [`MinDivByTen`](#).
"""
abstract type LREstimator end

"""
    estimatelr(::LREstimator, losses, lrs)

Estimate the optimal learning rate using `losses` and `lrs`.
"""
function estimatelr end

"""
    Steepest <: LREstimator

Estimate the optimal learning rate to be where the gradient of the loss
is the steepest, i.e. the decrease is largest.
"""
struct Steepest <: LREstimator
    beta
end
Steepest() = Steepest(0.98)

function estimatelr(est::Steepest, losses, lrs)
    slosses = smoothvalues(losses, est.beta)
    grads = (slosses[2:end] .- slosses[1:end-1]) ./ log.(lrs[2:end] .- lrs[1:end-1])
    i = length(lrs) ÷ 3
    lr = lrs[i:end][argmax(grads[i:end])]
    return lr
end

"""
    MinDivByTen <: LREstimator

Estimate the optimal learning rate to be value at the minimum loss divided by 10.
"""
struct MinDivByTen <: LREstimator
    beta
end
MinDivByTen() = MinDivByTen(0.98)

function estimatelr(est::MinDivByTen, losses, lrs)
    i = length(lrs) ÷ 3
    lr = lrs[i:end][argmin(smoothvalues(losses, est.beta)[i:end])] / 10
    return lr
end

LRFinderResult(losses, lrs) = LRFinderResult(losses, lrs, [Steepest(), MinDivByTen()])
LRFinderResult(losses, lrs, estimators) = LRFinderResult(
    losses,
    lrs,
    estimators,
    [estimatelr(est, losses, lrs) for est in estimators]
)

# Printing and plotting

function Base.show(io::IO, result::LRFinderResult)
    lrfindtextplot!(io, result)

    names = [typeof(est).name.name for est in result.estimators]
    println(io)
    println(io)
    pretty_table(
        io,
        hcat(names, result.estimates),
        header=["Estimator", "Suggestion"], tf=tf_borderless,
        alignment=[:r, :l])
end

lrfindtextplot(result::LRFinderResult) = lrfindtextplot!(stdout, result)
function lrfindtextplot!(io, result::LRFinderResult)
    p = UnicodePlots.lineplot(
        result.lrs, result.losses,
        height=10, xscale=:log10, width=displaysize(io)[2]-15)
    UnicodePlots.title!(p, "Learning rate finder result")

    for (i, estimate) in enumerate(result.estimates)
        UnicodePlots.lines!(p, estimate, maximum(result.losses), estimate, minimum(result.losses), :red)
        UnicodePlots.annotate!(p, estimate, maximum(result.losses), string(round(estimate; sigdigits=4)))
    end
    show(io, p)
    return p
end

InlineTest.@testset "LRFinderResult" begin
    res = LRFinderResult(100.:-1:1, 1:100.)
    InlineTest.@test_nowarn show(Base.DevNull(), res)
end

# Utilities

"""
    smoothvalues(xs, β)

Apply exponential smoothing with parameter `β` to vector `xs`.
"""
function smoothvalues(xs, β)
    res = similar(xs)
    val = 0
    for (i, x) in enumerate(xs)
        val = (val * β) + (x * (1 - β))
        res[i] = val/((1)-β^i)
    end
    res
end


text:
# ------------------------------------------------------------------------------------------
# # Julia's Type System
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# Julia programs manipulate *values*, and every value has two parts: a *type* part and a
# data part. The type part answers the question "what kind of thing is this?", and the data
# part distinguishes one thing of a certain kind from every other thing of that kind.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ## Part 1. DataType
# ------------------------------------------------------------------------------------------

typeof(3)

# ------------------------------------------------------------------------------------------
# In this case the type is `Int64` and the data part is the bits `...0011`.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# In Julia types are also values:
# ------------------------------------------------------------------------------------------

typeof(Int64)

typeof(DataType)

# ------------------------------------------------------------------------------------------
# In fact, the identity `typeof(typeof(x)) === DataType` holds for all values in Julia.
# `DataType` is the backbone of the entire system. It does many jobs, which can be
# identified by looking inside a `DataType` object:
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ### DataType Job 1: A symbolic description
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# This consists of a name (which is mostly a string), and a vector of sub-components:
# ------------------------------------------------------------------------------------------

T = typeof(1+2im)

T.name

T.parameters

# ------------------------------------------------------------------------------------------
# ### DataType Job 2: A nominal hierarchy of types
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# DataTypes form a tree of declared type relationships ("an x is-a y"):
# ------------------------------------------------------------------------------------------

T.super

T.super.super.super.super  # `Any` is the built-in top of the hierarchy.

# ------------------------------------------------------------------------------------------
# ### DataType Job 3: Describe the representation
# ------------------------------------------------------------------------------------------

T.types

T.name.names

T.size

T.mutable   # whether this was declared with `type` (vs. `immutable`)

T.abstract  # whether this was declared with `abstract`

T.ninitialized

T.layout

# ------------------------------------------------------------------------------------------
# ### Defining struct types
#
# ------------------------------------------------------------------------------------------

struct Point
    x::Float64
    y::Float64
end

Point(1,2)

# ------------------------------------------------------------------------------------------
# ### Abstract vs. Concrete
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# `abstract` types can have declared subtypes, while concrete types can have instances.
# These are separated because if an `X` IS-A `Y`, and `Y` specifies a representation, then
# `X` had better have the same representation.
#
# "car is-a vehicle" is correct because "vehicle" is an abstract concept that doesn't commit
# to any specifics. But if I tell you I'm giving you a Porsche, it needs to look like a
# Porsche.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# A type `T` is concrete if there could be some value `x` such that `typeof(x) === T`. This
# is also sometimes called a "leaf type".
# ------------------------------------------------------------------------------------------

abstract type PointLike end

# struct Point <: PointLike

# ------------------------------------------------------------------------------------------
# ## Part 2. Type parameters
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# Type parameters can be completely or partially specified:
# ------------------------------------------------------------------------------------------

Array{Int}

[1] isa Array

Array{Int,2}

# ------------------------------------------------------------------------------------------
# A type is concrete (can have instances) if
#     1. it is not declared `abstract`
#     2. all parameters are specified
# ------------------------------------------------------------------------------------------

[1] isa Array{Int,1}

[1] isa Array{Int}

[1] isa Array{Number}

Int <: Number

# ------------------------------------------------------------------------------------------
# Types with different *specified* parameters are just different, and have no subtype
# relationship. This is called *invariance*.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ### Defining types with parameters
# ------------------------------------------------------------------------------------------

struct GenericPoint{T<:Real}
    x::T
    y::T
end

GenericPoint(1,2)

GenericPoint(1.0,2.0)

GenericPoint(1,2.0)

# ------------------------------------------------------------------------------------------
# ### Tuple types
# ------------------------------------------------------------------------------------------

typeof((1,2.0))

# ------------------------------------------------------------------------------------------
# Very similar to other DataTypes, except
#     1. Have no field names, only indices
#     2. `T.parameters == T.types`
#     3. Are always immutable
#     4. Can have any number of fields
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# These factors conspire to make Tuples the only *covariant* types in Julia:
# ------------------------------------------------------------------------------------------

Tuple{Int} <: Tuple{Number}

# ------------------------------------------------------------------------------------------
# A Tuple type is concrete iff all its field types are.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# Tuple types can be abstract with respect to the number of elements. These are called
# variadic tuple types, or vararg types.
# ------------------------------------------------------------------------------------------

Tuple{Int,Vararg{Int}}

# ------------------------------------------------------------------------------------------
# Note that `Vararg` refers to the tail of a tuple type, and as such is not a first-class
# type itself. It only makes sense inside a Tuple type. This is a bit unfortunate.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# The second parameter to `Vararg` is a length, which can also be either unspecified (as
# above), or specified:
# ------------------------------------------------------------------------------------------

Tuple{Int,Vararg{Int,2}}

# ------------------------------------------------------------------------------------------
# ## Part 3. Larger type domains
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ### Union types
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# A type can be thought of as a set of possible values. A type expresses *uncertainty* about
# which value we have. You can do set operations on them.
# ------------------------------------------------------------------------------------------

Union{Int64,Float64}

1 isa Union{Int64,Float64}

Int64 <: Number

Int64 <: Union{Int64,Float64}

Union{Int,String} <: Union{Int,String,Float32}

typeintersect(Union{Int,String}, Union{Int,String,Float32})

# ------------------------------------------------------------------------------------------
# Union types naturally lend themselves to missing data.
# ------------------------------------------------------------------------------------------

data = [1.1, missing, 3.2, missing, 5.7, 0.4]

# ------------------------------------------------------------------------------------------
# ### UnionAll types
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# This is an *iterated union* of types.
#
# `Array{T,1} where T<:Integer`
#
# Means "the union of all types of the form Array{T,1} where T is a subtype of Integer".
#
# This expresses uncertainty about the value of a parameter.
#
# This concept exists in all versions of Julia, but does not have syntax or fully correct
# support within the system until upcoming v0.6.0 (currently on branch jb/subtype).
#
# * Since this kind of type introduces *variables*, its expressive power is (probably)
# equivalent to quantified boolean formulas.
# * Requires a quantified-SAT solver in the compiler.
# * Under common assumptions, harder than NP-complete.
# ------------------------------------------------------------------------------------------

# this definition is in the Base library
Vector = Array{T,1} where T

# ------------------------------------------------------------------------------------------
# These are used to express "unspecified parameters".
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# These also describe methods with "method parameters" (or "static parameters"):
# ------------------------------------------------------------------------------------------

func(a::Array{T,1}) where {T<:Integer} = T

func([0x00])

func([1.0])

# ------------------------------------------------------------------------------------------
# #### Question
#
# What is the difference between
#
# `Vector{Vector{T}} where T`
#
# and
#
# `Vector{Vector{T} where T}`?
#
# Is one a subtype of the other?
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ### Exercise
#
# Define a `UnitPoint{<:Real}` parametric struct type which has `x` and `y` fields of type
# `T` and which has an inner constructor that normalizes its arguments by diving them by
# `hypot(x, y)` upon construction, guaranteeing that the resulting point is on the unit
# circle.
# ------------------------------------------------------------------------------------------

max_stars_repo_name: grenkoca/JuliaTutorials

max_issues_repo_name: grenkoca/JuliaTutorials

max_forks_repo_name: grenkoca/JuliaTutorials


text:
module TSMLTransformers

using MLDataUtils
using Dates
using DataFrames
using Statistics
using Random
using CSV

export fit!,transform!

export Transformer,TSLearner
export Imputer,Pipeline,SKLLearner,OneHotEncoder,Wrapper

export Matrifier,Dateifier
export DateValizer,DateValgator,DateValNNer
export CSVDateValReader, CSVDateValWriter

using TSML.TSMLTypes
import TSML.TSMLTypes.fit! # to overload
import TSML.TSMLTypes.transform! # to overload
using TSML.Utils

include("valdate.jl")

# Transforms instances with nominal features into one-hot form
# and coerces the instance matrix to be of element type Float64.
mutable struct OneHotEncoder <: Transformer
  model
  args

  function OneHotEncoder(args=Dict())
    default_args = Dict(
      # Nominal columns
      :nominal_columns => nothing,
      # Nominal column values map. Key is column index, value is list of
      # possible values for that column.
      :nominal_column_values_map => nothing
    )
    new(nothing, mergedict(default_args, args))
  end
end

function fit!(ohe::OneHotEncoder, features::T, labels::Vector) where {T<:Union{Vector,Matrix,DataFrame}}
  instances=convert(Matrix,features)
  # Obtain nominal columns
  nominal_columns = ohe.args[:nominal_columns]
  if nominal_columns == nothing
    nominal_columns = find_nominal_columns(instances)
  end

  # Obtain unique values for each nominal column
  nominal_column_values_map = ohe.args[:nominal_column_values_map]
  if nominal_column_values_map == nothing
    nominal_column_values_map = Dict{Int, Any}()
    for column in nominal_columns
      nominal_column_values_map[column] = unique(instances[:, column])
    end
  end

  # Create model
  ohe.model = Dict(
    :nominal_columns => nominal_columns,
    :nominal_column_values_map => nominal_column_values_map
  )
end

function transform!(ohe::OneHotEncoder, features::T) where {T<:Union{Vector,Matrix,DataFrame}}
  instances=convert(Matrix,features)
  nominal_columns = ohe.model[:nominal_columns]
  nominal_column_values_map = ohe.model[:nominal_column_values_map]

  # Create new transformed instance matrix of type Float64
  num_rows = size(instances, 1)
  num_columns = (size(instances, 2) - length(nominal_columns))
  if !isempty(nominal_column_values_map)
    num_columns += sum(map(x -> length(x), values(nominal_column_values_map)))
  end
  transformed_instances = zeros(Float64, num_rows, num_columns)

  # Fill transformed instance matrix
  col_start_index = 1
  for column in 1:size(instances, 2)
    if !in(column, nominal_columns)
      transformed_instances[:, col_start_index] = instances[:, column]
      col_start_index += 1
    else
      col_values = nominal_column_values_map[column]
      for row in 1:size(instances, 1)
        entry_value = instances[row, column]
        entry_value_index = findfirst(isequal(entry_value),col_values)
        if entry_value_index == 0
          warn("Unseen value found in OneHotEncoder,
                for entry ($row, $column) = $(entry_value).
                Patching value to $(col_values[1]).")
          entry_value_index = 1
        end
        entry_column = (col_start_index - 1) + entry_value_index
        transformed_instances[row, entry_column] = 1
      end
      col_start_index += length(nominal_column_values_map[column])
    end
  end

  return transformed_instances
end

# Finds all nominal columns.
#
# Nominal columns are those that do not have Real type nor
# do all their elements correspond to Real.
function find_nominal_columns(features::T) where {T<:Union{Vector,Matrix,DataFrame}}
  instances=convert(Matrix,features)
  nominal_columns = Int[]
  for column in 1:size(instances, 2)
    col_eltype = infer_eltype(instances[:, column])
    if !<:(col_eltype, Real)
      push!(nominal_columns, column)
    end
  end
  return nominal_columns
end

# Imputes NaN values from Float64 features.
mutable struct Imputer <: Transformer
  model
  args

  function Imputer(args=Dict())
    default_args = Dict(
      # Imputation strategy.
      # Statistic that takes a vector such as mean or median.
      :strategy => mean
    )
    new(nothing, mergedict(default_args, args))
  end
end

function fit!(imp::Imputer, instances::T, labels::Vector) where {T<:Union{Vector,Matrix,DataFrame}}
  imp.model = imp.args
end

function transform!(imp::Imputer, features::T)  where {T<:Union{Vector,Matrix,DataFrame}}
  instances=convert(Matrix,features)
  new_instances = copy(instances)
  strategy = imp.model[:strategy]

  for column in 1:size(instances, 2)
    column_values = instances[:, column]
    col_eltype = infer_eltype(column_values)

    if <:(col_eltype, Real)
      na_rows = map(x -> isnan(x), column_values)
      if any(na_rows)
        fill_value = strategy(column_values[.!na_rows])
        new_instances[na_rows, column] .= fill_value
      end
    end
  end

  return new_instances
end

# Chains multiple transformers in sequence.
mutable struct Pipeline <: Transformer
  model
  args

  function Pipeline(args=Dict())
    default_args = Dict(
      # Transformers as list to chain in sequence.
      :transformers => [OneHotEncoder(), Imputer()],
      # Transformer args as list applied to same index transformer.
      :transformer_args => nothing
    )
    new(nothing, mergedict(default_args, args))
  end
end

function fit!(pipe::Pipeline, features::T=[], labels::Vector=[]) where {T<:Union{Vector,Matrix,DataFrame}}
  #instances=convert(Matrix,features)
  instances=copy(features)
  transformers = pipe.args[:transformers]
  transformer_args = pipe.args[:transformer_args]

  current_instances = instances
  new_transformers = Transformer[]
  for t_index in 1:length(transformers)
    transformer = createtransformer(transformers[t_index], transformer_args)
    push!(new_transformers, transformer)
    fit!(transformer, current_instances, labels)
    current_instances = transform!(transformer, current_instances)
  end

  pipe.model = Dict(
      :transformers => new_transformers,
      :transformer_args => transformer_args
  )
end

function transform!(pipe::Pipeline, features::T=[]) where {T<:Union{Vector,Matrix,DataFrame}}
  #instances = convert(Matrix,features)
  instances = copy(features)
  transformers = pipe.model[:transformers]

  current_instances = instances
  for t_index in 1:length(transformers)
    transformer = transformers[t_index]
    current_instances = transform!(transformer, current_instances)
  end

  return current_instances
end

# Wraps around an CombineML transformer.
mutable struct Wrapper <: Transformer
  model
  args

  function Wrapper(args=Dict())
    default_args = Dict(
      # Transformer to call.
      :transformer => OneHotEncoder(),
      # Transformer args.
      :transformer_args => nothing
    )
    new(nothing, mergedict(default_args, args))
  end
end

function fit!(wrapper::Wrapper, features::T, labels::Vector) where {T<:Union{Vector,Matrix,DataFrame}}
  instances=convert(Matrix,features)
  transformer_args = wrapper.args[:transformer_args]
  transformer = createtransformer(
    wrapper.args[:transformer],
    transformer_args
  )

  if transformer_args != nothing
    transformer_args = mergedict(transformer.args, transformer_args)
  end
  fit!(transformer, instances, labels)

  wrapper.model = Dict(
    :transformer => transformer,
    :transformer_args => transformer_args
  )
end

function transform!(wrapper::Wrapper, instances::T) where {T<:Union{Vector,Matrix,DataFrame}}
  transformer = wrapper.model[:transformer]
  return transform!(transformer, instances)
end

# Create transformer
#
# @param prototype Prototype transformer to base new transformer on.
# @param options Additional options to override prototype's options.
# @return New transformer.
function createtransformer(prototype::Transformer, args=nothing)
  new_args = copy(prototype.args)
  if args != nothing
    new_args = mergedict(new_args, args)
  end

  prototype_type = typeof(prototype)
  return prototype_type(new_args)
end

end


text:
@testset "Independent Bounds" begin
    m = BabyPOMDP()
    planner = solve(AdaOPSSolver(), m)
    belief = WPFBelief([true], [1.0], 1, 1, AdaOPSTree(planner, initialstate(m)), true)
    b = IndependentBounds(0.0, -1e-5)
    @test bounds(b, m, belief, 90, true) == (0.0, -1e-5)
    b = IndependentBounds(0.0, -1e-5, consistency_fix_thresh=1e-5)
    @test bounds(b, m, belief, 90, true) == (0.0, 0.0)
    b = IndependentBounds(0.0, -1e-4, consistency_fix_thresh=1e-5)
    @test bounds(b, m, belief, 90, true) == (0.0, -1e-4)
end


text:
using Pkg
Pkg.activate(joinpath(@__DIR__, "../.."))
Pkg.instantiate()

# ## visualizer
vis = Visualizer()
open(vis)

# ## setup
using Dojo
using IterativeLQR
using LinearAlgebra
using FiniteDiff
using DojoEnvironments
using JLD2
using Plots

# ## scripts
include(joinpath(module_dir(), "examples/policy/methods/continuation.jl"))
include(joinpath(module_dir(), "examples/policy/methods/tvlqr.jl"))
include(joinpath(module_dir(), "DojoEnvironments/src",
    "quadruped/methods/template.jl"))

################################################################################
# ## system
################################################################################
gravity = -9.81
timestep = 0.02
friction_coefficient = 0.8
damper = 0.5
spring = 1.0
env = get_environment(:quadruped,
    representation=:minimal,
    timestep=timestep,
    contact_body=false,
    gravity=gravity,
    friction_coefficient=friction_coefficient,
    damper=damper,
    spring=spring,
    infeasible_control=true,
    vis=vis)

# ## dimensions
nx = env.num_states
nu = env.num_inputs
nu_infeasible = 6

################################################################################
# Load trajectory
################################################################################
linear_velocities = [0, 0.02, 0.04, 0.07, 0.10, 0.13, 0.16, 0.19, 0.22, 0.25]
rotational_velocities = [0, 0.25]
unfiltered_velocities = []
for v in linear_velocities
    for ω in rotational_velocities
        push!(unfiltered_velocities, [v, ω])
    end
end
unfiltered_velocities

x_lib = Dict()
u_lib = Dict()

for velocity in unfiltered_velocities
    v = velocity[1]
    ω = velocity[2]
    try file = JLD2.jldopen(joinpath(@__DIR__, "../data/planar_v_" *
        string(round(v, digits=2)) * "_ω_" *
        string(round(ω, digits=2)) *
        ".jld2"))
        x_lib[velocity] = file["x"]
        u_lib[velocity] = file["u"]
        JLD2.close(file)
    catch e
    end
end
velocities = collect(keys(x_lib))
T = length(first(x_lib)[2])

################################################################################
# TVLQR
################################################################################
K_lib = Dict()
q_tracking=[0.3; 0.05; 0.05;
    5e-1 * ones(3);
    1e-6 * ones(3);
    1e-6 * ones(3);
    fill([2, 1e-6], 12)...]
r_tracking = env.mechanism.timestep * 30 * ones(nu)
for velocity in keys(x_lib)
    K_lib[velocity] = tvlqr(x_lib[velocity], u_lib[velocity], env;
        q_tracking=q_tracking, r_tracking=r_tracking)[1]
end

################################################################################
# Target trajectory
################################################################################
# divide timestep
S = 3
v_target = [0.20 * (1 - cos(2π * t/7))/2 for t=0:timestep/S:30]
ω_target = [0.20 * (1 - cos(2π * t/7))/2 for t=0:timestep/S:30]
velocity_target = [[v_target[i], ω_target[i]] for i=1:length(v_target)]
plot(0:timestep/S:30, v_target)

################################################################################
# Simulate accurately the TVLQR policy
################################################################################
mech_sim = get_mechanism(:quadruped,
    contact_body=true,
    timestep=timestep/S,
    gravity=gravity,
    friction_coefficient=friction_coefficient,
    damper=damper,
    spring=spring)

function find_closest(velocities, v_target)
    Δ = norm.([v - v_target for v in velocities], Inf)
    velocity = velocities[findmin(Δ)[2]]
    return velocity
end

# controller
function ctrl!(mechanism, k; velocities=velocities, x_lib=x_lib, u_lib=u_lib, K_lib=K_lib)
    @show k
    # choose reference based on target velocity
    velocity = find_closest(velocities, velocity_target[k])
    x_ref = x_lib[velocity]
    u_ref = u_lib[velocity]
    K_ref = K_lib[velocity]

    # find index in the loops
    N = length(u_ref)
    nu = input_dimension(mechanism)
    ind = (Int(floor((k-1)/S))) % N + 1

    # compute policy
    x = get_minimal_state(mechanism)
    Δx = (x_ref[ind] - x)
    Δx[[1,2]] .= 0.0
    u = u_ref[ind] + K_ref[ind] * Δx
    u = clamp.(u, -0.2, +0.2)
    u = [zeros(6); u[7:end]] / S
    set_input!(mechanism, SVector{nu}(u))
end

controller!(m, k) = ctrl!(m, k; velocities=velocities, x_lib=x_lib, u_lib=u_lib, K_lib=K_lib)

x1_roll = deepcopy(x_lib[[0, 0]][1])
x1_roll[3] += 0.0
x1_roll[4] += 0.0
set_minimal_state!(mech_sim, x1_roll)
Main.@elapsed storage = simulate!(mech_sim, 25.0, controller!,
    record=true,
    opts=SolverOptions(rtol=1e-5, btol=1e-4, undercut=5.0),
    )

vis, anim = Dojo.visualize(mech_sim, storage, vis=vis, build=true, name=:robot)
settransform!(vis[:tuned], Translation(0,0.5,0))
convert_frames_to_video_and_gif("quadruped_motion_matching_planar")


text:
using GeometricIntegrators
using ChargedParticleDynamics.GuidingCenter4d.TokamakFastSurface
using GeometricExamples

const Δt     = 10.
const ntime  = 2500
const nx     = 200
const ny     = 200
const nsave  = 5
const nplot  = 5
const run_id = "poincare_2nd_fast_firk_dt10"

include("guiding_center_4d_settings_firk.jl")

tableau_list = get_tableau_list_firk()

pinv = guiding_center_4d_ode_poincare_invariant_2nd(Δt, nx, ny, ntime, nsave)

include("guiding_center_4d_fast_poincare_invariant_2nd.jl")


text:
using CSV
using DataFrames
using Plots
using JLD2
using DiscretisedFluidQueues
using GLM
import Distributions
include("/Users/anguslewis/Documents/SFFMProject/DFQExamples.jl/plot_utils.jl")
include("/Users/anguslewis/Documents/SFFMProject/DFQExamples.jl/error_metrics.jl")

@halfwidth_plot_defaults()

plotlyjs()
function p()
ic = "exp"
ks_data = CSV.read((@__DIR__)*"/../data/kolmogorov_smirnov.csv",DataFrame)
ks_data.DG_lin_lim = CSV.read("/Users/anguslewis/Documents/SFFMProject/DFQExamples.jl/hitting_times_linearDGwLimiter/reflecting_model/data/"*ic*"/data/errors/ks_error.csv",DataFrame).DG_limit

ks_data_lwr = CSV.read((@__DIR__)*"/../data/kolmogorov_smirnov_ci_lwr.csv",DataFrame)
ks_data_lwr.DG_lin_lim = CSV.read("/Users/anguslewis/Documents/SFFMProject/DFQExamples.jl/hitting_times_linearDGwLimiter/reflecting_model/data/"*ic*"/data/errors/ks_error_lwr.csv",DataFrame).DG_limit

ks_data_upr = CSV.read((@__DIR__)*"/../data/kolmogorov_smirnovci_upr.csv",DataFrame)
ks_data_upr.DG_lin_lim = CSV.read("/Users/anguslewis/Documents/SFFMProject/DFQExamples.jl/hitting_times_linearDGwLimiter/reflecting_model/data/"*ic*"/data/errors/ks_error_upr.csv",DataFrame).DG_limit

lwr = [log10(ks_data[i,j]) .- log10(ks_data_lwr[i,j]) for i in 1:size(ks_data,1), j in 1:size(ks_data,2)]
upr = [log10(ks_data_upr[i,j]) .- log10(ks_data[i,j]) for i in 1:size(ks_data,1), j in 1:size(ks_data,2)]
ks_data = log10.(ks_data)
plot()
linestyles_vec = [:solid,:dash,:dashdot, :dot,:dash]
markerstyles_vec = [:false,:false,:false,:false,:false]
_names = ["DG"; "DG limit"; "Unif"; "QBD-RAP"; "DG lin lim"]
colours = [1;3;2;4;5]
for (c,col) in enumerate(names(ks_data))
    os = collect(1:2:21)
    (c==5)&&(os.+=1)
    plot!(log10.(os),ks_data[:,col], ribbon=(lwr[:,c],upr[:,c]),
        label=_names[c],linestyle=linestyles_vec[c],
        marker=markerstyles_vec[c],
        linewidth=2,color=colours[c])
end
plot!(xlabel="Dimension")
plot!(ylabel="Error"); error_ticks!(plot!())
plot!(title="KS error")
# plot!(legend=(0.8,1.0))
plot!()
@add_lines!(ks_data,("Unif","QBDRAP"),@__DIR__)

savefig((@__DIR__)*"/ks_error_formatted.pdf")
end
p()

@halfwidth_plot_defaults()
function p()
ic = "exp"
L1_cell_probs_data = CSV.read((@__DIR__)*"/../data/L1_cell_probs.csv",DataFrame)
L1_cell_probs_data.DG_lin_lim = log10.(CSV.read("/Users/anguslewis/Documents/SFFMProject/DFQExamples.jl/hitting_times_linearDGwLimiter/reflecting_model/data/"*ic*"/data/errors/l1_cell_probs_error.csv",DataFrame).DG_limit)

L1_cell_probs_data_lwr = CSV.read((@__DIR__)*"/../data/L1_cell_probs_ci_lwr.csv",DataFrame)
L1_cell_probs_data_lwr.DG_lin_lim = log10.(CSV.read("/Users/anguslewis/Documents/SFFMProject/DFQExamples.jl/hitting_times_linearDGwLimiter/reflecting_model/data/"*ic*"/data/errors/l1_cell_probs_error_lwr.csv",DataFrame).DG_limit)

L1_cell_probs_data_upr = CSV.read((@__DIR__)*"/../data/L1_cell_probsci_upr.csv",DataFrame)
L1_cell_probs_data_upr.DG_lin_lim = log10.(CSV.read("/Users/anguslewis/Documents/SFFMProject/DFQExamples.jl/hitting_times_linearDGwLimiter/reflecting_model/data/"*ic*"/data/errors/l1_cell_probs_error_upr.csv",DataFrame).DG_limit)

lwr = [(L1_cell_probs_data[i,j]) .- (L1_cell_probs_data_lwr[i,j]) for i in 1:size(L1_cell_probs_data,1), j in 1:size(L1_cell_probs_data,2)]
upr = [(L1_cell_probs_data_upr[i,j]) .- (L1_cell_probs_data[i,j]) for i in 1:size(L1_cell_probs_data,1), j in 1:size(L1_cell_probs_data,2)]
plot()
linestyles_vec = [:solid,:dash,:dashdot,:dot,:dash]
markerstyles_vec = [:false,:false,:false,:false,:false]
_names = ["DG"; "DG limit"; "Unif"; "QBD-RAP"; "DG lin lim"]
colours = [1;3;2;4;5]
for (c,col) in enumerate(names(L1_cell_probs_data))
    os = collect(1:2:21)
    (c==5)&&(os.+=1)
    plot!(log10.(os),L1_cell_probs_data[:,col], ribbon=(lwr[:,c],upr[:,c]),
        label=_names[c],linestyle=linestyles_vec[c],
        marker=markerstyles_vec[c],
        linewidth=2,color=colours[c])
end
plot!(xlabel="Dimension")
plot!(ylabel="Error"); error_ticks!(plot!())
plot!(title="Cell-wise error")
# plot!(legend=(0.8,1.02))
plot!()
@add_lines!(L1_cell_probs_data,("Unif","QBDRAP"),@__DIR__)

savefig((@__DIR__)*"/L1_cell_probs_error_formatted.pdf")
end
p()

# function p()
# ic = "exp"
# l2_cdf_data = CSV.read((@__DIR__)*"/../data/L2.csv",DataFrame)
# l2_cdf_data.DG_lin_lim = CSV.read("/Users/anguslewis/Documents/SFFMProject/DFQExamples.jl/hitting_times_linearDGwLimiter/reflecting_model/data/"*ic*"/data/errors/l1_cell_probs_error.csv",DataFrame).DG_limit

# l2_cdf_data_lwr = CSV.read((@__DIR__)*"/../data/L2_ci_lwr.csv",DataFrame)
# l2_cdf_data_upr = CSV.read((@__DIR__)*"/../data/L2ci_upr.csv",DataFrame)
# markerstyles_vec = [:false,:false,:false,:false]
# _names = ["DG"; "DG limit"; "Unif"; "QBD-RAP"]
# colours = [1;3;2;4]
# lwr = [(l2_cdf_data[i,j]) .- (l2_cdf_data_lwr[i,j]) for i in 1:size(l2_cdf_data,1), j in 1:size(l2_cdf_data,2)]
#  upr = [(l2_cdf_data_upr[i,j]) .- (l2_cdf_data[i,j]) for i in 1:size(l2_cdf_data,1), j in 1:size(l2_cdf_data,2)]
# plot()
# linestyles_vec = [:solid,:dash,:dashdot, :dot]
# for (c,col) in enumerate(names(l2_cdf_data))
#     plot!(log10.(1:2:21),l2_cdf_data[:,col],ribbon=(lwr[:,c],upr[:,c]),
#         label=_names[c],linestyle=linestyles_vec[c],
#         marker=markerstyles_vec[c],
#         linewidth=2,color=colours[c])
# end
# plot!(xlabel="Dimension")
# plot!(ylabel="Error"); error_ticks!(plot!())
# plot!(title="L² error - CDF")
# plot!(legend=(0.8,1.0))
# plot!()
# @add_lines!(l2_cdf_data,("Unif","QBDRAP"),@__DIR__)

# savefig((@__DIR__)*"/l2_cdf_error_formatted.pdf")
# end
# p()

function p()
_names = ["DG"; "DG limit"; "Unif"; "QBD-RAP"; "DG lin lim"]
ic = "exp"
l1_cdf_data = CSV.read((@__DIR__)*"/../data/L1.csv",DataFrame)
l1_cdf_data.DG_lin_lim = log10.(CSV.read("/Users/anguslewis/Documents/SFFMProject/DFQExamples.jl/hitting_times_linearDGwLimiter/reflecting_model/data/"*ic*"/data/errors/l1_error.csv",DataFrame).DG_limit)

l1_cdf_data_lwr = CSV.read((@__DIR__)*"/../data/L1_ci_lwr.csv",DataFrame)
l1_cdf_data_lwr.DG_lin_lim = log10.(CSV.read("/Users/anguslewis/Documents/SFFMProject/DFQExamples.jl/hitting_times_linearDGwLimiter/reflecting_model/data/"*ic*"/data/errors/l1_error_lwr.csv",DataFrame).DG_limit)

l1_cdf_data_upr = CSV.read((@__DIR__)*"/../data/L1ci_upr.csv",DataFrame)
l1_cdf_data_upr.DG_lin_lim = log10.(CSV.read("/Users/anguslewis/Documents/SFFMProject/DFQExamples.jl/hitting_times_linearDGwLimiter/reflecting_model/data/"*ic*"/data/errors/l1_error_upr.csv",DataFrame).DG_limit)

markerstyles_vec = [:false,:false,:false,:false,:false]
lwr = [(l1_cdf_data[i,j]) .- (l1_cdf_data_lwr[i,j]) for i in 1:size(l1_cdf_data,1), j in 1:size(l1_cdf_data,2)]
upr = [(l1_cdf_data_upr[i,j]) .- (l1_cdf_data[i,j]) for i in 1:size(l1_cdf_data,1), j in 1:size(l1_cdf_data,2)]
colours = [1;3;2;4;5]

@halfwidth_plot_defaults()
plot()
linestyles_vec = [:solid,:dash,:dashdot,:dot,:dash]
for (c,col) in enumerate(names(l1_cdf_data))
    os = collect(1:2:21)
    (c==5)&&(os.+=1)
    plot!(log10.(os),l1_cdf_data[:,col],ribbon=(lwr[:,c],upr[:,c]),
        label=_names[c],linestyle=linestyles_vec[c],
        marker=markerstyles_vec[c],
        linewidth=2,color=colours[c])
end
plot!(xlabel="Dimension")
plot!(ylabel="Error"); error_ticks!(plot!())
plot!(title="L¹ error - CDF")
# plot!(legend=(0.8,1.0))
plot!()
@add_lines!(l1_cdf_data,("Unif","QBDRAP"),@__DIR__)

savefig((@__DIR__)*"/l1_cdf_error_formatted.pdf")
end
p()


text:
export Todo, Todos

using Channels

type Todo <: AbstractModel
  ### internals
  _table_name::String
  _id::String

  id::Nullable{SearchLight.DbId}
  title::String
  completed::Bool
  created_at::DateTime
  updated_at::DateTime

  validator::ModelValidator

  on_dehydration::Function
  on_hydration::Function

  after_save::Function

  scopes::Dict{Symbol,Vector{SearchLight.SQLWhereEntity}}

  ### constructor
  Todo(;
    id = Nullable{SearchLight.DbId}(),
    title = "",
    completed = false,
    created_at = Dates.now(),
    updated_at = Dates.now(),

    validator = ModelValidator([
      (:title, Validation.TodoValidator.not_empty)
    ]),

    on_dehydration = (m::Todo, field::Symbol, value::Any) -> begin
      if field == :updated_at
        Dates.now()
      end
    end,
    on_hydration = (m::Todo, field::Symbol, value::Any) -> begin
      if field == :updated_at

      end
    end,

    after_save = (m::Todo) -> begin
      Channels.broadcast("""{ "channel":"todos", "message":"update", "payload": { "todo": { "id": $(m.id |> Base.get) } } }""")
    end,

    scopes = Dict(:active => [SQLWhereExpression("completed = FALSE")])

  ) = new("todos", "id",
          id,
          title,
          completed,
          created_at,
          updated_at,

          validator,

          on_dehydration,
          on_hydration,

          after_save,

          scopes
          )
end

module Todos

using App, Faker

function random()
  todo = Todo()
  todo.title = Faker.sentence()

  todo
end

end

max_stars_repo_name: essenciary/genie-todo-mvc

max_issues_repo_name: essenciary/genie-todo-mvc

max_forks_repo_name: essenciary/genie-todo-mvc


text:
# ------------------------------------------ #
#                  Constants                 #
# ------------------------------------------ #
Pr = 13.5
k = 0.57  # thermal conductivity of water, W / (K*m)
nu = mu/rhow # kin. viscosity

# ------------------------------------------ #
#       Nusselt number parameterisations     #
# ------------------------------------------ #

# Where A, a and b are defined: Dittus-Boelter correlation A * Pr^a * Re^b
coefs = (
    standard = (A=0.023, a=2/5, b=0.8),
    sommers = (A=0.0025, a=0.4, b=0.95),
    lunardini = (A=0.0078, a=1/3, b=0.927),
    vincent = (A=0.332, a=1/3, b=0.74),
    ogier = (A=1.78, a=1/3, b=0.58),
    gnielinski = ()
)

# ------------------------------------------ #
#                  Formulas                  #
# ------------------------------------------ #

r(S) = sqrt(S/pi)
Re(Q,S) = Q/S * 2 * r(S) / nu
Nu_Dittus_Boelter(Re, p) = p.A * Pr^p.a .* Re.^p.b
Nu_Gnielinski(f, Re) = f./8 .* (Re .- 1000) .* Pr ./
                       (1 .+ 12.7 .* (f./8).^0.5 .* (Pr^(2/3) .- 1))
# select between the two depending on parameter:
Nu_fn(f, Re, p::Tuple{}) = Nu_Gnielinski(f, Re)
Nu_fn(f, Re, p) = Nu_Dittus_Boelter(Re, p)

melt_energy_eq(Q, dphi_dz, dp_dz) = -Q * ( dphi_dz + rhow*cw*ct*dp_dz ) # Energy

# from Sommers and RAJARAM 2020, Eq. 29
z_eq_fn(Q, Nu) = rhow * cw * Q ./ (pi * k * Nu) # e-folding length, m
tau_eq_fn(Q, dphi_dz, dp_dz, Nu) = melt_energy_eq(Q, dphi_dz, dp_dz) / (k * pi * Nu ) # equilibrium offset-temperature: T_eq - T_melt, °C

# same paper, modification of Eq. 21
tau_w_diff(z_eq, dtauw_dz, tau_eq) = tau_eq - z_eq * dtauw_dz # water offset-temperature, T_w - T_melt, °C, calculated from differential equation
tau_w_sol(z, tau_eq, t0, z_eq) = tau_eq + (t0-tau_eq)*exp(-z/z_eq) # water offset-temperature, T_w - T_melt, °C, calculated from solution of differential equation as a function of z

# ------------------------------------------ #
#          Initialise dictionaries           #
# ------------------------------------------ #

z_eq = Dict()
tau_eq = Dict()
tau_w = Dict()
tau_diff = Dict() # Tw - Teq = tau_w - tau_eq
T_surf = Dict()
tauw_measured = Dict()
dtauw_dz = Dict()
Nus = Dict()

# ------------------------------------------ #
#    Calculate heat transfer parameters      #
# ------------------------------------------ #

days = ["0908", "2108"]

for date in days

    # Quantities derived from field experiments
    Qs = mid_309_265[date][:Q]
    Ss = mid_309_265[date][:S]
    f = mid_309_265[date][:f]
    dphi_dzs = mid_309_265[date][:dphi_dz]
    dp_dzs = mid_309_265[date][:dp] ./ mid_309_265[date][:dz]
    T_ct = ct * mid_309_265[date][:p_mean]

    # Data using output of free-gradient model: dtauw_dz = dTw_dz - dTmelt_dz, the first term comes from the MCMC runs
    local dtauw_dz = bootstrap(Particles(model_runs[date][:dTdz_MCMC]), n_partcl) .-  model_runs[date][:dTdz_ct]  # first term has to be converted from array to particles

    # Temperature measurements
    i = indexin(DateTime.(mid_309_265[date][:t_inj]), DateTime.(ctd309[date][:t]))  # extract the indices where there was a tracer experiment
    tauw_measured[date] = Particles(n_partcl, Normal(                             # compute the mean value over the date since it does not vary much
            mean(
            [mean(ctd309[date][:temp_PT][i]),
             mean(ctd265[date][:temp_PT][i])]), e_T)
            ) .- mean(T_ct, dims=1)

    # mean depth of test-section
    d = mean([ctd309[date][:depth], ctd265[date][:depth]], dims=1)

    # compute Reynolds number
    Res = Re.(Qs, Ss)

    # e-folding length
    z_eq[date] = Dict()
    for n = 1:length(coefs)
        c = coefs[n]
        z_eq[date][keys(coefs)[n]] = z_eq_fn.(Qs, Nu_fn(f, Res, c))
    end

    # equilibrium offset-temperature and Nusselt number
    tau_eq[date] = Dict()
    Nus[date] = Dict()
    for n = 1:length(coefs)
        c = coefs[n]
        tau_eq[date][keys(coefs)[n]] = tau_eq_fn.(Qs, dphi_dzs, dp_dzs, Nu_fn(f, Res, c))
        Nus[date][keys(coefs)[n]] = Nu_fn(f, Res, c)
    end

    # water offset-temperature, its difference to tau_eq and the estimated surface temperature
    tau_w[date] = Dict() # Tw - Tmelt
    tau_diff[date] = Dict()
    T_surf[date] = Dict()
    for k = keys(coefs)
        c = coefs[k]
        tau_w[date][k] = tau_w_diff.(z_eq[date][k], dtauw_dz, tau_eq[date][k])
        tau_diff[date][k] = tau_w[date][k] .- tau_eq[date][k]
        T_surf[date][k] = tau_w_sol.(-d, tau_eq[date][k], tau_w[date][k], z_eq[date][k])
    end
end


text:
function matrixnetwork_test()
    A = load_matrix_network("dfs_example")
    M = MatrixNetwork(A)
    B = sparse(M)
    @test A == B
    C = sparse_transpose(M)
    @test A' == C

    load_matrix_network_all("minnesota")
    load_matrix_network_all("U3A")

    #@show issym(sparse([0 1; 0 0]'))
    @test is_undirected(MatrixNetwork(sparse([0 1; 1 0]))) == true
    @test is_undirected(MatrixNetwork(sparse([0 1; 0 1]))) == false
    @test is_undirected(MatrixNetwork(load_matrix_network("dfs_example"))) == false

    @test is_connected(havel_hakimi_graph([1,1,1,1])) == false
    @test is_connected(havel_hakimi_graph([2,2,2])) == true
    @test is_connected(sparse([0 1 0; 0 0 1; 1 0 0])) == true
    @test is_connected(spzeros(1,1)) == true

    return true

end


text:
using ModelingToolkit
using Base.Test

struct Lorenz <: AbstractComponent
  x::Variable
  y::Variable
  z::Variable
  σ::Variable
  ρ::Variable
  β::Variable
  eqs::Vector{Expression}
end
function generate_lorenz_eqs(t,x,y,z,σ,ρ,β)
  D = Differential(t)
  [D*x ~ σ*(y-x)
   D*y ~ x*(ρ-z)-y
   D*z ~ x*y - β*z]
end
Lorenz(t) = Lorenz(first(@DVar(x(t))),first(@DVar(y(t))),first(@DVar(z(t))),first(@Param(σ)),first(@Param(ρ)),first(@Param(β)),generate_lorenz_eqs(t,x,y,z,σ,ρ,β))

@IVar t
lz1 = Lorenz(t)
lz2 = Lorenz(t)
Expression[lz1.x ~ lz2.x
           lz1
           lz2]


text:
using StatsBase
using LinearAlgebra

function ridge(x::Matrix,y::Vector,λ=0)
    @assert λ ≥ 0
    X = copy(x)
    n,p = size(X)
    X̄, σ = mean_and_std(X,1) #説明変数の平均と標準偏差を記録
    ȳ = mean(y) #目的変数については平均値を取得
    zscore!(X,X̄,σ) #標準化処理
    y .-= ȳ
    β = (X'X + n*λ*I(p)) \ X'*y
    β ./= vec(σ) #X'Xで割っているので、βをσで割っておく
    β₀ = ȳ .- X̄*β
    Dict("β" => β, "β₀" => β₀) #vcatでつなげた方が良いかも
end

soft_th(λ::Number,x::Number) = sign(x) * max(abs(x)-λ,0)

function lasso(x::Matrix,y::Vector,λ=0)
    @assert λ ≥ 0
    X = copy(x)
    n,p = size(X)
    X̄, σ = mean_and_std(X,1) #説明変数の平均と標準偏差を記録
    ȳ = mean(y) #目的変数については平均値を取得
    zscore!(X,X̄,σ) #標準化処理
    y .-= ȳ

    #ここまでridgeと同じ
    ϵ = 1.0
    β = zeros(p); β_old = zeros(p)
    while ϵ > 0.001
        for j in 1:p
            index = setdiff(1:p,j)
            r = y - X[:,index]*β[index]
            β[j] = soft_th(λ,r'*X[:,j]/n)
        end
        ϵ = abs.(β-β_old) |> maximum
        β_old = β
    end
    β ./= vec(σ) #X'Xで割っているので、βをσで割っておく
    β₀ = ȳ .- X̄*β
    Dict("β" => β, "β₀" => β₀) #vcatでつなげた方が良いかも
end


text:
const PT{T, N} = Union{PredictiveAsymmetryTest{T, N},
    NormalisedPredictiveAsymmetryTest{T, N}} where {T, N}
const SPT{T, N} = SubsampledDataCausalityTest{<:PT{T, N}} where {T, N}

"""
    summarise(f::Function, ::CausalAnalysis{MetaCausalityTest}, args...; kwargs...)

Summarise a causality analysis performed by repeated subsampling of the input data using
the statistic `f`. The statistic is apply on a per-subsample basis, so that one value
for the summary statistic is returned for each evaluation of the causality test.
"""
function summarise(f::Function, analysis::CausalAnalysis{MetaCausalityTest},
        args...; kwargs...)
    f.(analysis.result, args...; kwargs...)
end

"""
For predictive tests, don't summarise on a per-subsample basis, but gather the
subsamples and summarise per prediction lag (i.e. over the `ηs`).
"""
function summarise(f::Function, analysis::CausalAnalysis{CT, DT, RT},
        args...; kwargs...) where {CT <: SPT{T, N}, DT, RT} where {T, N}

    M = hcat(analysis.result...)
    η_results = [M[i, :] for i = 1:N]
    f.(η_results, args...; kwargs...)
end


text:
"""
Reverse function for `asinh`.
"""
function asinh_rev(y::Interval,x::Interval)
    x = x ∩ sinh(y)

    return y, x
end

"""
Reverse function for `acosh`.
"""
function acosh_rev(y::Interval,x::Interval)
    y_new = y ∩ Interval(0.0,∞)
    x = x ∩ cosh(y)

    return y_new, x
end

"""
Reverse function for `atanh`.
"""
function atanh_rev(y::Interval,x::Interval)
    x = x ∩ tanh(y)

    return y, x
end


text:
# Autogenerated wrapper script for rtlsdr_jll for i686-linux-musl
export librtlsdr

using libusb_jll
JLLWrappers.@generate_wrapper_header("rtlsdr")
JLLWrappers.@declare_library_product(librtlsdr, "librtlsdr.so.0")
function __init__()
    JLLWrappers.@generate_init_header(libusb_jll)
    JLLWrappers.@init_library_product(
        librtlsdr,
        "lib/librtlsdr.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()