text:
using BOMBS
using Test
using CSV
using DataFrames
using LinearAlgebra
using Dates
using Distributions
using Random
using CmdStan

# -------------------------------------------------------- MODEL DEFINITION TESTS
model_def = Dict();
model_def["NameF"] = [];
model_def["nStat"] = [];
model_def["nPar"] = [];
model_def["nInp"] = [];
model_def["stName"] = [];
model_def["parName"] = [];
model_def["inpName"] = [];
model_def["eqns"] = [];
model_def["Y0eqs"] = [];
model_def["Y0Sim"] = [];
model_def["tols"] = [];
model_def["solver"] = [];

model_def2 = defModStruct();
model_def2["NameF"] = "Test";
model_def2["nStat"] = 2;
model_def2["nPar"] = 4;
model_def2["nInp"] = 1;
model_def2["stName"] = ["A", "B"];
model_def2["parName"] = ["k1", "k2", "k3", "k4"];
model_def2["inpName"] = ["inp1"];
model_def2["eqns"] = ["dA = k1*A - k2*A*inp1", "dB = k3*(A^2) - k4*B"];
model_def2["Y0eqs"] = ["A = 1", "B = k3*(A^2)/k4"];
model_def2["Y0Sim"] = false;
model_def2["tols"] = [1e-5,1e-5];
model_def2["solver"] = "CVODE_BDF";

# -------------------------------------------------------- SIMULATION TESTS
simul_def = Dict()
simul_def["Nexp"] = [];
simul_def["finalTime"] = [];
simul_def["switchT"] = [];
simul_def["y0"] = [];
simul_def["preInd"] = [];
simul_def["uInd"] = [];
simul_def["theta"] = [];
simul_def["tsamps"] = [];
simul_def["plot"] = [];
simul_def["flag"] = [];

simul_def2 = defSimulStruct();
simul_def2["Nexp"] = 2;
simul_def2["finalTime"] = [40, 40];
simul_def2["switchT"] = [[0,20,40], [0,20,40]];
simul_def2["y0"] = [[0,0], [0,0]];
simul_def2["preInd"] = [[0.1], [0.1]];
simul_def2["uInd"] = [[1,1], [1,1]];
simul_def2["theta"] = [0.1 0.2 0.2 0.02; 0.2 0.1 0.2 0.01; 1 1 1 0.1];
simul_def2["tsamps"] = [[0,5,10,15,20,25,30,35,40], [0,3,5,10,12,15,20,25,27, 30,33, 35,40]];
simul_def2["plot"] = false;
simul_def2["flag"] = "testsim";

simul_def3 = Dict()
simul_def3["ObservablesFile"] = [];
simul_def3["EventInputsFile"] = [];
simul_def3["theta"] = [];
simul_def3["MainDir"] = [];
simul_def3["plot"] = [];
simul_def3["flag"] = [];

simul_defCSV1 = Dict();
simul_defCSV1["ObservablesFile"] = ["TestSimulCSVObs.csv"];
simul_defCSV1["EventInputsFile"] = ["TestSimulCSVInps.csv"];
simul_defCSV1["theta"] = ["Theta1.csv"];
simul_defCSV1["MainDir"] = [];
simul_defCSV1["plot"] = false;
simul_defCSV1["flag"] = "testsimCSV";

simul_defCSV2 = Dict();
simul_defCSV2["Nexp"] = 1;
simul_defCSV2["finalTime"] = [40];
simul_defCSV2["switchT"] = [[0.0,20.0,40.0]];
simul_defCSV2["y0"] = [[0,0]];
simul_defCSV2["preInd"] = [[0.1]];
simul_defCSV2["uInd"] = [reshape([1, 0.5], 2,1)];
simul_defCSV2["theta"] = [0.1 0.2 0.2 0.002];
simul_defCSV2["tsamps"] = [[0,5,10,15,20,25,30,35,40]];
simul_defCSV2["plot"] = false;
simul_defCSV2["flag"] = "testsimCSV";

# -------------------------------------------------------- PSEUDO-DATA TESTS

pseudo_def = Dict();
pseudo_def["Nexp"] = [];
pseudo_def["finalTime"] = [];
pseudo_def["switchT"] = [];
pseudo_def["y0"] = [];
pseudo_def["preInd"] = [];
pseudo_def["uInd"] = [];
pseudo_def["theta"] = [];
pseudo_def["tsamps"] = [];
pseudo_def["plot"] = [];
pseudo_def["flag"] = [];
pseudo_def["Obs"] = [];
pseudo_def["Noise"] = [];

pseudo_def2 = defPseudoDatStruct();
pseudo_def2["Nexp"] = 2;
pseudo_def2["finalTime"] = [40, 40];
pseudo_def2["switchT"] = [[0,20,40], [0,20,40]];
pseudo_def2["y0"] = [[0,0], [0,0]];
pseudo_def2["preInd"] = [[0.1], [0.1]];
pseudo_def2["uInd"] = [[1,1], [1,1]];
pseudo_def2["theta"] = [0.1 0.2 0.2 0.02; 0.2 0.1 0.2 0.01; 1 1 1 0.1];
pseudo_def2["tsamps"] = [[0,5,10,15,20,25,30,35,40], [0,3,5,10,12,15,20,25,27, 30,33, 35,40]];
pseudo_def2["plot"] = false;
pseudo_def2["flag"] = "testPD";
pseudo_def2["Obs"] = ["B*2"];
pseudo_def2["Noise"] = [0.1];

pseudo_def3 = Dict()
pseudo_def3["ObservablesFile"] = [];
pseudo_def3["EventInputsFile"] = [];
pseudo_def3["theta"] = [];
pseudo_def3["MainDir"] = [];
pseudo_def3["plot"] = [];
pseudo_def3["flag"] = [];
pseudo_def3["Obs"] = [];
pseudo_def3["Noise"] = [];

pseudo_defCSV1 = Dict();
pseudo_defCSV1["ObservablesFile"] = ["TestSimulCSVObsPD.csv"];
pseudo_defCSV1["EventInputsFile"] = ["TestSimulCSVInpsPD.csv"];
pseudo_defCSV1["theta"] = ["Theta1PD.csv"];
pseudo_defCSV1["MainDir"] = [];
pseudo_defCSV1["plot"] = false;
pseudo_defCSV1["flag"] = "testPDCSV";
pseudo_defCSV1["Obs"] = ["B*2"];
pseudo_defCSV1["Noise"] = [0.1];

pseudo_defCSV2 = Dict();
pseudo_defCSV2["Nexp"] = 1;
pseudo_defCSV2["finalTime"] = [40];
pseudo_defCSV2["switchT"] = [[0.0,20.0,40.0]];
pseudo_defCSV2["y0"] = [[0,0]];
pseudo_defCSV2["preInd"] = [[0.1]];
pseudo_defCSV2["uInd"] = [reshape([1, 0.5], 2,1)];
pseudo_defCSV2["theta"] = [0.1 0.2 0.2 0.002];
pseudo_defCSV2["tsamps"] = [[0,5,10,15,20,25,30,35,40]];
pseudo_defCSV2["plot"] = false;
pseudo_defCSV2["flag"] = "testPDCSV";
pseudo_defCSV2["Obs"] = ["B*2"];
pseudo_defCSV2["Noise"] = [0.1];

# -------------------------------------------------------- MLE TESTS
mle_def = Dict();
mle_def["Nexp"] = [];
mle_def["finalTime"] = [];
mle_def["switchT"] = [];
mle_def["y0"] = [];
mle_def["preInd"] = [];
mle_def["uInd"] = [];
mle_def["tsamps"] = [];
mle_def["plot"] = [];
mle_def["flag"] = [];
mle_def["thetaMAX"] = [];
mle_def["thetaMIN"] = [];
mle_def["runs"] = [];
mle_def["parallel"] = [];
mle_def["DataMean"] = [];
mle_def["DataError"] = [];
mle_def["Obs"] = [];
mle_def["OPTsolver"] = [];
mle_def["MaxTime"] = [];
mle_def["MaxFuncEvals"] = [];

mle_def2 = Dict();
mle_def2["Nexp"] = 1;
mle_def2["finalTime"] = [40];
mle_def2["switchT"] = [[0,20,40]];
mle_def2["y0"] = [[0,0]];
mle_def2["preInd"] = [[0.1]];
mle_def2["uInd"] = [[1,1]];
mle_def2["tsamps"] = [[0,5,10,15,20,25,30,35,40]];
mle_def2["plot"] = false;
mle_def2["flag"] = "testmle";
mle_def2["thetaMAX"] = [0.15 0.25 0.25 0.025];
mle_def2["thetaMIN"] = [0.1 0.2 0.2 0.02];
mle_def2["runs"] = 2;
mle_def2["parallel"] = false;
mle_def2["DataMean"] = [[0 1 2 3 4 5 6 7 8; 0 1 2 3 4 5 6 7 8]'];
mle_def2["DataError"] = [[[0,1,2,3,4,5,6,7,8], [0,1,2,3,4,5,6,7,8]]];
mle_def2["Obs"] = ["A", "B"];
mle_def2["OPTsolver"] = "adaptive_de_rand_1_bin_radiuslimited";
mle_def2["MaxTime"] = [];
mle_def2["MaxFuncEvals"] = 10;

cvmle_def = Dict();
cvmle_def["Nexp"] = [];
cvmle_def["finalTime"] = [];
cvmle_def["switchT"] = [];
cvmle_def["y0"] = [];
cvmle_def["preInd"] = [];
cvmle_def["uInd"] = [];
cvmle_def["theta"] = [];
cvmle_def["tsamps"] = [];
cvmle_def["plot"] = [];
cvmle_def["flag"] = [];
cvmle_def["DataMean"] = [];
cvmle_def["DataError"] = [];
cvmle_def["Obs"] = [];

cvmle_def2 = Dict();
cvmle_def2["Nexp"] = 1;
cvmle_def2["finalTime"] = [40];
cvmle_def2["switchT"] = [[0,20,40]];
cvmle_def2["y0"] = [[0,0]];
cvmle_def2["preInd"] = [[0.1]];
cvmle_def2["uInd"] = [[1,1]];
cvmle_def2["theta"] = convert(Array, [0.1 0.2 0.2 0.02;0.12 0.22 0.22 0.022]');
cvmle_def2["tsamps"] = [[0,5,10,15,20,25,30,35,40]];
cvmle_def2["plot"] = false;
cvmle_def2["flag"] = "cvmletest";
cvmle_def2["DataMean"] = [[0 1 2 3 4 5 6 7 8]'];
cvmle_def2["DataError"] = [[[0,1,2,3,4,5,6,7,8]]];;
cvmle_def2["Obs"] = ["B"];

testdict1 = Dict();
testdict2 = Dict();
testdict1["a"] = "a";
testdict1["b"] = "b";
testdict2["a"] = [];
testdict2["b"] = 1;

@testset "ModelGenTests" begin

    # In case any modification to the structure is made, this will be the test to see if I have forgot some entry
    @test model_def == defModStruct()

    # Check function does some modifications, so it is good to see that it does what it should
    @test model_def2 == checkStruct(model_def2)

    # Check that a file has been generated (the functions containing the ODEs)
    GenerateModel(model_def2)
    @test isfile(string(pwd(),"\\ModelsFunctions\\", model_def2["NameF"], "_Model.jl"))

    rm(string(pwd(),"\\ModelsFunctions\\", model_def2["NameF"], "_Model.jl"))
    # rm(string(pwd(),"\\ModelsFunctions"))

end

@testset "ModelSimulationTests" begin

    CSV.write("Theta1.csv", DataFrame([0.1 0.2 0.2 0.002]))
    CSV.write("TestSimulCSVInps.csv", DataFrame([0 40 0.1 1; 20 40 0.1 0.5]))
    CSV.write("TestSimulCSVObs.csv", DataFrame([0 0 0; 5 0 0; 10 0 0; 15 0 0; 20 0 0; 25 0 0; 30 0 0; 35 0 0; 40 0 0]))

    # In case any modification to the structure is made, this will be the test to see if I have forgot some entry
    @test simul_def == defSimulStruct()

    # Check function does some modifications, so it is good to see that it does what it should
    @test simul_def2 == checkStructSimul(model_def2, simul_def2)

    # In case any modification to the structure is made, this will be the test to see if I have forgot some entry
    @test simul_def3 == defSimulStructFiles()

    # Might need to add a test for this, but then I need some example CSV files
    @test simul_defCSV2 == extractSimulCSV(model_def2, simul_defCSV1)

    rm("Theta1.csv")
    rm("TestSimulCSVInps.csv")
    rm("TestSimulCSVObs.csv")

    GenerateModel(model_def2);
    simuls, ~, ~ = simulateODEs(model_def2, simul_def2);

    @test length(simuls) == simul_def2["Nexp"]
    for i in 1:2
        a,b,c = size(simuls[string("Exp_", i)])
        @test a == length(simul_def2["tsamps"][i])
        @test b == (model_def2["nStat"])
        @test c == length(simul_def2["theta"])/model_def2["nPar"]
        @test sum(simuls[string("Exp_", i)]) != 0
    end

    plotSimsODE(simuls,model_def2,simul_def2)
    for i in 1:2
        @test isfile(string(simul_def2["savepath"], "\\PlotSimulation_Exp", i,"_", simul_def2["flag"], ".png"))
        rm(string(simul_def2["savepath"], "\\PlotSimulation_Exp", i,"_", simul_def2["flag"], ".png"))
    end

    rm(string(pwd(),"\\ModelsFunctions\\", model_def2["NameF"], "_Model.jl"))
    rm(string(simul_def2["savepath"], "\\", simul_def2["savename"]))
    rm(string(simul_def2["savepath"]))

end

@testset "PseudoDataGenerationTests" begin

    CSV.write("Theta1PD.csv", DataFrame([0.1 0.2 0.2 0.002]))
    CSV.write("TestSimulCSVInpsPD.csv", DataFrame([0 40 0.1 1; 20 40 0.1 0.5]))
    CSV.write("TestSimulCSVObsPD.csv", DataFrame([0 0 0; 5 0 0; 10 0 0; 15 0 0; 20 0 0; 25 0 0; 30 0 0; 35 0 0; 40 0 0]))

    # In case any modification to the structure is made, this will be the test to see if I have forgot some entry
    @test pseudo_def == defPseudoDatStruct()

    # Check function does some modifications, so it is good to see that it does what it should
    @test pseudo_def2 == checkStructPseudoDat(model_def2, pseudo_def2)

    # In case any modification to the structure is made, this will be the test to see if I have forgot some entry
    @test pseudo_def3 == defPseudoDatStructFiles()

    # Might need to add a test for this, but then I need some example CSV files
    @test pseudo_defCSV2 == extractPseudoDatCSV(model_def2, pseudo_defCSV1)

    rm("Theta1PD.csv")
    rm("TestSimulCSVInpsPD.csv")
    rm("TestSimulCSVObsPD.csv")

    GenerateModel(model_def2)
    pdat, ~, ~ = GenPseudoDat(model_def2, pseudo_def2);

    @test length(pdat["SimsObs"]) == pseudo_def2["Nexp"]
    @test length(pdat["Sims"]) == pseudo_def2["Nexp"]
    @test length(pdat["PData"]) == pseudo_def2["Nexp"]
    @test length(pdat["PError"]) == pseudo_def2["Nexp"]

    for i in 1:2
        a,b,c = size(pdat["Sims"][string("Exp_", i)])
        @test a == length(pseudo_def2["tsamps"][i])
        @test b == (model_def2["nStat"])
        @test c == length(pseudo_def2["theta"])/model_def2["nPar"]
        @test sum(pdat["Sims"][string("Exp_", i)]) != 0

        a,b,c = size(pdat["SimsObs"][string("PDExp_", i)])
        @test b == length(pseudo_def2["Obs"]);
        @test pdat["SimsObs"][string("PDExp_", i)][:,1,:] == pdat["Sims"][string("Exp_", i)][:,2,:].*2

        a,b,c = size(pdat["PData"][string("PDExp_", i)])
        @test b == length(pseudo_def2["Obs"]);

        @test length(pdat["PData"][string("PDExp_", i)]) == length(pdat["PError"][string("PDExp_", i)])
    end

    plotPseudoDatODE(pdat,model_def2,pseudo_def2)
    for i in 1:2
        @test isfile(string(pseudo_def2["savepath"], "\\PlotPseudoDat_Exp", i,"_", pseudo_def2["flag"], ".png"))
        rm(string(pseudo_def2["savepath"], "\\PlotPseudoDat_Exp", i,"_", pseudo_def2["flag"], ".png"))
    end

    rm(string(pwd(),"\\ModelsFunctions\\", model_def2["NameF"], "_Model.jl"))
    rm(string(pseudo_def2["savepath"], "\\", pseudo_def2["savename"]))

    for i in 1:2
        rm(string(pseudo_def2["savepath"], "\\PseudoDataFiles\\", model_def2["NameF"], "_EXP", i, "_", pseudo_def2["flag"], "_Events_Inputs.csv"))
        rm(string(pseudo_def2["savepath"], "\\PseudoDataFiles\\", model_def2["NameF"], "_EXP", i, "_", pseudo_def2["flag"], "_Observables.csv"))
        rm(string(pseudo_def2["savepath"], "\\PseudoDataFiles\\", model_def2["NameF"], "_EXP", i, "_", pseudo_def2["flag"], "_Simulations.csv"))
    end
    rm(string(pseudo_def2["savepath"], "\\PseudoDataFiles"))
    rm(string(pseudo_def2["savepath"]))
end

@testset "MLESeriesTests" begin

    # First check that both structure calls give what it is supposed to
    @test mle_def == defMLEStruct();
    @test cvmle_def == defCrossValMLEStruct();

    # Check functionality of helpping functions
    testdict3 = SimToMle(testdict1, testdict2);
    @test testdict3 == testdict2;

    simst = reshape([1 2 3 4; 1 2 3 4; 1 2 3 4], 3,4,1);
    obst = ["c*2"];
    stnamest = ["a","b","c","d"];
    @test sum(selectObsSim_te(simst, obst, stnamest)) == sum(simst[:,3,:]*2);
    @test (selectObsSim_te(simst, [2], stnamest)) == reshape([2,2,2], 3,1,1);

    ttess1 = Dict();
    ttess1["nInp"] = 2;
    ttess2 = Dict();
    ttess2["uInd"] = [[1 1 1; 2 2 2]'];
    @test restructInputs_te(ttess1, ttess2, 1) == [1,2,1,2,1,2];
    ttess1["nInp"] = 3;
    ttess2["uInd"] = [[1 1 1; 2 2 2; 3 3 3]'];
    @test restructInputs_te(ttess1, ttess2, 1) == [1,2,3,1,2,3,1,2,3];

    # More tests might be needed?
    m=10;
    s=1;
    d=9;
    @test (-0.5)*(log(2*pi) + 0 + (1)) == UVloglike(d, m, s);
    m=[10, 10, 10, 10];
    s=[1,1,1,1];
    d=[9, 11, 9, 11];
    @test ((-0.5)*(log(2*pi) + 0 + (1)))*4 == UVloglike(d, m, s);
    m=[10, 15, 10, 15];
    s=[0.5,0.5,0.5,0.5];
    d=[12, 12, 12, 12];
    @test ((-0.5)*(log(2*pi) + log(0.25) + (2^2)/0.25)) + ((-0.5)*(log(2*pi) + log(0.25) + (2^2)/0.25)) +
    ((-0.5)*(log(2*pi) + log(0.25) + (3^2)/0.25)) + ((-0.5)*(log(2*pi) + log(0.25) + (3^2)/0.25)) == UVloglike(d, m, s);

    m=[10, 10];
    s=[1 1; 1 1]; # Diagonal element of 0.1 is added to avoid Infs due to non positive definite covariances.
    d=[9, 11];
    correcmat = Diagonal(ones(length(d))).*0.1;
    @test [MVloglike(d, m, s)] == (-0.5) * ((length(d)*log(2*pi)) .+ log(det(s.+correcmat)) .+ ([-1 1]*inv(s.+correcmat)*[-1,1]));

    # Check structure check functions (with and without CSVs)
    @test mle_def2 == checkStructMLE(model_def2, mle_def2)
    @test cvmle_def2 == checkStructCrossValMLE(model_def2, cvmle_def2)

    # Check MLE
    GenerateModel(model_def2);
    mle_res, model_def3, mle_def3 = MLEtheta(model_def2, mle_def2);

    @test !isempty(mle_res);
    @test typeof(mle_res["StanDict"]) <: Array;
    @test typeof(mle_res["StanDict"][1]) <: Dict;
    @test typeof(mle_res["StanDict"][2]) <: Dict;
    @test symdiff(["k1","k2","k3","k4"],keys(mle_res["StanDict"][1])) == [];
    @test symdiff(["k1","k2","k3","k4"],keys(mle_res["StanDict"][2])) == [];

    @test length(mle_res["Theta"]) == 4*2;
    @test size(mle_res["Theta"]) == (4,2);

    @test length(mle_res["convCurv"]) == 2;

    @test (typeof(mle_res["convCurv"][1][1]) == Tuple{Int,Float64} || typeof(mle_res["convCurv"][1][1]) == Tuple{Int,Float32});
    @test (typeof(mle_res["convCurv"][1][2]) == Tuple{Int,Float64} || typeof(mle_res["convCurv"][1][2]) == Tuple{Int,Float32});
    @test (typeof(mle_res["convCurv"][2][1]) == Tuple{Int,Float64} || typeof(mle_res["convCurv"][2][1]) == Tuple{Int,Float32});
    @test (typeof(mle_res["convCurv"][2][2]) == Tuple{Int,Float64} || typeof(mle_res["convCurv"][2][2]) == Tuple{Int,Float32});

    @test length(mle_res["BestTheta"]) == 4;
    @test length(mle_res["BestCFV"]) == 2;

    #----------------------------- the parallel stuff!!!!
    mle_def2["parallel"] = true;
    mle_res2, model_def3, mle_def3 = MLEtheta(model_def2, mle_def2);
    @test isempty(mle_res2);
    @test isfile(string(mle_def3["savepath"], "\\MLEScripts\\", model_def3["NameF"], "_MLE.jl"));
    mle_res3, ~, ~ = finishMLEres(mle_res, model_def2, mle_def2);
    @test !isempty(mle_res3);
    mle_def2["parallel"] = false;

    # Check CV for MLE results
    cvmle_res, ~, ~ = CrossValMLE(model_def2, cvmle_def2);
    @test !isempty(cvmle_res);
    @test !isempty(cvmle_res["BestSimulations"][1]);
    @test size(cvmle_res["BestSimulations"][1])[2] == 2;
    @test length(cvmle_res["Costs"]) == 2;
    @test !isempty(cvmle_res["BestSimObservables"][1]);
    @test size(cvmle_res["BestSimObservables"][1])[2] == 1;
    @test length(cvmle_res["BestTheta"]) == 4;
    @test typeof(cvmle_res["SimObservables"]) <: Dict;
    @test size(cvmle_res["SimObservables"]["ExpObs_1"])[2] == 1;
    @test typeof(cvmle_res["Simulations"]) <: Dict;
    @test size(cvmle_res["Simulations"]["Exp_1"])[2] == 2;

    # Test Plotting
    plotMLEResults(mle_res,model_def2,mle_def2)
    @test isfile(string(mle_def2["savepath"], "\\PlotMLEResults_Exp", 1,"_", mle_def2["flag"], ".png"))
    @test isfile(string(mle_def2["savepath"], "\\Plot_MLEConvergence", "_", mle_def2["flag"], ".png"))
    rm(string(mle_def2["savepath"], "\\PlotMLEResults_Exp", 1,"_", mle_def2["flag"], ".png"))
    rm(string(mle_def2["savepath"], "\\Plot_MLEConvergence", "_", mle_def2["flag"], ".png"))

    rm(string(pwd(),"\\ModelsFunctions\\", model_def2["NameF"], "_Model.jl"))
    rm(string(mle_def2["savepath"], "\\", mle_def2["savename"]))

    rm(string(mle_def2["savepath"], "\\MLEScripts\\", model_def2["NameF"], "_MLE.jl"))
    rm(string(mle_def2["savepath"], "\\MLEScripts"))

    rm(string(cvmle_def2["savepath"], "\\", cvmle_def2["savename"]))
    rm(string(mle_def2["savepath"], "\\",model_def2["NameF"], "_", today(), "_SimulationResults_MLEsimulations1.jld"))
    rm(string(mle_def2["savepath"], "\\",model_def2["NameF"], "_", today(), "_SimulationResults_MLEsimulations2.jld"))
    rm(string(mle_def2["savepath"]))

end

@testset "StanInferenceTests" begin

    # Tests on structure definitions
    bayinf_def = defBayInfStruct();
    @test typeof(bayinf_def) <: Dict;
    entries1 = ["Priors", "Data", "StanSettings", "flag", "plot", "runInf", "MultiNormFit"];
    @test isempty(symdiff(entries1,keys(bayinf_def)));
    datinf_def = defBayInfDataStruct();
    @test typeof(datinf_def) <: Dict;
    entries2 = ["Nexp", "finalTime", "switchT", "y0", "preInd", "uInd", "tsamps", "Obs", "DataMean", "DataError"];
    @test isempty(symdiff(entries2,keys(datinf_def)));
    csvinf_def = defBayInfDataFromFilesStruct();
    @test typeof(csvinf_def) <: Dict;
    entries3 = ["Obs", "Observables", "Inputs", "y0"];
    @test isempty(symdiff(entries3,keys(csvinf_def)));
    stainf_def = defBasicStanSettingsStruct();
    @test typeof(stainf_def) <: Dict;
    entries4 = ["cmdstan_home", "nchains", "nsamples", "nwarmup", "printsummary", "init", "maxdepth", "adaptdelta", "jitter"];
    @test isempty(symdiff(entries4,keys(stainf_def)));

    # Check structures
    stainf_def["cmdstan_home"] = "C:/Users/David/.cmdstanpy/cmdstan-2.20.0";
    stainf_def["nchains"] = 2;
    stainf_def["nsamples"] = 20;
    stainf_def["nwarmup"] = 20;
    stainf_def["printsummary"] = false;
    stainf_def["init"] = [];
    stainf_def["maxdepth"] = 13;
    stainf_def["adaptdelta"] = 0.95;
    stainf_def["jitter"] = 0.5;
    @test stainf_def == checkStructBayInfStanSettings(model_def2, stainf_def);

    datinf_def["Nexp"] = 1;
    datinf_def["finalTime"] = [40];
    datinf_def["switchT"] = [[0,20,40]];
    datinf_def["y0"] = [[0,0]];
    datinf_def["preInd"] = [[0.1]];
    datinf_def["uInd"] = [[1,1]];
    datinf_def["tsamps"] = [[0,5,10,15,20,25,30,35,40]];
    datinf_def["DataMean"] = [[0 1 2 3 4 5 6 7 8; 0 1 2 3 4 5 6 7 8]'];
    datinf_def["DataError"] = [[[0,1,2,3,4,5,6,7,8], [0,1,2,3,4,5,6,7,8]]];
    datinf_def["Obs"] = ["A", "B"];
    @test datinf_def == checkStructBayInfData(model_def2, datinf_def);

    bayinf_def["Priors"] = [0.15 0.25 0.25 0.025; 0.1 0.2 0.2 0.02]'';
    bayinf_def["Data"] = datinf_def;
    bayinf_def["StanSettings"] = stainf_def;
    bayinf_def["flag"] = "stantest";
    bayinf_def["plot"] = false;
    bayinf_def["runInf"] = false;
    bayinf_def["MultiNormFit"] = false;
    bayinf_def2 =  checkStructBayInf(model_def2, bayinf_def);
    @test bayinf_def["Data"] == bayinf_def2["Data"]
    @test bayinf_def["StanSettings"] == bayinf_def2["StanSettings"]
    @test bayinf_def["flag"] == bayinf_def2["flag"]
    @test bayinf_def["plot"] == bayinf_def2["plot"]
    @test bayinf_def["runInf"] == bayinf_def2["runInf"]
    @test bayinf_def["MultiNormFit"] == bayinf_def2["MultiNormFit"]
    @test length(bayinf_def2["Priors"]) == 3;
    @test typeof(bayinf_def2["Priors"]) <: Dict
    # checkStructBayInfDataFiles(model_def, data_def)

    # Helpper functions tests
    s1 = [0,1,3,9,10]
    @test convertBoundTo2(s1, 10, 20) == s1.+10;
    @test convertBoundTo2(s1.+10, 0, 10) == s1;

    samp1 = rand(Normal(10,1), 80000,4);
    fipri1 = fitPriorSamps(samp1, model_def2);
    @test length(fipri1["pars"]) == 4;
    @test length(fipri1["transpars"]) == 5;
    @test length(fipri1["pridis"]) == 4;
    fipri2 = fitPriorSampsMultiNorm(samp1, model_def2);
    @test length(fipri2["pars"]) < 4;
    @test length(fipri2["transpars"]) == 5;
    @test length(fipri2["pridis"]) < 4;
    kee = ["transpars","mera","cora","pars","pridis","numN"]
    @test isempty(symdiff(kee,keys(fipri2)));

    dis1 = genStanInitDict(samp1', model_def2["parName"], 10);
    @test typeof(dis1[1]) <: Dict;
    @test isempty(symdiff(model_def2["parName"],keys(dis1[1])));
    @test length(dis1) == 10;

    # reparamDictStan(dis1, bayinf_def)

    # Stan
    genStanModel(model_def2, bayinf_def)
    @test isfile(string(pwd(),"\\ModelsFunctions\\", model_def2["NameF"], "_StanModel.stan"))
    sdat = restructureDataInference(model_def2, bayinf_def);
    datks = ["elm","tml","ts","tsl","tsmax","Nsp","inputs","evnT","m","stsl","stslm","sts","obser","obSta","nindu","preInd","Means","Y0us","Erros"]
    for i in 1:length(datks)
        @test datks[i] in keys(sdat)
        @test !isempty(sdat[datks[i]])
    end
    modelpath, Model, StanModel, inferdata, init, ~, ~ = getStanInferenceElements(model_def2, bayinf_def);
    @test isfile(modelpath);
    rm(modelpath)
    @test typeof(Model) == String;
    @test typeof(StanModel) == CmdStan.Stanmodel;
    @test inferdata == sdat;
    @test isempty(init);

    # saveStanResults(rc, chns, cnames, model_def, bayinf_def)
    # runStanInference(model_def, bayinf_def)

    stin, ~, ~ = StanInfer(model_def2, bayinf_def);
    kssm = ["init", "StanModel", "inferdata", "Model", "modelpath"]
    @test isempty(symdiff(kssm,keys(stin)));

    # Entropy tests
    ps = genSamplesPrior(model_def2, bayinf_def, 100);
    @test size(ps) == (100, 4);

    w = [1];
    E = [[1 0; 0 1]];
    MU = reshape([10, 10], 2, 1);
    x = reshape([10, 10], 2, 1);
    @test round(BOMBS.H_Upper(w,E)) == 3;
    @test round(BOMBS.mvGauss(x, MU, E[1])[1], digits=2) == 0.16;
    @test round(BOMBS.H_Lower(w, E, MU'), digits = 1) == 2.5;
    @test round(BOMBS.GaussMix(x[:,1], convert(Array, MU'), E, w), digits=2) == 0.16;
    @test round(BOMBS.ZOTSE(MU', E, w)) == 2;
    # Due to the use of global variables I cannot test the other functions... But I can test that the main script runs.
    # Do not know why, but I cannot make ScikitLearn work in the test file (works in scripts and jupyter)??? Need to have a look at it.
    # computeH(ps, model_def2, "test")

    # Plots
    # plotStanResults(staninf_res, model_def, bayinf_def)

    rm(string(pwd(),"\\ModelsFunctions\\", model_def2["NameF"], "_StanModel.stan"));
    rm(string(pwd(),"\\tmp\\", model_def2["NameF"], "_Stan_",bayinf_def["flag"],".stan"));
end

@testset "OEDModelSelectionTests" begin

    # structure definition
    oedms_def = defODEModelSelectStruct();
    @test typeof(oedms_def) <: Dict;
    entries1 = ["Model_1", "Model_2", "Obs", "Theta_M1", "Theta_M2", "y0_M1", "y0_M2", "preInd_M1", "preInd_M2",
                "finalTime", "switchT", "tsamps", "equalStep",
                "fixedInp", "fixedStep", "plot", "flag", "uUpper", "uLower", "maxiter"];
    @test isempty(symdiff(entries1,keys(oedms_def)));

    #structure check
    model_def3 = defModStruct();
    model_def3["NameF"] = "Test2";
    model_def3["nStat"] = 2;
    model_def3["nPar"] = 4;
    model_def3["nInp"] = 1;
    model_def3["stName"] = ["A", "B"];
    model_def3["parName"] = ["k1", "k2", "k3", "k4"];
    model_def3["inpName"] = ["inp1"];
    model_def3["eqns"] = ["dA = k1*A - k2*A*inp1", "dB = k3*(A^2) - k4*B"];
    model_def3["Y0eqs"] = ["A = 1", "B = k3*(A^2)/k4"];
    model_def3["Y0Sim"] = false;
    model_def3["tols"] = [1e-5,1e-5];
    model_def3["solver"] = "CVODE_BDF";

    oedms_def = Dict()
    oedms_def["Model_1"] = model_def2;
    oedms_def["Model_2"] = model_def3;
    oedms_def["Obs"] = ["B"];
    oedms_def["Theta_M1"] = [0.1 0.2 0.2 0.02; 0.11 0.21 0.21 0.021; 0.12 0.22 0.22 0.022];
    oedms_def["Theta_M2"] = [0.13 0.23 0.23 0.023; 0.14 0.24 0.24 0.024; 0.15 0.25 0.25 0.025];
    oedms_def["y0_M1"] = [0,0];
    oedms_def["y0_M2"] = [0,0];
    oedms_def["preInd_M1"] = [0.1];
    oedms_def["preInd_M2"] = [0.1];
    oedms_def["finalTime"] = [40];
    oedms_def["switchT"] = [0,20,30,40];
    oedms_def["tsamps"] = [0,5,10,15,20,25,30,35,40];
    oedms_def["fixedInp"] = [];
    oedms_def["fixedStep"] = [(1,[0])];
    oedms_def["equalStep"] = [[2,3]];
    oedms_def["plot"] = false;
    oedms_def["flag"] = "testoedms";
    oedms_def["uUpper"] = [1];
    oedms_def["uLower"] = [0];
    oedms_def["maxiter"] = 5;
    @test oedms_def == checkStructOEDMS(oedms_def);

    # Distance functions
    m1 = [10, 10];
    s1 = [1 0; 0 1];
    m2 = [11, 11];
    s2 = [1 0; 0 1];
    m3 = [15, 15];
    s3 = [1 0; 0 1];
    @test BhattacharyyaDist(m1, m1, s1, s1) == 0;
    @test BhattacharyyaDist(m1, m2, s1, s2) == 0.25;
    @test BhattacharyyaDist(m1, m3, s1, s3) == 6.25;

    @test EuclideanDist(m1, m1) == 0;
    @test EuclideanDist(m1, m2) == sqrt(2);
    @test EuclideanDist(m1, m3) == sqrt(50);

    # Generate utility script
    oedms_def = genOptimMSFuncts(oedms_def);
    @test isfile(string(pwd(), "\\ModelsFunctions\\", model_def2["NameF"], "_Model.jl"));
    @test isfile(string(pwd(), "\\ModelsFunctions\\", model_def3["NameF"], "_Model.jl"));
    @test isfile(string(pwd(), "\\Results\\", model_def2["NameF"], "_VS_", model_def3["NameF"], "_", today(),
                        "\\OEDModelSelectionScripts\\", model_def2["NameF"], "_VS_", model_def3["NameF"], "_OEDMS.jl"));

    # Bayes settings
    # Need to figure out if there is any test I can do here
    # opt = settingsBayesOpt(oedms_def);

    # Main function
    oedms_res, oedms_def = mainOEDMS(oedms_def);
    @test !isempty(oedms_res);
    @test typeof(oedms_res) <: Dict;
    @test isfile(string(pwd(), "\\Results\\", model_def2["NameF"], "_VS_", model_def3["NameF"], "_", today(),"\\OEDModelSelectResults_", oedms_def["flag"], ".jld"));
    @test oedms_res["uInpOpt"]["inp1"][1] == 0;
    @test oedms_res["uInpOpt"]["inp1"][2] == oedms_res["uInpOpt"]["inp1"][3];
    @test length(oedms_res["BestUtil"]) == 1;
    @test size(oedms_res["ConvCurv"])[1] == 5;
    @test size(oedms_res["Simul_M1"]) == size(oedms_res["Simul_M2"]);
    @test size(oedms_res["SimulObs_M1"])[2] == size(oedms_res["Simul_M1"])[2]-1;
    @test size(oedms_res["SimulObs_M2"])[2] == size(oedms_res["Simul_M2"])[2]-1;

    # Plot results
    plotOEDMSResults(oedms_res, oedms_def);
    @test isfile(string(oedms_def["savepath"], "\\Plot_OEDMSConvergence_", oedms_def["flag"], ".png"));
    rm(string(oedms_def["savepath"], "\\Plot_OEDMSConvergence_", oedms_def["flag"], ".png"));
    @test isfile(string(oedms_def["savepath"], "\\PlotOEDMSResults_Exp1_", oedms_def["flag"], ".png"));
    rm(string(oedms_def["savepath"], "\\PlotOEDMSResults_Exp1_", oedms_def["flag"], ".png"));

    rm(string(pwd(), "\\Results\\", model_def2["NameF"], "_VS_", model_def3["NameF"], "_", today(),"\\OEDModelSelectResults_", oedms_def["flag"], ".jld"));
    rm(string(pwd(), "\\ModelsFunctions\\", model_def2["NameF"], "_Model.jl"))
    rm(string(pwd(), "\\ModelsFunctions\\", model_def3["NameF"], "_Model.jl"))
    rm(string(pwd(), "\\Results\\", model_def2["NameF"], "_VS_", model_def3["NameF"], "_", today(),
                        "\\OEDModelSelectionScripts\\", model_def2["NameF"], "_VS_", model_def3["NameF"], "_OEDMS.jl"));
    rm(string(pwd(), "\\Results\\", model_def2["NameF"], "_VS_", model_def3["NameF"], "_", today(),
                        "\\OEDModelSelectionScripts"));
    rm(string(pwd(), "\\Results\\", model_def2["NameF"], "_VS_", model_def3["NameF"], "_", today()));

end

@testset "OEDModelCalibrationTests" begin

    # structure definition
    oedmc_def = defODEModelCalibrStruct();
    @test typeof(oedmc_def) <: Dict;
    entries1 = ["Model", "Obs", "Theta", "y0", "preInd", "finalTime", "switchT", "tsamps", "equalStep",
                "fixedInp", "fixedStep", "plot", "flag", "uUpper", "uLower", "maxiter", "util"];
    @test isempty(symdiff(entries1,keys(oedmc_def)));

    oedmc_def = Dict()
    oedmc_def["Model"] = model_def2;
    oedmc_def["Obs"] = ["B"];
    oedmc_def["Theta"] = [0.1 0.2 0.2 0.02; 0.11 0.21 0.21 0.021; 0.12 0.22 0.22 0.022];
    oedmc_def["y0"] = [0,0];
    oedmc_def["preInd"] = [0.1];
    oedmc_def["finalTime"] = [40];
    oedmc_def["switchT"] = [0,20,30,40];
    oedmc_def["tsamps"] = [0,5,10,15,20,25,30,35,40];
    oedmc_def["fixedInp"] = [];
    oedmc_def["fixedStep"] = [(1,[0])];
    oedmc_def["equalStep"] = [[2,3]];
    oedmc_def["plot"] = false;
    oedmc_def["flag"] = "testoedmc";
    oedmc_def["uUpper"] = [1];
    oedmc_def["uLower"] = [0];
    oedmc_def["maxiter"] = 5;
    oedmc_def["util"] = "perc";
    @test oedmc_def == checkStructOEDMC(oedmc_def);

    # Generate utility script
    oedmc_def = genOptimMCFuncts(oedmc_def);
    @test isfile(string(pwd(), "\\ModelsFunctions\\", model_def2["NameF"], "_Model.jl"));
    @test isfile(string(pwd(), "\\Results\\", model_def2["NameF"], "_", today(),
                        "\\OEDModelCalibrationScripts\\", model_def2["NameF"], "_OEDMC.jl"));

    # Main function
    oedmc_res, oedmc_def = mainOEDMC(oedmc_def);
    @test !isempty(oedmc_res);
    @test typeof(oedmc_res) <: Dict;
    @test isfile(string(pwd(), "\\Results\\", model_def2["NameF"], "_", today(),"\\OEDModelCalibrationResults_", oedmc_def["flag"], ".jld"));
    @test oedmc_res["uInpOpt"]["inp1"][1] == 0;
    @test oedmc_res["uInpOpt"]["inp1"][2] == oedmc_res["uInpOpt"]["inp1"][3];
    @test length(oedmc_res["BestUtil"]) == 1;
    @test size(oedmc_res["ConvCurv"])[1] == 5;
    @test size(oedmc_res["SimulObs_MC"])[2] == size(oedmc_res["Simul_MC"])[2]-1;

    # Plot results
    plotOEDMCResults(oedmc_res, oedmc_def);
    @test isfile(string(oedmc_def["savepath"], "\\Plot_OEDMCConvergence_", oedmc_def["flag"], ".png"));
    rm(string(oedmc_def["savepath"], "\\Plot_OEDMCConvergence_", oedmc_def["flag"], ".png"));
    @test isfile(string(oedmc_def["savepath"], "\\PlotOEDMCResults_Exp1_", oedmc_def["flag"], ".png"));
    rm(string(oedmc_def["savepath"], "\\PlotOEDMCResults_Exp1_", oedmc_def["flag"], ".png"));

    rm(string(pwd(), "\\ModelsFunctions\\", model_def2["NameF"], "_Model.jl"));
    rm(string(pwd(), "\\Results\\", model_def2["NameF"], "_", today(),
                        "\\OEDModelCalibrationScripts\\", model_def2["NameF"], "_OEDMC.jl"));
    rm(string(pwd(), "\\Results\\", model_def2["NameF"], "_", today(),"\\OEDModelCalibrationResults_", oedmc_def["flag"], ".jld"));

    rm(string(pwd(), "\\Results\\", model_def2["NameF"],  "_", today(),
                        "\\OEDModelCalibrationScripts"));
    rm(string(pwd(), "\\Results\\", model_def2["NameF"], "_", today()));
end

# rm(string(pwd(), "\\ModelsFunctions"))
# rm(string(pwd(), "\\Results"))
# rm(string(pwd(), "\\tmp"))

max_stars_repo_name: DavidGomezC/BOMBSbeta

max_issues_repo_name: DavidGomezC/BOMBSbeta

max_forks_repo_name: DavidGomezC/BOMBSbeta


text:
#prelude.jl

#some common imports used by almost all scripts:
using Base.MPFR: ROUNDING_MODE, big_ln2
using Base.Math: @horner, libm, nan_dom_err

ext: jl


text:
struct Graph
    n :: Int32 # |V|
    m :: Int32 # |E|
    V :: Array{Int32, 1} # V[i] = Real Index of node i
    E :: Array{Tuple{Int32, Int32, Int32, Float64}, 1} # (ID, u, v, w) in Edge Set
end

function readGraph(fileName, graphType) # Read graph from a file | graphType = [weighted, unweighted]
    # Initialized
    n = 0
    origin = Dict{Int32, Int32}()
    label = Dict{Int32, Int32}()
    edge = Set{Tuple{Int32, Int32, Float64}}()

    getid(x :: Int32) = haskey(label, x) ? label[x] : label[x] = n += 1

    open(fileName) do f1
        for line in eachline(f1)
            # Read origin data from file
            buf = split(line)
            u = parse(Int32, buf[1])
            v = parse(Int32, buf[2])
            if graphType == "weighted"
                w = parse(Float64, buf[3])
            else
                w = 1.0
            end
            if u == v
                continue
            end
            # Label the node
            u1 = getid(u)
            v1 = getid(v)
            origin[u1] = u
            origin[v1] = v
            # Store the edge
            if u1 > v1
                u1, v1 = v1, u1
            end
            push!(edge, (u1, v1, w))
        end
    end

    # Store data into the struct Graph
    m = length(edge)
    V = Array{Int32, 1}(undef, n)
    E = Array{Tuple{Int32, Int32, Int32, Float64}, 1}(undef, m)

    for i = 1 : n
        V[i] = origin[i]
    end

    ID = 0
    for (u, v, w) in edge
        ID = ID + 1
        E[ID] = (ID, u, v, w)
    end

    return Graph(n, m, V, E)
end

function getConnectedComponents(G)
    F = zeros(Int, G.n)
    foreach(i -> F[i] = i, 1 : G.n)

    find(x) = begin
        if F[x] != x
            F[x] = find(F[x])
        end
        return F[x]
    end

    for (ID, u, v, w) in G.E
        p = find(u)
        q = find(v)
        (p != q) ? F[p] = q : nothing
    end

    CC = zeros(Int, G.n)
    foreach(i -> CC[i] = find(i), 1 : G.n)

    return CC
end

function getBiconnectedComponents(G)
    dfn = zeros(Int, G.n)
    low = zeros(Int, G.n)
    cnt = 0
    cutPoint = zeros(Bool, G.n)
    bccid = zeros(Int, G.n)
    edgeC = zeros(Int, G.m)
    S = zeros(Int, 2*G.m)
    top = 0

    nbcc = 0
    bcc = Array{Array{Int, 1}, 1}(undef, G.n)
    g = Array{Array{Int, 1}, 1}(undef, G.n)
    gID = Array{Array{Int, 1}, 1}(undef, G.n)
    for i = 1 : G.n
        g[i] = []
        gID[i] = []
        bcc[i] = []
    end
    for (ID, u, v, w) in G.E
        push!(g[u], v)
        push!(gID[u], ID)
        push!(g[v], u)
        push!(gID[v], ID)
    end

    DFS(u, fa) = begin
        cnt += 1
        dfn[u] = cnt
        low[u] = cnt
        son = 0
        for i = 1 : size(g[u], 1)
            v = g[u][i]
            if dfn[v] == 0
                top += 1
                S[top] = gID[u][i]
                son += 1
                DFS(v, u)
                low[u] = min(low[u], low[v])
                if dfn[u] <= low[v]
                    cutPoint[u] = true
                    nbcc += 1
                    while true
                        x = S[top]
                        top -= 1
                        if bccid[G.E[x][2]] != nbcc
                            bccid[G.E[x][2]] = nbcc
                            push!(bcc[nbcc], G.E[x][2])
                        end
                        if bccid[G.E[x][3]] != nbcc
                            bccid[G.E[x][3]] = nbcc
                            push!(bcc[nbcc], G.E[x][3])
                        end
                        if x == gID[u][i]
                            break
                        end
                    end
                end
            elseif ((v != fa) && (dfn[v] < dfn[u]))
                top += 1
                S[top] = gID[u][i]
                low[u] = min(low[u], dfn[v])
            end
        end
        if (fa == 0) && (son <= 1)
            cutPoint[u] = false
        end
    end

    for i = 1 : G.n
        if dfn[i] == 0
            DFS(i, 0)
        end
    end

    timeStamp = 0
    fill!(bccid, 0)
    for i = 1 : nbcc
        timeStamp += 1
        for j = 1 : size(bcc[i], 1)
            bccid[bcc[i][j]] = timeStamp
        end
        for u in bcc[i]
            for j = 1 : size(g[u], 1)
                v = g[u][j]
                ei = gID[u][j]
                (bccid[v] == timeStamp) ? edgeC[ei] = timeStamp : nothing
            end
        end
    end
    for i = 1 : G.m
        if edgeC[i] == 0
            timeStamp += 1
            edgeC[i] = timeStamp
        end
    end

    cutList = Array{Int, 1}()
    for i = 1 : G.n
        if cutPoint[i] == true
            push!(cutList, i)
        end
    end

    return cutList, edgeC
end

function IsBiconnected(G)
    cutList, edgeC = getBiconnectedComponents(G)
    return size(cutList, 1) == 0
end

function getAdjacentList(G)
    g = Array{Array{Int, 1}, 1}(undef, G.n)
    foreach(i -> g[i] = [], 1 : G.n)
    for (ID, u, v, w) in G.E
        push!(g[u], v)
        push!(g[v], u)
    end
    return g
end

function SubGraph(G, Selected, EC)
    n = 0
    label = Dict{Int32, Int32}()
    nodeList = Array{Int32, 1}()

    for i = 1 : G.n
        if Selected[i]
            n += 1
            push!(nodeList, i)
            label[i] = n
        end
    end

    m = 0
    edges = Array{Tuple{Int32, Int32, Int32, Float64}, 1}()
    for (ID, u, v, w) in G.E
        if haskey(label, u) && haskey(label, v)
            m += 1
            push!(edges, (m, label[u], label[v], EC[ID]))
        end
    end

    V = zeros(Int, n)
    foreach(i -> V[i] = G.V[nodeList[i]], 1 : n)

    return Graph(n, m, V, edges)
end


text:
function Sockets.connect(fun::Function, args...)
    client = connect(args...)
    res = fun(client)
    close(client)
    res
end

@enum Magic::Int NEW_WORKER NEW_TASK RESULT ALIVE FINISHED

struct Tasker{XT,YT,Sampler<:AbstractSampler{XT,YT},
              Server,ActiveTasks,
              LifeSigns,AliveThreshold,Mutex}
    sampler::Sampler
    server::Server
    active_tasks::ActiveTasks
    life_signs::LifeSigns
    alive_threshold::AliveThreshold
    mutex::Mutex
end

Tasker(sampler::AbstractSampler, server, alive_threshold::Number) =
    Tasker(sampler, server, Int[],
           Vector{typeof(now())}(), alive_threshold,
           Threads.SpinLock())

function Base.show(io::IO, t::Tasker)
    s = t.sampler
    n = length(s)
    write(io, "$(n) ($(count(s.done)) done) samples Tasker listening on ")
    show(io, t.server)
end

function Base.show(io::IO, ::MIME"text/plain", t::Tasker)
    show(io, t)
    println(io)

    tn = now()
    pretty_table(io,
                 hcat(eachindex(t.active_tasks), t.active_tasks,
                      tn .- t.life_signs),
                 header=["Worker id", "Current task", "Alive"],
                 hlines=[1], vlines=[])
    println(io)
end

function serve_tasks!(tasker::Tasker{XT,YT};
                      plot_fun::Union{Function,Nothing}=nothing) where {XT,YT}
    server = tasker.server
    s = tasker.sampler
    load_samples!(s)

    pf = () -> begin
        if !isnothing(plot_fun)
            sel = done(s)
            plot_fun(s.x[sel], s.y[sel])
        end
    end
    pf()

    # At the moment, we run the task server synchronously, i.e. we do
    # not allow any workers in the same Julia process.
    while !isdone(s) || any(!iszero, tasker.active_tasks)
        sock = accept(server)

        magic = read(sock, Magic)
        tasks_left = !isdone(s)
        @info "Got connected" sock magic tasks_left islocked(tasker.mutex)
        if magic == NEW_WORKER
            if tasks_left
                lock(tasker.mutex) do
                    push!(tasker.active_tasks, 0)
                    push!(tasker.life_signs, now())
                end
                worker_id = length(tasker.active_tasks)
                @info "Assigning new worker id #$(worker_id)"
                write(sock, NEW_WORKER, worker_id)
            else
                write(sock, FINISHED)
            end
        elseif magic == NEW_TASK
            worker_id = read(sock, Int)
            if tasks_left
                lock(tasker.mutex) do
                    nd = filter(∉(tasker.active_tasks), not_done(s))
                    if isempty(nd)
                        @warn "Could not find next sample, weird"
                        tasker.active_tasks[worker_id] = 0
                        write(sock, FINISHED)
                    else
                        i = first(nd)
                        tasker.active_tasks[worker_id] = i
                        tasker.life_signs[worker_id] = now()
                        x = get_sample!(s, i)
                        write(sock, NEW_TASK, i, x)
                        @info "Worker #$(worker_id) asks for work, gets it" i x
                    end
                end
            else
                tasker.active_tasks[worker_id] = 0
                write(sock, FINISHED)
                @info "Worker #$(worker_id) asks for work, none left"
            end
        elseif magic == RESULT
            worker_id = read(sock, Int)
            x = read(sock, XT)
            y = read(sock, YT)
            lock(tasker.mutex) do
                i = tasker.active_tasks[worker_id]
                @info "Worker #$(worker_id) tells us the following result" i x y
                s[i] = (x,y)
                tasker.life_signs[worker_id] = now()
                save_samples!(s)
            end
            pf()
        elseif magic == ALIVE
            worker_id = read(sock, Int)
            @info "Worker #$(worker_id) is still alive, nice"
            lock(tasker.mutex) do
                tasker.life_signs[worker_id] = now()
            end
        end
        horizontal_line(color=:red)
        display(tasker)
        horizontal_line()
        horizontal_line(color=:yellow)
    end
    close(server)
end

function task_farm(fun::Function, sampler::AbstractSampler{XT,YT};
                   host=localhost, port=2000,
                   alive_sleep=60, alive_threshold=3alive_sleep,
                   kwargs...) where {XT,YT}
    server = try
        if host == localhost || hostname() == host
            listen(host == localhost ? localhost : IPv4(0), port)
        else
            nothing
        end
    catch IOError
        nothing
    end

    if !isnothing(server)
        @info "We are the server, yay!"
        tasker = Tasker(sampler, server, alive_threshold)
        serve_tasks!(tasker; kwargs...)
    else
        @info "We are a measly worker"
        worker_id = connect(host, port) do client
            write(client, NEW_WORKER)
            response = read(client, Magic)
            if response == NEW_WORKER
                read(client, Int)
            elseif response == FINISHED
                @info "There's not even anything to do"
            end
        end

        if !isnothing(worker_id)
            @info "We are measly worker #$(worker_id)"

            @sync begin
                # @async while isopen(client)
                #     lock(mutex) do
                #         write(client, ALIVE, worker_id)
                #     end
                #     sleep(alive_sleep)
                # end

                while true
                    horizontal_line(color=:green)
                    ix = connect(host, port) do client
                        write(client, NEW_TASK, worker_id)
                        response = read(client, Magic)
                        @info "Requested work" response
                        if response == NEW_TASK
                            read(client, Int), read(client, XT)
                        elseif response == FINISHED
                            @info "We must go home"
                        end
                    end
                    isnothing(ix) && break

                    i,x = ix
                    @info "We are asked to work" i x
                    y = fun(i, x)

                    connect(host, port) do client
                        write(client, RESULT, worker_id, x, y)
                    end
                end
            end
        end
    end
    horizontal_line(color=:blue)
    @info "We are done"
end

export task_farm


text:
# High-level tests

makearray(D, val) = fill(val, ntuple(d -> 1, D))

if comm_rank == 0
    const dirname2 = Filesystem.mktempdir(; cleanup=true)
    const filename2 = "$dirname2/test.bp"

    @testset "High-level write tests " begin
        file = adios_open_serial(filename2, mode_write)

        adios_define_attribute(file, "a1", float(π))
        adios_define_attribute(file, "a2", [float(π)])
        adios_define_attribute(file, "a3", [float(π), 0])

        adios_put!(file, "v1", float(ℯ))
        adios_put!(file, "v3", makearray(1, float(ℯ)))
        adios_put!(file, "v4", makearray(2, float(ℯ)))
        adios_put!(file, "v5", makearray(3, float(ℯ)))
        adios_put!(file, "g1/v6", makearray(4, float(ℯ)))
        adios_put!(file, "g1/g2/v7", makearray(5, float(ℯ)))

        @test shapeid(inquire_variable(file.io, "v1")) == shapeid_local_value
        @test shapeid(inquire_variable(file.io, "v3")) == shapeid_local_array
        @test shapeid(inquire_variable(file.io, "v4")) == shapeid_local_array
        @test shapeid(inquire_variable(file.io, "v5")) == shapeid_local_array
        @test shapeid(inquire_variable(file.io, "g1/v6")) == shapeid_local_array
        @test shapeid(inquire_variable(file.io, "g1/g2/v7")) ==
              shapeid_local_array

        adios_define_attribute(file, "v4/a4", float(π))
        adios_define_attribute(file, "v5", "a5", [float(π)])
        adios_define_attribute(file, "g1/v6", "a6", [float(π), 0])

        adios_perform_puts!(file)
        close(file)
    end
    # run(`/Users/eschnett/src/CarpetX/Cactus/view/bin/bpls -aD $filename2`)

    @testset "High-level read tests " begin
        file = adios_open_serial(filename2, mode_read)

        @test Set(adios_subgroup_names(file, "")) == Set(["g1"])
        @test_broken Set(adios_subgroup_names(file, "g1")) == Set(["g2"])
        @test Set(adios_subgroup_names(file, "g1")) == Set(["/g2"]) # don't want this

        @test Set(adios_all_attribute_names(file)) ==
              Set(["a1", "a2", "a3", "v4/a4", "v5/a5", "g1/v6/a6"])
        @test Set(adios_group_attribute_names(file, "g1")) == Set()
        @test Set(adios_group_attribute_names(file, "g1/v6")) ==
              Set(["g1/v6/a6"])

        @test adios_attribute_data(file, "a1") == float(π)
        @test adios_attribute_data(file, "a2") == float(π)
        @test adios_attribute_data(file, "a3") == [float(π), 0]
        @test adios_attribute_data(file, "v4", "a4") == float(π)
        @test adios_attribute_data(file, "v5/a5") == float(π)
        @test adios_attribute_data(file, "g1/v6", "a6") == [float(π), 0]

        @test Set(adios_all_variable_names(file)) ==
              Set(["v1", "v3", "v4", "v5", "g1/v6", "g1/g2/v7"])
        @test Set(adios_group_variable_names(file, "g1")) == Set(["g1/v6"])
        @test Set(adios_group_variable_names(file, "g1/g2")) ==
              Set(["g1/g2/v7"])

        # Local values are converted to global arrays
        @test shapeid(inquire_variable(file.io, "v1")) == shapeid_global_array
        @test shapeid(inquire_variable(file.io, "v3")) == shapeid_local_array
        @test shapeid(inquire_variable(file.io, "v4")) == shapeid_local_array
        @test shapeid(inquire_variable(file.io, "v5")) == shapeid_local_array
        @test shapeid(inquire_variable(file.io, "g1/v6")) == shapeid_local_array
        @test shapeid(inquire_variable(file.io, "g1/g2/v7")) ==
              shapeid_local_array

        @test ndims(inquire_variable(file.io, "v1")) == 1
        @test ndims(inquire_variable(file.io, "v3")) == 1
        @test ndims(inquire_variable(file.io, "v4")) == 2
        @test ndims(inquire_variable(file.io, "v5")) == 3
        @test ndims(inquire_variable(file.io, "g1/v6")) == 4
        @test ndims(inquire_variable(file.io, "g1/g2/v7")) == 5

        @test shape(inquire_variable(file.io, "v1")) == (1,)
        @test shape(inquire_variable(file.io, "v3")) ≡ nothing
        @test shape(inquire_variable(file.io, "v4")) ≡ nothing
        @test shape(inquire_variable(file.io, "v5")) ≡ nothing
        @test shape(inquire_variable(file.io, "g1/v6")) ≡ nothing
        @test shape(inquire_variable(file.io, "g1/g2/v7")) ≡ nothing

        @test start(inquire_variable(file.io, "v1")) == (0,)
        @test start(inquire_variable(file.io, "v3")) ≡ nothing
        @test start(inquire_variable(file.io, "v4")) ≡ nothing
        @test start(inquire_variable(file.io, "v5")) ≡ nothing
        @test start(inquire_variable(file.io, "g1/v6")) ≡ nothing
        @test start(inquire_variable(file.io, "g1/g2/v7")) ≡ nothing

        @test count(inquire_variable(file.io, "v1")) == (1,)
        @test count(inquire_variable(file.io, "v3")) == (1,)
        @test count(inquire_variable(file.io, "v4")) == (1, 1)
        @test count(inquire_variable(file.io, "v5")) == (1, 1, 1)
        @test count(inquire_variable(file.io, "g1/v6")) == (1, 1, 1, 1)
        @test count(inquire_variable(file.io, "g1/g2/v7")) == (1, 1, 1, 1, 1)

        v1 = adios_get(file, "v1")
        @test !isready(v1)
        @test fetch(v1) == fill(float(ℯ), 1)
        @test isready(v1)
        v3 = adios_get(file, "v3")
        v4 = adios_get(file, "v4")
        @test !isready(v3)
        @test fetch(v3) == makearray(1, float(ℯ))
        @test fetch(v4) == makearray(2, float(ℯ))
        @test isready(v3)
        v5 = adios_get(file, "v5")
        v6 = adios_get(file, "g1/v6")
        v7 = adios_get(file, "g1/g2/v7")
        @test !isready(v5)
        adios_perform_gets(file)
        @test isready(v5)
        @test fetch(v5) == makearray(3, float(ℯ))
        @test fetch(v6) == makearray(4, float(ℯ))
        @test fetch(v7) == makearray(5, float(ℯ))
        close(file)
    end
end


text:
module DecimalNumbers

export Decimal

importall Base.Operators

# https://books.google.com/books?id=pNdiJMvPoZMC&pg=PA11&lpg=PA11&dq=IEEE+BID+format&source=bl&ots=-GbBYEh6Sa&sig=7_LN575lYEa8443zTGtFEOuLv1o&hl=en&sa=X&ved=0ahUKEwiqudHy0fLSAhUP62MKHbigDkwQ6AEIPjAI#v=onepage&q=IEEE%20BID%20format&f=false

immutable Decimal{T <: Integer}
    value::T
    exp::T
end

Decimal{T, T2}(value::T, exp::T2) = Decimal(promote(value, exp)...)

# normalize
function normalize{T <: Integer}(x::T)
    value = x
    exp = Z = zero(T)
    for d in digits(x)
        if d == Z
            value = div(value, T(10))
            exp += T(1)
        else
            break
        end
    end
    return value, exp
end

function Decimal{T}(d::Decimal{T})
    v, e = normalize(d.value)
    return Decimal(v, e + d.exp)
end

Base.convert{T, T2}(::Type{Decimal{T}}, d::Decimal{T2}) = Decimal(T(d.value), T(d.exp))

# from int
Base.convert{T}(::Type{Decimal{T}}, x) = convert(Decimal, x)
function Base.convert{T <: Integer}(::Type{Decimal}, x::T)
    x == T(0) && return Decimal(T(0), T(0))
    return Decimal(normalize(x)...)
end

# to int
function Base.convert{T <: Integer}(::Type{T}, d::Decimal)
    d.exp < 0 && throw(InexactError())
    return T(d.value * T(10)^d.exp)
end

Base.zero{T}(::Union{Decimal{T}, Type{Decimal{T}}}) = Decimal(zero(T))
Base.one{T}(::Union{Decimal{T}, Type{Decimal{T}}}) = Decimal(one(T))

float2int(::Type{BigFloat}) = BigInt
float2int(::Type{Float64}) = Int64
float2int(::Type{Float32}) = Int32
float2int(::Type{Float16}) = Int16

# from float
function Base.convert{T <: AbstractFloat}(::Type{Decimal}, x::T)
    # easy if float is int
    trunc(x) == x && return Decimal(float2int(T)(x))
    # otherwise, go string route for now
    return parse(Decimal{float2int(T)}, string(x))
end

# to float
Base.convert{T <: AbstractFloat}(::Type{T}, d::Decimal) = T(d.value * exp10(d.exp))

function ==(a::Decimal, b::Decimal)
    aa, bb = _scale(a, b)
    return aa.value == bb.value && aa.exp == bb.exp
end

=={T}(a::Decimal, b::T) = ==(promote(a, b)...)
=={T}(a::T, b::Decimal) = ==(promote(a, b)...)

function <(a::Decimal, b::Decimal)
    aa, bb = _scale(a, b)
    return aa.value < bb.value
end

<{T}(a::Decimal, b::T) = <(promote(a, b)...)
<{T}(a::T, b::Decimal) = <(promote(a, b)...)

const ZERO = UInt8('0')
const DOT = UInt8('.')
const MINUS = UInt8('-')
const PLUS = UInt8('+')

"""
Parse a Decimal from a string. Supports decimals of the following form:

  * "101"
  * "101."
  * "101.0"
  * "1.01"
  * ".101"
  * "0.101"
  * "0.0101"
"""
function Base.parse{T}(::Type{Decimal{T}}, str::String)
    str = strip(str)
    bytes = Vector{UInt8}(str)
    value = exp = zero(T)
    frac = neg = false
    for i = 1:length(bytes)
        b = bytes[i]
        if b == MINUS
            neg = true
        elseif b == PLUS
            continue
        elseif b == DOT
            frac = true
            continue
        else
            value *= T(10)
            value += T(b - ZERO)
            frac && (exp -= T(1))
        end
    end
    for i = length(bytes):-1:1
        b = bytes[i]
        if b == ZERO
            if exp < 0
                value = div(value, T(10))
            end
            exp += 1
        else
            break
        end
    end
    return Decimal(ifelse(neg, T(-1), T(1)) * value, exp)
end

function Base.show(io::IO, d::Decimal)
    print(io, "dec\"")
    if d.value == 0
        print(io, "0.0")
    else
        sn = sign(d.value) < 0 ? "-" : ""
        str = string(abs(d.value))
        if d.exp == 0
            print(io, sn, str, ".0")
        else
            if d.exp > 0
                print(io, sn, rpad(str, length(str) + d.exp, '0'), ".0")
            else
                d = length(str) + d.exp
                if d == 0
                    print(io, sn, "0.", str)
                elseif d > 0
                    print(io, sn, str[1:d], ".", str[d+1:end])
                else
                    print(io, sn, "0.", "0"^abs(d), str)
                end
            end
        end
    end
    print(io, '"')
    return
end

# math
-(d::Decimal) = Decimal(-d.value, d.exp)
Base.abs(d::Decimal) = Decimal(abs(d.value), d.exp)

# 10, 100
# (1, 1), (1, 2)
# (1, 1), (10, 1)

# 1.1, 0.001
# (11, -1), (1, -3)
# (1100, -3), (1, -3)

# 10, 0.1
# (1, 1), (1, -1)
# (100, -1), (1, -1)

# scales two decimals to the same exp
_scale{T, T2}(a::Decimal{T}, b::Decimal{T2}) = _scale(promote(a, b)...)
function _scale{T}(a::Decimal{T}, b::Decimal{T})
    a.exp == b.exp && return a, b
    if a.exp < b.exp
        return a, Decimal(b.value * T(10)^(abs(b.exp - a.exp)), a.exp)
    else
        return Decimal(a.value * T(10)^(abs(a.exp - b.exp)), b.exp), b
    end
end

function +(a::Decimal, b::Decimal)
    a2, b2 = _scale(a, b)
    return Decimal(Decimal(a2.value + b2.value, a2.exp))
end

function -(a::Decimal, b::Decimal)
    a2, b2 = _scale(a, b)
    return Decimal(Decimal(a2.value - b2.value, a2.exp))
end

function *(a::Decimal, b::Decimal)
    exp = a.exp + b.exp
    val = a.value * b.value
    return Decimal(Decimal(val, exp))
end

maxprec{T <: Integer}(::Type{T}) = T(length(string(typemax(T))))

function /{T}(a::Decimal{T}, b::Decimal{T})
    b.value == 0 && throw(DivideError())
    # scale num up to max precision
    scale = maxprec(T) - T(length(string(a.value)))
    aa = a.value * (widen(T)(10) ^ scale)
    # simulate division
    q, r = divrem(aa, b.value)
    # return scaled results
    return Decimal(Decimal(T(q + ifelse(div(r * 10, b.value) > 5, 1, 0)), a.exp - b.exp - scale))
end

Base.promote_rule{T, T2}(::Type{Decimal{T}}, ::Type{Decimal{T2}}) = Decimal{promote_type(T, T2)}
Base.promote_rule{T, TI <: Integer}(::Type{Decimal{T}}, ::Type{TI}) = Decimal{promote_type(T, TI)}
Base.promote_rule{T, TF <: AbstractFloat}(::Type{Decimal{T}}, ::Type{TF}) = Decimal{float2int(promote_type(T, TF))}

# TODO:
 # rounding, trunc, floor, ceil
 # maybe:
   # equality: isapprox?
   # ranges
   # fld, mod, rem, divrem, divmod, mod1,
   # fma, muladd
   # shifts?
   # trig functions
   # log, log2, log10
   # exp, ldexp, modf
   # sqrt
   # special functions

end # module

max_stars_repo_name: quinnj/DecimalNumbers

max_issues_repo_name: quinnj/DecimalNumbers

max_forks_repo_name: quinnj/DecimalNumbers


text:
using HDF5,JLD,KUparser,Base.Test
@date d = load("conll07.tst.jld4")
@show ndeps = length(d["deprel"])
@date corpus = d["corpus"]
@date feats = Flist.acl11eager

# We try each arctype with rparse and oparse
# They will complain if movecosts or nmoves is buggy
for pt in (ArcEager13, ArcEagerR1, ArcHybrid13, ArcHybridR1)
    @show pt
    @date r = rparse(pt, corpus, ndeps)
    @show evalparse(r, corpus)
    @date (p,x,y) = oparse(pt, corpus, ndeps, feats)
    @show evalparse(p, corpus)
    p1 = pt(1, ndeps)
    @test @show size(x) == KUparser.xsize(p1, corpus, feats)
    @test @show size(y) == KUparser.ysize(p1, corpus)
end


text:
@testset "reified fixed to 1 or constraint violated or solved?" begin
    m = Model(optimizer_with_attributes(CS.Optimizer, "no_prune" => true, "logging" => []))
    @variable(m, 0 <= x[1:2] <= 5, Int)
    @constraint(m, x in CS.AllDifferent() || sum(x) > 2 || x[1] > 1)
    optimize!(m)
    com = CS.get_inner_model(m)

    variables = com.search_space
    constraint = com.constraints[1]
    @test CS.fix!(com, variables[constraint.indices[1]], 1; check_feasibility = false)
    @test CS.fix!(com, variables[constraint.indices[2]], 1; check_feasibility = false)
    @test CS.is_constraint_violated(com, constraint, constraint.fct, constraint.set)

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

    m = Model(optimizer_with_attributes(CS.Optimizer, "no_prune" => true, "logging" => []))
    @variable(m, b >= 1, Bin)
    @variable(m, 0 <= x[1:2] <= 5, Int)
    @constraint(m, b := {x in CS.AllDifferent() || sum(x) <= 2 || 2x[1]+2x[2] > 8})
    optimize!(m)
    com = CS.get_inner_model(m)

    variables = com.search_space
    constraint = com.constraints[1]
    @test CS.fix!(com, variables[constraint.indices[2]], 2; check_feasibility = false)
    @test CS.fix!(com, variables[constraint.indices[3]], 2; check_feasibility = false)
    @test CS.is_constraint_violated(com, constraint, constraint.fct, constraint.set)

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

    m = Model(optimizer_with_attributes(CS.Optimizer, "no_prune" => true, "logging" => []))
    @variable(m, b >= 1, Bin)
    @variable(m, 0 <= x[1:2] <= 5, Int)
    @constraint(m, b := {x in CS.AllDifferent() || sum(x) > 2})
    optimize!(m)
    com = CS.get_inner_model(m)

    variables = com.search_space
    constraint = com.constraints[1]
    or_constraint = com.constraints[1].inner_constraint
    @test CS.fix!(com, variables[or_constraint.indices[1]], 1; check_feasibility = false)
    @test CS.fix!(com, variables[or_constraint.indices[2]], 1; check_feasibility = false)
    @test CS.is_constraint_violated(com, or_constraint, or_constraint.fct, or_constraint.set)

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

    m = Model(optimizer_with_attributes(CS.Optimizer, "no_prune" => true, "logging" => []))
    @variable(m, b >= 1, Bin)
    @variable(m, 0 <= x[1:2] <= 5, Int)
    @constraint(m, b := {sum(x) > 2 || x in CS.AllDifferent()})
    optimize!(m)
    com = CS.get_inner_model(m)

    variables = com.search_space
    constraint = com.constraints[1]
    or_constraint = com.constraints[1].inner_constraint
    @test CS.fix!(com, variables[or_constraint.indices[1]], 1; check_feasibility = false)
    @test CS.fix!(com, variables[or_constraint.indices[2]], 1; check_feasibility = false)
    @test CS.is_constraint_violated(com, or_constraint, or_constraint.fct, or_constraint.set)

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

    m = Model(optimizer_with_attributes(CS.Optimizer, "no_prune" => true, "logging" => []))
    @variable(m, b >= 1, Bin)
    @variable(m, 0 <= x[1:2] <= 5, Int)
    @constraint(m, b := {x in CS.AllDifferent() || sum(x) <= 2})
    optimize!(m)
    com = CS.get_inner_model(m)

    variables = com.search_space
    constraint = com.constraints[1]
    @test CS.fix!(com, variables[constraint.indices[2]], 1; check_feasibility = false)
    @test CS.fix!(com, variables[constraint.indices[3]], 1; check_feasibility = false)
    @test CS.is_constraint_solved(com, constraint, constraint.fct, constraint.set)
end

@testset "Indicator fixed to 1 or constraint violated or solved?" begin
    m = Model(optimizer_with_attributes(CS.Optimizer, "no_prune" => true, "logging" => []))
    @variable(m, b >= 1, Bin)
    @variable(m, 0 <= x[1:2] <= 5, Int)
    @constraint(m, b => {x in CS.AllDifferent() || sum(x) > 2 || x[1] > 1 })
    optimize!(m)
    com = CS.get_inner_model(m)

    variables = com.search_space
    constraint = com.constraints[1]
    @test CS.fix!(com, variables[constraint.indices[2]], 1; check_feasibility = false)
    @test CS.fix!(com, variables[constraint.indices[3]], 1; check_feasibility = false)
    @test CS.is_constraint_violated(com, constraint, constraint.fct, constraint.set)

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

    m = Model(optimizer_with_attributes(CS.Optimizer, "no_prune" => true, "logging" => []))
    @variable(m, b >= 1, Bin)
    @variable(m, 0 <= x[1:2] <= 5, Int)
    @constraint(m, b => {x in CS.AllDifferent() || sum(x) <= 2 || 2x[1]+2x[2] > 8})
    optimize!(m)
    com = CS.get_inner_model(m)

    variables = com.search_space
    constraint = com.constraints[1]
    @test CS.fix!(com, variables[constraint.indices[2]], 2; check_feasibility = false)
    @test CS.fix!(com, variables[constraint.indices[3]], 2; check_feasibility = false)
    @test CS.is_constraint_violated(com, constraint, constraint.fct, constraint.set)

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

    m = Model(optimizer_with_attributes(CS.Optimizer, "no_prune" => true, "logging" => []))
    @variable(m, b >= 1, Bin)
    @variable(m, 0 <= x[1:2] <= 5, Int)
    @constraint(m, b => {x in CS.AllDifferent() || sum(x) > 2})
    optimize!(m)
    com = CS.get_inner_model(m)

    variables = com.search_space
    constraint = com.constraints[1]
    or_constraint = com.constraints[1].inner_constraint
    @test CS.fix!(com, variables[or_constraint.indices[1]], 1; check_feasibility = false)
    @test CS.fix!(com, variables[or_constraint.indices[2]], 1; check_feasibility = false)
    @test CS.is_constraint_violated(com, or_constraint, or_constraint.fct, or_constraint.set)

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

    m = Model(optimizer_with_attributes(CS.Optimizer, "no_prune" => true, "logging" => []))
    @variable(m, b >= 1, Bin)
    @variable(m, 0 <= x[1:2] <= 5, Int)
    @constraint(m, b => {sum(x) > 2 || x in CS.AllDifferent()})
    optimize!(m)
    com = CS.get_inner_model(m)

    variables = com.search_space
    constraint = com.constraints[1]
    or_constraint = com.constraints[1].inner_constraint
    @test CS.fix!(com, variables[or_constraint.indices[1]], 1; check_feasibility = false)
    @test CS.fix!(com, variables[or_constraint.indices[2]], 1; check_feasibility = false)
    @test CS.is_constraint_violated(com, or_constraint, or_constraint.fct, or_constraint.set)

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

    m = Model(optimizer_with_attributes(CS.Optimizer, "no_prune" => true, "logging" => []))
    @variable(m, b >= 1, Bin)
    @variable(m, 0 <= x[1:2] <= 5, Int)
    @constraint(m, b => {x in CS.AllDifferent() || sum(x) <= 2})
    optimize!(m)
    com = CS.get_inner_model(m)

    variables = com.search_space
    constraint = com.constraints[1]
    @test CS.fix!(com, variables[constraint.indices[2]], 1; check_feasibility = false)
    @test CS.fix!(com, variables[constraint.indices[3]], 1; check_feasibility = false)
    @test CS.is_constraint_solved(com, constraint, constraint.fct, constraint.set)
end

@testset "reified or constraint prune_constraint!" begin
    m = Model(optimizer_with_attributes(CS.Optimizer, "no_prune" => true, "logging" => []))
    @variable(m, b >= 1, Bin)
    @variable(m, 0 <= x[1:2] <= 5, Int)
    @constraint(m, b := {2x[1]+x[2] >= 3 || x[1] > 1})
    optimize!(m)

    com = CS.get_inner_model(m)

    variables = com.search_space
    constraint = com.constraints[1]

    constr_indices = constraint.indices
    @test CS.fix!(com, variables[constraint.indices[2]], 0; check_feasibility = false)
    CS.changed!(com, constraint, constraint.fct, constraint.set)
    @test CS.prune_constraint!(com, constraint, constraint.fct, constraint.set)
    for v in 0:2
        @test !CS.has(variables[constraint.indices[3]], v)
    end

    # prune rhs
    m = Model(optimizer_with_attributes(CS.Optimizer, "no_prune" => true, "logging" => []))
    @variable(m, b >= 1, Bin)
    @variable(m, 0 <= x[1:2] <= 5, Int)
    @constraint(m, b := {x[1] > 1 || 2x[1]+x[2] >= 3})
    optimize!(m)

    com = CS.get_inner_model(m)

    variables = com.search_space
    constraint = com.constraints[1]

    constr_indices = constraint.indices
    @test CS.fix!(com, variables[constraint.indices[2]], 0; check_feasibility = false)
    CS.changed!(com, constraint, constraint.fct, constraint.set)
    @test CS.prune_constraint!(com, constraint, constraint.fct, constraint.set)
    for v in 0:2
        @test !CS.has(variables[constraint.indices[4]], v)
    end
end

@testset "reified or constraint still_feasible" begin
    m = Model(optimizer_with_attributes(CS.Optimizer, "no_prune" => true, "logging" => []))
    @variable(m, b >= 1, Bin)
    @variable(m, 0 <= x[1:2] <= 5, Int)
    @variable(m, 0 <= y <= 5, Int)
    @constraint(m, b := { x in CS.AllDifferent() || y < 1})
    optimize!(m)

    com = CS.get_inner_model(m)

    variables = com.search_space
    constraint = com.constraints[1]

    or_constraint =constraint.inner_constraint

    constr_indices = or_constraint.indices
    @test CS.fix!(com, variables[constraint.indices[2]], 0; check_feasibility = false)
    @test CS.fix!(com, variables[constraint.indices[3]], 0; check_feasibility = false)
    @test !CS.still_feasible(com, or_constraint, or_constraint.fct, or_constraint.set, constr_indices[3], 1)

    # swap lhs and rhs
    m = Model(optimizer_with_attributes(CS.Optimizer, "no_prune" => true, "logging" => []))
    @variable(m, b >= 1, Bin)
    @variable(m, 0 <= x[1:2] <= 5, Int)
    @variable(m, 0 <= y <= 5, Int)
    @constraint(m, b := { y < 1 || x in CS.AllDifferent()})
    optimize!(m)

    com = CS.get_inner_model(m)

    variables = com.search_space
    constraint = com.constraints[1]

    or_constraint =constraint.inner_constraint

    constr_indices = or_constraint.indices
    @test CS.fix!(com, variables[constr_indices[2]], 0; check_feasibility = false)
    @test CS.fix!(com, variables[constr_indices[3]], 0; check_feasibility = false)
    @test !CS.still_feasible(com, or_constraint, or_constraint.fct, or_constraint.set, constr_indices[1], 1)
end


text:
using Test, DataFrames, CSV, BDisposal
println("Testing BDisposal...")

println("Testing BDisposal on js-data...")
# Byung M. Jeon, Robin C. Sickles, "The Role of Environmental Factors in
# Growth Accounting: a Nonparametric Analysis", Journal of Applied
# Econometrics, Vol. 19, No. 5, 2004, pp. 567-591.

data = CSV.read(joinpath(dirname(pathof(BDisposal)),"..","test","data","js-data","oecd.txt"),DataFrame; delim=' ',ignorerepeated=true,copycols=true,header=false)
rename!(data,[:ccid,:year,:gdp,:co2,:capital,:labour,:energy])
goodInputsLabels  = ["capital","labour"]
badInputsLabels   = ["energy"]
goodOutputsLabels = ["gdp"]
badOutputsLabels  = ["co2"]
sort!(data, [:year, :ccid]) # sort data by period and dmu
periods = unique(data.year)
dmus    = unique(data.ccid)

nGI, nBI, nGO, nBO, nPer, nDMUs,  = length(goodInputsLabels), length(badInputsLabels), length(goodOutputsLabels), length(badOutputsLabels), length(periods),length(dmus)

gI = Array{Float64}(undef, (nDMUs,nGI,nPer))
bI = Array{Float64}(undef, (nDMUs,nBI,nPer))
gO = Array{Float64}(undef, (nDMUs,nGO,nPer))
bO = Array{Float64}(undef, (nDMUs,nBO,nPer))

for (p,period) in enumerate(periods)
    periodData = data[data.year .== period,:]
    gI[:,:,p] = Matrix{Float64}(periodData[:,goodInputsLabels])
    if nBI > 0
        bI[:,:,p] =
         Matrix{Float64}(periodData[:,badInputsLabels])
    end
    gO[:,:,p] = Matrix{Float64}(periodData[:,goodOutputsLabels])
    bO[:,:,p] = Matrix{Float64}(periodData[:,badOutputsLabels])
end

################################################################################
# Testing prodIndex

oecdAnalysis  = prodIndex(gI,gO,bO,bI;
                   retToScale="variable",prodStructure="multiplicative",convexAssumption=true)

@test oecdAnalysis.prodIndexes[1,10] == 1.02891196267105

@test isapprox(oecdAnalysis.prodIndexes_G .* oecdAnalysis.prodIndexes_B, oecdAnalysis.prodIndexes, atol=0.000001)
decomp = isapprox.(oecdAnalysis.prodIndexes_T .* oecdAnalysis.prodIndexes_E .* oecdAnalysis.prodIndexes_S, oecdAnalysis.prodIndexes, atol=0.000001)
@test all(skipmissing(decomp)) == true

oecdAnalysisA  = prodIndex(gI,gO,bO,bI;
                   retToScale="variable",prodStructure="addittive",convexAssumption=true)

@test oecdAnalysisA.prodIndexes[1,10] == 0.02111291275337509
@test isapprox(oecdAnalysisA.prodIndexes_G .+ oecdAnalysisA.prodIndexes_B, oecdAnalysisA.prodIndexes, atol=0.000001)
decomp =  isapprox.(oecdAnalysisA.prodIndexes_T .+ oecdAnalysisA.prodIndexes_E .+ oecdAnalysisA.prodIndexes_S, oecdAnalysisA.prodIndexes, atol=0.000001)
@test all(skipmissing(decomp)) == true

oecdAnalysis_nc  = prodIndex(gI,gO,bO,bI;
                   retToScale="variable",prodStructure="multiplicative",convexAssumption=false)

@test isapprox(oecdAnalysis_nc.prodIndexes_G .* oecdAnalysis_nc.prodIndexes_B, oecdAnalysis_nc.prodIndexes, atol=0.000001)

oecdAnalysis_ncA  = prodIndex(gI,gO,bO,bI;
                   retToScale="variable",prodStructure="additive",convexAssumption=false)

@test isapprox(oecdAnalysis_ncA.prodIndexes_G .+ oecdAnalysis_ncA.prodIndexes_B, oecdAnalysis_ncA.prodIndexes, atol=0.000001)

@test isapprox(log.(oecdAnalysis_nc.prodIndexes), oecdAnalysis_ncA.prodIndexes, atol=0.05)

################################################################################
# Testing efficiencyScore using the Airports dataset..

println("Testing BDisposal on airport data...")
airportData = CSV.read(joinpath(dirname(pathof(BDisposal)),"..","test","data","airports.csv"),DataFrame; delim=';',copycols=true)
airportGoodInputs  = ["employees","totalCosts"]
airportBadInputs   = []
airportGoodOutputs = ["passengers"]
airportBadOutputs  = ["co2emissions"]
sort!(airportData, [:period, :dmu]) # sort data by period and dmu
periods = unique(airportData.period)
dmus    = unique(airportData.dmu)

nGI, nBI, nGO, nBO, nPer, nDMUs,  = length(airportGoodInputs), length(airportBadInputs), length(airportGoodOutputs), length(airportBadOutputs), length(periods),length(dmus)

# Setting empty containers for our data
# Each of them is a 3D matrix where the first dimension is the decision units, the second one is the individual input or output item and the third dimension is the period to which the data refer
gI = Array{Float64}(undef, (nDMUs,nGI,nPer)) # Good inputs
bI = Array{Float64}(undef, (nDMUs,nBI,nPer)) # Bad inputs (optional)
gO = Array{Float64}(undef, (nDMUs,nGO,nPer)) # Good outputs, aka "desiderable" outputs
bO = Array{Float64}(undef, (nDMUs,nBO,nPer)) # Bad outputs, aka "undesiderable" outputs

for (p,period) in enumerate(periods)
    periodData = airportData[airportData.period .== period,:]
    gI[:,:,p] = Matrix{Float64}(periodData[:,airportGoodInputs])
    if nBI > 0
         bI[:,:,p] = Matrix{Float64}(periodData[:,airportBadInputs])
    end
    gO[:,:,p] = Matrix{Float64}(periodData[:,airportGoodOutputs])
    bO[:,:,p] = Matrix{Float64}(periodData[:,airportBadOutputs])
end

(λ_crs, λ_convex_crs, λ_nonconvex_crs, nonConvTest_crs, nonConvTest_value_crs) = efficiencyScores(
  gI,gO,bO,bI,retToScale="variable", dirGI=0,dirBI=0,dirGO=1,dirBO=0, prodStructure="multiplicative")
@test nonConvTest_value_crs[3,2]  ≈ 0.8842007631310966

################################################################################
# Basic testing of dmuEfficiency

println("Testing of the basic dmuEfficiency function..")

I  = [10 2; 8 4; 12 1.5; 24 3]
O =  [100;80;120;120]
I₀ = [24 3]
O₀ = [120]

results = dmuEfficiency(I₀,O₀,I,O)

I = [
  3	5
  2.5	4.5
  4	6
  6	7
  2.3	3.5
  4	6.5
  7	10
  4.4	6.4
  3	5
  5	7
  5	7
  2	4
  5	7
  4	4
  2	3
  3	6
  7	11
  4	6
  3	4
  5	6 ]
O = [
  40	55	30
  45	50	40
  55	45	30
  48	20	60
  28	50	25
  48	20	65
  80	65	57
  25	48	30
  45	64	42
  70	65	48
  45	65	40
  45	40	44
  65	25	35
  38	18	64
  20	50	15
  38	20	60
  68	64	54
  25	38	20
  45	67	32
  57	60	40]

nDMU = size(I,1)
efficiencies = [dmuEfficiency(I[d,:],O[d,:],I,O)[:obj] for  d in 1:nDMU]
efficiencies = hcat(1:nDMU,efficiencies)
efficiencies = efficiencies[sortperm(efficiencies[:, 2],rev=true), :]

# Test treating last output as bad
O2 = convert(Array{Float64,2},copy(O))
O2[:,3] = 1 ./ O2[:,3] # third is a bad output
efficiencies = [dmuEfficiency(I[d,:],O2[d,:],I,O2)[:obj] for  d in 1:nDMU]
efficiencies = hcat(1:nDMU,efficiencies)
efficiencies = efficiencies[sortperm(efficiencies[:, 2],rev=true), :]

#Table 13.1 Cooper 2006
I = ones(Float64,9)
O = [
  1 1
  2 1
  6 2
  8 4
  9 7
  5 2
  4 3
  6 4
  4 6
]
nDMU = size(I,1)
O2 = convert(Array{Float64,2},copy(O))
O2[:,2] = 1 ./ O2[:,2] # third is a bad output

efficiencies = [dmuEfficiency(I[d,:],O2[d,:],I,O2)[:obj] for  d in 1:nDMU]
efficiencies = hcat(1:nDMU,efficiencies)
efficiencies = efficiencies[sortperm(efficiencies[:, 2],rev=true), :]

I = [
4	140
5	90
6	36
10	300
11	66
8	36
9	12
5	210
5.5	33
8	288
10	80
8	8
]

O =[
2	28
1	22.5
6	12
8	60
7	16.5
6	12
7	6
3	30
4.4	5.5
4	72
2	20
1	4
]

nDMU = size(I,1)
efficiencies = [dmuEfficiency(I[d,:],O[d,:],I,O)[:obj] for  d in 1:nDMU]
efficiencies = hcat(1:nDMU,efficiencies)
efficiencies = efficiencies[sortperm(efficiencies[:, 2],rev=true), :]

I = [
 1
 1
 1
 1
 1
 1
 1
 1]

O = [
 1 7
 2 7
 4 6
 6 4
 7 2
 7 1
 2 2
 5.3 5.3
]
nDMU = size(I,1)
efficiencies = [dmuEfficiency(I[d,:],O[d,:],I,O)[:obj] for  d in 1:nDMU]
#efficiencies = hcat(1:nDMU,efficiencies)
#scatter(O[:,1],O[:,2])
@test efficiencies == [1.0,1.0,1.0,1.0,1.0,1.0,0.37735849056603776,1.0]

# From Example 2.2 of Cooper et oth. 2006 "Data envelopment Analysis. A Comprehensive Text with...."

X = [4 7 8 4 2 10; 3 3 1 2 4 1]
Y = [1 1 1 1 1 1]

nDMU = size(X,2)
out = [dmuEfficiency(X[:,d],Y[:,d],X',Y') for  d in 1:nDMU]
@test out[1].obj ≈ 0.8571428571428571
@test out[1].wI ≈ [0.14285714285714285, 0.14285714285714285]
@test out[1].wO ≈ [0.8571428571428571]
@test collect(keys(out[1].refSet)) ≈ [5,4] || collect(keys(out[1].refSet)) ≈ [4,5]
@test collect(values(out[1].refSet)) ≈ [0.2857142857142857,0.7142857142857143] ||  collect(values(out[1].refSet)) ≈ [0.7142857142857143,0.2857142857142857]
@test [i[:eff] for i in out] == [false,false,true,true,true,false]

outDual = [dmuEfficiencyDual(X[:,d],Y[:,d],X',Y') for  d in 1:nDMU]
@test [i[:eff] for i in outDual] == [false,false,true,true,true,false]

# From Example 2.2 of Cooper et oth. 2006 "Data envelopment Analysis. A Comprehensive Text with...."

X = [4 7 8 4 2 10; 3 3 1 2 4 1]
Y = [1 1 1 1 1 1]

nDMU = size(X,2)
out = [dmuEfficiency(X[:,d],Y[:,d],X',Y') for  d in 1:nDMU]
@test out[1].obj ≈ 0.8571428571428571
@test out[1].wI ≈ [0.14285714285714285, 0.14285714285714285]
@test out[1].wO ≈ [0.8571428571428571]
@test collect(keys(out[1].refSet)) ≈ [5,4] || collect(keys(out[1].refSet)) ≈ [4,5]
@test collect(values(out[1].refSet)) ≈ [0.2857142857142857,0.7142857142857143] ||  collect(values(out[1].refSet)) ≈ [0.7142857142857143,0.2857142857142857]
@test [i[:eff] for i in out] == [false,false,true,true,true,false]

outDual = [dmuEfficiencyDual(X[:,d],Y[:,d],X',Y') for  d in 1:nDMU]
@test [i[:eff] for i in outDual] == [false,false,true,true,true,false]


text:
############################################################################################
# Import External Packages
using Test
using Random: Random, AbstractRNG, seed!
using UnPack: UnPack, @unpack, @pack!
using Distributions

############################################################################################
# Import Baytes Packages
using BaytesCore, ModelWrappers, BaytesMCMC

############################################################################################
# Include Files
include("TestHelper.jl")

############################################################################################
# Run Tests
@testset "All tests" begin
    include("test-construction.jl")
    #include("test-nuts.jl")
end


text:
function spdiags(B,d,m,n)
  d = d[:]
  p = length(d)

  len = zeros(p+1,1)
  for k = 1:p
      len[k+1] = int(len[k]+length(Base.max(1,1-d[k]): Base.min(m,n-d[k])))
  end
  a = zeros(int(len[p+1]),3)
  for k = 1:p
      # Append new d[k]-th diagonal to compact form
      i = Base.max(1,1-d[k]):Base.min(m,n-d[k])
      a[(int(len[k])+1):int(len[k+1]),:] = [i i+d[k] B[i+(m>=n)*d[k],k]]
  end

  A = sparse(int(a[:,1]),int(a[:,2]),a[:,3],m,n)

  return A
end

y = readcsv("snp_500.txt")
n = length(y)

e = ones(n, 1)
D = spdiags([e -2*e e], 0:2, n-2, n)

lambda = 50
start = time()
D*x
println(time() - start)

x = Variable(n)
p = minimize(lambda*sum(D*x))
println(time() - start)
solve!(p)
p.optval
println(time() - start)


text:
abstract type AbstractJones{T} <: AbstractMatrix{T} end

struct Gain{S, T<:Number} <: AbstractJones{T}
    g1::T
    g2::T
end

struct DTerm{S, T<:Number} <: AbstractJones{T}
    d1::T
    d2::T
end


text:
@testset "Heat balance for single line" begin
    I_lim   = 732.147 # A
    V_base  = 138.0e3 # V
    emm     = 0.7  # emissivity, [0.23, 0.91]
    T_s     = 70.0 # conductor surface temperature, C
    T_a     = 35.0 # ambient temperature, C
    phi     = 90.0 # wind/line angle in degrees: wind perpendicular to line.
    H_e     = 61.0 # height above sea level in m. Avg PJM elevation.
    V_w     = 0.61 # wind speed in m/s
    alpha   = 0.9  # solar absorptivity, [0.23, 0.91]
    lat     = 40.0 # latitude in deg
    N       = 161  # day of year: June 10
    Z_l     = 90.0 # line azimuth: West-to-East
    hours_from_noon = 0.0 # time: noon

    # R and bundle are not necessary here
    a = acsr_interpolation(I_lim, V_base)
    D, Al_m, St_m, R, bundle, label = a.D, a.Al_m, a.St_m, a.R, a.bundle, a.label

    A_prime = D # conductor area, m^2 per linear m
    eta_r   = eq_eta_r(D, emm)
    T_film  = eq6_T_film(T_s, T_a)
    k_f     = eq15a_k_f(T_film)
    K_angle = eq4a_K_angle(phi)
    p_f     = eq14a_p_f(H_e, T_film)
    mu_f    = eq13a_mu_f(T_film)
    N_Re    = eq2c_N_Re(D, p_f, V_w, mu_f)
    eta_c   = eq_eta_c(k_f, K_angle, N_Re)

    omega   = eq_omega(hours_from_noon)
    delta   = eq16b_delta(N)
    H_c     = eq16a_H_c(lat, delta, omega)
    chi     = eq17b_chi(omega, lat, delta)
    C       = eqtable2_C(omega, chi)
    Z_c     = eq17a_Z_c(C, chi)
    theta   = eq9_theta(H_c, Z_c, Z_l)
    Q_s     = eq18_Q_s(H_c)
    K_solar = eq20_K_solar(H_e)
    Q_se    = eq19_Q_se(K_solar, Q_s)
    eta_s   = eq8_q_s(alpha, Q_se, theta, A_prime)

    @test D ≈ 0.0232156 atol=atol
    @test Al_m ≈ 0.779797 atol=atol
    @test St_m ≈ 0.285727 atol=atol
    @test R ≈ 6.167979e-5 atol=atol
    @test bundle == 1
    @test label == "Parakeet"

    @test eta_r ≈ 0.289266376 atol=atol
    @test T_film == 52.5
    @test k_f ≈ 0.02815328 atol=atol
    @test K_angle ≈ 1.2311481927818961 atol=atol
    @test p_f ≈ 1.0763378424625318 atol=atol
    @test mu_f ≈ 1.9642517300084803e-5 atol=atol
    @test N_Re ≈ 775.999091387987 atol=atol
    @test eta_c ≈ 1.5240287413704434 atol=atol

    @test omega == 0.0
    @test delta ≈ 23.021449792571953 atol=atol
    @test H_c ≈ 73.02144979257197 atol=atol
    @test chi == 0.0
    @test C == 180.0
    @test Z_c == 180.0
    @test theta == 90.0
    @test Q_s ≈ 1025.3693031729497 atol=atol
    @test K_solar ≈ 1.00696157132 atol=atol
    @test Q_se ≈ 1032.5074847063267 atol=atol
    @test eta_s ≈ 21.57325268575338 atol=atol
end

@testset "ACSR interpolation" begin
    acsr = acsr_interpolation(400.0, 138e3)
    @test acsr isa ACSRSpecsMetric
    @test acsr.label == "Penguin"
    @test acsr.R ≈ 1.9520997e-4 atol=atol

    acsr = acsr_interpolation(400.0, 138e3; metric=false)
    @test acsr isa ACSRSpecsEnglish
    @test acsr.label == "Penguin"
    @test acsr.R ≈ 0.119 atol=atol

    acsr = acsr_interpolation(1200.0, 138e3)
    @test acsr.label == "Lark"
    @test acsr.bundle == 2

    acsr = acsr_interpolation(10.0, 138e3)
    @test acsr.label == "Turkey"

    acsr = acsr_interpolation(3e3, 338e3)
    @test acsr.label == "Chukar"
    @test acsr.bundle == 2

    acsr = acsr_interpolation(200.0, 20e3)
    @test acsr.label == "Sparrow"
    @test acsr.bundle == 1

    @test_logs (:warn, "No ACSR match for line with I_lim = 6.0e6 A and V_base = 20000.0 kV") acsr_interpolation(6e6, 20e3)
    acsr = acsr_interpolation(6e6, 20e3)
    @test acsr.label == ""
    @test acsr.bundle == 1000
    @test isnan(acsr.D)
end

@testset "Line lengths" begin
    I_lim = 400.0
    V_base = 138e3
    acsr = acsr_interpolation(I_lim, V_base)

    S_base = 100e6
    R_pu = 0.001
    l = estimate_length(S_base, V_base, R_pu, acsr.R, acsr.bundle)

    l_test = 975.5649076

    @test l ≈ l_test atol=atol

    l = estimate_length(S_base, V_base, R_pu * 10, acsr.R, acsr.bundle)
    @test l ≈ l_test * 10 atol=atol

    l = estimate_length(S_base, V_base, R_pu, acsr.R, acsr.bundle * 2)
    @test l ≈ l_test * 2 atol=atol

    l = estimate_length(S_base, V_base / 10, R_pu, acsr.R, acsr.bundle * 2)
    @test l ≈ 19.511298 atol=atol
end


text:
# save current configuration to EEPROM.
# power supply will power up in this state.
# EEP
export psp_savepreset

const SAVE_PRESET = [0x45, 0x45, 0x50, 0x0d]

"""
    psp_savepreset(io)

Save current configuration to EEPROM.
"""
function psp_savepreset(io_psp::IO)
  write(io_psp, SAVE_PRESET)
  return nothing
end


text:
function metronome(bpm::Real=72, bpb::Int=4)
    s = 60.0 / bpm
    counter = 0
    while true
        counter += 1
        if counter % bpb != 0
            println("tick")
        else
            println("TICK")
        end
        sleep(s)
    end
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
include("bimodal_cauchy_distribution.jl")

ADE_DefaultOptions = merge(DE_DefaultOptions, {
  # Distributions we will use to generate new F and CR values.
  "fdistr" => bimodal_cauchy(0.65, 0.1, 1.0, 0.1),
  "crdistr" => bimodal_cauchy(0.1, 0.1, 0.95, 0.1),
})

# An Adaptive DE typically change parameters of the search dynamically. This is
# typically done in the tell! function when we know if the trial vector
# was better than the target vector.
type AdaptConstantsDiffEvoOpt <: DifferentialEvolutionOpt
  name::ASCIIString

  # A population is a matrix of floats.
  population::Array{Float64, 2}

  search_space::SearchSpace

  # Options
  options

  # Set of functions that together define a specific DE strategy.
  sample::Function
  mutate::Function
  crossover::Function
  bound::Function

  # Specific data and functions for adaptation
  fs::Vector{Float64}   # One f value per individual in population
  crs::Vector{Float64}  # One cr value per individual in population

  function AdaptConstantsDiffEvoOpt(name, pop, ss, options, sample, mutate, crossover, bound)
    popsize = size(pop, 1)
    fs = [sample_bimodal_cauchy(options["fdistr"]; truncateBelow0 = false) for i in 1:popsize]
    crs = [sample_bimodal_cauchy(options["crdistr"]) for i in 1:popsize]
    new(name, pop, ss, merge(DE_DefaultOptions, options),
      sample, mutate, crossover, bound, fs, crs)
  end
end

# To get the constants for an adaptive DE we access the vectors of constants.
fconst(ade::AdaptConstantsDiffEvoOpt, i) = ade.fs[i]
crconst(ade::AdaptConstantsDiffEvoOpt, i) = ade.crs[i]

# To sample we use the distribution given as options
sample_f(ade::AdaptConstantsDiffEvoOpt) = sample_bimodal_cauchy(ade.options["fdistr"]; truncateBelow0 = false)
sample_cr(ade::AdaptConstantsDiffEvoOpt) = sample_bimodal_cauchy(ade.options["crdistr"]; truncateBelow0 = false)

# Tell the optimizer about the ranking of candidates. Returns the number of
# better candidates that were inserted into the population.
function tell!(de::AdaptConstantsDiffEvoOpt,
  # archive::Archive, # Skip for now
  rankedCandidates)
  num_candidates = length(rankedCandidates)
  num_better = 0
  for i in 1:div(num_candidates, 2)
    candidate, index = rankedCandidates[i]
    if candidate != de.population[index, :]
      num_better += 1
      old = de.population[index,:]
      de.population[index,:] = candidate
      # Since the trial vector was better we keep the f and cr values for this target.
    else
      # The trial vector for this target was not better so we change the f and cr constants.
      de.fs[index] = sample_f(de)
      de.crs[index] = sample_cr(de)
    end
  end
  num_better
end

function adaptive_de_rand_1_bin(parameters = Dict())
  params = Parameters(parameters, ADE_DefaultOptions)
  ss = get(params, :SearchSpace, BlackBoxOptim.symmetric_search_space(1))
  population = get(params, :Population, BlackBoxOptim.rand_individuals_lhs(ss, 50))
  AdaptConstantsDiffEvoOpt("AdaptiveDE/rand/1/bin", population, ss, params,
    random_sampler,
    de_mutation_rand_1,
    de_crossover_binomial,
    rand_bound_from_target!)
end

function adaptive_de_rand_1_bin_radiuslimited(parameters = Dict())
  params = Parameters(parameters, ADE_DefaultOptions)
  ss = get(params, :SearchSpace, BlackBoxOptim.symmetric_search_space(1))
  population = get(params, :Population, BlackBoxOptim.rand_individuals_lhs(ss, 50))
  AdaptConstantsDiffEvoOpt("AdaptiveDE/rand/1/bin/radiuslimited", population, ss, params,
    radius_limited_sampler,
    de_mutation_rand_1,
    de_crossover_binomial,
    rand_bound_from_target!)
end


text:
@testset "SMeasure" begin
    x, y = rand(100), rand(100)
    X, Y = Dataset(rand(100, 3)), Dataset(rand(100, 2))
    Z, W = Dataset(rand(110, 2)), Dataset(rand(90, 4))
    dx, τx = 2, 1
    dy, τy = 2, 1

    @test s_measure(x, y) isa Float64
    @test s_measure(x, Y, dx = dx, τx = τx) isa Float64
    @test s_measure(X, y, dy = dy, τy = τy) isa Float64
    @test s_measure(X, Y) isa Float64
    # test that multivariate datasets are being length-matched
    @test s_measure(X, Z) isa Float64
    @test s_measure(W, X) isa Float64
end


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

"""
	The number of GPU or other accelerators found by the CVortex library.

It is possible that accelerators may be listed multiple times if they can
be used by more than one installed platform.
To know how many are in use see number_of_enabled_accelerators()
"""
function number_of_accelerators()
	res = ccall(("cvtx_num_accelerators", libcvortex),
		Cint, ())
	return res
end

"""
	The number of accelerators that CVortex has been directed to use.

If no accelerators are in use this is zero, and the CPU is used for
all computation.

Find which accelerators are enabled using accelerator_enabled,
and enable and disable with accelerator_enable and
accelerator_disable.
"""
function number_of_enabled_accelerators()
	# int cvtx_num_enabled_accelerators();
	res = ccall(("cvtx_num_enabled_accelerators", libcvortex), Cint, ())
	return res
end

"""
	The name of an accelerator.

Input is an integer in the range 1:number_of_accelerators(). Returns
the name of the accelerator as a string.
"""
function accelerator_name(accelerator_id :: Int)
	@assert(accelerator_id >= 1, "Minimum accelerator id is 1.")
	@assert(accelerator_id <= number_of_accelerators(),
		"accelerator_id is higher than the number of accelerators found.")
	# char* cvtx_accelerator_name(int accelerator_id);
	res = ccall(("cvtx_accelerator_name", libcvortex),
		Cstring, (Cint,), accelerator_id-1)
	return unsafe_string(res)
end

"""
	States whether CVortex is in use.
"""
function accelerator_enabled(accelerator_id :: Int)
	@assert(accelerator_id >= 1, "Minimum accelerator id is 1.")
	@assert(accelerator_id <= number_of_accelerators(),
		"accelerator_id is higher than the number of accelerators found.")
	# int cvtx_accelerator_enabled(int accelerator_id);
	res = ccall(("cvtx_accelerator_enabled", libcvortex),
		Cint, (Cint,), accelerator_id-1)
	return res
end

"""
	Allows CVortex to use an accelerator.
"""
function accelerator_enable(accelerator_id :: Int)
	@assert(accelerator_id >= 1, "Minimum accelerator id is 1.")
	@assert(accelerator_id <= number_of_accelerators(),
		"accelerator_id is higher than the number of accelerators found.")
	# void cvtx_accelerator_enable(int accelerator_id);
	ccall(("cvtx_accelerator_enable", libcvortex),
		Cvoid, (Cint,), accelerator_id-1)
	return
end

"""
	Stops CVortex using an accelerator.
"""
function accelerator_disable(accelerator_id :: Int)
	@assert(accelerator_id >= 1, "Minimum accelerator id is 1.")
	@assert(accelerator_id <= number_of_accelerators(),
		"accelerator_id is higher than the number of accelerators found.")
	# void cvtx_accelerator_disable(int accelerator_id);
	ccall(("cvtx_accelerator_disable", libcvortex),
		Cvoid, (Cint,), accelerator_id-1)
	return
end


text:
#=
  https://dtai.cs.kuleuven.be/problog/tutorial/basic/03_dice.html
  """
  The following example illustrates the use of a logic program with recursion and lists.
  We start by rolling the first die. Every roll determines the next die to roll, but we stop
  if we have used that die before. We query for the possible sequences of rolled dice.
  We use three-sided dice instead of the regular six-sided ones simply to restrict the number
  of possible outcomes (and thus inference time).
  """

  Cf ~/blog/rolling_dice5.blog
     ~/webppl/rolling_dice5.wppl
=#

using Turing, StatsPlots, DataFrames
include("jl_utils.jl")

@model function rolling_dice5(n=3)

    function roll(a)
        # t ~ DiscreteUniform(1,n)
        t = rand(DiscreteUniform(1,n))
        if t in a
            return a
        else
            return roll(vcat(a,t))
        end
    end

    a = roll([])
    len ~ DiscreteUniform(1,n)
    len ~ Dirac(length(a))

    # In the example, the first roll is 1.
    true ~ Dirac(a[1] == 1)

end

function run_model(len=0)
    model = rolling_dice5(len)

    num_chns = 4

    # chns = sample(model, Prior(), 10_000)

    # chns = sample(model, MH(), 10_000)
    chns = sample(model, MH(), 100_000)

    # chns = sample(model, PG(15), MCMCThreads(), 1_000, num_chns)

    # chns = sample(model, SMC(1000), MCMCThreads(), 10_000, num_chns)
    # chns = sample(model, SMC(1000), 10_000)
    # chns = sample(model, IS(), 10_000)
    #
    display(chns)
    show_var_dist_pct(chns,:len,1000)

end

genq = undef
for val in 1:6
    global genq
    println("\nval:$val")
    @time run_model(val)
end

println("\nval=10")
# println("mean:", mean(run_model(10)))
run_model(10)
println("\nval=100")
# println("mean:", mean(run_model(100)))
run_model(100)

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
# # Traversal
#
# This code demonstrates how to use and extend advanced optics.
# The examples are taken from the README.md of the specter clojure library.
# Many of the features here are experimental, please consult the docstrings of the
# involved optics.

using Test
using Accessors

import Accessors: modify, OpticStyle
using Accessors: ModifyBased, SetBased, setindex

# ### Increment all even numbers
# We have the following data and the goal is to increment all nested even numbers.
data = (a = [(aa=1, bb=2), (cc=3,)], b = [(dd=4,)])
# To acomplish this, we define a new optic `Vals`.
function mapvals(f, d)
    Dict(k => f(v) for (k,v) in pairs(d))
end

mapvals(f, nt::NamedTuple) = map(f, nt)

struct Vals end
OpticStyle(::Type{Vals}) = ModifyBased()
modify(f, obj, ::Vals) = mapvals(f, obj)

# Now we can increment as follows:
out = @set data |> Vals() |> Elements() |> Vals() |> If(iseven) += 1

@test out == (a = [(aa = 1, bb = 3), (cc = 3,)], b = [(dd = 5,)])

struct Filter{F}
    keep_condition::F
end
OpticStyle(::Type{<:Filter}) = ModifyBased()
(o::Filter)(x) = filter(o.keep_condition, x)
function modify(f, obj, optic::Filter)
    I = eltype(eachindex(obj))
    inds = I[]
    for i in eachindex(obj)
        x = obj[i]
        if optic.keep_condition(x)
            push!(inds, i)
        end
    end
    vals = f(obj[inds])
    setindex(obj, vals, inds)
end

# ### Append to nested vector
data = (a = 1:3,)
out = @modify(v -> vcat(v, [4,5]), data.a)

@test out == (a = [1,2,3,4,5],)

# ### Increment last odd number in a sequence

data = 1:4
out = @set data |> Filter(isodd) |> last += 1
@test out == [1,2,4,4]

### Map over a sequence

data = 1:3
out = @set data |> Elements() += 1
@test out == [2,3,4]

# ### Increment all values in a nested Dict

data = Dict(:a => Dict(:aa =>1), :b => Dict(:ba => -1, :bb => 2))
out = @set data |> Vals() |> Vals() += 1
@test out == Dict(:a => Dict(:aa => 2),:b => Dict(:bb => 3,:ba => 0))

# ### Increment all the even values for :a keys in a sequence of maps

data = [Dict(:a => 1), Dict(:a => 2), Dict(:a => 4), Dict(:a => 3)]
out = @set data |> Elements() |> _[:a] += 1
@test out == [Dict(:a => 2), Dict(:a => 3), Dict(:a => 5), Dict(:a => 4)]

# ### Retrieve every number divisible by 3 out of a sequence of sequences

function getall(obj, optic)
    out = Any[]
    modify(obj, optic) do val
        push!(out, val)
    end
    out
end

data = [[1,2,3,4],[], [5,3,2,18],[2,4,6], [12]]
optic = @optic _ |> Elements() |> Elements() |> If(x -> mod(x, 3) == 0)
out = getall(data, optic)
@test out == [3, 3, 18, 6, 12]
@test_broken eltype(out) == Int

# ### Increment the last odd number in a sequence

data = [2, 1, 3, 6, 9, 4, 8]
out = @set data |> Filter(isodd) |> _[end] += 1
@test out == [2, 1, 3, 6, 10, 4, 8]
@test_broken eltype(out) == Int

# ### Remove nils from a nested sequence

data = (a = [1,2,missing, 3, missing],)
optic = @optic _.a |> Filter(!ismissing)
out = optic(data)
@test out == [1,2,3]


text:
using SplitApplyCombine

"""
    normalize_year(fun, fams)

Compute function `fun` for each family in `fams` and divide the result
by the average for all families with the same earliest filing date.
"""
function normalize_year(fun, fams::Vector{Family})
    year = Dates.year.(earliest_filing.(fams))
    vals = fun.(fams)
    dict = group(year, vals)
    refs = mean.(dict)
    map(zip(year, vals)) do (y, v)
        r = refs[y]
        r == 0 ? 0 : v / r
    end
end

"""
    normalize_year(fun, g, fams)

Compute function `fun` for all vertices in `g` and divide the result
by the average for all families with the same earliest filing date.
"""
function normalize_year(fun, g::AbstractGraph, fams::Vector{Family})
    year = Dates.year.(earliest_filing.(fams))
    vals = fun(g)
    dict = group(year, vals)
    refs = mean.(dict)
    map(zip(year, vals)) do (y, v)
        r = refs[y]
        r == 0.0 ? 0.0 : v / r
    end
end

function component(g::AbstractGraph, which=1)
    cc = weakly_connected_components(g)
    idx = sortperm(length.(cc), rev=true)
    vids = reduce(vcat, cc[idx[which]])
    induced_subgraph(g, vids)
end

function top_applicants(families, k)
    apps = reduce(vcat, applicants.(families))
    c = StatsBase.countmap(apps)
    ks = collect(keys(c)); vs = collect(values(c))
    idx = sortperm(vs, rev=true)
    NamedTuple(Symbol.(first(ks[idx], k)) .=> first(vs[idx], k))
end

function maxcentralization_degree(g)
	N = nv(g)
	if is_directed(g)
		return (N-1)*(N-1)
	else
		return (N-1)*(N-2)
	end
end

function centralization(g, measure)
	c = measure(g)
	@assert measure in [degree, indegree, outdegree]
	val, ind = findmax(c)
	return sum(val .- c) / maxcentralization_degree(g)
end

using MainPaths

function component(mp::MainPaths.MainPathResult, which=1)
    g, cs = Patents.component(mp.mainpath, which)
    vs = mp.vertices[cs]
    s = Set(mp.vertices[mp.start])
    start = findall(v -> v in s, vs)
    MainPaths.MainPathResult(g, vs, start)
end

function get_stats(mp, segment, families, g)
    idxs = mp.vertices[segment.intermediates]
    fams = families[idxs]
    ef = Dates.year.(earliest_filing.(fams))
    (;
        size = length(fams),
        span = maximum(ef) - minimum(ef),
        density = (maximum(ef) - minimum(ef)) / length(fams),
        weight = MainPaths.meanweight(mp, segment, SPCEdge(:log)(g)),
        familysize = mean(f.size for f in fams),
        citedinternal = mean(outdegree(g, idxs)),
        citedexternal = mean(citedby_count.(fams)),
        citinginternal = mean(indegree(g, idxs)),
        citingexternal = mean(cites_count.(fams)),
        classdiversity = Patents.diversity(fams, subclass),
        jurdiversity = Patents.jurdiversity(fams),
        applicant_homogeneity = Patents.share_same_applicant(fams),
        topapplicant = Patents.top_applicants(fams, 2)
    )
end

function findstart(g, fams, k; period=(2000, 2015), cpc=["Y", "B", "C"], levelfun=section)
    idx = map(fams) do f
        y = Dates.year(earliest_filing(f))
        c = levelfun.(classification(f))
        (period[1] .<= y .<= period[2]) && any(s in cpc for s in c)
    end |> findall
    cit = Patents.normalize_year(outdegree, g, fams)
    ck = sort(cit[idx], rev=true)[k]
    intersect(idx, findall(cit .>= ck))
end