text:
# Tests for Radial DEA Models
@testset "AdditiveDEAModel" begin

    ## Test Additive DEA Models with FLS Book data
    # Test against results in R
    X = [5 13; 16 12; 16 26; 17 15; 18 14; 23 6; 25 10; 27 22; 37 14; 42 25; 5 17]
    Y = [12; 14; 25; 26; 8; 9; 27; 30; 31; 26; 12]

    # Additive CRS
    deaaddcrs1 = deaadd(X, Y, :Custom, rhoX = ones(size(X)), rhoY = ones(size(Y)), rts = :CRS)

    @test typeof(deaaddcrs1) == AdditiveDEAModel

    @test nobs(deaaddcrs1) == 11
    @test ninputs(deaaddcrs1) == 2
    @test noutputs(deaaddcrs1) == 1
    @test efficiency(deaaddcrs1) ≈ [0.0000000000;
                               10.76923077;
                               10.83783784;
                               0.0000000000;
                               22.15384615;
                               17.92307692;
                               0.0000000000;
                               12.07692308;
                               11.61379310;
                               35.000000000;
                               4.000000000]

    @test convert(Matrix, peers(deaaddcrs1)) ≈
    [ 1.0000000000  0  0 0.0000000000  0  0 0.0000000000  0  0   0   0
      0.0000000000  0  0 0.5384615385  0  0 0.0000000000  0  0   0   0
      0.1216216216  0  0 0.9054054054  0  0 0.0000000000  0  0   0   0
      0.0000000000  0  0 1.0000000000  0  0 0.0000000000  0  0   0   0
      0.0000000000  0  0 0.3076923077  0  0 0.0000000000  0  0   0   0
      0.0000000000  0  0 0.3461538462  0  0 0.0000000000  0  0   0   0
      0.0000000000  0  0 0.0000000000  0  0 1.0000000000  0  0   0   0
      0.0000000000  0  0 1.1538461538  0  0 0.0000000000  0  0   0   0
      0.0000000000  0  0 0.4689655172  0  0 0.6965517241  0  0   0   0
      0.0000000000  0  0 1.0000000000  0  0 0.0000000000  0  0   0   0
      1.0000000000  0  0 0.0000000000  0  0 0.0000000000  0  0   0   0]

    # Additive VRS model
    deaaddvrs1 = deaadd(X, Y, :Custom, rhoX = ones(size(X)), rhoY = ones(size(Y)), rts = :VRS)

    @test nobs(deaaddvrs1) == 11
    @test ninputs(deaaddvrs1) == 2
    @test noutputs(deaaddvrs1) == 1
    @test efficiency(deaaddvrs1) ≈ [0.0000000000;
                               7.333333333;
                               0.0000000000;
                               0.0000000000;
                               18.000000000;
                               0.0000000000;
                               0.0000000000;
                               0.0000000000;
                               0.0000000000;
                               35.000000000;
                               4.000000000]

    @test convert(Matrix, peers(deaaddvrs1)) ≈
    [1.000000000  0  0  0  0  0 0.00000000000  0  0   0   0;
     0.6666666667 0  0  0  0  0 0.3333333333   0  0   0   0;
     0.000000000  0  1  0  0  0 0.00000000000  0  0   0   0;
     0.000000000  0  0  1  0  0 0.00000000000  0  0   0   0;
     1.000000000  0  0  0  0  0 0.00000000000  0  0   0   0;
     0.000000000  0  0  0  0  1 0.00000000000  0  0   0   0;
     0.000000000  0  0  0  0  0 1.00000000000  0  0   0   0;
     0.000000000  0  0  0  0  0 0.00000000000  1  0   0   0;
     0.000000000  0  0  0  0  0 0.00000000000  0  1   0   0;
     0.000000000  0  0  1  0  0 0.00000000000  0  0   0   0;
     1.000000000  0  0  0  0  0 0.00000000000  0  0   0   0]

    # Check weights equal to one equals default model
    deaadddefaultcrs = deaadd(X, Y, rts = :CRS)
    @test efficiency(deaaddcrs1) ≈ efficiency(deaadddefaultcrs)

    deaadddefaultvrs = deaadd(X, Y)
    @test efficiency(deaaddvrs1) ≈ efficiency(deaadddefaultvrs)

    # Test model :Ones equals model with all weights equal to 1
    deaaddonescrs = deaadd(X, Y, :Ones, rts = :CRS)
    @test deaaddonescrs.weights == :Ones
    @test efficiency(deaaddonescrs) ≈ efficiency(deaadddefaultcrs)

    @test efficiency(deaadd(targets(deaaddonescrs, :X), targets(deaaddonescrs, :Y), :Ones)) ≈ zeros(11) atol=1e-13

    deaaddonesvrs = deaadd(X, Y, :Ones)
    @test deaaddonesvrs.weights == :Ones
    @test efficiency(deaaddonesvrs) ≈ efficiency(deaadddefaultvrs)

    @test efficiency(deaadd(targets(deaaddonesvrs, :X), targets(deaaddonesvrs, :Y), :Ones, rts = :VRS)) ≈ zeros(11) atol=1e-12

    # MIP CRS
    deaaddmipcrs = deaadd(X, Y, :MIP, rts = :CRS)
    @test deaaddmipcrs.weights == :MIP
    @test efficiency(deaaddmipcrs) ≈ [0.0000000000;
                               0.7577160494;
                               0.4168399168;
                               0.0000000000;
                               2.2219512195;
                               1.1478260870;
                               0.0000000000;
                               0.4867909868;
                               0.4041184041;
                               1.0726153846;
                               0.2352941176]
    @test deaaddmipcrs.slackX ≈ [0 0;
                                 3.037037037 6.814814815;
                                 0 10.837837838;
                                 0 0;
                                 0 0;
                                 8.000000000  0;
                                 0 0;
                                 7.384615385 4.692307692;
                                 8.296296296 2.518518519;
                                 0.000000000  8.200000000;
                                 0 4.000000000]

    @test deaaddmipcrs.slackY ≈ [0;
                                  0;
                                  0;
                                  0;
                                  17.77560976 ;
                                  7.20000000;
                                  0;
                                  0;
                                  0;
                                  19.36000000;
                                  0]

    # MIP VRS
    deaaddmipvrs = deaadd(X, Y, :MIP, rts = :VRS)
    @test deaaddmipvrs.weights == :MIP
    @test efficiency(deaaddmipvrs) ≈ [0.0000000000;
                               0.5075187970;
                               0.0000000000;
                               0.0000000000;
                               2.2039473684;
                               0.0000000000;
                               0.0000000000;
                               0.0000000000;
                               0.0000000000;
                               1.0432234432;
                               0.2352941176]
    @test deaaddmipvrs.slackX ≈ [0 0;
                                 0 0;
                                 0 0;
                                 0 0;
                                 0 0;
                                 0 0;
                                 0 0;
                                 0 0;
                                 0 0;
                                 17 15;
                                 0 4]
    @test deaaddmipvrs.slackY ≈ [0;
                                 7.105263158;
                                 0;
                                 0;
                                 17.631578947;
                                 0;
                                 0;
                                 0;
                                 0;
                                 1.000000000;
                                 0]

    # Normalized CRS
    deaaddnormalizedcrs = deaadd(X, Y, :Normalized, rts = :CRS)
    @test deaaddnormalizedcrs.weights == :Normalized
    @test efficiency(deaaddnormalizedcrs) ≈ [0.0000000000;
                                 1.3569256615;
                                 1.7407259078;
                                 0.0000000000;
                                 2.6877810368;
                                 1.6953153107;
                                 0.0000000000;
                                 1.6540884810;
                                 1.1212042237;
                                 4.0172855145;
                                 0.6424624298]

    @test deaaddnormalizedcrs.slackX ≈ [0.000000000  0.000000000;
                                 3.037037037  6.814814815;
                                 0.000000000 10.837837838;
                                 0.000000000  0.000000000;
                                 10.592592593 11.037037037;
                                 14.666666667  2.666666667;
                                 0.000000000  0.000000000;
                                 0.000000000 10.298429319;
                                 8.296296296  2.518518519;
                                 17.925925926 15.370370370;
                                 0.000000000  4.000000000]

    @test deaaddnormalizedcrs.slackY ≈ [0;
                                  0;
                                  0;
                                  0;
                                  0 ;
                                  0;
                                  0;
                                  0;
                                  0;
                                  0;
                                  0]

    # Normalized VRS
    deaaddnormalizedvrs = deaadd(X, Y, :Normalized, rts = :VRS)
    @test deaaddnormalizedvrs.weights == :Normalized
    @test efficiency(deaaddnormalizedvrs) ≈ [0.0000000000;
                                 0.8049248943;
                                 0.0000000000;
                                 0.0000000000;
                                 2.0149704908;
                                 0.0000000000;
                                 0.0000000000;
                                 0.0000000000;
                                 0.0000000000;
                                 3.9898943952;
                                 0.6424624298]

    @test deaaddnormalizedvrs.slackX ≈ [0  0.00;
                                 0  0.65;
                                 0  0.00;
                                 0  0.00;
                                 0  2.95;
                                 0  0.00;
                                 0  0.00;
                                 0  0.00;
                                 0  0.00;
                                 17 15.00;
                                 0  4.00]

    @test deaaddnormalizedvrs.slackY ≈ [0.00;
                                6.25;
                                0.00;
                                0.00;
                                13.75;
                                0.00;
                                0.00;
                                0.00;
                                0.00;
                                1.00;
                                0.00]

    # RAM CRS
    deaaddramcrs = deaadd(X, Y, :RAM, rts = :CRS)
    @test deaaddramcrs.weights == :RAM
    @test efficiency(deaaddramcrs) ≈ [0.0000000000;
                                0.14094094094;
                                0.18063063063;
                                0.0000000000;
                                0.27937937938;
                                0.17657657658;
                                0.0000000000;
                                0.17164048866;
                                0.11671671672;
                                0.41766766767;
                                0.06666666667]

    @test deaaddramcrs.slackX ≈ [ 0.000000000  0.000000000;
                                3.037037037  6.814814815;
                                0.000000000  10.837837838;
                                0.000000000  0.000000000;
                                10.592592593 11.037037037;
                                14.666666667 2.666666667;
                                0.000000000  0.000000000;
                                0.000000000  10.298429319;
                                8.296296296  2.518518519;
                                17.925925926 15.370370370;
                                0.000000000  4.000000000]

    @test deaaddramcrs.slackY ≈ [0;
                                   0;
                                   0;
                                   0;
                                   0 ;
                                   0;
                                   0;
                                   0;
                                   0;
                                   0;
                                   0]

    # RAM VRS
    deaaddramvrs = deaadd(X, Y, :RAM, rts = :VRS)
    @test deaaddramvrs.weights == :RAM
    @test efficiency(deaaddramvrs) ≈ [0.0000000000;
                              0.10297482838;
                              0.0000000000;
                              0.0000000000;
                              0.25553012967;
                              0.0000000000;
                              0.0000000000;
                              0.0000000000;
                              0.0000000000;
                              0.41764590678 ;
                              0.06666666667]

    @test deaaddramvrs.slackX ≈ [0 0;
                                0 0;
                                0 0;
                                0 0;
                                0 0;
                                0 0;
                                0 0;
                                0 0;
                                0 0;
                                17 15;
                                0 4]
    @test deaaddramvrs.slackY ≈ [0;
                                7.105263158;
                                0;
                                0;
                                17.631578947;
                                0;
                                0;
                                0;
                                0;
                                1.000000000;
                                0]

    # BAM CRS
    deaaddbamcrs = deaadd(X, Y, :BAM, rts = :CRS)
    @test deaaddbamcrs.weights == :BAM
    @test efficiency(deaaddbamcrs) ≈ [0.0000000000;
                                0.4578892372;
                                0.3051750381;
                                0.0000000000;
                                0.6732181854;
                                0.3239568683;
                                0.0000000000;
                                0.5255235602;
                                0.1913580247;
                                0.6902867780;
                                0.1212121212]

    @test deaaddbamcrs.slackX ≈ [0.000000000  0.000000000;
                                4.110344828  6.000000000;
                                0.000000000  0.000000000;
                                0.000000000  0.000000000;
                                13.000000000  8.000000000;
                                17.493670886  0.000000000;
                                0.000000000  0.000000000;
                                0.000000000  9.225130890;
                                8.296296296  2.518518519;
                                13.296296296 13.518518519;
                                0.000000000  4.000000000]
    @test deaaddbamcrs.slackY ≈ [0;
                                0;
                                5.4931506849;
                                0;
                                0.4520547945 ;
                                0;
                                0;
                                1;
                                0;
                                5;
                                0]

    # BAM VRS
    deaaddbamvrs = deaadd(X, Y, :BAM, rts = :VRS)
    @test deaaddbamvrs.weights == :BAM
    @test efficiency(deaaddbamvrs) ≈ [0.0000000000;
                              0.1998936736;
                              0.0000000000;
                              0.0000000000;
                              0.4329710145;
                              0.0000000000;
                              0.0000000000;
                              0.0000000000;
                              0.0000000000;
                              0.5713608345;
                              0.1212121212]

    @test deaaddbamvrs.slackX ≈ [0 0;
                                6.596491228 0;
                                0 0;
                                0 0;
                                13.000000000 1;
                                0 0;
                                0 0;
                                0 0;
                                0 0;
                                5 11;
                                0 4]
    @test deaaddbamvrs.slackY ≈ [0;
                                0;
                                0;
                                0;
                                4;
                                0;
                                0;
                                0;
                                0;
                                5;
                                0]

    # Test model with Custom weights
    deaddcustomcrs = deaadd(X, Y, rhoX = 1 ./ X, rhoY = 1 ./ Y, rts = :CRS)
    @test deaddcustomcrs.weights == :Custom
    @test efficiency(deaddcustomcrs) ≈ efficiency(deaaddmipcrs)

    deaddcustomvrs = deaadd(X, Y, rhoX = 1 ./ X, rhoY = 1 ./ Y, rts = :VRS)
    @test deaddcustomvrs.weights == :Custom
    @test efficiency(deaddcustomvrs) ≈ efficiency(deaaddmipvrs)

    # Test that slacks are zero when weights are zero
    @test all(slacks(deaadd(X, Y, rhoX = zeros(size(X)), rhoY = ones(size(Y))) , :X) .== 0)
    @test all(slacks(deaadd(X, Y, rhoX = ones(size(X)) , rhoY = zeros(size(Y))), :Y) .== 0)

    ## Test if one-by-one DEA using evaluation and reference sets match initial results
    deaadddefaultcrs_ref_eff = zeros(size(X, 1))
    deaadddefaultvrs_ref_eff = zeros(size(X, 1))

    deaaddcustomcrs_ref_eff = zeros(size(X, 1))
    deaaddcustomvrs_ref_eff = zeros(size(X, 1))

    deaaddonescrs_ref_eff = zeros(size(X, 1))
    deaaddonesvrs_ref_eff = zeros(size(X, 1))

    deaaddmipcrs_ref_eff = zeros(size(X, 1))
    deaaddmipvrs_ref_eff = zeros(size(X, 1))

    Xref = X[:,:]
    Yref = Y[:,:]

    for i = 1:size(X, 1)
        Xeval = X[i:i,:]
        Xeval = Xeval[:,:]
        Yeval = Y[i:i,:]
        Yeval = Yeval[:,:]

        deaadddefaultcrs_ref_eff[i] = efficiency(deaadd(Xeval, Yeval, rts = :CRS, Xref = Xref, Yref = Yref))[1]
        deaadddefaultvrs_ref_eff[i] = efficiency(deaadd(Xeval, Yeval, rts = :VRS, Xref = Xref, Yref = Yref))[1]

        deaaddcustomcrs_ref_eff[i] = efficiency(deaadd(Xeval, Yeval, :Custom, rts = :CRS, rhoX = 1 ./ Xeval, rhoY = 1 ./ Yeval, Xref = Xref, Yref = Yref))[1]
        deaaddcustomvrs_ref_eff[i] = efficiency(deaadd(Xeval, Yeval, :Custom, rts = :VRS, rhoX = 1 ./ Xeval, rhoY = 1 ./ Yeval, Xref = Xref, Yref = Yref))[1]

        deaaddonescrs_ref_eff[i] = efficiency(deaadd(Xeval, Yeval, :Ones, rts = :CRS, Xref = Xref, Yref = Yref))[1]
        deaaddonesvrs_ref_eff[i] = efficiency(deaadd(Xeval, Yeval, :Ones, rts = :VRS, Xref = Xref, Yref = Yref))[1]

        deaaddmipcrs_ref_eff[i] = efficiency(deaadd(Xeval, Yeval, :MIP, rts = :CRS, Xref = Xref, Yref = Yref))[1]
        deaaddmipvrs_ref_eff[i] = efficiency(deaadd(Xeval, Yeval, :MIP, rts = :VRS, Xref = Xref, Yref = Yref))[1]

    end

    @test deaadddefaultcrs_ref_eff ≈ efficiency(deaadddefaultcrs)
    @test deaadddefaultvrs_ref_eff ≈ efficiency(deaadddefaultvrs)

    @test deaaddcustomcrs_ref_eff ≈ efficiency(deaaddmipcrs)
    @test deaaddcustomvrs_ref_eff ≈ efficiency(deaaddmipvrs)

    @test deaaddonescrs_ref_eff ≈ efficiency(deaaddcrs1)
    @test deaaddonesvrs_ref_eff ≈ efficiency(deaaddvrs1)

    @test deaaddmipcrs_ref_eff ≈ efficiency(deaaddmipcrs)
    @test deaaddmipvrs_ref_eff ≈ efficiency(deaaddmipvrs)

    # Print
    show(IOBuffer(), deaaddcrs1)

    # Test errors
    @test_throws DimensionMismatch deaadd([1; 2 ; 3], [4 ; 5], :Ones) #  Different number of observations
    @test_throws DimensionMismatch deaadd([1; 2], [4 ; 5], Xref = [1; 2; 3; 4], :Ones) # Different number of observations in reference sets
    @test_throws DimensionMismatch deaadd([1 1; 2 2], [4 4; 5 5], Xref = [1 1 1; 2 2 2], :Ones) # Different number of inputs
    @test_throws DimensionMismatch deaadd([1 1; 2 2], [4 4; 5 5], Yref = [4 4 4; 5 5 5], :Ones) # Different number of inputs
    @test_throws ArgumentError deaadd([1; 2; 3], [4; 5; 6], :Ones, rts = :Error) # Invalid returns to scale
    @test_throws ArgumentError deaadd([1; 2; 3], [4; 5; 6], :Error) # Invalid model
    @test_throws ArgumentError deaadd([1; 2; 3], [4; 5; 6], :Ones, orient = :Error) # Invalid orientation
    @test_throws ArgumentError deaadd([1; 2; 3], [4; 5; 6], :Custom, rhoX = [1; 1; 1], rhoY = [1; 1; 1], orient = :Error) # Invalid orientation with custom weights

    @test_throws DimensionMismatch deaadd([1; 2 ; 3], [4 ; 5]) #  Different number of observations
    @test_throws DimensionMismatch deaadd([1; 2], [4 ; 5], Xref = [1; 2; 3; 4]) # Different number of observations in reference sets
    @test_throws DimensionMismatch deaadd([1 1; 2 2], [4 4; 5 5], Xref = [1 1 1; 2 2 2]) # Different number of inputs
    @test_throws DimensionMismatch deaadd([1 1; 2 2], [4 4; 5 5], Yref = [4 4 4; 5 5 5]) # Different number of inputs
    @test_throws ArgumentError deaadd([1; 2; 3], [4; 5; 6], rts = :Error) # Invalid returns to scale
    @test_throws DimensionMismatch deaadd([1; 2; 3], [4; 5; 6], :Custom, rhoX = [1; 2; 3; 4], rhoY = [1; 1; 1]) # Different size of weights
    @test_throws DimensionMismatch deaadd([1; 2; 3], [4; 5; 6], :Custom, rhoX = [1; 1; 1], rhoY = [4; 5; 6; 7]) # Different size of weights

    @test_throws ArgumentError deaadd([1; 2; 3], [4; 5; 6], :Ones, rhoX = [1; 1; 1]) # Weights not allowed if model != :Custom

    # ------------------
    # Orientation Tests
    # ------------------

    X = [1; 2; 3; 2; 4]
    Y = [2; 3; 4; 1; 3]

    # Test default is graph
    @test efficiency(deaadd(X, Y)) == efficiency(deaadd(X, Y, orient = :Graph))

    # Graph orientaiton
    deaaddgraph = deaadd(X, Y, orient = :Graph)
    @test efficiency(deaaddgraph) ≈ [0; 0; 0; 2.0; 2.0]
    @test slacks(deaaddgraph, :X) ≈ [0; 0; 0; 0; 1.0]
    @test slacks(deaaddgraph, :Y) ≈ [0; 0; 0; 2.0; 1.0]

    # Input orientation
    deaaddinput = deaadd(X, Y, orient = :Input)
    @test efficiency(deaaddinput) ≈ [0; 0; 0; 1.0; 2.0]
    @test slacks(deaaddinput, :X) ≈ [0; 0; 0; 1.0; 2.0]
    @test slacks(deaaddinput, :Y) ≈ [0; 0; 0; 1.0; 0.0]

    # Output orientation
    deaaddoutput = deaadd(X, Y, orient = :Output)
    @test efficiency(deaaddoutput) ≈ [0; 0; 0; 2.0; 1.0]
    @test slacks(deaaddoutput, :X) ≈ [0; 0; 0; 0.0; 1.0]
    @test slacks(deaaddoutput, :Y) ≈ [0; 0; 0; 2.0; 1.0]

    # Input orientation with Weak disposability in outputs
    deaaddinputweak = deaadd(X, Y, orient = :Input, disposY = :Weak)
    @test efficiency(deaaddinputweak) ≈ [0; 0; 0; 0.0; 2.0]
    @test slacks(deaaddinputweak, :X) ≈ [0; 0; 0; 0.0; 2.0]
    @test slacks(deaaddinputweak, :Y) ≈ [0; 0; 0; 0.0; 0.0]

    # Output orientation with Weak disposability in inputs
    deaaddoutputweak = deaadd(X, Y, orient = :Output, disposX = :Weak)
    @test efficiency(deaaddoutputweak) ≈ [0; 0; 0; 2.0; 0.0]
    @test slacks(deaaddoutputweak, :X) ≈ [0; 0; 0; 0.0; 0.0]
    @test slacks(deaaddoutputweak, :Y) ≈ [0; 0; 0; 2.0; 0.0]

    # Test errors with orientation and disposability
    @test_throws ArgumentError deaadd([1; 2 ; 3], [4 ; 5; 6], orient = :Graph, disposX = :Error)  # Invalid inputs disposability
    @test_throws ArgumentError deaadd([1; 2 ; 3], [4 ; 5; 6], orient = :Graph, disposY = :Error)  # Invalid output disposability
    @test_throws ArgumentError deaadd([1; 2 ; 3], [4 ; 5; 6], orient = :Graph, disposX = :Weak)  # Weak disposability not possible in graph oriented model
    @test_throws ArgumentError deaadd([1; 2 ; 3], [4 ; 5; 6], orient = :Graph, disposX = :Weak)  # Weak disposability not possible in graph oriented model

    @test_throws ArgumentError deaadd([1; 2 ; 3], [4 ; 5; 6], orient = :Input, disposX = :Weak)  # Weak input disposability not possible in input oriented model

    @test_throws ArgumentError deaadd([1; 2 ; 3], [4 ; 5; 6], orient = :Output, disposY = :Weak)  # Weak output disposability not possible in output oriented model

    # ------------------
    # Test Vector and Matrix inputs and outputs
    # ------------------
    # Tests against results in R

    # Inputs is Matrix, Outputs is Vector
    X = [2 2; 1 4; 4 1; 4 3; 5 5; 6 1; 2 5; 1.6	8]
    Y = [1; 1; 1; 1; 1; 1; 1; 1]

    @test efficiency(deaadd(X, Y, :Ones)) ≈ [0.0; 0.0; 0.0; 3.0; 6.0; 2.0; 3.0; 5.2]

    # Inputs is Vector, Output is Matrix
    X = [1; 1; 1; 1; 1; 1; 1; 1]
    Y = [7 7; 4 8; 8 4; 3 5; 3 3; 8 2; 6 4; 1.5 5]

    @test efficiency(deaadd(X, Y, :Ones)) ≈ [0.0; 0.0; 0.0; 6.0; 8.0; 2.0; 4.0; 7.5]

    # Inputs is Vector, Output is Vector
    X = [2; 4; 8; 12; 6; 14; 14; 9.412]
    Y = [1; 5; 8; 9; 3; 7; 9; 2.353]

    @test efficiency(deaadd(X, Y, :Ones)) ≈ [0; 0; 0; 0;4; 7.333333333; 2; 8.059]

    # ------------------
    # RAM and BAM with orientation
    # ------------------

    X = [1; 2; 3; 2; 4]
    Y = [2; 3; 4; 1; 3]

    @test efficiency(deaadd(X, Y, :RAM, orient = :Graph)) ≈ [0.0; 0.0; 0.0; 1/3; 1/3]
    @test efficiency(deaadd(X, Y, :RAM, orient = :Input)) ≈ [0.0; 0.0; 0.0; 1/3; 2/3]
    @test efficiency(deaadd(X, Y, :RAM, orient = :Output)) ≈ [0.0; 0.0; 0.0; 2/3; 1/3]

    @test efficiency(deaadd(X, Y, :BAM, orient = :Graph)) ≈ [0.0; 0.0; 0.0; 2/3; 2/3]
    @test efficiency(deaadd(X, Y, :BAM, orient = :Input)) ≈ [0.0; 0.0; 0.0; 1; 2/3]
    @test efficiency(deaadd(X, Y, :BAM, orient = :Output)) ≈ [0.0; 0.0; 0.0; 2/3; 1]

end


text:
import RoME
include(joinpath(pkgdir(RoME), "test", "runtests.jl"))

ext: jl


text:
module RandomMonad

export Distribution

export Bernoulli, Bind, Binomial, Categorical, CloseClose, CloseOpen, counts,
       Counts, discretize!, Exponential, Fill, Filter, FisherYates, Iterate,
       Iterated, Join, Keep, Lift, Map,
       MixtureModel, Multinomial, Normal, OpenClose, OpenOpen, Pack,
       pmf, Poisson, Pure,
       Reduce, SelfAvoid, Shuffle, SubIter, SubSeq, support, Thunk, Uniform,
       Unique, variate_size, Zip,
       normalize!, denormalize!

using Random: AbstractRNG, gentype, randexp, randn, Random, randsubseq!,
              Repetition

using Statistics: mean

import Random: rand, rand!, Sampler

"""
    Distribution{T}

An instance of a subtype of `Distribution{T}` is an object able
to produce random values of type `T` via `rand`-related functions.
"""
abstract type Distribution{T} end

"""
    Base.eltype(::Type{<:Distribution{T}}) where {T}

Return the type parameter `T` of a distribution, i.e. the type of
values drawn from it.
"""
Base.eltype(::Type{<:Distribution{T}}) where {T} = T

"""
    reset!(sp::Random.Sampler, [n::Integer]) -> sp

Prepare the sampler state for the generation of `n` new random objects, or of
an arbitrary number when `n` is not specified, and return `sp`. This is useful
for samplers which have "memory", where an internal state has to be
re-initialized. The default implementation has no side effects.

When the sampler type for a given distribution depends on the `Repetition`
argument `rep` passed to its constructor, `reset!` generally makes sense only
for the case `rep == Val(Inf)`.

This function should be called in each function generating a collection of
objects generated from a given sampler `sp`. As a notable exception,
`rand!(::AbstractRNG, ::AbstractArray, ::Sampler)` does not follow this rule,
as it is defined in `Random`. Therefore, samplers having a non-trivial `reset!`
should not assume that this function will be called (one easy way to handle
that is to call `reset!` at construction time).
However, a similar method `rand!(::AbstractRNG, ::AbstractArray, ::SamplerReset)`
is defined, which does call `reset!`. It is recommended that samplers with a
non-trivial `reset!` inherit from `SamplerReset`.
"""
reset!(sp::Sampler, _ = 0) = sp

## extended rand!

"""
    rand!([rng::AbstractRNG], A::AbstractArray, X::Distribution, ::Val{1})

Populate the array `A` with random values, according to the distribution specified
by `X`, whose `eltype` must be compatible with the type of `A`.
`A` will be resized as needed when possible.

!!! note
    The standard `rand!(A, X)` syntax is used to populate `A` according to
    distribution `X` which describes how to generate individual elements. On
    the other hand, in `rand!(A, X, Val(1))`, `X` describes how to generate a
    full array.

# Examples
```julia
julia> rand!(zeros(Int, 3), Fill(1:3, 4), Val(1))
4-element Array{Int64,1}:
 2
 3
 2
 1

julia> rand!(zeros(Int, 2, 2), Fill(1:3, 4, 1), Val(1))
ERROR: ArgumentError: can not resize destination array
Stacktrace:
[...]
```
"""
rand!(rng::AbstractRNG, A, X::Distribution, ::Val{N}) where {N} =
    rand!(rng, A, Sampler(rng, X), Val(N))

rand!(A, X::Distribution, ::Val{N}) where {N} =
    rand!(Random.GLOBAL_RNG, A, X, Val(N))

## SamplerReset

"""
    abstract type SamplerReset{T} <: Sampler{T} end

`SamplerReset` is meant to be as general as `Sampler`, but with the explicit
indication that `reset!` might be non-trivial.
In particular, `rand!([rng], a::AbstractArray, sp::SamplerReset)` is specialized
(including a call to `reset!(sp, length(a))` before populating `a`).
"""
abstract type SamplerReset{T} <: Sampler{T} end

function rand!(rng::AbstractRNG, A::AbstractArray, sp::SamplerReset)
    reset!(sp, length(A))
    for i in eachindex(A)
        @inbounds A[i] = rand(rng, sp)
    end
    A
end

## misc

if isdefined(Base, :require_one_based_indexing) # VERSION >= v"1.2"
    using Base: require_one_based_indexing
else
    require_one_based_indexing(A...) =
        !Base.has_offset_axes(A...) ||
        throw(ArgumentError("offset arrays are not supported but got an array with index other than 1"))
end

## includes

include("samplerhelpers.jl")
include("pmf.jl")
include("wrap.jl")
include("distributions.jl")
include("floatintervals.jl")
include("containers.jl")
include("adapters.jl")
include("sequences.jl")
include("util.jl")

end # module


text:
using HOODESolver
using DifferentialEquations
using LinearAlgebra
using Plots
using Random

function getmindif(tab::Vector{Vector{BigFloat}})
    nmmin = Inf
    ret = (0, 0)
    for i=1:size(tab,1)
        for j=(i+1):size(tab,1)
            nm = norm(tab[i]-tab[j], Inf)
            if nm < nmmin
                nmmin = nm
                ret = i, j
            end
        end
    end
    return ret, nmmin
end
function fct(u, p, t)
    s1, c1 = sincos(u[1]/2)
    s2, c2 = sincos(u[2])
    s3, c3 = sincos(u[3])
    return [0, 0, u[6], c1*s2*s3/2, s1*c2*s3, s1*s2*c3]
end

function fctmain(n_tau, prec)
    setprecision(512)
    Random.seed!(561909)
    u0 = 2rand(BigFloat,6)-ones(BigFloat,6)
    setprecision(prec)
    u0 = BigFloat.(u0)
    println("u0=$u0")
    t_max = big"1.0"
    epsilon=big"1e-7"
    println("epsilon=$epsilon")
    nbmaxtest=15
    ordmax=9
    debord=3
    pasord=1
    y = ones(Float64, nbmaxtest, div(ordmax-debord,pasord)+1 )
    res_gen = Array{ Array{BigFloat,1}, 2}(undef, nbmaxtest, div(ordmax-debord,pasord)+1 )
    x=zeros(Float64,nbmaxtest)
    A = [0 0 0 1 0 0; 0 0 0 0 1 0;0 0 0 0 0 0; 0 0 0 0 1 0; 0 0 0 -1 0 0; 0 0 0 0 0 0]
    t_0=big"0.0"
    t_max=big"1.0"
    prob = HOODEProblem(fct, u0, (t_0, t_max), missing, A, epsilon)
    nm = NaN
    ordmax += 1
    ordprep=undef
    nb = 10*2^(nbmaxtest)
    solref=undef
    while isnan(nm)
        ordmax -= 1
        @time sol = solve(prob, nb_t=nb, order=ordmax, getprecision=false, nb_tau=n_tau, dense=false)
        solref = sol[end]
	    nm = norm(solref, Inf)
    end

    tabsol = Array{Array{BigFloat,1},1}(undef,1)

    tabsol[1] = solref

    println("ordmax=$ordmax")

    solref_gen = solref
    indref = 1

    ind=1
    for order=debord:pasord:ordmax
        ordprep=order+2
        println("eps=$epsilon solRef=$solref order=$order")
        nb = 10
        indc = 1
        labels=Array{String,2}(undef, 1, order-debord+1)
        resnorm=0
        resnormprec=1
        sol =undef
        par_u0=missing
        println("preparation ordre $order + 2")
        while indc <= nbmaxtest
            @time solall = solve(prob, nb_t=nb, order=order, getprecision=false, nb_tau=n_tau, par_u0=par_u0, dense=false)
            par_u0 = solall.par_u0
            sol = solall[end]
            push!(tabsol, sol)
            res_gen[indc,ind] = sol
            (a, b), nm = getmindif(tabsol)
            if a != indref
                println("New solref !!!! a=$a, b=$b nm=$nm")
                indref = a
                solref = tabsol[a]
                for i=1:ind
                    borne = (i <ind) ? size(res_gen,1) : indc
                    for j = 1:borne
                        nm2 = min( norm(res_gen[j,i] - solref, Inf), 1.1)
                        y[j,i] = nm2 == 0 ? nm : nm2
                   end
                end
            else
                diff=solref-sol
                y[indc,ind] = min(norm(diff,Inf), 1.1)
            end
            x[indc] = 1.0/nb
            println("result=$y")
            println("log2(y)=$(log2.(y))")
            nb *= 2
            indc += 1
        end
        for i=debord:pasord:order
            labels[1,(i-debord)÷pasord+1] = " eps,order=$(convert(Float32,epsilon)),$i "
        end
        p=Plots.plot(
                        x,
                        view(y,:,1:ind),
                        xlabel="delta t",
                        xaxis=:log,
                        ylabel="error",
                        yaxis=:log,
                        legend=:topleft,
                        label=labels,
                        marker=2
                    )
        prec_v = precision(BigFloat)
        eps_v = convert(Float32,epsilon)
        Plots.savefig(p,"out/r4_$(prec_v)_$(eps_v)_$(order)_$(n_tau)_particle.pdf")
        if resnorm > resnormprec
            break
        end
        resnormprec = resnorm
        ind+= 1
    end
end
fctmain(32,512)


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

using Test
using Example
using BigProject

@test BigProject.f() == 1

# #306 Pkg.jl
cd("LibFoo.jl") do
    run(`$(Base.julia_cmd()) --project test/runtests.jl`)
end

max_stars_repo_name: greimel/julia

max_issues_repo_name: greimel/julia

max_forks_repo_name: greimel/julia


text:
struct FranklinParserException <: Exception
    kind::Symbol
    msg::String
end

@inline parser_exception(k::Symbol, m::String) = throw(FranklinParserException(k, m))


text:
export SpikedWishart

struct SpikedWishart <: ContinuousMatrixDistribution
    beta::Integer
    n::Integer
    p::Integer
    spikes::Array{Float64, 1}
end

SpikedWishart(beta, n, p) = SpikedWishart(beta, n, p, [])

# SAMPLERS

function randreduced(d::SpikedWishart)
    if length(d.spikes) <= 1
        return randtridiagonal(d)
    else
        return randbanded(d)
    end
end

function randtridiagonal(d::SpikedWishart)
    a = d.beta * d.n/2

    # diagonal and superdiagonal of B
    Bdv = [rand(Chi(2*a - d.beta*k)) for k in 0:(d.p-1)]
    Bev = [rand(Chi(d.beta*(d.p-k))) for k in 1:(d.p-1)]

    if length(d.spikes) == 1
        Bdv[1] *= sqrt(1 + d.spikes[1])
    end

    # diagonal and off-diagonals of BB'
    dv = Bdv.^2 + [0; Bev.^2]
    ev = Bdv[1:end-1] .* Bev

    SymTridiagonal(dv, ev)
end

# Multi-spike version only implemented for beta=1
function randbanded(d::SpikedWishart)
    @assert d.beta == 1
    r = length(d.spikes)

    U = BandedMatrix{Float64}(undef, (d.p, d.p), (0, r))

    dv = [rand(Chi(d.n - k + 1)) for k in 1:d.p]
    @. dv[1:r] *= sqrt(1 + d.spikes)

    U[band(0)] .= dv

    for k = 1:(r-1)
        ev = randn(d.p - k)
        @. ev[1:(r-k)] *= sqrt(1 + d.spikes[(k+1):end])

        U[band(k)] .= ev
    end

    U[band(r)] .= [rand(Chi(d.p - k + 1)) for k in 1:(d.p-r)]

    Symmetric(U' * U)
end


text:
using CSV, DataFrames, wgregseq, CairoMakie

wgregseq.plotting_style.default_makie!()
##
# Open Gene information
dir = @__DIR__
home_dir = joinpath(split(dir, '/')[1:end-2]...)

## Metadata for samples
df_samples = CSV.read("/$home_dir/data/iModulon/sample_table.csv", DataFrames.DataFrame)
names(df_samples)

# Fill missing values
df_samples[!, ["Carbon Source (g/L)"]] = coalesce.(df_samples[!, ["Carbon Source (g/L)"]], "none")
df_samples[!, ["Strain Description"]] = coalesce.(df_samples[!, ["Strain Description"]], "none")

# Take only MG1655
df_samples_MG = df_samples[
    Array(df_samples[!, "Strain Description"]) .== "Escherichia coli K-12 MG1655",
     :]

sample_IDs = df_samples_MG.Column1
df_samples[!, "Strain Description"] |> unique
df_samples_MG[!, ["Carbon Source (g/L)"]] |> unique |> println
df_samples_MG[!, ["Nitrogen Source (g/L)"]] |> unique  |> println
df_samples_MG[!, ["Electron Acceptor"]] |> unique |> println
df_samples_MG[!, ["Trace Element Mixture"]] |> unique |> println
## iModulon metadata
df_im = CSV.read("/$home_dir/data/iModulon/iM_table.csv", DataFrames.DataFrame)
names(df_im)
df_im.category |> unique
df_im.n_genes |> sum
df_im[43, :]
println(df_im[df_im.category .== "Carbon Metabolism", ["name", "regulator_readable", "n_genes", ]])

println(filter(x-> occursin("lac", x), df_im.name))

## Gene - iModulon list
df_gene_im_list = CSV.read("/$home_dir/data/iModulon/gene_presence_list.csv", DataFrame)
df_gene_im_list.Gene |> unique |> length
## Gene - iModulon boolean table
df_gene_im_table = CSV.read("/$home_dir/data/iModulon/gene_presence_matrix.csv", DataFrame)

## Raw Counts
df_counts = CSV.read("/$home_dir/data/iModulon/Precise2dot0counts.csv", DataFrame)

## Normalized reads
df_tpm = CSV.read("/$home_dir/data/iModulon/log_tpm.csv", DataFrame)
df_tpm = df_tpm[!, vcat(["Column1"], sample_IDs)]

##
df_M = CSV.read("/$home_dir/data/iModulon/M.csv", DataFrame)

# Import gene information
df_genes = CSV.read("/$home_dir/data/all_genes_table.csv", DataFrame)

## Test lacI
lacZ_acc = df_genes[df_genes.gene .== "lacZ", "accession"][1]
lacZ_iM = df_gene_im_list[df_gene_im_list.Gene .== lacZ_acc, "iModulon"][1]
lacZ_iM_genes = df_gene_im_list[df_gene_im_list.iModulon .== lacZ_iM, "Gene"]
df_genes[map(x -> x in lacZ_iM_genes, df_genes.accession), "gene"]

## Normalized reads
df = CSV.read("/$home_dir/data/iModulon/A.csv", DataFrames.DataFrame)
Array(df_samples[!, "Base Media"]) .== "M9"
df_samples = df_samples[
    (Array(df_samples[!, "Strain Description"]) .== "Escherichia coli K-12 MG1655") .&
    (Array(df_samples[!, "Base Media"]) .== "M9"), :]
df = df[!, df_samples.Column1]

x = df[24, names(df)[2:end]]
gal_samples = df_samples[df_samples[!, "Carbon Source (g/L)"] .== "galactose(4)", "Column1"]

sig_ind = Array(x) .> 2.5
names(x[sig_ind])
df_samples.Column1
_df = DataFrame()
for y in names(x[sig_ind])

    append!(_df, df_samples[(df_samples.Column1 .== y), :])
end
_df

fig = Figure()
ax = Axis(fig[1, 1])
scatter!(Array(x))
fig

## Knockouts

gdf = groupby(df_im, "category")
gene_cat_df = DataFrame()
for _df in gdf
    gene_list = []
    for i in _df.k
        genes = df_gene_im_list[df_gene_im_list.iModulon .== i, "Gene"]
        gene_names = df_genes[map(x -> x in genes, df_genes.accession), "gene"]
        push!(gene_list, gene_names...)
    end
    append!(gene_cat_df, DataFrame(genes=gene_list, category=fill(unique(_df.category)[1], length(gene_list))))
end
gene_cat_df[gene_cat_df.genes .== "bhsA", :]

gene_cat_df

cat_per_gene = sort(combine(groupby(unique(gene_cat_df), "genes"), nrow), "nrow", rev=true)
mean(cat_per_gene.nrow)

gene_cat_df[map(x -> occursin("lac", x), gene_cat_df.genes), :]

##
unique_genes = df_gene_im_list.Gene |> unique
genes_precise = df_genes[map(x -> x in unique_genes, df_genes.accession), "gene"]
df_heinemann = DataFrames.DataFrame(Pandas.read_excel("data/heinemann_data.xlsx", sheet_name="Table S6", header=[2]))
genes_heinemann = df_heinemann.Gene |> unique

sum(map(x -> x in genes_heinemann, genes_precise))
# compare synomyms

## Some genes may have the wrong synomym
all_gene_list = CSV.read(
    "/$home_dir/data/all_genes_table.csv",
    DataFrame,
    types=Dict(
        "ID"=>String,
        "gene"=>String,
        "gene_position"=>Float64,
        "direction"=>String
    )
)

all_gene_list.synonyms = parse.(Vector{String}, all_gene_list.synonyms)
all_gene_list.transcription_units = parse.(Vector{String}, all_gene_list.transcription_units)

for l in [genes_precise, genes_heinemann]
    for i in 1:length(l)
        name = l[i]
        if name ∉ all_gene_list.gene
            if length(all_gene_list[map(x -> name in x, all_gene_list.synonyms), :gene]) == 0
                println(name)
            else
                syn = all_gene_list[map(x -> name in x, all_gene_list.synonyms), :gene][1]
                l[i] = syn
            end
        end
    end
end

sum(map(x -> x in genes_heinemann, genes_precise))


text:
mutable struct UpdateTrigger
    execution_wait_count::Int
    current_count::Int
end

function set_execution_wait_count!(trigger::UpdateTrigger, val::Int)
    trigger.execution_wait_count = val
    return
end

function update_count!(trigger::UpdateTrigger)
    trigger.current_count += 1
    return
end

function trigger_update(trigger::UpdateTrigger)
    return trigger.current_count == trigger.execution_wait_count
end

function reset_trigger_count!(trigger::UpdateTrigger)
    trigger.current_count = 0
    return
end

function initialize_trigger_count!(trigger::UpdateTrigger)
    trigger.current_count = trigger.execution_wait_count
    return
end

function get_execution_wait_count(trigger::UpdateTrigger)
    return trigger.execution_wait_count
end

############################ Chronologies For FeedForward ###################################
@doc raw"""
    Synchronize(periods::Int)
Defines the co-ordination of time between Two problems.

# Arguments
- `periods::Int`: Number of time periods to grab data from
"""
mutable struct Synchronize <: FeedForwardChronology
    periods::Int
    current::Int
    trigger::UpdateTrigger
    function Synchronize(; periods)
        new(periods, 0, UpdateTrigger(-1, -1))
    end
end
# TODO: Add DocString
"""
    RecedingHorizon(period::Int)
"""
mutable struct RecedingHorizon <: FeedForwardChronology
    periods::Int
    trigger::UpdateTrigger
    function RecedingHorizon(; periods::Int = 1)
        new(periods, UpdateTrigger(-1, -1))
    end
end

mutable struct Consecutive <: FeedForwardChronology
    trigger::UpdateTrigger
    function Consecutive()
        new(UpdateTrigger(-1, -1))
    end
end

mutable struct FullHorizon <: FeedForwardChronology
    trigger::UpdateTrigger
    function FullHorizon()
        new(UpdateTrigger(-1, -1))
    end
end

mutable struct Range <: FeedForwardChronology
    range::UnitRange{Int}
    trigger::UpdateTrigger
    function Range(; range::UnitRange{Int})
        new(range, UpdateTrigger(-1, -1))
    end
end


text:
Base.displayable(d::ElectronDisplayType, ::MIME{Symbol("application/vnd.vegalite.v2+json")}) = true
Base.displayable(d::ElectronDisplayType, ::MIME{Symbol("application/vnd.vegalite.v3+json")}) = true
Base.displayable(d::ElectronDisplayType, ::MIME{Symbol("application/vnd.vegalite.v4+json")}) = true
Base.displayable(d::ElectronDisplayType, ::MIME{Symbol("application/vnd.vega.v3+json")}) = true
Base.displayable(d::ElectronDisplayType, ::MIME{Symbol("application/vnd.vega.v4+json")}) = true
Base.displayable(d::ElectronDisplayType, ::MIME{Symbol("application/vnd.vega.v5+json")}) = true

function _maybe_fallback_to(mime, d, x, payload)
    if sizeof(payload) > d.config.max_json_bytes
        @warn string(
            "The size of JSON representation (", sizeof(payload), ")",
            " exceeds `max_json_bytes = ", d.config.max_json_bytes, "`.",
            " Falling back to $mime."
        )
        return display(d, mime, x)
    end
    return nothing
end

function _display_vegalite(d, major_version_vegalite, major_version_vega, x)
    payload = stringmime(MIME("application/vnd.vegalite.v$major_version_vegalite+json"), x)
    ans = _maybe_fallback_to(MIME"image/png"(), d, x, payload)
    ans === nothing || return ans
    displayplot(d, "vega-lite", JSON.JSONText(payload), options=Dict("webPreferences" => Dict("webSecurity" => false)))
end

function _display_vega(d, major_version, x)
    payload = stringmime(MIME("application/vnd.vega.v$major_version+json"), x)
    ans = _maybe_fallback_to(MIME"image/png"(), d, x, payload)
    ans === nothing || return ans
    displayplot(d, "vega", JSON.JSONText(payload), options=Dict("webPreferences" => Dict("webSecurity" => false)))
end

function Base.display(d::ElectronDisplayType, ::MIME{Symbol("application/vnd.vegalite.v2+json")}, x)
    return _display_vegalite(d, "2", "3", x)
end

function Base.display(d::ElectronDisplayType, ::MIME{Symbol("application/vnd.vegalite.v3+json")}, x)
    return _display_vegalite(d, "3", "5", x)
end

function Base.display(d::ElectronDisplayType, ::MIME{Symbol("application/vnd.vegalite.v4+json")}, x)
    return _display_vegalite(d, "4", "5", x)
end

function Base.display(d::ElectronDisplayType, ::MIME{Symbol("application/vnd.vega.v3+json")}, x)
    return _display_vega(d, "3", x)
end

function Base.display(d::ElectronDisplayType, ::MIME{Symbol("application/vnd.vega.v4+json")}, x)
    return _display_vega(d, "4", x)
end

function Base.display(d::ElectronDisplayType, ::MIME{Symbol("application/vnd.vega.v5+json")}, x)
    return _display_vega(d, "5", x)
end


text:
#function for segmentation in grids based on thresholding and connectivity
#
# based on functions from Images module

function im_labeling(tf::Matrix{U}; connectivity=8) where U <: Integer
    # create connectivity array
    if connectivity == 8
        connect_array = trues(3,3)
    elseif connectivity == 4
        connect_array = BitMatrix([0 1 0; 1 1 1; 0 1 0])
    end

    # label areas
    segments = label_components(tf, connect_array)
end

# copy of function above, not sure how to catch both types
function im_labeling(tf::BitMatrix; connectivity=8)
    # create connectivity array
    if connectivity == 8
        connect_array = trues(3,3)
    elseif connectivity == 4
        connect_array = BitMatrix([0 1 0; 1 1 1; 0 1 0])
    end

    # label areas
    segments = label_components(tf, connect_array)
end

"""cleanup binary images based on minimal area of segments
Input
segments        labeled array
min_area        min area of one segment [m2]
res             resolution of 1 pixel [m2]
connectivity    all surrounding pixels (=8) or only N-E-S-W pixels (=4) (default = 8)

output
updated label array
"""
function im_segments_cleanup(segments::Array{U, 2}, min_area::Real;
        res=1, connectivity=8) where U <: Integer
    # remove noise
    segments = copy(segments)
    segment_areas = component_lengths(segments)
    segment_idx = component_indices(segments)
    for lab in 1:maximum(segments) # skip background == 0
        area = segment_areas[lab+1]
        if area < (min_area / (res*res))
            for idx in segment_idx[lab+1]
                segments[idx] = 0
            end
        end
    end
    # re-label segments
    im_labeling(segments; connectivity=connectivity)
end

function split_segments(segments::Array{U, 2}, block_size, min_area;
    res=1, connectivity=8) where U <: Integer

    # split segments
    segments = copy(segments)
    segment_idx = component_subscripts(segments)
    segment_bbox = component_boxes(segments)
    block_size = Float64(block_size) / res
    nsegments = maximum(segments)
    for lab in 1:nsegments # skip background == 0
        (ymin, xmin), (ymax, xmax) = segment_bbox[lab+1]
        dx, dy = xmax-xmin+1, ymax-ymin+1
        if ((dx > 1.5*block_size) || (dy > 1.5*block_size))
            nrow, ncol = min(Int(cld(dy, block_size)),1), min(Int(cld(dx,block_size)),1)
            for (y, x) in segment_idx[lab+1]
                # cells include bottom and left boundaries
                col = Int(fld(x - xmin, block_size)+1)
                row = Int(fld(ymax - y, block_size)+1)
                iblock = sub2ind((nrow, ncol), row, col)
                segments[y,x] = iblock
            end
        end
    end
    # re-label segments
    segments = im_labeling(segments; connectivity=connectivity)

    # merge segments with too small area
    segment_areas = component_lengths(segments)
    segment_idx = component_indices(segments)
    for lab in 1:maximum(segments) # skip background == 0
        area = segment_areas[lab+1]
        if area < (min_area / (res*res))
            b = falses(segments)
            for idx in segment_idx[lab+1]
                b[idx] = true
            end
            bd = dilate(b)
            neighboring_cells = bd .& broadcast(~, b)
            new_lab = maximum(segments[neighboring_cells]) # find label neighboring segment
            new_lab == 0 && continue # no neighboring segments
            for idx in segment_idx[lab+1]
                segments[idx] = new_lab
            end
        end
    end
    # re-label segments
    im_labeling(segments; connectivity=connectivity)
end

"""reduce pixel in 'A' within a connected segment in binary grid 'tf' using function 'func'

Input
segments        labeled array
A               input image with to sample from can [Array{2,Real}]

output
image with reduced value in all segment pixels"""
function im_segments_reduce(segments::Array{U, 2}, A::Array{T, 2};
                            connectivity=8,
                            func=x -> maximum(x)-minimum(x),
                            nodata=-9999.0) where U <: Integer where T <: AbstractFloat
    @assert size(segments) == size(A) "the size of image A and binary grid should match"

    # get indices
    initiate = true
    B = 0

    for idx in component_indices(segments)
        if initiate # skip background, but initiate output
            # initiate output array with output type at first loop
            vals = func(T[1, 2])
            @assert isa(vals, Real) "incorrect reducer function"
            reduced_type = typeof(vals)
            B = fill(reduced_type(nodata), size(A))
            initiate = false
        else
            vals = [xi for xi=A[idx] if xi != nodata]
            if !isempty(vals)
                B[idx] = func(vals)
            end
        end
    end

    B
end


text:
"""

Stochastic average gradient method. Combining with a first-order optimization method results
in a variance-reduced method for finite-sum optimization.
"""
module SAG

export StochasticGradient, update!, initialized_fraction

"""
    StochasticGradient{T,N} <: AbstractArray{T,N}

Stochastic average gradient (SAG) method for tracking a gradient made up of the sum of several
sub-gradients. Combining SAG with an iterative method for first-order optimization method (e.g.,
gradient descent) results in a variance-reduced method for finite-sum optimization, e.g., empirical
risk minimization. Because it's a variance-reduced method, the iterate converges to the optimum
even when the gradients used to update the iterate are stochastic.

External links:

* [Minimizing Finite Sums with the Stochastic Average Gradient](https://arxiv.org/abs/1309.2388)

```julia

# initialization
n = 3                                       # number of sub-gradients
dims = (10,)                                # gradient dimension
sg = StochasticGradient(Float32, n, dims)

# compute one sub-gradient at a time
i = 1                                       # sub-gradient index
∇i = randn(dims)                            # replace with actual sub-gradient computation
update!(sg, i, ∇i)                          # writes ∇i into the gradient
sg ≈ ∇i                                     # = true; sg acts as an array

# compute another sub-gradient
j = 2
∇j = randn(dims)                            # replace with actual sub-gradient computation
update!(sg, j, ∇j)
sg ≈ ∇i + ∇j                                # = true; sg acts as an array

# using StochasticGradient for variance-reduced stochastic gradient descent
w = randn(dims)                             # iterate, e.g., a machine learning model
stepsize = 0.1                              # step size used for model updates
f = initialized_fraction(sg)                # fraction of sub-gradients that have been initialized; 2/3 in our case
w .-= (stepsize / f) .* sg                  # gradient descent step

```

"""
mutable struct StochasticGradient{T,N} <: AbstractArray{T,N}
    ∇s::Vector{Array{T,N}}          # sub-gradients, or components
    ∇::Array{T,N}                   # gradient, i.e., the sum of ∇s
    isinitialized::Vector{Bool}     # indicates which sub-gradients that have been initialized
    ninitialized::Int               # number of sub-gradients that have been initialized
end

function StochasticGradient(∇s::Vector{T}, ∇::T; overwrite=true) where T
    n = length(∇s)
    0 < n || throw(ArgumentError("n must be positive"))
    for i in 2:n
        size(∇s[i]) == size(∇s[1]) || throw(DimensionMismatch("The $(i)-th component has dimensions $(size(∇s[i])), but the first component has dimensions $(size(∇s[1]))"))
    end
    if overwrite
        ∇ .= ∇s[1]
        for i in 2:n
            ∇ .+= ∇s[i]
        end
    end
    StochasticGradient(∇s, ∇, zeros(Bool, n), 0)
end

function StochasticGradient(∇s::AbstractVector)
    length(∇s) > 0 || throw(ArgumentError("There must be at least 1 component"))
    StochasticGradient(∇s, similar(∇s[1]))
end

function StochasticGradient(∇::AbstractArray, n::Integer)
    0 < n || throw(ArgumentError("n must be positive"))
    ∇s = [similar(∇) for _ in 1:n]
    for ∇i in ∇s
        ∇i .= ∇ ./ n
    end
    StochasticGradient(∇s, ∇)
end

function StochasticGradient(T::DataType, n::Integer, dims::Tuple)
    0 < n || throw(ArgumentError("n must be positive"))
    StochasticGradient([zeros(T, dims) for _ in 1:n], zeros(T, dims))
end

StochasticGradient(n::Integer, dims::Tuple) = StochasticGradient(Float64, n, dims)

Base.length(sg::StochasticGradient) = length(sg.∇)
Base.size(sg::StochasticGradient, args...) = size(sg.∇, args...)
Base.getindex(sg::StochasticGradient, args...) = getindex(sg.∇, args...)
Base.view(sg::StochasticGradient, args...) = view(sg.∇, args...)
Base.eltype(sg::StochasticGradient) = eltype(sg.∇)
Base.iterate(sg::StochasticGradient, args...; kwargs...) = iterate(sg.∇, args...; kwargs...)
Base.similar(sg::StochasticGradient, args...) = StochasticGradient(similar(sg.∇, args...), length(sg.∇s))
Base.show(io::IO, sg::StochasticGradient) = write(io, "SAG{eltype: $(eltype(sg)), n: $(length(sg.∇s)), dims: $(size(sg.∇))}")

"""
    update!(sg::StochasticGradient, i::Integer, ∇i::AbstractArray)

Set the `i`-th component to `∇i`, and update the overall sum `∇` accordingly.
"""
function update!(sg::StochasticGradient, i::Integer, ∇i::AbstractArray)
    0 < i <= length(sg.∇s) || throw(ArgumentError("i is $i, but there are $(length(sg.∇s)) components"))
    size(∇i) == size(sg.∇) || throw(DimensionMismatch("∇i has dimensions $(size(∇i)), but ∇ has dimensions $(size(sg.∇))"))
    sg.∇ .+= ∇i .- sg.∇s[i]
    sg.∇s[i] .= ∇i
    if !sg.isinitialized[i]
        sg.ninitialized += 1
        sg.isinitialized[i] = true
    end
    return ∇i
end

"""
    update!(sg::StochasticGradient, pair::Pair{<:Integer, <:AbstractArray})

"""
update!(sg::StochasticGradient, pair::Pair{<:Integer, <:AbstractArray}) = update!(sg, pair[1], pair[2])

"""
    update!(sg::StochasticGradient, pairs)

Perform an update for each element of `pairs`, which must be an iterable of collections of length 2
(e.g., pairs or tuples), where the second element is a component and the first is the index of that
component.
"""
function update!(sg::StochasticGradient, pairs)
    if 2*length(pairs) <= length(sg.∇s)
        for (i, ∇i) in pairs
            update!(sg, i, ∇i)
        end
    else
        for (i, ∇i) in pairs
            sg.∇s[i] .= ∇i
            if !sg.isinitialized[i]
                sg.ninitialized += 1
                sg.isinitialized[i] = true
            end
        end
        sg.∇ .= sg.∇s[1]
        for i in 2:length(sg.∇s)
            sg.∇ .+= sg.∇s[i]
        end
    end
    return pairs
end

initialized_fraction(::Any) = 1.0

"""

Return the fraction of components of `sg` that have been initialized.
"""
initialized_fraction(sg::StochasticGradient) = sg.ninitialized / length(sg.∇s)

end


text:
using DReal
# using Gadfly

## Function Approximation
## ======================
f(x) = sin(x) + 2x

npowers = 5
θ = [DReal.Var(Float64, -10.0, 10.0) for i = 1:npowers]
g(x) = sum([θ[i]*(x^i) for i = 1:(npowers - 1)]) + θ[npowers]

x = ForallVar(global_context(), Float64, -1000.0, 1000.0)
constraint = abs(g(x) - f(x)) < 1.0

add!(constraint)
result = model(θ...)
@show result
# pop_ctx!()


text:
@inline maybe_particles(x) = x
@inline maybe_particles(p::AbstractParticles) = p.particles
@inline maybe_logweights(x) = 0
@inline maybe_logweights(p::WeightedParticles) = p.logweights

function register_primitive(ff, eval=eval)
    register_primitive_multi(ff, eval)
    register_primitive_single(ff, eval)
end

function register_primitive_multi(ff, eval=eval)
    f = nameof(ff)
    m = Base.parentmodule(ff)
    for PT in (:Particles, :StaticParticles)
        eval(quote
            function ($m.$f)(p::$PT{T,N},a::Real...) where {T,N}
                $PT{T,N}(($m.$f).(p.particles, MonteCarloMeasurements.maybe_particles.(a)...)) # maybe_particles introduced to handle >2 arg operators
            end
            function ($m.$f)(a::Real,p::$PT{T,N}) where {T,N}
                $PT{T,N}(map(x->($m.$f)(a,x), p.particles))
            end
            function ($m.$f)(p1::$PT{T,N},p2::$PT{T,N}) where {T,N}
                $PT{T,N}(map(($m.$f), p1.particles, p2.particles))
            end
            function ($m.$f)(p1::$PT{T,N},p2::$PT{S,N}) where {T,S,N} # Needed for particles of different float types :/
                $PT{promote_type(T,S),N}(map(($m.$f), p1.particles, p2.particles))
            end
        end)
    end
    # The code below is resolving some method ambiguities
    eval(quote
        function ($m.$f)(p1::StaticParticles{T,N},p2::Particles{T,N}) where {T,N}
            StaticParticles{T,N}(map(($m.$f), p1.particles, p2.particles))
        end
        function ($m.$f)(p1::StaticParticles{T,N},p2::Particles{S,N}) where {T,S,N} # Needed for particles of different float types :/
            StaticParticles{promote_type(T,S),N}(map(($m.$f), p1.particles, p2.particles))
        end

        function ($m.$f)(p1::Particles{T,N},p2::StaticParticles{T,N}) where {T,N}
            StaticParticles{T,N}(map(($m.$f), p1.particles, p2.particles))
        end
        function ($m.$f)(p1::Particles{T,N},p2::StaticParticles{S,N}) where {T,S,N} # Needed for particles of different float types :/
            StaticParticles{promote_type(T,S),N}(map(($m.$f), p1.particles, p2.particles))
        end
    end)
end

function register_primitive_single(ff, eval=eval)
    f = nameof(ff)
    m = Base.parentmodule(ff)
    for PT in (:Particles, :StaticParticles)
        eval(quote
            function ($m.$f)(p::$PT)
                $PT(map(($m.$f), p.particles))
            end
        end)
    end
end


text:
pwd() == joinpath(@__DIR__, "bin") && cd(@__DIR__) # allow starting app from bin/ dir

using Klingelbeutel
push!(Base.modules_warned_for, Base.PkgId(Klingelbeutel))
Klingelbeutel.main()

max_stars_repo_name: Kompendix/Klingelbeutel

max_issues_repo_name: Kompendix/Klingelbeutel

max_forks_repo_name: Kompendix/Klingelbeutel


text:
using DashBootstrapComponents, DashHtmlComponents

card = dbc_card([
    dbc_cardheader(
        dbc_tabs(
            [
                dbc_tab(label = "Tab 1", tab_id = "tab-1"),
                dbc_tab(label = "Tab 2", tab_id = "tab-2"),
            ],
            id = "card-tabs",
            active_tab = "tab-1",
        ),
    ),
    dbc_cardbody(html_p(id = "card-content", className = "card-text")),
]);

callback!(
    app,
    Output("card-content", "children"),
    Input("card-tabs", "active_tab"),
) do active_tab
    return "This is tab $active_tab"
end;

max_stars_repo_name: glsdown/dash-bootstrap-components

max_issues_repo_name: glsdown/dash-bootstrap-components

max_forks_repo_name: glsdown/dash-bootstrap-components


text:
using Statistics

function format(x)
  if abs(x) < 999
    return @sprintf("%12.7f",x)
  else
    return @sprintf("%12.5e",x)
  end
end

import Base.write
"""

```
write(R::ComplexMixtures.Result, filename::String, solute::Selection, solvent::Selection)
```

Function to write the final results to output files as simple tables
that are human-readable and easy to analyze with other software

If the solute and solvent selections are provides, pass on the atom
names.

"""
write(R::Result, filename::String,
      solute::Selection,
      solvent::Selection) =
  write(R,filename,solute_names=solute.names,solvent_names=solvent.names)

"""

```
write(R::ComplexMixtures.Result, filename::String;
      solute_names::Vector{String} = ["nothing"],
      solvent_names::Vector{String} = ["nothing"])
```

Optional passing of atom names.

"""
function write(R::Result, filename::String;
               solute_names::Vector{String} = ["nothing"],
               solvent_names::Vector{String} = ["nothing"])

  # Names of output files containing atomic contibutions
  atom_contrib_solvent = FileOperations.remove_extension(filename)*
                                 "-ATOM_CONTRIB_SOLVENT."*
                            FileOperations.file_extension(filename)
  atom_contrib_solute = FileOperations.remove_extension(filename)*
                            "-ATOM_CONTRIB_SOLUTE."*
                            FileOperations.file_extension(filename)

  #
  # GMD computed with minimum distance
  #

  # Conversion factor for volumes (as KB integrals), from A^3 to cm^3/mol
  mole = 6.022140857e23
  convert = mole / 1.e24

  output = open(filename,"w")
  println(output,@sprintf("#"))
  println(output,@sprintf("# Output of ComplexMixtures - MDDF"))
  println(output,@sprintf("# Trajectory files and weights:"))
  for i in 1:length(R.files)
    println(output,"#  $(R.files[i]) - w = $(R.weights[i])")
  end
  println(output,@sprintf("#"))
  println(output,@sprintf("# Density of solvent in simulation box (sites/A^3): %15.8f",R.density.solvent))
  println(output,@sprintf("# Density of solvent in bulk (estimated) (sites/A^3): %15.8f",R.density.solvent_bulk))
  println(output,@sprintf("# Molar volume of solvent in simulation (cc/mol): %15.8f",convert/R.density.solvent))
  println(output,@sprintf("# Molar volume of solvent in bulk (estimated) (cc/mol): %15.8f",convert/R.density.solvent_bulk))
  println(output,@sprintf("#"))
  println(output,@sprintf("# Number of atoms solute: %i9",size(R.solute_atom,2)))
  println(output,@sprintf("# Number of atoms of the solvent: %i9",size(R.solvent_atom,2)))
  println(output,@sprintf("#"))

  if R.options.usecutoff
    ibulk = setbin(R.options.dbulk,R.options.binstep)
  else
    ibulk = round(Int,R.nbins-1/R.options.binstep)
  end
  bulkerror = Statistics.mean( R.mddf[ibulk:R.nbins] )
  sdbulkerror = Statistics.std( R.mddf[ibulk:R.nbins] )
  println(output,"#")
  println(output,@sprintf("# Average and standard deviation of bulk-gmd: %12.5f +/- %12.5f",bulkerror,sdbulkerror))
  println(output,"#")

  # Output table

  println(output,"# COLUMNS CORRESPOND TO: ")
  println(output,"#       1  Minimum distance to solute (dmin)")
  println(output,"#       2  GMD distribution (md count normalized by md count of random-solute distribution)")
  println(output,"#       3  Kirwood-Buff integral (cc/mol) computed [(1/bulkdensity)*(col(6)-col(7))].")
  println(output,"#       4  Minimum distance site count for each dmin.")
  println(output,"#       5  Minimum distance site count for each dmin for random solute distribution.")
  println(output,"#       6  Cumulative number of molecules within dmin in the simulation.")
  println(output,"#       7  Cumulative number of molecules within dmin for random solute distribution.")
  println(output,"#       8  Volume of the shell of distance dmin and width binstep.")
  println(output,"#")
  println(output,"#   1-DISTANCE         2-GMD      3-KB INT    4-MD COUNT  5-COUNT RAND      6-SUM MD    7-SUM RAND   8-SHELL VOL")

  for i in 1:R.nbins
    line = "  "*format(R.d[i])                                 #  1-DISTANCE
    line = line*"  "*format(R.mddf[i])                         #  2-GMD
    line = line*"  "*format(R.kb[i])                           #  3-KB INT
    line = line*"  "*format(R.md_count[i])                     #  4-MD COUNT
    line = line*"  "*format(R.md_count_random[i])              #  5-COUNT RAND
    line = line*"  "*format(R.sum_md_count[i])                 #  6-SUM MD
    line = line*"  "*format(R.sum_md_count_random[i])          #  7-SUM RAND
    line = line*"  "*format(R.volume.shell[i])                 #  8-SHELL VOL
    println(output,line)
  end
  close(output)

  # Writting gmd per atom contributions for the solvent

  output = open(atom_contrib_solvent,"w")
  println(output,"#")
  println(output,"# Solvent atomic contributions to total MDDF. ")
  println(output,"#")
  println(output,"# Trajectory files: ", R.files)
  println(output,"#")
  println(output,"# Atoms: ")
  for i in 1:size(R.solvent_atom,2)
    if solvent_names[1] == "nothing"
      println(output,@sprintf("# %9i",i))
    else
      println(output,@sprintf("# %9i %5s",i, solvent_names[i]))
    end
  end
  println(output,"#")
  string = "#   DISTANCE     GMD TOTAL"
  for i in 1:size(R.solvent_atom,2)
    string = string*@sprintf("  %12i",i)
  end
  println(output,string)
  for i in 1:R.nbins
    string = format(R.d[i])
    string = string*"  "*format(R.mddf[i])
    for j in 1:size(R.solvent_atom,2)
      string = string*"  "*format(R.solvent_atom[i,j])
    end
    println(output,string)
  end
  close(output)

  # Writting gmd per atom contributions for the solute

  output = open(atom_contrib_solute,"w")
  println(output,"#")
  println(output,"# Solute atomic contributions to total MDDF. ")
  println(output,"#")
  println(output,"# Trajectory files: ", R.files)
  println(output,"#")
  println(output,"# Atoms:")
  for i in 1:size(R.solute_atom,2)
    if solute_names[1] == "nothing"
      println(output,@sprintf("# %9i",i))
    else
      println(output,@sprintf("# %9i %5s",i, solute_names[i]))
    end
  end
  println(output,"#")
  string = "#   DISTANCE      GMD TOTAL"
  for i in 1:size(R.solute_atom,2)
    string = string*@sprintf("  %12i",i)
  end
  println(output,string)
  for i in 1:R.nbins
    string = format(R.d[i])*"  "*format(R.mddf[i])
    for j in 1:size(R.solute_atom,2)
      string = string*"  "*format(R.solute_atom[i,j])
    end
    println(output,string)
  end
  close(output)

  # Write final messages with names of output files and their content

  println()
  println(" OUTPUT FILES: ")
  println()
  println(" Wrote solvent atomic GMD contributions to file: ", atom_contrib_solvent)
  println(" Wrote solute atomic GMD contributions to file: ", atom_contrib_solute)
  println()
  println(" Wrote main output file: ", filename)
  println()

  return nothing
end


text:
@rule AR(:γ, Marginalisation) (q_y_x::MultivariateNormalDistributionsFamily, q_θ::NormalDistributionsFamily, meta::ARMeta) = begin
    order = getorder(meta)
    F     = getvform(meta)

    y_x_mean, y_x_cov = mean_cov(q_y_x)
    θ_mean, θ_cov = mean_cov(q_θ)

    mA, Vθ = as_companion_matrix(θ_mean), θ_cov
    my, Vy = ar_slice(F, y_x_mean, 1:order), ar_slice(F, y_x_cov, 1:order, 1:order)
    mx, Vx = ar_slice(F, y_x_mean, (order + 1):2order), ar_slice(F, y_x_cov, (order + 1):2order, (order + 1):2order)
    Vyx    = ar_slice(F, y_x_cov, (order + 1):2order, 1:order)

    C = rank1update(Vx, mx)
    R = rank1update(Vy, my)
    L = rank1update(Vyx, mx, my)

    B = first(R) - 2 * first(mA * L) + first(mA * C * mA') + mul_trace(Vθ, C)

    return GammaShapeRate(convert(eltype(B), 3//2), B / 2)
end


text:
export assessfitness!,
    PopulationModel

mutable struct PopulationModel{G, F<:AbstractFitness{>:G}} <: AbstractPopulationModel{G}
    population::Population{G}
    fitness::F
    fittest::Individual{G}

    PopulationModel(population::Population{G}, fitness) where {G} =
        new{G, typeof(fitness)}(population, fitness)
end

assessfitness!(model::PopulationModel) = model.fittest = maximumby(i -> assess!(i, model.fitness),
                                                                   model.population)


text:
using OrdinaryDiffEq: phi, phi, phiv, phiv_timestep, expv, expv_timestep, arnoldi, getH, getV
using OrdinaryDiffEq: WOperator, set_gamma!, calc_W!
using OrdinaryDiffEq, LinearAlgebra, SparseArrays, Random, Test, DiffEqOperators

@testset "Exponential Utilities" begin
  # Scalar phi
  K = 4
  z = 0.1
  P = fill(0., K+1); P[1] = exp(z)
  for i = 1:K
    P[i+1] = (P[i] - 1/factorial(i-1))/z
  end
  @test phi(z, K) ≈ P

  # Matrix phi
  A = [0.1 0.2; 0.3 0.4]
  P = Vector{Matrix{Float64}}(undef, K+1); P[1] = exp(A)
  for i = 1:K
    P[i+1] = (P[i] - 1/factorial(i-1)*I) / A
  end
  @test phi(A, K) ≈ P

  # Krylov
  n = 20; m = 5
  Random.seed!(0)
  A = randn(n, n)
  t = 1e-2
  b = randn(n)
  @test exp(t * A) * b ≈ expv(t, A, b; m=m)
  P = phi(t * A, K)
  W = fill(0., n, K+1)
  for i = 1:K+1
    W[:,i] = P[i] * b
  end
  Ks = arnoldi(A, b; m=m)
  W_approx = phiv(t, Ks, K)
  @test W ≈ W_approx

  # Happy-breakdown in Krylov
  v = normalize(randn(n))
  A = v * v' # A is Idempotent
  Ks = arnoldi(A, b)
  @test Ks.m == 2

  # Arnoldi vs Lanczos
  A = Hermitian(randn(n, n))
  Aperm = A + 1e-10 * randn(n, n) # no longer Hermitian
  w = expv(t, A, b; m=m)
  wperm = expv(t, Aperm, b; m=m)
  @test w ≈ wperm

  # Internal time-stepping for Krylov (with adaptation)
  n = 100
  K = 4
  t = 5.0
  tol = 1e-7
  A = spdiagm(-1=>ones(n-1), 0=>-2*ones(n), 1=>ones(n-1))
  B = randn(n, K+1)
  Phi_half = phi(t/2 * A, K)
  Phi = phi(t * A, K)
  uhalf_exact = sum((t/2)^i * Phi_half[i+1] * B[:,i+1] for i = 0:K)
  u_exact = sum(t^i * Phi[i+1] * B[:,i+1] for i = 0:K)
  U = phiv_timestep([t/2, t], A, B; adaptive=true, tol=tol)
  @test norm(U[:,1] - uhalf_exact) / norm(uhalf_exact) < tol
  @test norm(U[:,2] - u_exact) / norm(u_exact) < tol
  # p = 0 special case (expv_timestep)
  u_exact = Phi[1] * B[:, 1]
  u = expv_timestep(t, A, B[:, 1]; adaptive=true, tol=tol)
  @test norm(u - u_exact) / norm(u_exact) < tol
end

@testset "Derivative Utilities" begin
  @testset "WOperator" begin
    Random.seed!(0); y = zeros(2); b = rand(2)
    mm = I; _J = rand(2,2)
    W = WOperator(mm, 1.0, DiffEqArrayOperator(_J))
    set_gamma!(W, 2.0)
    _W = mm - 2.0 * _J
    @test convert(AbstractMatrix,W) ≈ _W
    @test W * b ≈ _W * b
    mul!(y, W, b); @test y ≈ _W * b
  end

  @testset "calc_W!" begin
    A = [-1.0 0.0; 0.0 -0.5]; mm = [2.0 0.0; 0.0 1.0]
    u0 = [1.0, 1.0]; tmp = zeros(2)
    tspan = (0.0,1.0); dt = 0.01; dtgamma = 0.5dt
    concrete_W = mm - dtgamma * A
    concrete_Wt = mm/dtgamma - A

    # Out-of-place
    fun = ODEFunction((u,p,t) -> A*u;
                      mass_matrix=mm,
                      jac=(u,p,t) -> A)
    integrator = init(ODEProblem(fun,u0,tspan), ImplicitEuler(); adaptive=false, dt=dt)
    W = calc_W!(integrator, integrator.cache, dtgamma, false)
    @test convert(AbstractMatrix, W) ≈ concrete_W
    @test W \ u0 ≈ concrete_W \ u0

    # In-place
    fun = ODEFunction((du,u,p,t) -> mul!(du,A,u);
                      mass_matrix=mm,
                      jac_prototype=DiffEqArrayOperator(A))
    integrator = init(ODEProblem(fun,u0,tspan), ImplicitEuler(); adaptive=false, dt=dt)
    calc_W!(integrator, integrator.cache, dtgamma, false)
    @test convert(AbstractMatrix, integrator.cache.W) ≈ concrete_W
    ldiv!(tmp, lu!(integrator.cache.W), u0); @test tmp ≈ concrete_W \ u0
  end

  @testset "Implicit solver with lazy W" begin
    A = sparse([-1.0 0.0; 0.0 -0.5])
    mm = sparse([2.0 0.0; 0.0 1.0])
    u0 = [1.0, 1.0]; tspan = (0.0,1.0)

    _f = (u,p,t) -> t*(A*u); _f_ip = (du,u,p,t) -> lmul!(t,mul!(du,A,u))
    fun1 = ODEFunction(_f; mass_matrix=mm)
    fun2 = ODEFunction(_f; mass_matrix=mm, jac=(u,p,t) -> t*A)
    fun1_ip = ODEFunction(_f_ip; mass_matrix=mm)
    fun2_ip = ODEFunction(_f_ip; mass_matrix=mm,
    jac_prototype=DiffEqArrayOperator(similar(A); update_func=(J,u,p,t) -> (J .= t .* A; J)))

    for Alg in [ImplicitEuler, Rosenbrock23]
      println(Alg)
      sol1 = solve(ODEProblem(fun1,u0,tspan), Alg(); adaptive=false, dt=0.01)
      sol2 = solve(ODEProblem(fun2,u0,tspan), Alg(); adaptive=false, dt=0.01)
      @test sol1(1.0) ≈ sol2(1.0)

      sol1_ip = solve(ODEProblem(fun1_ip,u0,tspan), Alg(); adaptive=false, dt=0.01)
      sol2_ip = solve(ODEProblem(fun2_ip,u0,tspan), Alg(linsolve=LinSolveFactorize(lu)); adaptive=false, dt=0.01)
      @test sol1_ip(1.0) ≈ sol2_ip(1.0)
    end
  end
end


text:
using ReachabilityAnalysis
Dict([
    "opts" => Dict(:T=>20.0, :alg=>GLGM06(δ=0.01)),
    "X0"   => BallInf(zeros(270), 1e-4),
    "info" => Dict(
        "name" => "iss",
        "dim" => 270,
        "linear" => true,
        "hybrid" => false
    )
])

max_stars_repo_name: mforets/ReachabilityLibrary

max_issues_repo_name: mforets/ReachabilityLibrary

max_forks_repo_name: mforets/ReachabilityLibrary


text:
function WelzlEncloser(arg0::jdouble, arg1::SupportBallGenerator)
    return WelzlEncloser((jdouble, SupportBallGenerator), arg0, arg1)
end

function enclose(obj::WelzlEncloser, arg0::Iterable)
    return jcall(obj, "enclose", EnclosingBall, (Iterable,), arg0)
end

function select_farthest(obj::WelzlEncloser, arg0::Iterable, arg1::EnclosingBall)
    return jcall(obj, "selectFarthest", Point, (Iterable, EnclosingBall), arg0, arg1)
end


text:
function asString(m::Dates.Minute)
	minutes = m.value % 60
	hours = m.value ÷ 60
	@sprintf("%02d:%02d", hours, minutes)
end

pm(fun::Function, x) = !ismissing(x) ? fun(x) : missing

function ellipsis(str::String, limit::Int)
	if limit <= 0 || limit >= length(str)
		return str
	end

	str[1:nextind(str, limit)] * "..."
end

function applyNormalizations(str::AbstractString, norm::Vector{Tuple{String, String}})
	for (oldval, newval) in norm
		if occursin(str, oldval)
			return newval
		end
	end
	return str
end

max_stars_repo_name: julianschick/flugbuch

max_issues_repo_name: julianschick/flugbuch

max_forks_repo_name: julianschick/flugbuch


text:
module ODBC

using DataFrames
using Datetime

export advancedconnect,
	   query, querymeta, @query, @sql_str,
	   Connection, Metadata, conn, Connections,
	   disconnect, listdrivers, listdsns

import Base: show, string

include("ODBC_Types.jl")
include("ODBC_API.jl")

#Metadata type holds metadata related to an executed query resultset
type Metadata
	querystring::String
	cols::Int
	rows::Int
	colnames::Array{UTF8String}
	coltypes::Array{(String,Int16)}
	colsizes::Array{Int}
	coldigits::Array{Int16}
	colnulls::Array{Int16}
end
function show(io::IO,meta::Metadata)
	if meta == null_meta
		print(io,"No metadata")
	else
		println(io,"Resultset metadata for executed query")
		println(io,"------------------------------------")
		print(io,"Query: $(meta.querystring)")
		println(io,"Columns: $(meta.cols)")
		println(io,"Rows: $(meta.rows)")
		println(io,DataFrame([meta.colnames meta.coltypes meta.colsizes meta.coldigits meta.colnulls], ["Column Names","Types","Sizes","Digits","Nullable"]))
	end
end
#Connection object that holds information related to each established connection and retrieved resultsets
type Connection
	dsn::String
	number::Int
	dbc_ptr::Ptr{Void}
	stmt_ptr::Ptr{Void}
	#Holding a reference to the last resultset is useful if the user runs several test queries just using `query()` or `sql"..."`
	#then realizes the last resultset should actually be saved to a variable (happended to me all the time in RODBC)
	resultset::Any
end
function show(io::IO,conn::Connection)
	if conn == null_connection
		print("Null ODBC Connection Object")
	else
		println("ODBC Connection Object")
		println("----------------------")
		println("Connection Data Source: $(conn.dsn)")
		println("$(conn.dsn) Connection Number: $(conn.number)")
		if isequal(conn.resultset,null_resultset)
			print("Contains resultset? No")
		else
			print("Contains resultset(s)? Yes (access by referencing the resultset field (e.g. conn.resultset))")
		end
	end
end
#There was a weird bug where Connections was showing each Connection 3 times, this seems to solve it
show(io::IO,conns::Array{Connection,1}) = map(show,conns)

typealias Output Union(DataType,String)
#Global module consts and variables
const null_resultset = DataFrame(0)
const null_connection = Connection("",0,C_NULL,C_NULL,null_resultset)
const null_meta = Metadata("",0,0,UTF8String[],Array((String,Int16),0),Int[],Int16[],Int16[])
env = C_NULL
Connections = Connection[] #For managing references to multiple connections
conn = null_connection #Create default connection = null
ret = ""

include("backend.jl")
include("userfacing.jl")

end #ODBC module


text:
module MultiTruncatedDistributions

using Distributions, IntervalSets
export TruncatedDistribution, cdf

prob(d::UnivariateDistribution, range::AbstractInterval) = cdf(d, range.left) - cdf(d, range.right)

struct TruncatedDistribution{S,T}
    d::S # distribution
    Z::T # truncated region
end

function Distributions.cdf(td::TruncatedDistribution, x)
    d, Z = td.d, td.Z

    length(Z) > 1 && begin
        function modify(Z)
            Z_new = eltype(Z)[Z[begin]]
            for I in Z[2:end]
                if abs(Z_new[end].right - I.left) < 1e-2 * std(d)
                    Z_new[end] = (Z_new[end].left)..(I.right)
                else
                    push!(Z_new, I)
                end
            end
            return Z_new
        end
        Z = modify(Z)
    end

    numerator = zero(prob(d, x..x)); denominator = numerator
    for I in Z
        denominator += prob(d, I)
        if I.right ≤ x
            numerator += prob(d, I)
        elseif x ∈ I
            numerator += prob(d, (I.left)..x)
        end
    end
    return numerator / denominator
end

end # module MultiTruncatedDistributions


text:
module PhysicalConstants

import Base.sqrt

export μB, ħ, m_ca40, c, e, α, ϵ₀, kB, ca40_qubit_transition_frequency

"""
    PhysicalConstant(x::Real)
Useful physical constants. Values are in SI units.
"""
struct PhysicalConstant <: Real
    x::Real
    units::String
end

# some useful constants, everything in SI units
"""`μB` = 9.27400994e-24 J⋅T⁻¹
 (Bohr Magneton)"""
const μB = PhysicalConstant(9.27400994e-24, "J⋅T⁻¹")
"""`ħ` = 1.0545718e-34 m²kg/s
 (Planck's constant / 2π)"""
const ħ = PhysicalConstant(1.0545718e-34, "m²kg/s")
"""`m_ca40` = 6.635943757345042e-26 kg
 (mass of 40Ca)"""
const m_ca40 = PhysicalConstant(6.635943757345042e-26, "kg")
"""`c` = 2.99792458e8 m/s
 (speed of light in vacuum)"""
const c = PhysicalConstant(2.99792458e8, "m/s")
"""`e` = 1.60217662e-19 C
 (charge of electron)"""
const e = PhysicalConstant(1.60217662e-19, "C")
"""`ϵ₀` = 8.85418782e-12 ``(s^4A^2) / (m^3 kg)``"""
const ϵ₀ = PhysicalConstant(8.85418782e-12, "(s^4A^2) / (m^3 kg)")
"""`α` = e²/4πϵ₀ħc``"""
const α = PhysicalConstant(0.007297352557920479, "")
"""`kB` = 1.38064852e-23 ``m^2kg/(s^2K)``"""
const kB = PhysicalConstant(1.38064852e-23, "m^2kg/(s^2K)")
"""`ca40_qubit_transition_frequency` = c / 729.147e-9 ``Hz`` """
const ca40_qubit_transition_frequency = PhysicalConstant(2.99792458e8 / 729.147e-9, "Hz")

Base.print(pc::PhysicalConstant) = print("$(pc.x) [$(pc.units)]")
Base.show(io::IO, pc::PhysicalConstant) = print(io, "$(pc.x) [$(pc.units)]")

Base.convert(::Type{<:Number}, x::PhysicalConstant) = map(x->x.x, x)
Base.promote_rule(::Number, ::PhysicalConstant) = Number
Base.Fix2(f, x::PhysicalConstant) = Base.Fix2(f, x.x)

Base.:*(x1::PhysicalConstant, x2::PhysicalConstant) = x1.x * x2.x
Base.:/(x1::PhysicalConstant, x2::PhysicalConstant) = x1.x / x2.x
Base.:+(x1::PhysicalConstant, x2::PhysicalConstant) = x1.x + x2.x
Base.:-(x1::PhysicalConstant, x2::PhysicalConstant) = x1.x - x2.x
Base.:^(x1::PhysicalConstant, x2::PhysicalConstant) = x1.x^x2.x
Base.:(>)(x1::PhysicalConstant, x2::PhysicalConstant) = Base.:(>)(x1.x, x2.x)
Base.:(<)(x1::PhysicalConstant, x2::PhysicalConstant) = Base.:(<)(x1.x, x2.x)
Base.:(≥)(x1::PhysicalConstant, x2::PhysicalConstant) = Base.:(≥)(x1.x, x2.x)
Base.:(≤)(x1::PhysicalConstant, x2::PhysicalConstant) = Base.:(≤)(x1.x, x2.x)

sqrt(x1::PhysicalConstant) = sqrt(x1.x)

end  # module

export x̂, ŷ, ẑ, ndot

const x̂ = (x=1, y=0, z=0)
const ŷ = (x=0, y=1, z=0)
const ẑ = (x=0, y=0, z=1)

function _print_axis(a::NamedTuple{(:x,:y,:z)})
    if a == x̂
        return "x̂"
    elseif a == ŷ
        return "ŷ"
    elseif a == ẑ
        return "ẑ"
    else
        return string(a)
    end
end

ndot(a::NamedTuple{(:x,:y,:z)}, b::NamedTuple{(:x,:y,:z)}) = a.x * b.x + a.y * b.y + a.z * b.z
function Base.:+(a::NamedTuple{(:x,:y,:z)}, b::NamedTuple{(:x,:y,:z)})
    (x=a.x + b.x, y=a.y + b.y, z=a.z + b.z)
end
function Base.:-(a::NamedTuple{(:x,:y,:z)}, b::NamedTuple{(:x,:y,:z)})
    (x=a.x - b.x, y=a.y - b.y, z=a.z - b.z)
end
Base.:/(a::NamedTuple{(:x,:y,:z)}, b::Number) = (x=a.x/b, y=a.y/b, z=a.z/b)
Base.:*(a::NamedTuple{(:x,:y,:z)}, b::Number) = (x=a.x*b, y=a.y*b, z=a.z*b)
Base.:*(b::Number, a::NamedTuple{(:x,:y,:z)}) = (x=a.x*b, y=a.y*b, z=a.z*b)


text:
macro visualtest(testfun, refname, popup=true, tol=0.02)
  esc(quote
    testFilename = tempname()*".png"
    $testfun(testFilename)
    @test test_images(testFilename, $refname, popup=$popup, tol=$tol) |> success
  end)
end

macro plottest(plotfun, refname, popup=true, tol=0.02)
  esc(quote
    testFilename = tempname()*".png"
    if $plotfun isa Function
      $plotfun()
    else # function body
      $plotfun
    end
    png(testFilename)
    @test test_images(testFilename, $refname, popup=$popup, tol=$tol) |> success
  end)
end


text:
module Mesoscale
using Pkg
export add_branch

function add_branch()
    Pkg.add(url = "https://github.com/CliMA/Oceananigans.jl.git", rev = "glw-ncc/overturning-channel-example")
end

end # end of module

max_stars_repo_name: sandreza/Mesoscale

max_issues_repo_name: sandreza/Mesoscale

max_forks_repo_name: sandreza/Mesoscale


text:
"Check if term is a predicate."
is_pred(term::Term, domain::Domain) =
    term.name in keys(get_predicates(domain))

"Check if term is a non-boolean fluent (i.e. function)."
is_func(term::Term, domain::Domain) =
    term.name in keys(get_functions(domain))

"Check if term is an external function attached to a domain."
is_attached_func(term::Term, domain::Domain) =
    term.name in keys(get_funcdefs(domain))

"Check if term is an external function (attached or global)."
is_external_func(term::Term, domain::Domain) =
    is_global_func(term) || is_attached_func(term, domain)

"Check if term is a derived predicate"
is_derived(term::Term, domain::Domain) =
    term.name in keys(get_axioms(domain))

"Check if term is a domain fluent."
is_fluent(term::Term, domain::Domain) =
    is_pred(term, domain) || is_func(term, domain)

"Check if term is a universal or existential quantifier."
is_quantifier(term::Term) =
    term.name == :forall || term.name == :exists

"Check if term is a global predicate (comparison, equality, etc.)."
is_global_pred(term::Term) =
    term.name in keys(GLOBAL_PREDICATES)

"Check if term is a global function, including global predicates."
is_global_func(term::Term) =
    term.name in keys(GLOBAL_FUNCTIONS)

"Check if term is a type predicate."
is_type(term::Term, domain::Domain) =
    term.name in get_types(domain)

"Check if term is a type predicate with subtypes."
has_subtypes(term::Term, domain::Domain) =
    !isempty(get_subtypes(domain, term.name))

"Check if `term` has a (sub-term with a) name in `names`."
has_name(term::Const, names) = term.name in names
has_name(term::Var, names) = false
has_name(term::Compound, names) =
    term.name in names || any(has_name(f, names) for f in term.args)

"Check if term contains a predicate name."
has_pred(term::Term, domain::Domain) =
    has_name(term, keys(get_predicates(domain)))

"Check if term contains the name of numeric fluent (i.e. function)."
has_func(term::Term, domain::Domain) =
    has_name(term, keys(get_functions(domain)))

"Check if term contains a derived predicate"
has_derived(term::Term, domain::Domain) =
    has_name(term, keys(get_axioms(domain)))

"Check if term contains a universal or existential quantifier."
has_quantifier(term::Term) =
    has_name(term, Set(Symbol[:forall, :exists]))

"Check if term contains a fluent name."
has_fluent(term::Term, domain::Domain) =
    has_pred(term, domain) || has_func(term, domain)

"Returns list of constituent fluents."
constituents(term::Const, domain::Domain) =
    is_fluent(term, domain) ? [term] : []
constituents(term::Var, domain::Domain) =
    false
constituents(term::Compound, domain) = is_fluent(term, domain) ?
    Term[term] : reduce(vcat, (constituents(a, domain) for a in term.args))