text:
# Tests for Directional Multiplier DEA Models
@testset "DirectionalMultiplierDEAModel" begin

    ## Test Directional DF DEA Models with FLS Book data
    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]

    # Observed CRS
    deaddfmobs = deaddfm(X, Y, Gx = X, Gy = Y, rts = :CRS)

    @test typeof(deaddfmobs) == DirectionalMultiplierDEAModel

    @test nobs(deaddfmobs) == 11
    @test ninputs(deaddfmobs) == 2
    @test noutputs(deaddfmobs) == 1
    @test efficiency(deaddfmobs) ≈ [
        0.00000000000;
        0.23282544774;
        0.09898790422;
        0.00000000000;
        0.52628538417;
        0.28571428571;
        0.00000000000;
        0.13787061050;
        0.09883720930;
        0.34177215190;
        0.00000000000] atol = 1e-5
    @test multipliers(deaddfmobs, :X) ≈ [
        0.0450913    0.0211187
        0.0193798    0.0255279
        0.0195014    0.00913355
        0.0136023    0.0179174
        0.0209417    0.0275852
        0.0          0.107143
        0.0130985    0.0172538
        0.0101633    0.0133875
        0.0          0.0392442
        0.00895338   0.0117938
        0.1         -1.11022e-16
    ] atol = 1e-5
    @test multipliers(deaddfmobs, :Y) ≈ [
        0.04166666666666667
        0.027399091152098366
        0.0180202419155764
        0.01923076923076923
        0.02960716348931253
        0.0396825396825397
        0.018518518518518517
        0.014368823158337424
        0.014534883720930238
        0.012658227848101266
        0.04166666666666667
    ] atol = 1e-5

    @test efficiency(deaddfm(X, Y, Gx = :Observed, Gy = :Observed, rts = :CRS)) == efficiency(deaddfmobs)

    @test efficiency(deaddfm(targets(deaddfmobs, :X), targets(deaddfmobs, :Y), Gx = :Observed, Gy = :Observed, rts = :CRS)) ≈ zeros(11) atol = 1e-5
    @test rts(deaddfmobs) ≈ zeros(11)

    # Observed VRS
    deaddfmobsvrs = deaddfm(X, Y, Gx = X, Gy = Y, rts = :VRS)

    @test nobs(deaddfmobsvrs) == 11
    @test ninputs(deaddfmobsvrs) == 2
    @test noutputs(deaddfmobsvrs) == 1
    @test efficiency(deaddfmobsvrs) ≈ [
        0.00000000000;
        0.1076130509;
        0.0000000000;
        0.0000000000;
        0.2883597884;
        0.0000000000;
        0.00000000000;
        0.0000000000;
        0.0000000000;
        0.1821192053;
        0.00000000000] atol = 1e-5
    @test multipliers(deaddfmobsvrs, :X) ≈ [
        0.0450913    0.0211187
        0.0163137    0.0472238
        0.0264496    0.000339098
        0.0136023    0.0179174
        0.0185185    0.047619
        0.0          0.125
        0.0130985    0.0172538
        0.00627716   0.00482859
        0.0          0.0222222
        0.00331126   0.0
        0.1         -1.11022e-16
    ] atol = 1e-5
    @test multipliers(deaddfmobsvrs, :Y) ≈ [
        0.04166666666666667
        0.012306811677160837
        0.02271956595456087
        0.01923076923076923
        0.0
        0.027777777777777773
        0.018518518518518517
        0.024142926122646062
        0.022222222222222206
        0.033112582781456984
        0.04166666666666667
    ] atol = 1e-5

    @test efficiency(deaddfm(X, Y, Gx = :Observed, Gy = :Observed, rts = :VRS)) == efficiency(deaddfmobsvrs)

    @test efficiency(deaddfm(targets(deaddfmobsvrs, :X), targets(deaddfmobsvrs, :Y), Gx = :Observed, Gy = :Observed, rts = :VRS)) ≈ zeros(11) atol = 1e-5
    @test rts(deaddfmobsvrs) ≈ [
        0.0
        -0.5477962220950202
         0.1359782977280435
         0.0
        -0.7116402116402119
        -0.5000000000000003
         0.0
         0.448575567358764
         0.37777777777777694
         0.9039735099337763
         0.0
    ] atol = 1e-5

    # Ones CRS
    deaddfmones = deaddfm(X, Y, Gx = ones(size(X)), Gy = ones(size(Y)), rts = :CRS)

    @test nobs(deaddfmones) == 11
    @test ninputs(deaddfmones) == 2
    @test noutputs(deaddfmones) == 1
    @test efficiency(deaddfmones) ≈ [
        0.00000000000;
        3.219963031;
        2.121693122;
        0.00000000000;
        6.735674677;
        1.945945946;
        0.00000000000;
        3.635859519;
        1.837837838;
        10.231053604;
        0.00000000000]
    @test multipliers(deaddfmones, :X) ≈ [
        0.705882  0.0
        0.268022  0.35305
        0.417989  0.195767
        0.417989  0.195767
        0.268022  0.35305
        0.0       0.72973
        0.268022  0.35305
        0.268022  0.35305
        0.0       0.72973
        0.268022  0.35305
        0.705882  0.0
    ] atol = 1e-5
    @test multipliers(deaddfmones, :Y) ≈ [
        0.29411764705882354
        0.3789279112754159
        0.3862433862433863
        0.3862433862433863
        0.3789279112754159
        0.27027027027027034
        0.3789279112754159
        0.3789279112754159
        0.27027027027027034
        0.3789279112754159
        0.29411764705882354
    ] atol = 1e-5

    @test efficiency(deaddfm(X, Y, Gx = :Ones, Gy = :Ones, rts = :CRS)) == efficiency(deaddfmones)
    @test rts(deaddfmobs) ≈ zeros(11)

    # Ones VRS
    deaddfmonesvrs = deaddfm(X, Y, Gx = ones(size(X)), Gy = ones(size(Y)), rts = :VRS)

    @test nobs(deaddfmonesvrs) == 11
    @test ninputs(deaddfmonesvrs) == 2
    @test noutputs(deaddfmonesvrs) == 1
    @test efficiency(deaddfmonesvrs) ≈ [
        0.000000000;
        1.418867925;
        0.000000000;
        0.000000000;
        4.067924528;
        0.000000000;
        0.000000000;
        0.000000000;
        0.000000000;
        5.000000000;
        0.000000000] atol = 1e-5
    @test multipliers(deaddfmonesvrs, :X) ≈ [
        0.705882   0.0
        0.215094   0.622642
        0.534247   0.00684932
        0.417989   0.195767
        0.215094   0.622642
        0.0        0.818182
        0.268022   0.35305
        0.178082   0.136986
        0.0        0.5
        0.0909091  0.0
        0.705882   0.0
    ] atol = 1e-5
    @test multipliers(deaddfmonesvrs, :Y) ≈ [
        0.29411764705882354
        0.16226415094339627
        0.4589041095890411
        0.3862433862433863
        0.16226415094339627
        0.18181818181818193
        0.3789279112754159
        0.6849315068493145
        0.4999999999999999
        0.9090909090909107
        0.29411764705882354
    ] atol = 1e-5

    @test efficiency(deaddfm(X, Y, Gx = :Ones, Gy = :Ones, rts = :VRS)) == efficiency(deaddfmonesvrs)
    @test rts(deaddfmonesvrs) ≈ [
        0.0
        -7.222641509433959
         2.746575342465754
         0.0
        -7.222641509433959
        -3.2727272727272703
         0.0
        12.72602739726025
         8.499999999999993
        24.81818181818191
         0.0
    ] atol = 1e-5

    # Only X CRS
    deaddfmonlyX = deaddfm(X, Y, Gx =  X, Gy = zeros(size(Y)), rts = :CRS)

    @test nobs(deaddfmonlyX) == 11
    @test ninputs(deaddfmonlyX) == 2
    @test noutputs(deaddfmonlyX) == 1
    @test efficiency(deaddfmonlyX) ≈ [
        0.0000000000;
        0.3777103209;
        0.1801437556;
        0.0000000000;
        0.6896290689;
        0.4444444444;
        0.0000000000;
        0.2423309104;
        0.1798941799;
        0.5094339623;
        0.0000000000] atol = 1e-5
    @test multipliers(deaddfmonlyX, :X) ≈ [
        0.0901826  0.0422374
        0.0314397  0.0414137
        0.0354897  0.0166217
        0.0416228  0.0194942
        0.0274413  0.0361469
        0.0        0.166667
        0.0261969  0.0345077
        0.0178637  0.0235309
        0.0        0.0714286
        0.0133456  0.0175794
        0.2        2.61229e-17
    ] atol = 1e-5
    @test multipliers(deaddfmonlyX, :Y) ≈ [
        0.08333333333333336
        0.04444926279271466
        0.032794249775381853
        0.038461538461538464
        0.03879636638909917
        0.06172839506172836
        0.03703703703703704
        0.02525563631883701
        0.026455026455026464
        0.018867924528301886
        0.0833333333333333
    ] atol = 1e-5

    @test efficiency(deaddfm(X, Y, Gx = :Observed, Gy = :Zeros, rts = :CRS)) == efficiency(deaddfmonlyX)

    # Only X VRS
    deaddfmonlyXvrs = deaddfm(X, Y, Gx = X, Gy = zeros(size(Y)), rts = :VRS)

    @test nobs(deaddfmonlyXvrs) == 11
    @test ninputs(deaddfmonlyXvrs) == 2
    @test noutputs(deaddfmonlyXvrs) == 1
    @test efficiency(deaddfmonlyXvrs) ≈ [
        0.0000000000;
        0.1300138313;
        0.0000000000;
        0.0000000000;
        0.2883597884;
        0.0000000000;
        0.0000000000;
        0.0000000000;
        0.0000000000;
        0.5068790731;
        0.0000000000] atol = 1e-5
    @test multipliers(deaddfmonlyXvrs, :X) ≈ [
        0.0901826  0.0422374
        0.0197095  0.0570539
        0.0612245  0.000784929
        0.0416228  0.0194942
        0.0185185  0.047619
        0.0        0.166667
        0.0261969  0.0345077
        0.0227671  0.0175131
        0.0        0.0714286
        0.013034   0.0181028
        0.2        2.61229e-17
    ] atol = 1e-5
    @test multipliers(deaddfmonlyXvrs, :Y) ≈ [
        0.08333333333333336
        0.014868603042876911
        0.05259026687598115
        0.038461538461538464
        0.0
        0.03703703703703698
        0.03703703703703704
        0.08756567425569173
        0.07142857142857138
        0.013758146270818254
        0.0833333333333333
    ] atol = 1e-5

    @test efficiency(deaddfm(X, Y, Gx = :Observed, Gy = :Zeros, rts = :VRS)) == efficiency(deaddfmonlyXvrs)

    # Only Y CRS
    deaddfmonlyY = deaddfm(X, Y, Gx = zeros(size(X)), Gy = Y, rts = :CRS)

    @test nobs(deaddfmonlyY) == 11
    @test ninputs(deaddfmonlyY) == 2
    @test noutputs(deaddfmonlyY) == 1
    @test efficiency(deaddfmonlyY) ≈ [
        0.0000000000;
        0.6069686411;
        0.2197260274;
        0.0000000000;
        2.2219512195;
        0.8000000000;
        0.0000000000;
        0.3198373984;
        0.2193548387;
        1.0384615385;
        0.0000000000] atol = 1e-5
    @test multipliers(deaddfmonlyY, :X) ≈ [
        0.0901826  0.0422374
        0.0505226  0.0665505
        0.0432877  0.020274
        0.0272045  0.0358349
        0.0884146  0.116463
        0.0        0.3
        0.0261969  0.0345077
        0.0235772  0.0310569
        0.0        0.0870968
        0.0272045  0.0358349
        0.2        0.0
    ] atol = 1e-5
    @test multipliers(deaddfmonlyY, :Y) ≈ [
        0.08333333333333329
        0.07142857142857142
        0.04000000000000001
        0.03846153846153846
        0.12499999999999999
        0.11111111111111112
        0.037037037037037035
        0.03333333333333334
        0.03225806451612903
        0.03846153846153846
        0.08333333333333318
    ] atol = 1e-5

    @test efficiency(deaddfm(X, Y, Gx = :Zeros, Gy = :Observed, rts = :CRS)) == efficiency(deaddfmonlyY)

    # Only Y VRS
    deaddfmonlyYvrs = deaddfm(X, Y, Gx = zeros(size(X)), Gy = Y, rts = :VRS)

    @test nobs(deaddfmonlyYvrs) == 11
    @test ninputs(deaddfmonlyYvrs) == 2
    @test noutputs(deaddfmonlyYvrs) == 1
    @test efficiency(deaddfmonlyYvrs) ≈ [
        0.0000000000;
        0.5075187970;
        0.0000000000;
        0.0000000000;
        2.2039473684;
        0.0000000000;
        0.0000000000;
        0.0000000000;
        0.0000000000;
        0.1923076923;
        0.0000000000] atol = 1e-5
    @test multipliers(deaddfmonlyYvrs, :X) ≈ [
        0.0901826   0.0422374
        0.0676692   0.093985
        0.0465672   0.000597015
        0.0272045   0.0358349
        0.118421    0.164474
        0.0         0.5
        0.0261969   0.0345077
        0.00866667  0.00666667
        0.0         0.0322581
        0.0         0.0
        0.2         0.0
    ] atol = 1e-5
    @test multipliers(deaddfmonlyYvrs, :Y) ≈ [
        0.08333333333333329
        0.07142857142857144
        0.040000000000000015
        0.03846153846153846
        0.125
        0.11111111111111109
        0.037037037037037035
        0.03333333333333333
        0.03225806451612903
        0.038461538461538464
        0.08333333333333318
    ] atol = 1e-5

    @test efficiency(deaddfm(X, Y, Gx = :Zeros, Gy = :Observed, rts = :VRS)) == efficiency(deaddfmonlyYvrs)

    # Mean CRS
    deaddfmmean = deaddfm(X, Y, Gx = :Mean, Gy = :Mean, rts = :CRS)

    @test efficiency(deaddfmmean) ≈ [0; 0.1713509018; 0.1082533683; 0; 0.3584401184; 0.1148158887; 0; 0.1934829069; 0.1084372282; 0.5444473258; 0] atol = 1e-5

    # Mean VRS
    deaddfmmeanvrs = deaddfm(X, Y, Gx = :Mean, Gy = :Mean, rts = :VRS)

    @test efficiency(deaddfmmeanvrs) ≈ [0; 0.08056567388; 0; 0;0.23098350118; 0; 0; 0; 0; 0.24886877828; 0] atol = 1e-5

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

    deamddfvrs_ref_eff = zeros(size(X, 1))

    deamddfvrs_ref_multX = zeros(size(X))
    deamddfvrs_ref_multY = zeros(size(Y))

    Gx = X[:,:]
    Gy = Y[:,:]

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

    for i = 1:size(X, 1)
        Xeval = X[i:i,:]
        Xeval = Xeval[:,:]
        Yeval = Y[i:i,:]
        Yeval = Yeval[:,:]
        Gxeval = Gx[i:i,:]
        Gxeval = Gxeval[:,:]
        Gyeval = Gy[i:i,:]
        Gyeval = Gyeval[:,:]

        deamddf_ref_eff[i] = efficiency(deaddfm(Xeval, Yeval, Gx = Gxeval, Gy = Gyeval, rts = :CRS, Xref = Xref, Yref = Yref))[1]

        deamddfvrs_ref_eff[i] = efficiency(deaddfm(Xeval, Yeval, Gx = Gxeval, Gy = Gyeval, rts = :VRS, Xref = Xref, Yref = Yref))[1]

        deamddfvrs_ref_multX[i,:] = multipliers(deaddfm(Xeval, Yeval, Gx = Gxeval, Gy = Gyeval, rts = :VRS, Xref = Xref, Yref = Yref), :X)
        deamddfvrs_ref_multY[i,:] = multipliers(deaddfm(Xeval, Yeval, Gx = Gxeval, Gy = Gyeval, rts = :VRS, Xref = Xref, Yref = Yref), :Y)
    end

    @test deamddf_ref_eff ≈ efficiency(deaddfmobs)

    @test deamddfvrs_ref_eff ≈ efficiency(deaddfmobsvrs)

    @test deamddfvrs_ref_multX ≈ multipliers(deaddfmobsvrs, :X) atol = 1e-5
    @test deamddfvrs_ref_multY ≈ multipliers(deaddfmobsvrs, :Y) atol = 1e-5

    # Print
    show(IOBuffer(), deaddfmobs)

    # Test errors
    @test_throws DimensionMismatch deaddfm([1; 2 ; 3], [4 ; 5], Gx = [1; 2 ; 3], Gy = [4 ; 5]) #  Different number of observations
    @test_throws DimensionMismatch deaddfm([1; 2], [4 ; 5], Gx = [1; 2], Gy = [4 ; 5], Xref = [1; 2; 3; 4]) # Different number of observations in reference sets
    @test_throws DimensionMismatch deaddfm([1 1; 2 2], [4 4; 5 5], Gx = [1 1; 2 2], Gy = [4 4; 5 5], Xref = [1 1 1; 2 2 2]) # Different number of inputs
    @test_throws DimensionMismatch deaddfm([1 1; 2 2], [4 4; 5 5], Gx = [1 1; 2 2], Gy = [4 4; 5 5], Yref = [4 4 4; 5 5 5]) # Different number of inputs
    @test_throws ArgumentError deaddfm([1; 2; 3], [4; 5; 6], Gx = [1; 2; 3], Gy = [4; 5; 6], rts = :Error) # Invalid returns to scale
    @test_throws DimensionMismatch deaddfm([1 1; 2 2; 3 3], [4; 5; 6], Gx = [1 1 1; 2 2 2; 3 3 3], Gy = [4; 5; 6]) # Different size of inputs direction
    @test_throws DimensionMismatch deaddfm([1; 2; 3], [4 4; 5 5; 6 6], Gx = [1; 2; 3], Gy = [4 4 4; 5 5 5; 6 6 6]) # Different size of inputs direction
    @test_throws ArgumentError deaddfm([1; 2; 3], [1; 2; 3], Gx = :Error, Gy = :Ones) # Invalid inuts direction
    @test_throws ArgumentError deaddfm([1; 2; 3], [1; 2; 3], Gx = :Ones, Gy = :Error) # Invalid outputs direction

    # ------------------
    # 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(deaddfm(X, Y, Gx = X, Gy = Y)) ≈ [0; 0; 0; 0.25; 0.4285714286; 0; 0.2; 0.2307692308] atol = 1e-5

    # 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(deaddfm(X, Y, Gx = X, Gy = Y)) ≈ [0; 0; 0; 0.2173913043; 0.4; 0; 0.12; 0.2307692308] atol = 1e-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(deaddfm(X, Y, Gx = X, Gy = Y)) ≈ [0.4285714286; 0; 0.1111111111; 0.25; 0.4285714286; 0.4285714286; 0.3207547170; 0.6666666667] atol = 1e-5

    # ------------------
    # Test model with zero Directions
    # ------------------
    deaddfmzeros = deaddfm(X, Y, Gx = zeros(size(X)), Gy = zeros(size(Y)), rts = :VRS)
    @test efficiency(deaddfmzeros) == zeros(8)

end


text:
# Run forward model demonstrating hand-adjusting some stuff
println("Time to load BITEModel.jl")
@time using BITEModel
const BM=BITEModel
if !isdefined(:plotyes)
    plotyes = true
end

# run forward model for ITMIX glaciers
glaciers = [:Academy, :Aqqutikitsoq, :Austfonna, :Brewster, :Columbia,
            :Devon, :Elbrus, :Freya, :Hellstugubreen, :Kesselwandferner, :Mocho,
            :NorthGlacier, :SouthGlacier, :Starbuck, :Synthetic1, :Synthetic2,
            :Synthetic3, :Tasman, :Unteraar, :Urumqi, :Washmawapta]
# glacier = glaciers[gln]
glacier = :Unteraar
use_glogem = false # use GloGEM bdot, if available, and use the SyntheticBenchLoader for dhdt
update_cache = false

##################
gid = BM.ITMIXGlacier(glacier, 2)

println("Time to to make load-paras")
@time pl = BM.LoadPara(gid, update_cache=update_cache, use_glogem=use_glogem)

# Now customize pl.  This essentially follows what is encoded in the
# functions init_forward and the custom LoadPara method.

# E.g. make sigma larger
pl.dataset_opts[BM.ThicknessData][:sigma]=40.0

if glacier==:Elbrus
    # why not some more down-sampling?
    pl.dataset_opts[BM.DataKind][:grid_downsample_step]=3
end

pl.dataset_opts[BM.ParaData][:pm][:bandsize] = 40 # bigger band

dt, pl = BM.make_datatable(gid,pl)
# swap out the bdot against one from the synthetic bench run
bdot_synth = BM.DataSet{BM.BdotData, BM.SyntheticBenchLoader}([""],
                                                              pl.dataset_opts[BM.BdotData],
                                                              (Date(),Date())
                                                              )
dt = BM.DataTable(dt, bdot=bdot_synth)

gl,pp,pm,pn,pl = load_glacier(dt,pl)
println("Time for loading glacier:")
@time gl,pp,pm,pn,pl = load_glacier(dt,pl)

Base.Test.@inferred BM.make_bands(gl, pp, pm, pn);
println("Time for making bands:")
@time gb,pm = BM.make_bands(gl, pp, pm, pn);

Base.Test.@inferred BM.fwdm(gb, pp, pm, pn, gb.bdot-gb.dhdt, gb.fsl, gb.temp)

println("Time to run forward model:")
gc()
@time fwdsol = BM.fwdm(gb, pp, pm, pn, gb.bdot-gb.dhdt, gb.fsl, gb.temp)

qtot1d, qd1d = postproc1d(gb, pp, pm)

if plotyes
    println("Time to plot:")
    @time display(BM.plot_run(fwdsol, reuse=false, plot2d=[:h,:iv,:ivlog][1]))
end

ext: jl


text:
"""
# module Earley

- Julia version:
- Author: bramb
- Date: 2019-02-13

based on https://github.com/tomerfiliba/tau/blob/master/earley3.py

# Examples

```jldoctest
julia>using HyperGraphTools.Earley
```
"""
module Earley

DEBUG = false

import Base.push!, Base.isequal, Base.hash, Base.iterate, Base.length, Base.getindex

export earley_parse,
    str,
    Grammar,
    Rule,
    Production

mutable struct Production
    terms::Vector{}

    Production(term::String) = new([term])

    Production(terms::Vector{}) = new(terms)
end

length(production::Production) = length(production.terms)

getindex(production::Production, i::Int64) = production.terms[i]

isequal(p1::Production, p2::Production) = p1.terms == p2.terms

mutable struct Rule
    name::String
    productions::Vector{Production}

    Rule(name::String) = new(name,[])

    Rule(name::String, production::Production) = new(name,[production])

    Rule(name::String, productions::Vector{Production}) = new(name,productions)
end

push!(rule::Rule, production::Production) = push!(rule.productions, production)

mutable struct Grammar
    rules_index::Dict{String,Rule}

    Grammar() = new(Dict{String,Rule}())
end

push!(grammar::Grammar, rule::Rule) = grammar.rules_index[rule.name] = rule

abstract type AState end

mutable struct Column
    index::Int64
    token::String
    states::Vector{AState}
    _unique::Set{AState}

    Column(index::Int64, token::String) = new(index,token,[],Set())
end

iterate(column::Column) = iterate(column.states)

iterate(column::Column, i::Int64) = iterate(column.states,i)

function Base.push!(column::Column, state::AState)
    if !(state in column._unique)
        push!(column._unique,state)
        state.end_column = column
        push!(column.states, state)
        return true
    else
        return false
    end
end

mutable struct State <: AState
    name::String
    production::Production
    dot_index::Int64
    start_column::Column
    end_column::Union{Column,Nothing}

    State(name::String, production::Production, dot_index::Int64, start_column::Column) =
      new(name,production,dot_index,start_column,nothing)
end

function Base.isequal(state1::State, state2::State)
    return (state1.name, state1.production, state1.dot_index, state1.start_column) == (state2.name, state2.production, state2.dot_index, state2.start_column)
end

hash(state::State) = hash(state.name) + hash(state.production) + hash(state.dot_index) + hash(state.start_column)

iscompleted(state::State) = state.dot_index - 1 >= length(state.production)

next_term(state::State) = (iscompleted(state)) ? nothing : state.production[state.dot_index]

GAMMA_RULE = "γ"
function earley_parse(tokens, grammar)
    table = [Column(i,token) for (i,token) in enumerate(vcat(["∅"],tokens))]
    start_rule = grammar.rules_index["S"]
    push!(table[1],State(GAMMA_RULE, Production([start_rule]), 1, table[1]))

    for (i, col) in enumerate(table)
        DEBUG && println("## $(i - 1) $(col.token)")
        for state in col.states
            DEBUG && println(str(state))
            if iscompleted(state)
                DEBUG && println("> state $(str(state)) is completed!")
                complete(col,state)
            else
                term = next_term(state)
                DEBUG && println("term = $(str(term))")
                if (isa(term,Rule))
                    predict(col,term)
                    DEBUG && println("after predict: col.states =\n\t", join([str(s) for s in col.states],"\n\t"))
                elseif (i < length(table))
                    scan(table[i + 1], state, term)
                    DEBUG && println("after scan: table[i + 1].states =\n\t", join([str(s) for s in table[i + 1].states],"\n\t"))
                end
            end
            DEBUG && println()
        end
    end
    # find gamma rule in last table column (otherwise fail)
#     gamma_states = filter(s->s.name == GAMMA_RULE, collect(table[end]))
    for state in table[end]
        DEBUG && println(str(state))
        if state.name == GAMMA_RULE && iscompleted(state)
            return state
        end
    end
#     !isempty(gamma_states) && iscompleted(gamma_states[1]) && return gammma_state[1]
    throw("parsing failed")
end

predict(col::Column, rule::Rule) = map(prod -> push!(col, State(rule.name, prod, 1, col)), rule.productions)

function scan(column::Column, state::State, token::String)
    if (token != column.token)
        return nothing
    end
    push!(column, State(state.name, state.production, state.dot_index + 1, state.start_column))
end

function complete(column::Column, state::State)
    if !iscompleted(state)
        return nothing
    end
    for st in state.start_column
        term = next_term(st)
        if !isa(term,Rule)
            continue
        end
        if (term.name == state.name)
            push!(column, State(st.name, st.production, st.dot_index + 1, st.start_column))
        end
    end
end

function str(grammar::Grammar)
    buf = IOBuffer()
    for r in values(grammar.rules_index)
        write(buf,str(r) * "\n")
    end
    return String(take!(buf))
end

function str(rule::Rule)
    lhs = rule.name
    rhs = join(map(str, rule.productions), " | ")
    return "$lhs -> $rhs"
end

function str(production::Production)
    parts = [isa(t,Rule) ? t.name : "\"$t\"" for t in production.terms ]
    return join(parts," ")
end

function str(state::State)
    terms = [isa(t,Rule) ? t.name : "\"$t\"" for t in state.production.terms]
    insert!(terms, state.dot_index, "•")
    _end = state.end_column == nothing ? -1 : state.end_column.index - 1
    return "$(state.name)\t->\t$(join(terms," "))\t\t[$(state.start_column.index - 1)-$(_end)]"
end

str(s::String) = "\"$s\""

# http://matt.might.net/articles/parsing-with-derivatives/
end

max_stars_repo_name: brambg/HyperGraphTools

max_issues_repo_name: brambg/HyperGraphTools

max_forks_repo_name: brambg/HyperGraphTools


text:
using GeometryTypes
using Makie, Colors

struct Ring
    movement::Float32
    radius::Float32
    speed::Float32
    spin::Vec2f0
end
function cartesian(ll::Point2)
    return Point3(
        cos(ll[1]) * sin(ll[2]),
        sin(ll[1]) * sin(ll[2]),
        cos(ll[2])
    )
end
fract(x) = x - floor(x)
function positions(rings, index, time, audio)
    position = Point3f0(0.0)
    precision = 0.2f0
    for ring in rings
        position += ring.radius * cartesian(
            precision *
            index *
            Point2f0(ring.spin + Point2f0(sin(time * ring.speed), cos(time * ring.speed)) * ring.movement)
        )
    end
    amplitude = audio[round(Int, clamp(fract(position[1] * 0.1), 0, 1) * (25000-1)) + 1]; # index * 0.002
    position *= 1.0 + amplitude * 0.5;
    position
end
scene = Scene()
rings = [Ring(0.1, 1.0, 0.00001, (0.2, 0.1)), Ring(0.1, 0.0, 0.0002, (0.052, 0.05))]
N = 25000
t_audio = sin.(linspace(0, 10pi, N)) .+ (cos.(linspace(-3, 7pi, N)) .* 0.6) .+ (rand(Float32, N) .* 0.1) ./ 2f0
start = time()
t = (time() - start) * 100
pos = positions.((rings,), 1:N, t, (t_audio,))
l = lines(pos, color = map(x-> RGBA{Float32}(x, 0.6), colormap("RdBu", N)), thickness = 0.6f0)
center!(scene)
io = VideoStream(scene, homedir() * "/Desktop/", "moire")
for i = 1:1000
    t = (time() - start) * 100
    pos .= positions.((rings,), 1:N, t, (t_audio,))
    l[:positions] = pos
    recordframe!(io)
    sleep(1/60)
end
finish(io, "mp4")


text:
import SpecialFunctions
const Γ = SpecialFunctions.gamma

# * Properties
@doc raw"""
    turning_point(η, ℓ)

Computes the turning point ``\rho_{\mathrm{TP}}`` of the Coulomb
functions, which separates the monotonic region
``\rho<\rho_{\mathrm{TP}}`` from the oscillatory region
``\rho>\rho_{\mathrm{TP}}``.

See Eq. (2) of

- Barnett, A. (1982). COULFG: Coulomb and Bessel Functions and Their
  Derivatives, for Real Arguments, by Steed's Method. Computer Physics
  Communications, 27(2),
  147–166. http://dx.doi.org/10.1016/0010-4655(82)90070-4

or

- Barnett, A. R. (1996). The Calculation of Spherical Bessel and
  Coulomb Functions. In (Eds.), Computational Atomic Physics
  (pp. 181–202). : Springer Berlin Heidelberg.
"""
turning_point(η, ℓ) = η + √(η^2 + ℓ*(ℓ+1))

# https://dlmf.nist.gov/33.2
coulomb_normalization(η, ℓ) = 2^ℓ*exp(-π*η/2)*abs(Γ(ℓ+1+im*η))/Γ(2ℓ+2)

# * Continued fractions

coulomb_fraction1(x::T, η::T, n::Integer; kwargs...) where T =
    lentz_thompson((n+1)/x+η/(n+1),
                   k -> -(one(T)+η^2/(n+k)^2),
                   k -> (2(n+1)+1)*(inv(x)+η/((n+k)^2+(n+k))); kwargs...)

function coulomb_fraction2(x::T, η::T, n::Integer, ω; kwargs...) where T
    imω = im*ω
    r = lentz_thompson(x-η,
                       k -> (imω*η - n - 1 + k)*(imω*η + n + k),
                       k -> 2*(x - η + imω*k); kwargs...)
    ((imω/x)*r[1],r[2:end]...)
end

# * Recurrences

function coulomb_downward_recurrence!(F, F′, x⁻¹::T, η, ℓ::UnitRange, cf1, s) where T
    Fₙ = s ? 1 : -1
    F′ₙ = cf1*Fₙ
    η² = η^2
    for (i,n) in zip(reverse(eachindex(ℓ)), reverse(ℓ))
        S = n*x⁻¹ + η/n
        R = √(1 + η²/n^2)

        Fₙ₋₁ = (S*Fₙ + F′ₙ)/R
        F′ₙ₋₁ = S*Fₙ₋₁ - R*Fₙ

        F[i] = Fₙ
        F′[i] = F′ₙ
        Fₙ = Fₙ₋₁
        F′ₙ = F′ₙ₋₁
    end
end

function coulomb_upward_recurrence!(G, G′, x⁻¹::T, η, ℓ::UnitRange) where T
    Gₙ = G[1]
    G′ₙ = G′[1]
    η² = η^2
    for i = 2:length(ℓ)
        n = ℓ[i]

        S = (n+1)*x⁻¹ + η/(n+1)
        R = √(1 + η²/(n+1)^2)

        G[i] = Gₙ₊₁ = (S*Gₙ - G′ₙ)/R
        G′[i] = G′ₙ₊₁ = R*Gₙ - S*Gₙ₊₁

        Gₙ = Gₙ₊₁
        G′ₙ = G′ₙ₊₁
    end
end

coulomb_upward_recurrence!(::Nothing, ::Nothing, args...; _...) = nothing

# * Driver

function coulombs!(F::FF, F′::FF, G::GG, G′::GG, x::T, η::T, ℓ::UnitRange; verbosity=0, kwargs...) where {FF,GG,T<:Number}
    verbose = get(kwargs, :verbose, false)
    ρtp = turning_point(η, first(ℓ))
    ρtp > x && verbosity > 0 &&
        @warn "Turning point ρ_TP = $(ρtp) > $(x) which may result in loss of accuracy, consider decreasing first ℓ below $(first(ℓ))."

    if iszero(x)
        F .= false
        for (i,ℓ) in enumerate(ℓ)
            C = coulomb_normalization(η, ℓ)
            F′[i] = (ℓ+1)*C*x^ℓ
            G[i] = 1/((2ℓ+1)*C) # This seems wrong for ℓ ≥ 1
            G′[i] = -T(Inf) # This is not necessarily true
        end
        return
    end

    verbosity > 1 && @show x, η, ℓ
    cf1,n,_,s,converged = coulomb_fraction1(x, η, ℓ[end]; verbosity=verbosity, kwargs...)
    converged || verbosity > 0 && @warn "Consider increasing ℓmax beyond $(ℓ[end])"

    x⁻¹ = inv(x)
    coulomb_downward_recurrence!(F, F′, x⁻¹, η, ℓ, cf1, s)

    cf2,n,_,s,converged = coulomb_fraction2(x, η, ℓ[1], 1; verbosity=verbosity, kwargs...)

    p,q = real(cf2),imag(cf2)
    Fₘ = F[1]
    F′ₘ = F′[1]
    Gₘ = (F′ₘ-p*Fₘ)/q
    G′ₘ = p*Gₘ - q*Fₘ

    ω₁ = 1/√(F′ₘ*Gₘ - Fₘ*G′ₘ)

    lmul!(ω₁, F)
    lmul!(ω₁, F′)

    if !isnothing(G)
        G[1] = ω₁*Gₘ
        G′[1] = ω₁*G′ₘ

        coulomb_upward_recurrence!(G, G′, x⁻¹, η, ℓ)
    end
end

# * Interface

"""
    coulombs!(F, F′, G, G′, x::AbstractVector; kwargs...)

Loop through all values of `x` and compute all Coulomb
functions, storing the results in the preallocated matrices `F`, `F′`,
`G`, `G′`.
"""
function coulombs!(F, F′, G, G′, x::AbstractVector, η, ℓ::AbstractVector; kwargs...)
    size(F,1) == size(F′,1) == size(G,1) == size(G′,1) == length(x) &&
        size(F,2) == size(F′,2) == size(G,2) == size(G′,2) == length(ℓ) ||
        throw(DimensionError("The dimension of the output arrays do not agree"))
    for (i,x) in enumerate(x)
        coulombs!(view(F, i, :),
                  view(F′, i, :),
                  view(G, i, :),
                  view(G′, i, :),
                  x, η, ℓ; kwargs...)
    end
end

"""
    coulombs(x::AbstractVector, nℓ; kwargs...)

Convenience wrapper around [`coulombs!`](@ref) that preallocates output
matrices of the appropriate dimensions.
"""
function coulombs(x::AbstractVector{T}, η, ℓ::AbstractVector; kwargs...) where T
    nx = length(x)
    nℓ = length(ℓ)
    F = zeros(T, nx, nℓ)
    F′ = zeros(T, nx, nℓ)
    G = zeros(T, nx, nℓ)
    G′ = zeros(T, nx, nℓ)
    coulombs!(F, F′, G, G′, x, η, ℓ; kwargs...)
    F, F′, G, G′
end

coulombs(x::AbstractVector, η, nℓ::Integer; kwargs...) =
    coulombs(x, η, 0:nℓ-1; kwargs...)

export coulombs!, coulombs


text:
module UW
using SeisIO, SeisIO.Quake
using Dates: DateTime
import SeisIO: BUF, KW, checkbuf!, checkbuf_8!, dtconst, fillx_i16_be!, fillx_i32_be!, sμ
import SeisIO.Quake: unsafe_convert
import SeisIO.Formats: formats, FmtVer, FormatDesc, HistVec

include("UW/uwdf.jl")
include("UW/uwpf.jl")
include("UW/uwevt.jl")
include("UW/desc.jl")

# exports
export readuwevt, uwdf, uwdf!, uwpf, uwpf!, formats

end


text:
#t = randn(100)
#y = randn(100)
using CSV
df = CSV.read("lc.csv")
using DataFrames
df2 = sort!(df)
t = convert(Array{Float64},df2[1])
y = convert(Array{Float64},df2[4])

N = length(t)
jPeak = 1
jDip = 1
ratio = 1.005
flag = "Dip"
Peak = ones(N,2) # 1 for t, 2 for y
Dip = ones(N,2)
PeakTmp = ones(N,2)
DipTmp = ones(N,2)
PeakTmp[jPeak,1:2] = [t[1],y[1]]
DipTmp[jDip,:] = [t[1],y[1]]
Peak[jPeak,:] = [t[1],y[1]]
for i in 2:N
    if y[i] > y[i-1]
        if y[i] > PeakTmp[jPeak,2]
            PeakTmp[jPeak,:] = [t[i],y[i]]
            println("PeakTmp[jPeak,:]", PeakTmp[jPeak,:], " y[i]: ", y[i], " y[i-1]: ", y[i-1])
        end
#        if flag == "Dip"
#            DipTmp[jDip,:] = [t[i],y[i]]
            #println(DipTmp[jDip,:])
#        end
        if flag == "Peak"
            println("Now in peak")
            if PeakTmp[jPeak,2]/ratio > DipTmp[jDip,2]
                Dip[jDip,1:2] = DipTmp[jDip,1:2]
                jDip = jDip + 1
                println("jDip:", jDip)
                flag = "Dip"
                DipTmp[jDip,:] = [t[i],y[i]]
            end
        end
    end
    if y[i] < y[i-1]
        if y[i] < DipTmp[jDip,2]
            DipTmp[jDip,:] = [t[i],y[i]]
            println("DipTmp[jDip,:] ", DipTmp[jDip,:], " y[i]: ", y[i], " y[i-1]: ", y[i-1])
        end
#        if flag == "Peak"
#            PeakTmp[jPeak,:] = [t[i],y[i]]
#            println("Now in Peak 2, PeakTmp:", PeakTmp[jPeak,:])
#        end
        if flag == "Dip"
            println(i, " Now in Dip 2, PeakTmp: ", PeakTmp[jPeak,2],", DipTmp:", DipTmp[jDip,2])
            if PeakTmp[jPeak,2]/ratio > DipTmp[jDip,2]
                Peak[jPeak,1:2] = PeakTmp[jPeak,1:2]
                jPeak = jPeak + 1
                println("jPeak:", jPeak)
                flag = "Peak"
                PeakTmp[jPeak,:] = [t[i],y[i]]
            end
        end
    end
end
Peak2 = ones(jPeak,2)
for i=1:jPeak
    Peak2[i,:] = Peak[i,:]
end
Peaks = Peak2


text:
module PeaceFounder

### Internal modules
#include("BraidChains/BraidChains.jl")
import DemeNet
import DemeNet: Plugin
const plugin = Plugin[@__MODULE__]

### The rest of the code

include("types.jl")
include("serialization.jl")
include("tools.jl")

import PeaceVote: record
record(config::PeaceFounderConfig,data) = record(config.braidchain,data)

import PeaceVote: sync!
sync!(chain::BraidChain,config::PeaceFounderConfig) = sync!(chain,config.braidchain.syncport)

using DemeNet: datadir

function updateconfig(deme::Deme)
    ddir = datadir(deme.spec.uuid)
    from = ddir * "/ledger/config/00000000"
    to = ddir * "/PeaceFounder.toml"
    ### Before copying one can verify that the config had been signed by the maintainer
    cp(from,to,force=true)
end

function init(deme::Deme,port::Port)
    ledger = Ledger(deme.spec.uuid)
    sync!(ledger,port)
    cp(chain.ledger * "/config/00000000",fname,force=true)
end

init(deme::Deme,port::Dict) = init(deme,Port(port))

function config(deme::Deme)
    cert = deserialize(deme,Certificate{PeaceFounderConfig})
    intent = Intent(cert,deme.notary)
    @assert intent.reference==deme.spec.maintainer
    return intent.document
end

### The initializer ###

using DemeNet: AbstractInitializer

struct Init <: AbstractInitializer
    deme::Deme
end

plugin(::Type{AbstractInitializer},deme::Deme) = Init(deme)
DemeNet.init(x::Init,port::Dict) = init(x.deme,port)
DemeNet.config(x::Init) = config(x.deme)
DemeNet.updateconfig(x::Init) = updateconfig(x.deme)

import PeaceVote: BraidChain
BraidChain(deme::Deme,config::PeaceFounderConfig) = BraidChain(config.braidchain,deme)
BraidChain(deme::Deme) = BraidChain(config(deme),deme) ### This is piracy.

using PeaceVote: AbstractChain
plugin(::Type{AbstractChain}, deme::Deme,config::PeaceFounderConfig) = BraidChain(deme,config)

# function certify(chain::BraidChain,signer::Signer)
#     @assert chain.deme.spec.maintainer==signer.id "You are not eligible to certify PeaceFounder.toml for this deme"
#     sc = deserialize(chain.deme,PeaceFounderConfig)
#     cert = Certificate(sc,signer)
#     writebootrecord(chain.ledger,cert)
#     updateconfig(chain.deme)
# end

function certify(deme::Deme,signer::Signer)
    @assert deme.spec.maintainer==signer.id "You are not eligible to certify PeaceFounder.toml for this deme"
    sc = deserialize(deme,PeaceFounderConfig)
    cert = Certificate(sc,signer)

    writebootrecord(deme,cert)
    updateconfig(deme)
end

end # module

max_stars_repo_name: akels/DemocracyCredit

max_issues_repo_name: akels/DemocracyCredit

max_forks_repo_name: akels/DemocracyCredit


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

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 4c21907e-580c-11eb-1b67-e9342e06a291
push!(LOAD_PATH, "/media/pi/RASPBIAN A/julia", "/media/pi/RASPBIAN A/julia/MyPackages", "/media/pi/RASPBIAN A/julia/MyProjects", "/media/pi/RASPBIAN A/julia/pluto_notebooks")

# ╔═╡ ac57f83e-587f-11eb-06f2-dfccc003eebf
LOAD_PATH

# ╔═╡ a9aa2f9e-59b5-11eb-2c3d-93f40c6b8626
using Makie

# ╔═╡ 23bdedca-568f-11eb-1dac-bda2b3839216
using ImageIO, ImageMagick, InteractiveUtils, PlutoUI, Markdown, Plots

# ╔═╡ f9219eac-5692-11eb-36ea-c3473d601641
begin
	import Pkg
	Pkg.activate(tempdir())
end

# ╔═╡ 87b7da98-5693-11eb-343a-df49ef084c17
begin
	Pkg.add(["Images", "ImageMagick"])
	using Images
end

# ╔═╡ 447d1ef4-569c-11eb-1643-b9b608892a84
begin
	Pkg.add("UrlDownload")
	using UrlDownload
end

# ╔═╡ 4e8fc752-599e-11eb-0fa8-0997b7982c7a
begin
	Pkg.add("FileIO")
	using FileIO
end

# ╔═╡ 75b77a10-59af-11eb-1535-477cc15569ae
begin
Pkg.add("VideoIO")
using VideoIO
end

# ╔═╡ fad639ec-59b4-11eb-0e6e-473f167c95c8
begin
Pkg.add("Makie")
end

# ╔═╡ f730134c-59af-11eb-0be7-23dd1b8c10cd
println(pwd())

# ╔═╡ 03631d6c-59b0-11eb-1694-adc77c3b2f98
vid1 = openvideo("./color.ts")

# ╔═╡ 13bb7778-59b4-11eb-26b4-a385bf49c16a
io = VideoIO.testvideo("annie_oakley")

# ╔═╡ 606cff8c-59b5-11eb-053b-73e603cd801a
VideoIO.playvideo(io)

# ╔═╡ 1841ffc2-59b4-11eb-3d34-31280d300fe6

# ╔═╡ 4d69fe20-5687-11eb-2b0c-e972689711fc
@bind p1 Slider(1:10)

# ╔═╡ 3a38fd4a-568e-11eb-29dd-fb556d0263bb
p1

# ╔═╡ 4defa406-568e-11eb-2f31-9b9c14208910
md"""`power` = $(@bind power Slider(1:10; default=3, show_value=true))"""

# ╔═╡ 203ed75e-5691-11eb-2ce1-319fcaf5005c
Total = power + p1

# ╔═╡ 6fd5e2f6-58dc-11eb-2452-6b81153c9206
@bind fname TextField((30, 7))

# ╔═╡ 1fe5e86a-5691-11eb-0c01-8fbbd89e3ff1
fname

# ╔═╡ a7e12faa-58de-11eb-1ec1-817c0b7ac06f
@bind go Button("Go")

# ╔═╡ a8c4ab7a-58ea-11eb-35cc-e552f47ef4e9
data = let
    go # we reference `go`, making this cell a dependency of the first cell
    rand(4)

end

# ╔═╡ 4e33c052-58df-11eb-32b3-c715c0cb45e2
@bind tf TimeField()

# ╔═╡ ab2ce008-58e3-11eb-09f2-791ce015a696
@bind NewFile FilePicker()

# ╔═╡ e6aa5ba6-58e0-11eb-345b-7d1dcf27b1eb
@bind nf PasswordField()

# ╔═╡ d261a068-58e1-11eb-3bb3-33d3a382378b
@bind cb CheckBox()

# ╔═╡ bf3fc830-58e5-11eb-0082-a7df477c6e39
@bind col ColorStringPicker()

# ╔═╡ 41092694-58e7-11eb-3575-e3bccc642b63
@bind cl Clock()

# ╔═╡ 3643da94-58ef-11eb-0aab-81ef964a1748
@bind dt DateField()

# ╔═╡ 40b790ec-58e7-11eb-284e-33ee25592d26
@bind rd Radio(["Male", "Female"])

# ╔═╡ dc40d2da-58ef-11eb-00c0-15a506975c3d
@bind location Select(["London", "New York"])

# ╔═╡ 86470950-58f0-11eb-0530-57a8c76c03b2
@bind occupation MultiSelect(["D","L","T"])

# ╔═╡ ec07b806-58f3-11eb-3221-09e35beb4fe9
Show(MIME"text/html"(),"""
Iam line one.

	I am line two.

I am line three.$nf

	""")

# ╔═╡ 3bc6e572-58f1-11eb-2e2f-77aa592f16b8
md"""media  MIME  resource """

# ╔═╡ d12381ac-594b-11eb-00d2-cd71a8369a62
audio1 = load("./andersen_zinnsoldat.mp3")

# ╔═╡ 0d203ffa-594d-11eb-0037-c14f0535cfd1
video1 = load("./color.ts")

# ╔═╡ f3d50782-5948-11eb-1b53-bdbbb80d9151
Show(MIME"audio/mp3"(), read("./andersen_zinnsoldat.mp3"))

# ╔═╡ 63b100a4-594b-11eb-210d-49da18818575

# ╔═╡ dbd04822-58df-11eb-188a-6db81486e9d1
string("PassWord-->", nf, " check box1-->", cb, " color-->",col, "   Clock-->", cl, "Radio-->", rd, "Date-->", dt )

# ╔═╡ b7952fbc-59a3-11eb-1b3d-37ee8d447e13
md"""
!!! hint
    If you're stuck, feel free to cheat, this is homework 0 after all 🙃
    Julia has a function called `sqrt`
"""

# ╔═╡ dc69a05e-5695-11eb-39b2-a5b1c77f1c6a
img = download("https://i.imgur.com/VGPeJ6s.jpg")

# ╔═╡ 4ae42cc2-5694-11eb-382a-731cd20be683
philip_file = download("https://i.imgur.com/VGPeJ6s.jpg")

# ╔═╡ a73fe2a8-58f2-11eb-31ed-7ddc409554fd
Show(MIME"image/png"(), read( philip_file))

# ╔═╡ ce97aa22-5698-11eb-07d4-0d734a140258
begin
    struct Wow
    filename
    end

    function Base.show(io::IO, ::MIME"image/png", w::Wow)
    write(io, read(w.filename))
    end
end

# ╔═╡ d348c3b2-5698-11eb-2e3a-777dd2e06997
Wow("/media/pi/RASPBIAN A/julia/Images/corgi.jpg")

# ╔═╡ 663f42d4-5699-11eb-1d0a-65e26b7e7f6f
Wow("/home/pi/in.png")

# ╔═╡ 2058ca04-5783-11eb-17b9-939bb1981c5c
typeof(Wow)

# ╔═╡ 6e7c9ec4-5786-11eb-1d20-598e54949837
RGBX(.5,.5,.5)

# ╔═╡ c11dccca-5786-11eb-1291-0bf70a4b494e
begin
#using Images
MyHald = load("/home/pi/in.png")
end

# ╔═╡ 37b430d8-5796-11eb-1018-f95eab7481d2
dog = load("/media/pi/RASPBIAN A/julia/Images/corgi.jpg")

# ╔═╡ 3a30109a-5788-11eb-3b59-275c20f7a1ce
typeof(MyHald)

# ╔═╡ 885ac5d0-5788-11eb-06fe-7176b0569f2c
size(MyHald)

# ╔═╡ 9fa37750-5788-11eb-3c99-29f77db69e21
#(h, w) = size(myHald)
#MyHald[(h÷2):h (w÷10): (9w÷10)]
# half way down,
MyHald[(27÷2):27, (27÷2):(9*27÷10)]

# ╔═╡ d55125ae-5789-11eb-16af-8f3725160a12

 MyHald[(27÷2):27, (27÷2): (7*27÷10)]

# ╔═╡ d8fdd574-578b-11eb-000a-ef83b71e25b2
[MyHald reverse(MyHald, dims=1)
MyHald        MyHald
MyHald        MyHald         ]

# ╔═╡ d4fc9a66-5789-11eb-2da1-d5ae1269e2f0
MyHaldCopy = copy(MyHald)

# ╔═╡ 1940094c-56f1-11eb-1812-553f21b3eb2b
var = 5

# ╔═╡ d3bc1140-5789-11eb-30f3-bfc1219d64a3
for i in 1:10
	for j in 1:15
MyHaldCopy[i,j] = RGBX(.9,.2,.2)
	end
end

# ╔═╡ d34a62a0-5789-11eb-0cbf-db00d11ae9fa
MyHaldCopy

# ╔═╡ c19e0c6a-578f-11eb-1b6d-3b15ff9529a1
begin
	MyHaldCopy2 = copy(MyHald)
	MyHaldCopy2[5:10, 8:20] .= RGBX(0,1,0)
	MyHaldCopy2
#   .=  look into aray and do it to each individually
end

# ╔═╡ 30e20880-5793-11eb-323c-a5ef82071356
function redify(c)
	return RGB(c.r,0,0)
end

# ╔═╡ 2d4b1d44-5791-11eb-09b5-290cdacfdba2
begin
	color = RGBX(.5, .9, .9)
	[color, redify(color)]
end

# ╔═╡ dde234aa-5799-11eb-075c-97b06ff5b4d3
fill(.25, (2,2))

# ╔═╡ 636159cc-5799-11eb-3633-eb7040c6fab1
#blur = [1:2 3:4]
function blur(n)
 return [.25 .25; .25 .25]
end

# ╔═╡ 64e58744-5796-11eb-108b-29f3c10d9834
begin
	convolve(dog, blur(2))
#blur()  not defined,  it is a 2x2 aray of floats   .25 .25
#                                                   .25 .25
end

# ╔═╡ 649062a0-5796-11eb-369e-072beec02205
varinfo(ImageMagick)

# ╔═╡ 64421e38-5796-11eb-357c-23cfc8b93e86

# ╔═╡ 04e2c54c-569d-11eb-1b5b-7997183e9c9d
pwd()

# ╔═╡ c553586e-56f0-11eb-024e-ff83b61885e3
varinfo(Main,r".*field.*")

# ╔═╡ 6d2cda82-5710-11eb-2e7e-d9094e5e8ad1
typeof(555555555555555555555555555555555555555555555555555555555555555)

# ╔═╡ ad2f732e-5710-11eb-1cc7-d712d3a78587
typeof("Ḧello world")

# ╔═╡ d399e060-5710-11eb-04dc-dfb5f59afd76
typeof(-2039415715155.262326)

# ╔═╡ d1faaf20-5710-11eb-38ba-ed675c6b2b02
"
    vms
Document for variable vms
use: this is the velosity
give in meters per second
"
const vms = 23

# ╔═╡ 040048d8-5716-11eb-29f7-f18b53b66ffe
"""
    bar(x[, y])

Compute the Bar index between `x` and `y`. If `y` is missing, compute
the Bar index between all pairs of columns of `x`.

# Examples
```julia-repl
julia> bar([1, 2], [1, 2])
1
```
"""
function bar(x, y)
end

# ╔═╡ 89fc8da8-5718-11eb-3384-03eb09e53581
?bar

# ╔═╡ a8a08162-571b-11eb-0018-ed11de6bac88
?vms

# ╔═╡ 6da694cc-5720-11eb-243c-d76ce6955f67
begin
	a = "Now is the time for apostrophy's"
	b = "Here is string #"
end

# ╔═╡ 8a4bf2b4-572f-11eb-3479-8d57725e4476
z = readlines("/home/pi/bc", keep=true)

# ╔═╡ d82df792-5735-11eb-1d98-a1a67422283e
varinfo(InteractiveUtils)

# ╔═╡ f3209632-5775-11eb-1751-53e997c0f0a2
clipboard(z)

# ╔═╡ aaf86596-587f-11eb-08d8-2356de69699b

# ╔═╡ 47e02b4a-587e-11eb-34b0-95ba8b809566

# ╔═╡ 03c96442-587e-11eb-35b5-2bd15578f118

# ╔═╡ 8affd2a2-5c0c-11eb-1eba-dbe7bc378ca9
md"""
!["A pretty tiger"](https://upload.wikimedia.org/wikipedia/commons/5/56/Tiger.50.jpg)
"""

# ╔═╡ 132dc3a4-5795-11eb-0145-656e5574331a
out = redify.(MyHald)

# ╔═╡ fe33c784-57ab-11eb-18d6-dda01ff3ab90
out = save("/home/pi/out.png")

# ╔═╡ Cell order:
# ╠═f9219eac-5692-11eb-36ea-c3473d601641
# ╠═87b7da98-5693-11eb-343a-df49ef084c17
# ╠═447d1ef4-569c-11eb-1643-b9b608892a84
# ╠═4e8fc752-599e-11eb-0fa8-0997b7982c7a
# ╠═75b77a10-59af-11eb-1535-477cc15569ae
# ╠═fad639ec-59b4-11eb-0e6e-473f167c95c8
# ╠═a9aa2f9e-59b5-11eb-2c3d-93f40c6b8626
# ╠═f730134c-59af-11eb-0be7-23dd1b8c10cd
# ╠═03631d6c-59b0-11eb-1694-adc77c3b2f98
# ╠═13bb7778-59b4-11eb-26b4-a385bf49c16a
# ╠═606cff8c-59b5-11eb-053b-73e603cd801a
# ╠═1841ffc2-59b4-11eb-3d34-31280d300fe6
# ╠═23bdedca-568f-11eb-1dac-bda2b3839216
# ╠═4d69fe20-5687-11eb-2b0c-e972689711fc
# ╠═3a38fd4a-568e-11eb-29dd-fb556d0263bb
# ╠═4defa406-568e-11eb-2f31-9b9c14208910
# ╠═203ed75e-5691-11eb-2ce1-319fcaf5005c
# ╠═6fd5e2f6-58dc-11eb-2452-6b81153c9206
# ╠═1fe5e86a-5691-11eb-0c01-8fbbd89e3ff1
# ╠═a7e12faa-58de-11eb-1ec1-817c0b7ac06f
# ╠═a8c4ab7a-58ea-11eb-35cc-e552f47ef4e9
# ╠═4e33c052-58df-11eb-32b3-c715c0cb45e2
# ╠═ab2ce008-58e3-11eb-09f2-791ce015a696
# ╠═e6aa5ba6-58e0-11eb-345b-7d1dcf27b1eb
# ╠═d261a068-58e1-11eb-3bb3-33d3a382378b
# ╠═bf3fc830-58e5-11eb-0082-a7df477c6e39
# ╠═41092694-58e7-11eb-3575-e3bccc642b63
# ╠═3643da94-58ef-11eb-0aab-81ef964a1748
# ╠═40b790ec-58e7-11eb-284e-33ee25592d26
# ╠═dc40d2da-58ef-11eb-00c0-15a506975c3d
# ╠═86470950-58f0-11eb-0530-57a8c76c03b2
# ╠═ec07b806-58f3-11eb-3221-09e35beb4fe9
# ╠═3bc6e572-58f1-11eb-2e2f-77aa592f16b8
# ╠═a73fe2a8-58f2-11eb-31ed-7ddc409554fd
# ╠═d12381ac-594b-11eb-00d2-cd71a8369a62
# ╠═0d203ffa-594d-11eb-0037-c14f0535cfd1
# ╠═f3d50782-5948-11eb-1b53-bdbbb80d9151
# ╠═63b100a4-594b-11eb-210d-49da18818575
# ╠═dbd04822-58df-11eb-188a-6db81486e9d1
# ╠═b7952fbc-59a3-11eb-1b3d-37ee8d447e13
# ╠═dc69a05e-5695-11eb-39b2-a5b1c77f1c6a
# ╠═4ae42cc2-5694-11eb-382a-731cd20be683
# ╠═ce97aa22-5698-11eb-07d4-0d734a140258
# ╠═d348c3b2-5698-11eb-2e3a-777dd2e06997
# ╠═663f42d4-5699-11eb-1d0a-65e26b7e7f6f
# ╠═2058ca04-5783-11eb-17b9-939bb1981c5c
# ╠═6e7c9ec4-5786-11eb-1d20-598e54949837
# ╠═c11dccca-5786-11eb-1291-0bf70a4b494e
# ╠═37b430d8-5796-11eb-1018-f95eab7481d2
# ╠═3a30109a-5788-11eb-3b59-275c20f7a1ce
# ╠═885ac5d0-5788-11eb-06fe-7176b0569f2c
# ╠═9fa37750-5788-11eb-3c99-29f77db69e21
# ╠═d55125ae-5789-11eb-16af-8f3725160a12
# ╠═d8fdd574-578b-11eb-000a-ef83b71e25b2
# ╠═d4fc9a66-5789-11eb-2da1-d5ae1269e2f0
# ╟─1940094c-56f1-11eb-1812-553f21b3eb2b
# ╠═d3bc1140-5789-11eb-30f3-bfc1219d64a3
# ╠═d34a62a0-5789-11eb-0cbf-db00d11ae9fa
# ╠═c19e0c6a-578f-11eb-1b6d-3b15ff9529a1
# ╠═30e20880-5793-11eb-323c-a5ef82071356
# ╠═2d4b1d44-5791-11eb-09b5-290cdacfdba2
# ╠═132dc3a4-5795-11eb-0145-656e5574331a
# ╠═fe33c784-57ab-11eb-18d6-dda01ff3ab90
# ╠═dde234aa-5799-11eb-075c-97b06ff5b4d3
# ╠═636159cc-5799-11eb-3633-eb7040c6fab1
# ╠═64e58744-5796-11eb-108b-29f3c10d9834
# ╠═649062a0-5796-11eb-369e-072beec02205
# ╠═64421e38-5796-11eb-357c-23cfc8b93e86
# ╠═04e2c54c-569d-11eb-1b5b-7997183e9c9d
# ╠═c553586e-56f0-11eb-024e-ff83b61885e3
# ╠═6d2cda82-5710-11eb-2e7e-d9094e5e8ad1
# ╠═ad2f732e-5710-11eb-1cc7-d712d3a78587
# ╠═d399e060-5710-11eb-04dc-dfb5f59afd76
# ╠═d1faaf20-5710-11eb-38ba-ed675c6b2b02
# ╠═040048d8-5716-11eb-29f7-f18b53b66ffe
# ╠═89fc8da8-5718-11eb-3384-03eb09e53581
# ╠═a8a08162-571b-11eb-0018-ed11de6bac88
# ╠═6da694cc-5720-11eb-243c-d76ce6955f67
# ╠═8a4bf2b4-572f-11eb-3479-8d57725e4476
# ╟─d82df792-5735-11eb-1d98-a1a67422283e
# ╠═f3209632-5775-11eb-1751-53e997c0f0a2
# ╠═4c21907e-580c-11eb-1b67-e9342e06a291
# ╠═ac57f83e-587f-11eb-06f2-dfccc003eebf
# ╟─aaf86596-587f-11eb-08d8-2356de69699b
# ╟─47e02b4a-587e-11eb-34b0-95ba8b809566
# ╠═03c96442-587e-11eb-35b5-2bd15578f118
# ╠═8affd2a2-5c0c-11eb-1eba-dbe7bc378ca9

max_stars_repo_name: paradocs/hello-world

max_issues_repo_name: paradocs/hello-world

max_forks_repo_name: paradocs/hello-world


text:
# Test utililties.jl

# Tests could be more comprehensive (of course).  In some cases all
# that is done is to execute the function in some way just to make
# sure no exceptions are thrown as a result of something breaking with
# a new version of Julia.

println("testing utilites")

using PerceptualColourMaps
using Test, Statistics

# bbspline
P = [0 1 2 3
     0 1 2 3]
k = 3
N = 10
S = bbspline(P, k, N)
#dx = Base.LinAlg.gradient(S[1,:],1)
#dy = Base.LinAlg.gradient(S[2,:],1)
# There seems to be an awkward method ambiguity between Base.LinAlg.gradient and
# Calculus.gradient, hence the simple replacement code below.
dx = S[1, 2:end] - S[1, 1:(end-1)]
dy = S[2, 2:end] - S[2, 1:(end-1)]

slope = dy./dx
@test maximum(slope .- 1) < 1e-3

# pbspline
P = [0 1 2 3 0
     0 1 2 3 0]
k = 3
N = 100
S = pbspline(P, k, N)

# gaussfilt1d
a = zeros(20,1)
a[10] = 1
sigma = 2
sa = gaussfilt1d(a, sigma)
@test abs(sa[10] - 1/(sigma*sqrt(2pi))) < 1e-3
@test abs(sum(sa) - 1) < 1e-3

# interp1
x = [0, 2, 7]
y = 2*x
xi = [1.1, 0.5, 6]
yi = interp1(x,y,xi)
@test maximum(abs.(yi - 2*xi) .< eps())

# normalise
a = rand(5,3)
an = normalise(a)
@test isapprox(minimum(an), 0)
@test isapprox(maximum(an), 1)

reqmean = 5
reqvar = 7
an = normalise(a, reqmean, reqvar)
@test isapprox(mean(an), reqmean)
@test isapprox(var(an), reqvar)

# histtruncate
rows = 100;
cols = 200;
img = rand(rows, cols)
lHistCut = 2
uHistCut = 4
htimg = histtruncate(img, lHistCut, uHistCut)
minval = minimum(htimg)
maxval = maximum(htimg)
# Check that the number of saturated pixels at each extreme is
# approximately lHistCut and uHistCut. Note some integer rounding has
# to occur.
v = findall(abs.(htimg.-minval).<eps())
@test abs(length(v)/(rows*cols) * 100 - lHistCut) < 1

v = findall(abs.(htimg.-maxval).<eps())
@test abs(length(v)/(rows*cols) * 100 - uHistCut) < 1

# sineramp
sr = sineramp()
sr = sineramp((256, 512), 12.5, 8, 2)

# circlesineramp
(img, alpha) = circlesineramp()
(img, alpha) = circlesineramp(512, pi/10, 8, 2, false)


text:
###############################################################################
#
# Leaf biological parameters
#
###############################################################################
"""
    mutable struct LeafBios{FT}

A struct of leaf biological parameters

# Fields
$(TYPEDFIELDS)
"""
Base.@kwdef mutable struct LeafBios{FT}
    # local storage of dimension information
    "Number of wave length"
    nWL ::Int = 10
    "Number of wave length for excitation"
    nWLE::Int = 10
    "Number of wave length for SIF"
    nWLF::Int = 10

    "Leaf structure parameter"
    N   ::FT = FT(1.4  )
    "Chlorophyll a+b content `[µg cm⁻²]`"
    Cab ::FT = FT(40.0 )
    "Carotenoid content `[µg cm⁻²]`"
    Car ::FT = FT(10.0 )
    "Anthocynanin content `[µg cm⁻²]`"
    Ant ::FT = FT(0.0  )
    "Senescent material fraction"
    Cs  ::FT = FT(0.0  )
    "Equivalent water thickness `[cm]`"
    Cw  ::FT = FT(0.009)
    "Dry matter content (dry leaf mass per unit area) `[g cm⁻²]`"
    Cm  ::FT = FT(0.012)
    "Fractionation between Zeaxanthin and Violaxanthin in Car (1=all Zeaxanthin) (-)"
    Cx  ::FT = FT(0.0  )
    "Leaf fluorescence efficiency (Fo standard)"
    fqe ::FT = FT(1.0  )
    "Broadband thermal reflectance (-)"
    ρ_LW::FT = FT(0.01 )
    "Broadband thermal transmission (-)"
    τ_LW::FT = FT(0.01 )

    "Shortwave leaf reflectance"
    ρ_SW       ::Array{FT,1} = zeros(FT, nWL)
    "Shortwave leaf transmission"
    τ_SW       ::Array{FT,1} = zeros(FT, nWL)
    "Shortwave absorption"
    α_SW       ::Vector{FT}  = zeros(FT, nWL)
    "Relative absorbtion by Chlorophyll+Car"
    kChlrel    ::Array{FT,1} = zeros(FT, nWL)
    "Relative absorbtion by Chlorophyll"
    kChlrel_old::Array{FT,1} = zeros(FT, nWL)
    "Fluorescence excitation matrix backwards"
    Mb         ::Array{FT,2} = zeros(FT,(nWLF,nWLE))
    "Fluorescence excitation matrix forwards"
    Mf         ::Array{FT,2} = zeros(FT,(nWLF,nWLE))
    "Doubling adding layers"
    ndub       ::Int = 10
end

max_stars_repo_name: climate-machine/Land


text:
abstract 動物

type 捕食者 <: 動物
  位置 :: Array{Int64, 1}
  捕食者(x, y) = new([y, x])
end

type 獲物 <: 動物
  位置 :: Array{Int64, 1}
  獲物(く, ゆ) = new([ゆ, く])
end

type 虚 <: 動物
  虚() = new()
end

type 天地
  地図 :: Array{動物, 2}
  天地(ん) = new(fill(虚(), ん, ん))
end

function 動物を加える{つ}(世界, く, ゆ, 物::Type{つ})
  if(typeof(世界.地図[ゆ, く]) != 虚)
    error("場所が取られた")
  end
  世界.地図[ゆ, く] = 物(く, ゆ)
end

function 距離(動物一::動物, 動物二::動物)
  return abs(動物一.位置[1] - 動物二.位置[1]) + abs(動物一.位置[1] - 動物二.位置[1])
end

function 除去(世界::天地, 動物一::動物)
  世界.地図[動物一.位置[1], 動物一.位置[2]] = 虚()
end

function 対話(世界::天地, 動物一::捕食者, 動物二::獲物)
  除去(世界, 動物二)
end

function 移転する(世界::天地, 動物一::獲物)
  for い in 1:size(世界.地図)[1]
    for や in 1:size(世界.地図)[2]
      if(typeof(世界.地図[い, や]) == 虚)
        世界.地図[い, や] = 動物一
        動物一.位置 = [い, や]
        return
      end
    end
  end
end

function 対話(世界::天地, 動物一::獲物, 動物二::捕食者)
  除去(世界, 動物一)
  移転する(世界::天地, 動物一::獲物)
end

function 対話(世界::天地, 動物一::捕食者, 動物二::捕食者)
  return "ぶるっっっ"
end

function 対話(世界::天地, 動物一::獲物, 動物二::獲物)
  return "べえええ"
end

# s =


text:
export Simulation

include("Grid.jl")

Base.@kwdef mutable struct Simulation{Dim}

		grid::Grid{Dim}

		λ₀::Float

		J_x::Array{CFloat,Dim} = zeros(CFloat,size(grid))
		J_y::Array{CFloat,Dim} = zeros(CFloat,size(grid))
		J_z::Array{CFloat,Dim} = zeros(CFloat,size(grid))

		Q::BitArray{Dim}	 = falses(size(grid))

		S_x::Array{CFloat,Dim} = ones(CFloat,size(grid))
		S_y::Array{CFloat,Dim} = ones(CFloat,size(grid))
		S_z::Array{CFloat,Dim} = ones(CFloat,size(grid))

		activate_tfsf::Bool = false

		ϵᵣ_xx::Array{CFloat,Dim}  = ones(CFloat,size(grid))
		ϵᵣ_yy::Array{CFloat,Dim}  = ones(CFloat,size(grid))
		ϵᵣ_zz::Array{CFloat,Dim}  = ones(CFloat,size(grid))

		μᵣ_xx::Array{CFloat,Dim}  = ones(CFloat,size(grid))
		μᵣ_yy::Array{CFloat,Dim}  = ones(CFloat,size(grid))
		μᵣ_zz::Array{CFloat,Dim}  = ones(CFloat,size(grid))

        e⁻ⁱᵏᴸ::NTuple{Dim,CFloat} = ntuple(i -> CFloat(1.0 + 0.0im) , dimension(grid))

		E::Vector{CFloat}  = zeros(CFloat,3 * ncells(grid))
	end

function Base.show(io::IO, sim::Simulation{Dim}) where {Dim}
        print(io, "A $(Dim)D Simulation is initialized.\n")
end

function (sim::Simulation)(sym::Symbol,dir::Direction, maskF::Function)
	mat = sim(sym,dir)
	mask = maskF.(Coordinates(sim.grid,p̂))
	mat[mask]
end

include("GridDerivatives.jl")
include("Set.jl")
include("Get.jl")
include("Add.jl")
include("Geometry.jl")
include("Material.jl")
include("TFSF.jl")
include("PML.jl")
include("Sources/Source.jl")
include("BC.jl")
include("SliceSimulation.jl")


text:
function run_simulation(p::PrecipParam)
    timer, mesh, soln, tsoln, prob = initialize_data(p)

    # initial brinkman
    mpbparam = BrinkmanMPParam(1.0,mpbFriction(soln.θf.u),soln.θf.u)
    sol = solve(prob.fluid,mesh.fMesh,mpbparam)
    tsoln.fluid.u = sol.u
    tsoln.fluid.v = sol.v
    tsoln.fluid.p = sol.p

    intializeTimeSteppingTimer!(timer)
    doPrinting!(timer,p,mesh,tsoln)
    while timer.simulationTime < timer.finalTime
        # reaction
        computeReaction!(p,tsoln)
        soln.θs.u = tsoln.θs.u
        soln.θf.u = tsoln.θf.u

        θfn = normalize_thetaf(tsoln.θf.u)
        peAnormalized = p.PeA./(θfn .+ 0.0 .*(1 .- θfn))        # 0.0 -> 0.5
        peBnormalized = p.PeB./(θfn .+ 0.0 .*(1 .- θfn))        # 0.0 -> 0.25
        ψAParam = AdvDiffParam(tsoln.fluid.u,tsoln.fluid.v,peAnormalized)
        ψBParam = AdvDiffParam(tsoln.fluid.u,tsoln.fluid.v,peBnormalized)
        ψCParam = AdvDiffParam(tsoln.fluid.u,tsoln.fluid.v,p.PeC./θfn)

        # advection-diffusion
        M = generateMassMatrixWScalar(mesh.sMesh,θfn)

        StiffA = GenerateSystem(mesh.sMesh,prob.ψA,ψAParam)
        ApplyBC!(StiffA,mesh.sMesh,prob.ψA,ψAParam,prob.ψA.OperatorType)

        StiffB = GenerateSystem(mesh.sMesh,prob.ψB,ψBParam)
        ApplyBC!(StiffB,mesh.sMesh,prob.ψB,ψBParam,prob.ψB.OperatorType)

        StiffC = GenerateSystem(mesh.sMesh,prob.ψC,ψCParam)
        ApplyBC!(StiffC,mesh.sMesh,prob.ψC,ψCParam,prob.ψC.OperatorType)

        if timer.stepIndex < 5
            soln.ψA.u = updateChemical(M,StiffA,tsoln.ψA.u,p;solver=:backwardeuler)
            soln.ψB.u = updateChemical(M,StiffB,tsoln.ψB.u,p;solver=:backwardeuler)
            soln.ψC.u = updateChemical(M,StiffC,tsoln.ψC.u,p;solver=:backwardeuler)
        else
            soln.ψA.u = updateChemical(M,StiffA,tsoln.ψA.u,p;solver=p.integrator)
            soln.ψB.u = updateChemical(M,StiffB,tsoln.ψB.u,p;solver=p.integrator)
            soln.ψC.u = updateChemical(M,StiffC,tsoln.ψC.u,p;solver=p.integrator)
        end

        scrub_scalars!(soln.ψA.u,soln.ψB.u,soln.ψC.u)

        mpbparam = BrinkmanMPParam(1.0,mpbFriction(θfn),θfn)
        sol = solve(prob.fluid,mesh.fMesh,mpbparam)
        soln.fluid.u = sol.u
        soln.fluid.v = sol.v
        soln.fluid.p = sol.p

        tsoln = deepcopy(soln)
        if isPrintIndex(timer)
            doPrinting!(timer,p,mesh,soln)
        end
        updateSimulationTime!(timer)
        updateStepIndex!(timer)
    end
    printSimulationFinishedMessage(timer)
    nothing
end

function scrub_scalars!(u1,u2,u3)
    scrub_scalar!(u1)
    scrub_scalar!(u2)
    scrub_scalar!(u3)
    nothing
end
function scrub_scalar!(u)
    for i=1:length(u); if u[i] < 0; u[i] = 0.0; end; end
    nothing
end

function initialize_data(p::precipParam)
    timer = initializeTimer(p)
    mesh = initialize_meshes(p)
    soln = initializeSolution(p)
    tsoln = initializeSolution(p)
    prob = defineProblems(p,mesh)

    return timer, mesh, soln, tsoln, prob
end

function initialize_meshes(p::precipParam)
    sMesh,fMesh = loadMesh(p)

    return precipMesh(sMesh,fMesh)
    if p.MeshFile == "sin_x2.msh"
        sin_x2_transform!(mesh,1.0)
    elseif p.MeshFile == "sin_wall.msh"
        sin_wall_transform!(mesh,1.0)
    end

    return precipMesh(sMesh,fMesh)
end


text:
using TimerOutputs, Printf

using Oceananigans
using Oceananigans.TurbulenceClosures

include("benchmark_utils.jl")

const timer = TimerOutput()

Ni = 2   # Number of iterations before benchmarking starts.
Nt = 10  # Number of iterations to use for benchmarking time stepping.

# Run benchmark across these parameters.
            Ns = [(32, 32, 32), (64, 64, 64)]
      closures = [ConstantIsotropicDiffusivity, ConstantAnisotropicDiffusivity, ConstantSmagorinsky]
   float_types = [Float32, Float64]     # Float types to benchmark.
         archs = [CPU()]                # Architectures to benchmark on.
@hascuda archs = [CPU(), GPU()]         # Benchmark GPU on systems with CUDA-enabled GPUs.

for arch in archs, float_type in float_types, N in Ns, Closure in closures
    Nx, Ny, Nz = N
    Lx, Ly, Lz = 1, 1, 1

    model = Model(N=(Nx, Ny, Nz), L=(Lx, Ly, Lz), arch=arch, float_type=float_type, closure=Closure(float_type))
    time_step!(model, Ni, 1)

    bn =  benchmark_name(N, string(Closure), arch, float_type; npad=2)
    @printf("Running benchmark: %s...\n", bn)
    for i in 1:Nt
        @timeit timer bn time_step!(model, 1, 1)
    end
end

print_timer(timer, title="Oceananigans.jl turbulence closure benchmarks")
println("")


text:
using MonteCarloMarkovKernels
using Test

import Random.seed!

@testset "AR1 test" begin
  @time include("ar1_test.jl")
end

@testset "RWM test" begin
  @time include("rwmh_test.jl")
end

@testset "AM test" begin
  @time include("am_test.jl")
  @time include("am_test2.jl")
end

@testset "Visualize test" begin
  @time include("visualize_test.jl")
end


text:
using ..Collocation: collocation_points

"""
    ApproxByInterpolation <: AbstractApproxMethod

Approximate function by a spline that passes through the given set of points.

The number of points must be equal to the length of the B-spline basis defining
the approximation space.

See belows for different ways of specifying the interpolation points.

---

    ApproxByInterpolation(xs::AbstractVector)

Specifies an approximation by interpolation at the given points `xs`.

---

    ApproxByInterpolation(B::AbstractBSplineBasis)

Specifies an approximation by interpolation at an automatically-determined set
of points, which are expected to be appropriate for the given basis.

The interpolation points are determined by calling [`collocation_points`](@ref).
"""
struct ApproxByInterpolation{Points <: AbstractVector} <: AbstractApproxMethod
    xs :: Points
end

function Base.show(io::IO, m::ApproxByInterpolation)
    let io = IOContext(io, :compact => true, :limit => true)
        print(io, "interpolation at ", m.xs)
    end
end

ApproxByInterpolation(B::AbstractBSplineBasis) =
    ApproxByInterpolation(collocation_points(B))

function approximate(f, B::AbstractBSplineBasis, m::ApproxByInterpolation)
    S = SplineInterpolation(undef, B, m.xs)
    A = SplineApproximation(m, spline(S), S)
    approximate!(f, A)
end

function _approximate!(f, A, m::ApproxByInterpolation)
    @assert method(A) === m
    S = data(A)
    xs = SplineInterpolations.interpolation_points(S)
    @assert xs === method(A).xs
    ys = coefficients(S)  # just to avoid allocating extra vector
    @inbounds for i in eachindex(xs, ys)
        ys[i] = f(xs[i])
    end
    interpolate!(S, ys)
    A
end

max_stars_repo_name: jipolanco/BasisSplines

max_issues_repo_name: jipolanco/BasisSplines

max_forks_repo_name: jipolanco/BasisSplines


text:
#=
This is a very crude first stab at the Tables.jl interface
https://github.com/JuliaData/Tables.jl
=#

using Tables

Tables.istable(::Type{<:KeyedArray}) = true

Tables.rowaccess(::Type{<:KeyedArray}) = true

function Tables.rows(A::Union{KeyedArray, NdaKa})
    L = hasnames(A) ? (dimnames(A)..., :value) :  # should gensym() if :value in dimnames(A)
        (ntuple(d -> Symbol(:dim_,d), ndims(A))..., :value)
    R = keys_or_axes(A)
    nt(inds) = NamedTuple{L}((map(getindex, R, inds)..., A[inds...]))
    # (nt(inds) for inds in Iterators.product(axes(A)...)) # should flatten?
    (nt(inds) for inds in Vectorator(Iterators.product(axes(A)...)))
end

#=
rr = wrapdims(rand(2,3), 11:12, 21:23)
nn = wrapdims(rand(2,3), a=11:12, b=21:23)

Tables.rows(rr) |> collect |> vec
Tables.rows(nn) |> collect |> vec

Tables.Schema(nn) # define a struct? Now below...

# No error if Tables.rows's generator has size,
# it uses Vectorator mostly to give Tables.Schema something to find.
=#

Tables.columnaccess(::Type{<:KeyedArray{T,N,AT}}) where {T,N,AT} =
    IndexStyle(AT) === IndexLinear()

function Tables.columns(A::Union{KeyedArray, NdaKa})
    L = hasnames(A) ? (dimnames(A)..., :value) :
        (ntuple(d -> Symbol(:dim_,d), ndims(A))..., :value)
    R = keys_or_axes(A)
    G = ntuple(ndims(A)) do d
        vec([rs[d] for rs in Iterators.product(R...)])
        # _vec(rs[d] for rs in Iterators.product(R...))
    end
    C = (G..., vec(parent(A)))
    NamedTuple{L}(C)
end

function Tables.Schema(nt::NamedTuple) # 🏴☠
    L = keys(nt)
    T = map(v -> typeof(first(v)), values(nt))
    Tables.Schema(L,T)
end

#=
Ah, iterators aren't allowed for columns, must be indexable:
https://github.com/JuliaData/Tables.jl/issues/101
They could be something like this, but seems overkill?
https://github.com/MichielStock/Kronecker.jl
https://github.com/JuliaArrays/LazyArrays.jl#kronecker-products

Tables.columns(nn)
map(collect, Tables.columns(nn))

using DataFrames

DataFrame(rand(2,3))
DataFrame(nn) # doesn't see Tables

dd1 = DataFrame(Tables.rows(nn))
dd2 = DataFrame(Tables.columns(nn))
=#

"""
    Vectorator(iter)

Wrapper for iterators which ensures they do not have an n-dimensional size.
Tries to ensure that `collect(Vectorator(iter)) == vec(collect(iter))`.
"""
struct Vectorator{T}
    iter::T
end

_vec(iter) = (x for x in Vectorator(iter))

Base.iterate(x::Vectorator, s...) = iterate(x.iter, s...)

Base.length(x::Vectorator) = length(x.iter)

Base.IteratorSize(::Type{Vectorator{T}}) where {T} =
    Base.IteratorSize(T) isa Base.HasShape ? Base.HasLength() : IteratorSize(T)

Base.IteratorEltype(::Type{Vectorator{T}}) where {T} = Base.IteratorEltype(T)

Base.eltype(::Type{Vectorator{T}}) where {T} = eltype(T)

function Tables.Schema(rows::Base.Generator{<:Vectorator})
    row = first(rows)
    Tables.Schema(keys(row), map(typeof, values(row)))
end

# struct OneKron{T, AT} <: AbstractVector{T}
#     data::AT
#     inner::Int
#     outer::Int
# end

# Tables.materializer(A::KeyedArray) = wrapdims

# function wrapdims(tab)
#     sch = Tables.Schema(tab)
#     for r in Tables.rows(tab)

#     end
# end


text:
using LightGraphs
using LightGraphs.SimpleGraphs
using LightGraphs.Experimental
using Test
using SparseArrays
using LinearAlgebra
using DelimitedFiles
using Base64
using Random
using Statistics: mean

const testdir = dirname(@__FILE__)

testgraphs(g) = [g, Graph{UInt8}(g), Graph{Int16}(g)]
testdigraphs(g) = [g, DiGraph{UInt8}(g), DiGraph{Int16}(g)]

# some operations will create a large graph from two smaller graphs. We
# might error out on very small eltypes.
testlargegraphs(g) = [g, Graph{UInt16}(g), Graph{Int32}(g)]
testlargedigraphs(g) = [g, DiGraph{UInt16}(g), DiGraph{Int32}(g)]

tests = [
    "simplegraphs/runtests",
    "linalg/runtests",
    "parallel/runtests",
    "interface",
    "core",
    "operators",
    "degeneracy",
    "distance",
    "digraph/transitivity",
    "cycles/hawick-james",
    "cycles/johnson",
    "cycles/karp",
    "cycles/basis",
    "cycles/limited_length",
    "edit_distance",
    "connectivity",
    "persistence/persistence",
    "shortestpaths/astar",
    "shortestpaths/bellman-ford",
    "shortestpaths/dijkstra",
    "shortestpaths/johnson",
    "shortestpaths/floyd-warshall",
    "shortestpaths/yen",
    "traversals/bfs",
    "traversals/bipartition",
    "traversals/greedy_color",
    "traversals/dfs",
    "traversals/maxadjvisit",
    "traversals/randomwalks",
    "traversals/diffusion",
    "community/cliques",
    "community/core-periphery",
    "community/label_propagation",
    "community/modularity",
    "community/clustering",
    "community/clique_percolation",
    "centrality/betweenness",
    "centrality/closeness",
    "centrality/degree",
    "centrality/katz",
    "centrality/pagerank",
    "centrality/eigenvector",
    "centrality/stress",
    "centrality/radiality",
    "utils",
    "spanningtrees/kruskal",
    "spanningtrees/prim",
    "biconnectivity/articulation",
    "biconnectivity/biconnect",
    "biconnectivity/bridge",
    "graphcut/normalized_cut",
    "graphcut/karger_min_cut",
    "dominatingset/degree_dom_set",
    "dominatingset/minimal_dom_set",
    "independentset/degree_ind_set",
    "independentset/maximal_ind_set",
    "vertexcover/degree_vertex_cover",
    "vertexcover/random_vertex_cover",
    "experimental/experimental"
]

@testset "LightGraphs" begin
    for t in tests
        tp = joinpath(testdir, "$(t).jl")
        include(tp)
    end
end