text:
using DynamicGrids, ModelParameters, Setfield, Test, StaticArrays,
      LinearAlgebra, CUDAKernels
import DynamicGrids: applyrule, applyrule!, maprule!, ruletype, extent, source, dest,
       _getreadgrids, _getwritegrids, _combinegrids, _readkeys, _writekeys,
       SimData, WritableGridData, Rule, Extent, CPUGPU

if CUDAKernels.CUDA.has_cuda_gpu()
    CUDAKernels.CUDA.allowscalar(false)
    # hardware = (SingleCPU(), ThreadedCPU(), CPUGPU(), CuGPU())
    hardware = (SingleCPU(), ThreadedCPU(), CPUGPU())
else
    hardware = (SingleCPU(), ThreadedCPU(), CPUGPU())
end

init  = [0 1 1 0
         0 1 1 0
         0 1 1 0
         0 1 1 0
         0 1 1 0]

proc = SingleCPU()
opt = NoOpt()

@testset "Generic rule constructors" begin
   rule1 = Cell{:a}(identity)
   @test rule1.f == identity
   @test_throws ArgumentError Cell()
   rule2 = Cell{:a,:a}(identity)
   @test rule1 == rule2
   #@test_throws ArgumentError Cell(identity, identity)
   rule1 = Neighbors{:a,:b}(identity, Moore(1))
   @test rule1.f == identity
   rule2 = Neighbors{:a,:b}(identity; neighborhood=Moore(1))
   @test rule1 == rule2
   @test typeof(rule1)  == Neighbors{:a,:b,typeof(identity),Moore{1,8,Nothing}}
   rule1 = Neighbors(identity, Moore(1))
   @test rule1.f == identity
   rule2 = Neighbors(identity; neighborhood=Moore(1))
   @test typeof(rule1)  == Neighbors{:_default_,:_default_,typeof(identity),Moore{1,8,Nothing}}
   @test rule1 == rule2
   @test_throws ArgumentError Neighbors()
   # @test_throws ArgumentError Neighbors(identity, identity, identity)
   rule1 = SetNeighbors{:a,:b}(identity, Moore(1))
   @test rule1.f == identity
   rule2 = SetNeighbors{:a,:b}(identity; neighborhood=Moore(1))
   @test rule1 == rule2
   @test typeof(rule1)  == SetNeighbors{:a,:b,typeof(identity),Moore{1,8,Nothing}}
   rule1 = SetNeighbors(identity, Moore(1))
   @test rule1.f == identity
   rule2 = SetNeighbors(identity; neighborhood=Moore(1))
   @test typeof(rule1)  == SetNeighbors{:_default_,:_default_,typeof(identity),Moore{1,8,Nothing}}
   @test rule1 == rule2
   @test_throws ArgumentError Neighbors()
   # @test_throws ArgumentError Neighbors(identity, identity, identity)
   rule1 = SetCell{:a,:b}(identity)
   @test rule1.f == identity
   @test_throws ArgumentError SetCell()
   # @test_throws ArgumentError SetCell(identity, identity)
end

@testset "Rulesets" begin
    rule1 = Cell(x -> 2x)
    rule2 = Cell(x -> 3x)
    rs1 = Ruleset((rule1, rule2); opt=NoOpt())
    rs2 = Ruleset(rule1, rule2; opt=NoOpt())
    @test rules(rs1) == rules(rs2)
    @test typeof(Ruleset(StaticRuleset(rs1))) == typeof(rs1)
    for fn in fieldnames(typeof(rs1))
        @test getfield(Ruleset(StaticRuleset(rs1)), fn) == getfield(rs1, fn)
    end
    @test typeof(Ruleset(StaticRuleset(rs1))) == typeof(rs1)
    ModelParameters.setparent!(rs2, (rule1,))
    @test rs2.rules == (rule1,)
end

@testset "Cell" begin
    rule = Cell((d, x, I) -> 2x)
    @test applyrule(nothing, rule, 1, (0, 0)) == 2
end

@testset "Neighbors" begin
    buf = [1 0 0; 0 0 1; 0 0 1]
    rule = Neighbors(VonNeumann(1, buf)) do data, hood, state, I
        sum(hood)
    end
    @test applyrule(nothing, rule, 0, (3, 3)) == 1
    rule = Neighbors(Moore{1}(buf)) do data, hood, state, I
        sum(hood)
    end
    @test applyrule(nothing, rule, 0, (3, 3)) == 3
    @test DynamicGrids._buffer(rule) === buf
end

@testset "Convolution" begin
    k = SA[1 0 1; 0 0 0; 1 0 1]
    @test Convolution{:a}(k) == Convolution{:a,:a}(; neighborhood=Kernel(Window(1), k))
    buf = SA[1 0 0; 0 0 1; 0 0 1]
    hood = Window{1,9,typeof(buf)}(buf)
    rule = Convolution{:a,:a}(; neighborhood=Kernel(hood, k))
    @test DynamicGrids.kernel(rule) === k
    @test applyrule(nothing, rule, 0, (3, 3)) == k ⋅ buf
    output = ArrayOutput((a=init,); tspan=1:2)
    sim!(output, rule)
end

@testset "SetNeighbors" begin
    init  = [0 1 0 0
             0 0 0 0
             0 0 0 0
             0 1 0 0
             0 0 1 0]
    @test_throws ArgumentError SetNeighbors()
    @testset "atomics" begin
        rule = SetNeighbors(VonNeumann(1)) do data, hood, state, I
            if state > 0
                for pos in positions(hood, I)
                    add!(data, 1, pos...)
                end
            end
        end
        output = ArrayOutput(init; tspan=1:2)
        for proc in hardware, opt in (NoOpt(), SparseOpt())
            @testset "$(nameof(typeof(opt))) $(nameof(typeof(proc)))" begin
                ref_output = [1 1 1 0
                              0 1 0 0
                              0 1 0 0
                              1 1 2 0
                              0 2 1 1]
                sim!(output, rule, proc=proc, opt=opt)
                @test output[2] == ref_output
            end
        end
    end

    @testset "setindex" begin
        rule = SetNeighbors(VonNeumann(1)) do data, hood, state, I
            state == 0 && return nothing
            for pos in positions(hood, I)
                data[pos...] = 1
            end
        end
        output = ArrayOutput(init; tspan=1:2)
        for proc in hardware, opt in (NoOpt(), SparseOpt())
            @testset "$(nameof(typeof(opt))) $(nameof(typeof(proc)))" begin
                ref_output = [1 1 1 0
                              0 1 0 0
                              0 1 0 0
                              1 1 1 0
                              0 1 1 1]
                sim!(output, rule, proc=proc, opt=opt)
                @test output[2] == ref_output
            end
        end
    end

end

@testset "SetCell" begin
    init  = [0 1 0 0
             0 0 0 0
             0 0 0 0
             0 1 0 0
             0 0 1 0]
    @testset "add!" begin
        output = ArrayOutput(init; tspan=1:2)
        rule = SetCell() do data, state, I
            if state > 0
                pos = I[1] - 2, I[2]
                isinbounds(data, pos) && add!(first(data), 1, pos...)
            end
        end
        for proc in hardware, opt in (NoOpt(), SparseOpt())
            @testset "$(nameof(typeof(opt))) $(nameof(typeof(proc)))" begin
                sim!(output, rule; proc=proc, opt=opt)
                ref_out = [0 1 0 0
                           0 1 0 0
                           0 0 1 0
                           0 1 0 0
                           0 0 1 0]
                @test output[2] == ref_out
            end
        end
    end
    @testset "setindex!" begin
        output = ArrayOutput(init; tspan=1:2)
        rule = SetCell() do data, state, I
            if state > 0
                pos = I[1] - 2, I[2]
                isinbounds(data, pos) && (data[pos...] = 5)
            end
        end
        sim!(output, rule)
        @test output[2] == [0 1 0 0
                            0 5 0 0
                            0 0 5 0
                            0 1 0 0
                            0 0 1 0]
    end
end

@testset "SetGrid" begin
    @test_throws ArgumentError SetGrid()
    rule = SetGrid() do r, w
        w .*= 2
    end
    init  = [0 1 0 0
             0 0 0 0
             0 0 0 0
             0 1 0 0
             0 0 1 0]
    output = ArrayOutput(init; tspan=1:2)
    for proc in hardware, opt in (NoOpt(), SparseOpt())
        @testset "$(nameof(typeof(opt))) $(nameof(typeof(proc)))" begin
            sim!(output, rule; proc=proc, opt=opt)
            @test output[2] == [0 2 0 0
                                0 0 0 0
                                0 0 0 0
                                0 2 0 0
                                0 0 2 0]
        end
    end
end

struct AddOneRule{R,W} <: CellRule{R,W} end
DynamicGrids.applyrule(data, ::AddOneRule, state, args...) = state + 1

@testset "Rulset mask ignores false cells" begin
    init = [0.0 4.0 0.0
            0.0 5.0 8.0
            3.0 6.0 0.0]
    mask = Bool[0 1 0
                0 1 1
                1 1 0]
    ruleset1 = Ruleset(AddOneRule{:_default_,:_default_}(); opt=NoOpt())
    ruleset2 = Ruleset(AddOneRule{:_default_,:_default_}(); opt=SparseOpt())
    output1 = ArrayOutput(init; tspan=1:3, mask=mask)
    output2 = ArrayOutput(init; tspan=1:3, mask=mask)
    sim!(output1, ruleset1)
    sim!(output2, ruleset2)
    @test output1[1] == output2[1] == [0.0 4.0 0.0
                                       0.0 5.0 8.0
                                       3.0 6.0 0.0]
    @test output1[2] == output2[2] == [0.0 5.0 0.0
                                       0.0 6.0 9.0
                                       4.0 7.0 0.0]
    @test output1[3] == output2[3] == [0.0 6.0 0.0
                                       0.0 7.0 10.0
                                       5.0 8.0 0.0]
end

# Single grid rules

struct TestRule{R,W} <: CellRule{R,W} end
applyrule(data, ::TestRule, state, index) = 0

@testset "A rule that returns zero gives zero outputs" begin
    final = [0 0 0 0
             0 0 0 0
             0 0 0 0
             0 0 0 0
             0 0 0 0]
    mask = nothing
    rule = TestRule{:a,:a}()

    for proc in (SingleCPU(), ThreadedCPU()), opt in (NoOpt(), SparseOpt())
        @testset "$(nameof(typeof(opt))) $(nameof(typeof(proc)))" begin
            ruleset = Ruleset(rule; opt=opt, proc=proc)
            @test DynamicGrids.boundary(ruleset) === Remove()
            @test DynamicGrids.opt(ruleset) === opt
            @test DynamicGrids.proc(ruleset) === proc
            @test DynamicGrids.cellsize(ruleset) === 1
            @test DynamicGrids.timestep(ruleset) === nothing
            @test DynamicGrids.ruleset(ruleset) === ruleset

            ext = Extent(; init=(a=init,), tspan=1:1)
            simdata = SimData(ext, ruleset)

            # Test maprules components
            rkeys, rgrids = _getreadgrids(rule, simdata)
            wkeys, wgrids = _getwritegrids(rule, simdata)
            @test rkeys == Val{:a}()
            @test wkeys == Val{:a}()
            newsimdata = @set simdata.grids = _combinegrids(rkeys, rgrids, wkeys, wgrids)
            @test newsimdata.grids[1] isa WritableGridData
            # Test type stability
            T = Val{DynamicGrids.ruletype(rule)}()
            @inferred maprule!(newsimdata, proc, opt, T, rule, rkeys, wkeys)

            resultdata = maprule!(simdata, rule)
            @test source(resultdata[:a]) == final
        end
    end
end

struct TestSetCell{R,W} <: SetCellRule{R,W} end
applyrule!(data, ::TestSetCell, state, index) = 0

@testset "A partial rule that returns zero does nothing" begin
    rule = TestSetCell()
    mask = nothing
    for proc in (SingleCPU(), ThreadedCPU()), opt in (NoOpt(), SparseOpt())
        @testset "$(nameof(typeof(opt))) $(nameof(typeof(proc)))" begin
            ruleset = Ruleset(rule; opt=NoOpt())
            # Test type stability
            ext = Extent(; init=(_default_=init,), tspan=1:1)
            simdata = SimData(ext, ruleset)
            rkeys, rgrids = _getreadgrids(rule, simdata)
            wkeys, wgrids = _getwritegrids(rule, simdata)
            newsimdata = @set simdata.grids = _combinegrids(wkeys, wgrids, rkeys, rgrids)
            T = Val{DynamicGrids.ruletype(rule)}()
            @inferred maprule!(newsimdata, proc, opt, T, rule, rkeys, wkeys)
            resultdata = maprule!(simdata, rule)
            @test source(resultdata[:_default_]) == init
        end
    end
end

struct TestSetCellWrite{R,W} <: SetCellRule{R,W} end
applyrule!(data, ::TestSetCellWrite{R,W}, state, index) where {R,W} = add!(data[W], 1, index[1], 2)

@testset "A partial rule that writes to dest affects output" begin
    init  = [0 1 1 0
             0 1 1 0
             0 1 1 0
             0 1 1 0
             0 1 1 0]
    final = [0 5 1 0;
             0 5 1 0;
             0 5 1 0;
             0 5 1 0;
             0 5 1 0]
    for proc in hardware, opt in (NoOpt(), SparseOpt())
        @testset "$(nameof(typeof(opt))) $(nameof(typeof(proc)))" begin
            rule = TestSetCellWrite()
            ruleset = Ruleset(rule; opt=opt, proc=proc)
            ext = Extent(; init=(_default_=init,), tspan=1:1)
            simdata = DynamicGrids._proc_setup(SimData(ext, ruleset));
            resultdata = maprule!(simdata, rule);
            @test Array(source(first(resultdata))) == final
        end
    end
end

struct TestCellTriple{R,W} <: CellRule{R,W} end
applyrule(data, ::TestCellTriple, state, index) = 3state

struct TestCellSquare{R,W} <: CellRule{R,W} end
applyrule(data, ::TestCellSquare, (state,), index) = state^2

@testset "Chained cell rules work" begin
    init  = [0 1 2 3;
             4 5 6 7]
    final = [0 9 36 81;
             144 225 324 441]
    for proc in hardware, opt in (NoOpt(), SparseOpt())
        @testset "$(nameof(typeof(opt))) $(nameof(typeof(proc)))" begin
            rule = Chain(TestCellTriple(), TestCellSquare())
            ruleset = Ruleset(rule; opt=opt, proc=proc)
            ext = Extent(; init=(_default_=init,), tspan=1:1)
            simdata = DynamicGrids._proc_setup(SimData(ext, ruleset))
            resultdata = maprule!(simdata, rule);
            @test Array(source(first(resultdata))) == final
        end
    end
end

struct PrecalcRule{R,W,P} <: CellRule{R,W}
    precalc::P
end
DynamicGrids.modifyrule(rule::PrecalcRule, simdata) = PrecalcRule(currenttime(simdata))
applyrule(data, rule::PrecalcRule, state, index) = rule.precalc[]

@testset "Rule precalculations work" begin
    init  = [1 1;
             1 1]
    out2  = [2 2;
             2 2]
    out3  = [3 3;
             3 3]
    for proc in hardware, opt in (NoOpt(), SparseOpt())
        @testset "$(nameof(typeof(opt))) $(nameof(typeof(proc)))" begin
            rule = PrecalcRule(1)
            ruleset = Ruleset(rule; proc=proc, opt=opt)
            output = ArrayOutput(init; tspan=1:3)
            sim!(output, ruleset)
            # Copy for GPU
            @test output[2] == out2
            @test output[3] == out3
        end
    end
end

# Multi grid rules

struct DoubleY{R,W} <: CellRule{R,W} end
applyrule(data, rule::DoubleY, (x, y), index) = y * 2

struct HalfX{R,W} <: CellRule{R,W} end
applyrule(data, rule::HalfX, x, index) = x, x * 0.5

struct Predation{R,W} <: CellRule{R,W} end
Predation(; prey=:prey, predator=:predator) =
    Predation{Tuple{predator,prey},Tuple{prey,predator}}()
applyrule(data, ::Predation, (predators, prey), index) = begin
    caught = 2predators
    # Output order is the reverse of input to test that can work
    prey - caught, predators + caught * 0.5
end

predation = Predation(; prey=:prey, predator=:predator)

@testset "Multi-grid keys are inferred" begin
    @test _writekeys(predation) == (:prey, :predator)
    @test _readkeys(predation) == (:predator, :prey)
    @test keys(predation) == (:prey, :predator)
    @inferred _writekeys(predation)
    @inferred _readkeys(predation)
    @inferred keys(predation)
end

@testset "Multi-grid rules work" begin
    init = (prey=[10. 10.], predator=[1. 0.])
    rules = DoubleY{Tuple{:predator,:prey},:prey}(), predation
    output = ArrayOutput(init; tspan=1:3)
    for proc in hardware, opt in (NoOpt(), SparseOpt())
        @testset "$(nameof(typeof(opt))) $(nameof(typeof(proc)))" begin
            sim!(output, rules; opt=opt, proc=proc)
            @test output[2] == (prey=[18. 20.], predator=[2. 0.])
            @test output[3] == (prey=[32. 40.], predator=[4. 0.])
        end
    end
end

@testset "Multi-grid rules work" begin
    init = (prey=[10. 10.], predator=[0. 0.])
    output = ArrayOutput(init; tspan=1:3)
    for proc in hardware, opt in (NoOpt(), SparseOpt())
        @testset "$(nameof(typeof(opt))) $(nameof(typeof(proc)))" begin
            ruleset = Ruleset((HalfX{:prey,Tuple{:prey,:predator}}(),); opt=opt, proc=proc)
            sim!(output, ruleset)
            @test output[2] == (prey=[10. 10.], predator=[5. 5.])
            @test output[3] == (prey=[10. 10.], predator=[5. 5.])
        end
    end
end

@testset "life with generic constructors" begin
    @test Life(Moore(1), (1, 1), (5, 5)) ==
        Life(; neighborhood=Moore(1), born=(1, 1), survive=(5, 5))
    @test Life{:a,:b}(Moore(1), (7, 1), (5, 3)) ==
          Life{:a,:b}(neighborhood=Moore(1), born=(7, 1), survive=(5, 3))
    # Defaults
    @test Life() == Life(
        Moore(1),
        Param(3, bounds=(0, 8)),
        (Param(2, bounds=(0, 8)), Param(3, bounds=(0, 8)))
    )
    @test Life{:a,:b}() == Life{:a,:b}(
         Moore(1),
         Param(3, bounds=(0, 8)),
         (Param(2, bounds=(0, 8)), Param(3, bounds=(0, 8)))
    )
end

@testset "generic ConstructionBase compatability" begin
    life = Life{:x,:y}(; neighborhood=Moore(2), born=(1, 1), survive=(2, 2))
    @set! life.born = (5, 6)
    @test life.born == (5, 6)
    @test life.survive == (2, 2)
    @test _readkeys(life) == :x
    @test _writekeys(life) == :y
    @test DynamicGrids.neighborhood(life) == Moore(2)
end


text:
########################################
## File Name: sac_control_test.jl
## Author: Haruki Nishimura (hnishimura@stanford.edu)
## Date Created: 2020/05/20
## Description: Test code for src/sac_controller.jl
########################################

using LinearAlgebra
using Random

@testset "SAC Control Test" begin

# Localization Task
dtc = 0.01;
dto = 0.2;
dtexec = dto*0.8;
Q = Matrix(1.0I, 2, 2);
Cu = Matrix(1.0I, 2, 2) #+ [0. 0.01; 0.01 0.];
transPos = TransModel_Pos();
ukfModel = UKFPosRange();
simPosUKF = SimulatePosRangeLocalization2D(transPos,ukfModel,dtc,dto,dtexec,Q,Cu);

rng = MersenneTwister(1234)
p = PhysPos(0.,zeros(2));
bVec = VecBelMvNormal([BelMvNormal(0.,rand(rng, 2),10*Matrix(1.0I, 2, 2)) for ii = 1:100]);
s = AugState(p,bVec);

UArray = [PosControl(round(dtc*t,digits=5),zeros(2)) for t = 0:Int64(2.0/dtc)-1];

coeffUKF = getControlCoeffs(simPosUKF,s,UArray,MersenneTwister(1234))

@test length(coeffUKF) == 201
@test isapprox(coeffUKF[1], [2298.44, 2561.91], rtol=1e-5)
@test isapprox(coeffUKF[end-1], [270.112, 308.623], rtol=1e-5)
@test coeffUKF[end] == [0.0, 0.0]

UCMatPosUKF = controlCoeffsExpected(simPosUKF,s,UArray,100,MersenneTwister(1234))

@test size(UCMatPosUKF) == (2, 201)
@test isapprox(UCMatPosUKF[1, 1], 2441.96, rtol=1e-5)
@test isapprox(UCMatPosUKF[2, 1], 2505.37, rtol=1e-5)
@test UCMatPosUKF[1, end] == 0.0
@test UCMatPosUKF[2, end] == 0.0
@test isapprox(UCMatPosUKF[1, end-1], 270.42, rtol=1e-5)
@test isapprox(UCMatPosUKF[2, end-1], 285.107, rtol=1e-5)

UOptArray, CostArray = optControlSchedule(simPosUKF, UArray, UCMatPosUKF, [-1., -1.], [1., 1.]);

@test all([u.vel == [-1., -1.] for u in UOptArray])
@test isapprox(CostArray[1], -4946.33, rtol=1e-5)
@test isapprox(CostArray[end], -554.527, rtol=1e-5)

UOpt,tcalc = determineControlTime(0.1,simPosUKF,UOptArray,CostArray)

@test UOpt.t == 0.27
@test UOpt.vel == [-1., -1.]

UArrayNew, act_time_init, act_time_final, tcalc = sacControlUpdate(simPosUKF,s,UArray,[-1.,-1.],[1.,1.],10,MersenneTwister(1234))

#@test act_time_init == 0.06
#@test act_time_final == 0.21
#@test UArrayNew[7].vel == [-1.0, -1.0]

# Manipulation Task
dtc = 0.02;
dto = 0.2;
dtexec = dto*0.8;
QMP = [Matrix(1.0I, 6, 6) zeros(6,5); zeros(5,6) 0.00*Matrix(1.0I, 5, 5)];
RMP = Matrix(1.0I, 9, 9);
CsMP = 10*QMP;
CuMP = 10*Matrix(1.0I, 3, 3) #+ [0. 0. 0.01; 0. 0. 0.; 0.01 0. 0.];
simMP = SimulateManipulate2D(dtc,dto,dtexec,QMP,RMP,CsMP,CuMP)

xMP = PhysManipulate2D(0.,[0.1,0.2],0.5,[0.3,0.4],0.6,3.,5.,[10.,11.],0.2);
uMP = MControl2D(0.,6.,7.,8.);
bMP = BelMvNormal(xMP.t,vec(xMP),diagm([10.,10.,10.,10.,10.,0.1,1.,1.,10.,10.,0.1]));

UArray = [MControl2D(round(dtc*t,digits=5),0.,0.0,0.0) for t = 0:Int64(2.0/dtc)-1];

cMP = getControlCoeffs(simMP,bMP,UArray,MersenneTwister(1234))

@test length(cMP) == 101
@test isapprox(cMP[1], [1929.07, -483.397, -134.458], rtol=1e-5)
@test isapprox(cMP[end], [-31.1805, 3.37386, -1.69659], rtol=1e-5)

cMPExpected = controlCoeffsExpected(simMP,bMP,UArray,10,MersenneTwister(1))

@test size(cMPExpected) == (3, 101)
@test isapprox(cMPExpected[1, 1], 6532.73, rtol=1e-5)
@test isapprox(cMPExpected[2, 1], -1597.58, rtol=1e-5)
@test isapprox(cMPExpected[3, 1], -453.342, rtol=1e-5)
@test isapprox(cMPExpected[1, end], 31.7051, rtol=1e-5)
@test isapprox(cMPExpected[2, end], 26.0962, rtol=1e-5)
@test isapprox(cMPExpected[3, end], 3.52962, rtol=1e-5)
@test isapprox(cMPExpected[1, end-1], 38.9432, rtol=1e-5)
@test isapprox(cMPExpected[2, end-1], 242.227, rtol=1e-5)
@test isapprox(cMPExpected[3, end-1], -8.07914, rtol=1e-5)

UOptArray, CostArray = optControlSchedule(simMP,UArray,cMPExpected, [-1.,-1.,-1.],[1.,1.,1])

@test length(UOptArray) == 100
@test isapprox(vec(UOptArray[1]), [-1.0, 1.0, 1.0], rtol=1e-5)
@test isapprox(vec(UOptArray[end]), [-1.0, -1.0, 0.807914], rtol=1e-5)

@test isapprox(CostArray[1], -8568.65, rtol=1e-5)
@test isapprox(CostArray[end], -274.434, rtol=1e-5)

UOpt,tcalc = determineControlTime(0.1,simMP,UOptArray,CostArray)

@test UOpt.t == 0.28
@test isapprox(vec(UOpt), [-1.0, -1.0, -0.9884261676986377], rtol=1e-5)

UArrayNew, act_time_init, act_time_final, tcalc = sacControlUpdate(simMP,bMP,UArray,[-1.,-1.,-1.],[1.,1.,1.],1,MersenneTwister(1234))

#@test act_time_init == 0.06
#@test act_time_final == 0.2

#@test isapprox(vec(UArrayNew[4]),  [-0.0158972, 0.544518, 0.0988366], rtol=1e-5)

end

ext: jl


text:
function update(fluid::WorkingFluid, inputpair::InputPair, input1::Real, input2::Real)
    CoolProp.AbstractState_update(fluid.handle, inputpair.v, input1, input2)
end

function keyed_output(fluid::WorkingFluid, output::Property)
    CoolProp.AbstractState_keyed_output(fluid.handle, output.v)
end

function first_partial_deriv(fluid::WorkingFluid,
                             of::Property,
                             with_respect_to::Property,
                             constant::Property )
    CoolProp.AbstractState_first_partial_deriv(
        fluid.handle,
        of.v,
        with_respect_to.v,
        constant.v)
end

function free(fluid::WorkingFluid)
    CoolProp.AbstractState_free(fluid.handle)
end


text:
# 1D model problem to test BSE code

# single particle code
# single particle problem
using LinearAlgebra, SparseArrays, ProgressMeter

#TODO: remove this type? just use it in setup of the bse problem?
abstract type SPProblem end

struct SPProblem1D <: SPProblem
    V_sp::Function
    l::Float64
    r_unit::Vector{Float64}
    r_super::Vector{Float64}
    k_bz::Vector{Float64}
end

struct SPSolution1D
    prob
    ev
    ef
end

function SPProblem1D(V_sp, l, N_unit, N_k)
    r_unit = range(0, stop = l, length = N_unit + 1)[1:(end-1)]
    r_super = range(-div(N_k, 2) * l, stop = div(N_k + 1, 2) * l, length = N_unit * N_k + 1)[1:(end-1)]
    k_bz = range(0.0, stop = 1.0, length = N_k + 1)[1:(end - 1)]

    return SPProblem1D(V_sp, l, r_unit, r_super, k_bz)
end

"""
function to compute the modified momentum operator \$-i \nabla + k\$ in 1D
used to assemble k-dependent Hamiltonians; uses periodic boundary conditions where the function values at r[1] and r[N + 1] are identified

input
r: vector; containing a 1D space grid with N + 1 points
k: number;
bc: 1 for periodic boundary conditions, 0 for homogeneous dirichlet boundary conditions

output
sparse matrix of size N x N
"""
function modified_momentum_operator(r, k) # with periodic boundary conditions
    N = length(r) - 1

    Δr_inv = 1 ./ vcat(diff(r))
    nabla = spdiagm(0 => -Δr_inv, 1 => Δr_inv[1:(end - 1)])
    nabla[end, 1] = Δr_inv[end]
    return -im * nabla + k * I
end

function modified_momentum_operator(r, l, k) # with periodic boundary conditions
    return modified_momentum_operator(vcat(r, l), k)
end

function potential_energy(r, V_fun)
    V = spdiagm(0 => V_fun.(r))
    return V
end

"""
assemble single particle hamiltonian with crystal momentum k
"""
function single_particle_hamiltonian(r, l, k, V_fun)
    momentum_k = modified_momentum_operator(r, l, k)
    H_k = 1 / 2 * momentum_k * momentum_k' + potential_energy(r, V_fun)
end

"""
solve single particle problem
"""
function solve(prob::SPProblem1D)
    V_sp = prob.V_sp
    l = prob.l
    r = prob.r_unit
    kc = 2 * pi / l * (mod.(prob.k_bz .+ 0.5, 1.0) .- 0.5)

    ev = []
    ef = []
    for k in kc
        H_k = Hermitian(Matrix(single_particle_hamiltonian(r, l, k, V_sp)))
        sol = eigen(H_k)
        push!(ev, sol.values)
        push!(ef, sol.vectors)
    end

    return SPSolution1D(prob, ev, ef)
end

# BSE problem

struct BSEProblem1D <: AbstractBSEProblem
    prob::SPProblem1D
    N_core::Int
    N_v::Int
    N_c::Int
    N_k::Int
    E_v::Array{Float64, 2}
    E_c::Array{Float64, 2}
    u_v::Array{Complex{Float64}, 3}
    u_c::Array{Complex{Float64}, 3}
    v_hat::Array{Complex{Float64}, 2}
    w_hat::Array{Complex{Float64}, 3}
end

function BSEProblem1D(sp_prob::SPProblem1D, N_core, N_v, N_c, V, W)
    sp_sol = solve(sp_prob)
    return BSEProblem1D(sp_sol, N_core, N_v, N_c, V, W)
end

function BSEProblem1D(sp_sol::SPSolution1D, N_core, N_v, N_c, V, W)
    sp_prob = sp_sol.prob
    l = sp_prob.l
    N_unit = length(sp_prob.r_unit)
    N_k = length(sp_prob.k_bz)

    E_v = [sp_sol.ev[ik][N_core + iv] for iv in 1:N_v, ik in 1:N_k]
    E_c = [sp_sol.ev[ik][N_core + N_v + ic] for ic in 1:N_c, ik in 1:N_k]
    u_v = [sqrt(N_unit / l) * sp_sol.ef[ik][ir, N_core + iv] for ir in 1:N_unit, iv in 1:N_v, ik in 1:N_k]
    u_c = [sqrt(N_unit / l) * sp_sol.ef[ik][ir, N_core + N_v + ic] for ir in 1:N_unit, ic in 1:N_c, ik in 1:N_k]

    # fourier representations of V and W
    v_hat = compute_hat_Gq(V, sp_prob.r_super, sp_prob.r_unit)
    w_hat = compute_hat_GGq(W, sp_prob.r_super, sp_prob.r_unit)

    return BSEProblem1D(sp_prob, N_core, N_v, N_c, N_k, E_v, E_c, u_v, u_c, v_hat, w_hat)
end

function compute_hat_Gq(V, r_super, r_unit)
    N_unit = length(r_unit)
    N_k = div(length(r_super), N_unit)
    Δr = r_unit[2] - r_unit[1]
    l = N_unit * Δr

    v = V.(r_super, 0.)
    v_hat_vec = (-1).^(0:(length(r_super) - 1)) .* fft(v) * (l / N_unit)
    v_hat = ifftshift(transpose(reshape(circshift(v_hat_vec, div(N_k, 2)), N_k, N_unit)), 2)

    return v_hat
end

function compute_hat_GGq(W, r_super, r_unit)
    N_unit = length(r_unit)
    N_k = div(length(r_super), N_unit)
    Δr = r_unit[2] - r_unit[1]
    l = N_unit * Δr
    L = N_k * l

    w = W.(r_super .+ r_unit', r_unit', l, L)
    w_hat_vec = (-1).^((0:(N_k * N_unit - 1))) .* fft(w) * (l / N_unit^2);

    A = (iG, iGp, iq) -> begin
        if iq > div(N_k, 2)
            iq_res = iq - N_k
        else
            iq_res = iq
        end
        if iG > div(N_unit, 2)
            iG_res = iG - N_unit
        else
            iG_res = iG
        end
        if iGp > div(N_unit, 2)
            iGp_res = iGp - N_unit
        else
            iGp_res = iGp
        end
        iqG_r = mod1(iq_res + N_k * (iG_res - 1), N_k * N_unit)
        iG_c = mod1(iG_res - iGp_res + 1, N_unit)

        return (iqG_r, iG_c)
    end

    w_hat = [w_hat_vec[A(iG, iGp, iq)...]
             for iG in 1:N_unit, iGp in 1:N_unit, iq in 1:N_k]

    return w_hat
end

function supercell_difference(r_1, r_2, L)
    r_1 - r_2 - round((r_1 - r_2) / L) * L
end

#

function size(prob::BSEProblem1D)
    return (prob.N_v, prob.N_c, prob.N_k)
end

function size_k(prob::BSEProblem1D)
    return (prob.N_k, 1, 1)
end

function size_r(prob::BSEProblem1D)
    return (length(prob.prob.r_unit), 1, 1)
end

function energies(prob::BSEProblem1D)
    return prob.E_v, prob.E_c
end

function orbitals(prob::BSEProblem1D)
    return prob.u_v, prob.u_c
end

function lattice_matrix(prob::BSEProblem1D)
    return [prob.prob.l 0 0; 0 1 0; 0 0 1]
end

function compute_v_hat(prob::BSEProblem1D)
    return prob.v_hat[:, 1]
end

function compute_w_hat(prob::BSEProblem1D)
    N_k = prob.N_k
    w_hat = [(prob.w_hat[:, :, iq], (0:127)' .* [1, 0, 0]) for iq in 1:N_k]
    q_2bz_ind = [1:N_k; 1:N_k]
    q_2bz_shift = hcat(zeros(Int, 3, div(N_k, 2)), ones(Int, div(N_k, 2))' .* [1, 0, 0], -ones(Int, div(N_k, 2))' .* [1, 0, 0], zeros(Int, 3, div(N_k, 2)))
    return w_hat, q_2bz_ind, q_2bz_shift
end

function find_r_μ(prob::BSEProblem1D, N_μ::Int) # copy of this function in exciting_interface.jl
    r_μ_indices = round.(Int, range(1, stop = size_r(prob)[1] + 1, length = N_μ + 1)[1:(end - 1)])
    return r_μ_indices
end

function optical_absorption_vector(prob::BSEProblem1D, direction)
    return optical_absorption_vector(prob.u_v, prob.u_c, prob.E_v, prob.E_c, prob.prob.r_unit, prob.prob.k_bz)
end

function optical_absorption_vector(u_v, u_c, E_v, E_c, r_unit, k_bz) # TODO: make general? only works in 1d atm
    N_unit = length(r_unit)
    N_v = size(u_v, 2)
    N_c = size(u_c, 2)
    N_k = length(k_bz)
    Δr = r_unit[2] - r_unit[1]
    l = N_unit * Δr

    d = vec([im * Δr * dot(u_c[:, ic, ik],
                  modified_momentum_operator(r_unit, l, 2 * pi / l * k_bz[ik]) *
                    u_v[:, iv, ik]) / (E_c[ic, ik] - E_v[iv, ik])
            for iv in 1:N_v, ic in 1:N_c, ik in 1:N_k])

    return d
end

# reference implementations for testing

function V_entry_realspace(V_func, iv, ic, ik, jv, jc, jk, u_v, u_c, r_super, r_unit)
    N_unit = size(u_v, 1)
    N_cells = div(length(r_super), N_unit)
    Δr = r_super[2] - r_super[1]
    L = N_cells * N_unit * Δr

    v = complex(0.)
    for (ir, r_1) in enumerate(r_super)
        for (jr, r_2) in enumerate(r_unit)
            diff_r = supercell_difference(r_1, r_2, L)
            v += (conj(u_c[mod1(ir, N_unit), ic, ik]) * u_v[mod1(ir, N_unit), iv, ik]) *
                V_func(diff_r, 0.) *
                (conj(u_v[jr, jv, jk]) * u_c[jr, jc, jk])
        end
    end
    return Δr^2 * v / N_cells
end
function assemble_exact_V_realspace(prob::BSEProblem1D, V_func)
    r_super = prob.prob.r_super
    r_unit = prob.prob.r_unit
    N_v = prob.N_v
    N_c = prob.N_c
    N_k = prob.N_k
    u_v = prob.u_v
    u_c = prob.u_c

    V_reshaped = zeros(Complex{Float64}, N_v, N_c, N_k, N_v, N_c, N_k)
    @showprogress 1 "Assemble exact V ..." for jk in 1:N_k, jc in 1:N_c, jv in 1:N_v, ik in 1:N_k, ic in 1:N_c, iv in 1:N_v
        V_reshaped[iv, ic, ik, jv, jc, jk] =
            V_entry_realspace(V_func, iv, ic, ik, jv, jc, jk, u_v, u_c, r_super, r_unit)
    end
    V = reshape(V_reshaped, N_v * N_c * N_k, N_v * N_c * N_k)
    for i in 1:(N_v * N_c * N_k)
        V[i, i] = real(V[i, i])
    end

    return V
end

function W_entry_fast(w_hat, iv, ic, ik, jv, jc, jk, u_v, u_c, r_super, r_unit, k_bz)
    N_unit = length(r_unit)
    N_cells = div(length(r_super), N_unit)
    N_k = length(k_bz)
    Δr = r_unit[2] - r_unit[1]
    l = N_unit * Δr
    L = N_cells * l

    ijk = mod1(ik - jk + 1, N_k)
    if ik - jk >= div(N_k, 2)
        G_shift = -1
    elseif ik - jk < -div(N_k, 2)
        G_shift = 1
    else
        G_shift = 0
    end

    U_c_hat_shift = circshift(fft(u_c[:, ic, ik] .* conj.(u_c[:, jc, jk])), -G_shift) * (l / N_unit)
    U_v_hat_shift = circshift(fft(u_v[:, iv, ik] .* conj.(u_v[:, jv, jk])), -G_shift) * (l / N_unit)

    w = 1 / L * U_c_hat_shift' * (@view(w_hat[:, :, ijk]) * U_v_hat_shift)

    return w
end
function W_entry_realspace(W_func, iv, ic, ik, jv, jc, jk, u_v, u_c, r_super, r_unit, k_bz)
    N_unit = size(u_v, 1)
    N_cells = div(length(r_super), N_unit)
    Δr = r_super[2] - r_super[1]
    l = N_unit * Δr
    L = N_cells * l

    w = complex(0.)
    for (ir, r_1) in enumerate(r_super)
        for (jr, r_2) in enumerate(r_unit)
            w += (conj(u_c[mod1(ir, N_unit), ic, ik]) * u_c[mod1(ir, N_unit), jc, jk]) *
                exp(-im * 2 * pi / l * (k_bz[ik] - k_bz[jk]) * (r_1 - r_2)) * W_func(r_1, r_2, l, L) *
                (conj(u_v[jr, jv, jk]) * u_v[jr, iv, ik])
        end
    end
    return Δr^2 * w / N_cells
end
function assemble_exact_W_realspace(prob::BSEProblem1D, W_func)
    r_super = prob.prob.r_super
    r_unit = prob.prob.r_unit
    k_bz = prob.prob.k_bz
    N_v = prob.N_v
    N_c = prob.N_c
    N_k = prob.N_k
    u_v = prob.u_v
    u_c = prob.u_c

    W_reshaped = zeros(Complex{Float64}, N_v, N_c, N_k, N_v, N_c, N_k)
    @showprogress 1 "Assemble exact W ..." for jk in 1:N_k, jc in 1:N_c, jv in 1:N_v, ik in 1:N_k, ic in 1:N_c, iv in 1:N_v
        W_reshaped[iv, ic, ik, jv, jc, jk] =
            W_entry_realspace(W_func, iv, ic, ik, jv, jc, jk, u_v, u_c, r_super, r_unit, k_bz)
    end
    W = reshape(W_reshaped, N_v * N_c * N_k, N_v * N_c * N_k)
    for i in 1:(N_v * N_c * N_k)
        W[i, i] = real(W[i, i])
    end

    return W
end

function assemble_exact_H_1d(prob)
    r_super = prob.prob.r_super
    r_unit = prob.prob.r_unit
    k_bz = prob.prob.k_bz
    N_v = prob.N_v
    N_c = prob.N_c
    N_k = prob.N_k
    E_v = prob.E_v
    E_c = prob.E_c
    u_v = prob.u_v
    u_c = prob.u_c
    v_hat = prob.v_hat
    w_hat = prob.w_hat

    H_reshaped = zeros(Complex{Float64}, N_v, N_c, N_k, N_v, N_c, N_k)
    for jk in 1:N_k, jc in 1:N_c, jv in 1:N_v, ik in 1:N_k, ic in 1:N_c, iv in 1:N_v
        if N_v * N_c * (ik - 1) + N_v * (ic - 1) + iv <= N_v * N_c * (jk - 1) + N_v * (jc - 1) + jv # exploit symmetry
            H_reshaped[iv, ic, ik, jv, jc, jk] =
                H_entry_fast(v_hat, w_hat, iv, ic, ik, jv, jc, jk, E_v, E_c, u_v, u_c, r_super, r_unit, k_bz)
        end
    end
    H = reshape(H_reshaped, N_v * N_c * N_k, N_v * N_c * N_k)
    for i in 1:(N_v * N_c * N_k)
        H[i, i] = real(H[i, i])
    end

    return Hermitian(H)
end


text:
module DataReader
  using ..DataParser
  function getDataSetInMemoryFromIndices(file,delimiter,indices,datatype=Any)
    size = DataParser.DataStream.getDataSetCount(file)
    #indices = datasetFunction(file,size,threshold)
    featureCount = DataParser.DataStream.getFeatureCount(file,delimiter);
    stream = DataParser.DataStream.getDataSet(file)
    dataset = cell(length(indices),featureCount)
    index=0;
    loc=1
    for line in eachline(stream)
      index+=1
      if in(index,indices)
        dataset[loc,:]= split(strip(line,'\n'),delimiter)
        loc+=1
      end

    end
    #println(size(testset))

    DataParser.DataStream.closeDataSet(stream)

    return dataset
  end

end

max_stars_repo_name: ajermaky/CustomMLLibrary

max_issues_repo_name: ajermaky/CustomMLLibrary

max_forks_repo_name: ajermaky/CustomMLLibrary


text:
# DOT Language
# ============

function nfa2dot(nfa::NFA)
    out = IOBuffer()
    println(out, "digraph {")
    println(out, "  graph [ rankdir = LR ];")
    println(out, "  0 -> 1;")
    println(out, "  0 [ shape = point ];")
    serials = Dict(s => i for (i, s) in enumerate(traverse(nfa.start)))
    for s in keys(serials)
        println(out, "  $(serials[s]) [ shape = $(s == nfa.final ? "doublecircle" : "circle") ];")
        for (e, t) in s.edges
            println(out, " $(serials[s]) -> $(serials[t]) [ label = \"$(edge2str(e))\" ];")
        end
    end
    println(out, "}")
    return String(take!(out))
end

function dfa2dot(dfa::DFA)
    out = IOBuffer()
    println(out, "digraph {")
    println(out, "  graph [ rankdir = LR ];")
    println(out, "  start -> 1;")
    println(out, "  start [ shape = point ];")
    serials = Dict(s => i for (i, s) in enumerate(traverse(dfa.start)))
    for s in keys(serials)
        println(out, "  $(serials[s]) [ shape = $(s.final ? "doublecircle" : "circle") ];")
        for (e, t) in s.edges
            println(out, "  $(serials[s]) -> $(serials[t]) [ label = \"$(edge2str(e))\" ];")
        end
        if !isempty(s.eof_actions)
            println(out, "  eof$(serials[s]) [ shape = point ];")
            println(out, "  $(serials[s]) -> eof$(serials[s]) [ label = \"$(eof_label(s.eof_actions))\", style = dashed ];")
        end
    end
    println(out, "}")
    return String(take!(out))
end

function machine2dot(machine::Machine)
    out = IOBuffer()
    println(out, "digraph {")
    println(out, "  graph [ rankdir = LR ];")
    println(out, "  start -> 1;")
    println(out, "  start [ shape = point ];")
    for s in traverse(machine.start)
        println(out, "  $(s.state) [ shape = $(s.state ∈ machine.final_states ? "doublecircle" : "circle") ];")
        for (e, t) in s.edges
            println(out, "  $(s.state) -> $(t.state) [ label = \"$(edge2str(e))\" ];")
        end
        if haskey(machine.eof_actions, s.state) && !isempty(machine.eof_actions[s.state])
            println(out, "  eof$(s.state) [ shape = point ];")
            println(out, "  $(s.state) -> eof$(s.state) [ label = \"$(eof_label(machine.eof_actions[s.state]))\", style = dashed ];")
        end
    end
    println(out, "}")
    return String(take!(out))
end

function edge2str(edge::Edge)
    out = IOBuffer()

    function printbyte(b, inrange)
        # TODO: does this work?
        if inrange && b == UInt8('-')
            print(out, "\\\\-")
        elseif inrange && b == UInt8(']')
            print(out, "\\\\]")
        else
            print(out, escape_string(b ≤ 0x7f ? escape_string(string(Char(b))) : @sprintf("\\x%x", b)))
        end
    end

    # output labels
    if isempty(edge.labels)
        print(out, 'ϵ')
    elseif length(edge.labels) == 1
        print(out, '\'')
        printbyte(first(edge.labels), false)
        print(out, '\'')
    else
        print(out, '[')
        for r in range_encode(edge.labels)
            if length(r) == 1
                printbyte(first(r), true)
            elseif length(r) == 2
                printbyte(first(r), true)
                printbyte(last(r), true)
            else
                @assert length(r) > 1
                printbyte(first(r), true)
                print(out, '-')
                printbyte(last(r), true)
            end
        end
        print(out, ']')
    end

    # output conditions
    if !isempty(edge.precond)
        print(out, '(')
        join(out, ((value == NONE ? "false" : value == TRUE ? name : value == FALSE ? string("!", name) : "true") for (name, value) in edge.precond), ',')
        print(out, ')')
    end

    # output actions
    if !isempty(edge.actions)
        print(out, '/')
        join(out, action_names(edge.actions), ',')
    end

    return String(take!(out))
end

function eof_label(actions::ActionList)
    out = IOBuffer()
    print(out, "EOF")
    if !isempty(actions)
        print(out, '/')
        join(out, action_names(actions), ',')
    end
    return String(take!(out))
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

if !isdefined(@__MODULE__, :DefsGroupId)
    const DefsGroupId = String
else
    @warn("Skipping redefinition of DefsGroupId to String")
end


text:
export getConvDiffMatrix

"""
function getConvDiffMatrix(M::RegularMesh,sig::Float64,v,iddir::Array{Int,1},idneu::Array{Int,1},
	                      iddirn::Array{Int,1},idneun::Array{Int,1},idint::Array{Int,1})

builds convection diffusion matrix on interior nodes and returns contribution of dirichlet and neuman boundary.
"""
function getConvDiffMatrix(M::RegularMesh,sig::Float64,v,iddir::Array{Int,1},idneu::Array{Int,1}, iddirn::Array{Int,1},idneun::Array{Int,1},idint::Array{Int,1})

	# biuld matrix on padded mesh
	Lap  = getLaplacianMatrix(M)
	Conv = getConvectionMatrix(M,v)
	Apad = sig*Lap + Conv

	# projection for gost point neighbors
	Id          = sparse(1.0I, prod(M.n.+2), prod(M.n.+2))[:,idint]
	Id[iddir,:] = -Id[iddirn,:]  # dirichlect boundary points
	Id[idneu,:] =  Id[idneun,:]  # neuman boundary points

	# matrix for interior nodes
	A    = Apad[idint,:]*Id

	# boundary nodes
	Adir = Apad[idint,iddir]
	Aneu = Apad[idint,idneu]

	return A, Adir, Aneu
end

max_stars_repo_name: JuliaInv/ConvDiffMIPDECO

max_issues_repo_name: JuliaInv/ConvDiffMIPDECO

max_forks_repo_name: JuliaInv/ConvDiffMIPDECO


text:
"""
    test_unbound_args(module::Module)

Test that all methods in `module` and its submodules do not have
unbound type parameter.
"""
function test_unbound_args(m::Module)
    # https://github.com/JuliaLang/julia/pull/31972
    # @test detect_unbound_args(m; recursive=true) == []
    @test detect_unbound_args_recursively(m) == []
end

function detect_unbound_args_recursively(m)
    methods = []
    walkmodules(m) do x
        append!(methods, detect_unbound_args(x))
    end
    return methods
end


text:
using MAT
using SparseArrays
using LinearAlgebra
include("../src/recursive_spectral.jl")
mat = matread("../data/Florida_Bay_Dataset.mat")
A125 = mat["A125"]                    # This is the processed dataset--a directed who-eats-who graph
n = size(A125,1)
C = round.(Int64,mat["C"])            # These are four diffrent ways to classify the species
label = C[:,3]                        # These is the labels of Li et al. NeurIPs 2017, ignoring "predator" and "no predator" tags

## Need to form different graphs based on different reductions
include("../src/Motif_Adj_Mat.jl")
AFan = ExtractFanMat(A125)      # So we can run Lambda-Louvain on the bi-fan hypergraph
ASim = turntounweighted(A125)   # So we can run Louvain on an unweighted version

# Scale the graphs so that their total volume is the same, for easier comparison
mSim = sum(ASim)
mFan = sum(AFan)
AFan = mSim/mFan*AFan
m = mSim
n = size(ASim,1)

numpoints = 20
mins = 10 .^range(log10(5),log10(75),length=numpoints)
times = zeros(20,2)

SpecFan = zeros(n,numpoints)
SpecSim = zeros(n,numpoints)
SpecFari = zeros(numpoints)
SpecSari = zeros(numpoints)

for i = 1:length(mins)
    minsize = mins[i]

    # Run recursive spectral on the original graph
    start = time()
    println("Minsize = $minsize")
    sets = recursive_spectral(ASim,minsize)
    c = set2clusvec(sets,n)
    times[i,1] = time()-start
    SpecSim[:,i] .= c
    SpecSari = ari(round.(Int64,c),label)

    # Run it on the motif adjacency matrix
    start = time()
    sets = recursive_spectral(AFan,minsize)
    c = set2clusvec(sets,n)
    times[i,2] = time()-start
    SpecFan[:,i] .= c
    SpecFari[i] = ari(round.(Int64,c),label)
end

matwrite("Output/RecursiveSpectral_Flbay.mat", Dict("times"=>times,
"SpecFan"=>SpecFan, "SpecSim"=>SpecSim))

max_stars_repo_name: nveldt/ParamCC

max_issues_repo_name: nveldt/ParamCC

max_forks_repo_name: nveldt/ParamCC


text:
function two_species_approx_wavenumber(ω::Number, medium, species)
  sp1 = species[1]
  sp2 = species[2]

  if outer_radius(sp1) > outer_radius(sp2)
    @warn("method two_species_approx was designed for species[1] to be the smallest.")
  end
  vol = volume_fraction(sp1)
  rhoeff =  medium.ρ*(medium.ρ + sp1.particle.medium.ρ - vol*(medium.ρ - sp1.particle.medium.ρ)) /
  (medium.ρ + sp1.particle.medium.ρ + vol*(medium.ρ - sp1.particle.medium.ρ))

  kTS = wavenumber_low_volumefraction(ω, medium, [sp1])
  mS = Acoustic(2; ρ=rhoeff, c= ω/kTS)
  wavenumber_low_volumefraction(ω, mS, [sp2])
end

function two_species_approx_wavenumber(ωs::AbstractArray, medium, species)
  sp1 = species[1]
  sp2 = species[2]

  if outer_radius(sp1) > outer_radius(sp2)
    @warn("method two_species_approx was designed for species[1] to be the smallest.")
  end
  vol = volume_fraction(sp1)
  rhoeff =  medium.ρ*(medium.ρ + sp1.particle.medium.ρ - vol*(medium.ρ - sp1.particle.medium.ρ))/
  (medium.ρ + sp1.particle.medium.ρ + vol*(medium.ρ - sp1.particle.medium.ρ))

  kTSs = wavenumber_low_volumefraction(ωs, medium, [sp1])
  kTLs = [
  begin
    mS = Acoustic(2; ρ=rhoeff, c= ωs[i]/kTSs[i])
    wavenumber_low_volumefraction(ωs[i], mS, [sp2])
  end
  for i in eachindex(ωs)]
end

max_stars_repo_name: arturgower/EffectiveWaves


text:
import Base.show

function show{R}(io::IO, w::Wire{R})
  show(io, typeof(w))
  print(io, "(0b")
  print(io, join([w.assigned[idx] ? (w.values[idx] ? "1" : "0") : "X" for idx = length(w):-1:1],""))
  print(io, ")")
end

function show{R}(io::IO, ::Type{Wire{R}})
  print(io, "Wire{",R.stop,":",R.start,"v}")
end


text:
module JavaCall

# Cuts runtests.jl time by one to two seconds
# https://github.com/JuliaLang/julia/pull/34896
@static if isdefined(Base, :Experimental) &&
           isdefined(Base.Experimental, Symbol("@optlevel"))
    Base.Experimental.@optlevel 1
end

export JavaObject, JavaMetaClass,
       jint, jlong, jbyte, jboolean, jchar, jshort, jfloat, jdouble, jvoid,
       JObject, JClass, JMethod, JConstructor, JField, JString,
       JavaRef, JavaLocalRef, JavaGlobalRef, JavaNullRef,
       @jimport, jcall, jfield, jlocalframe, isnull,
       getname, getclass, listmethods, getreturntype, getparametertypes, classforname,
       listfields, gettype,
       narrow

# using Compat, Compat.Dates

# using Sys: iswindows, islinux, isunix, isapple

import DataStructures: OrderedSet
using Dates

@static if Sys.iswindows()
    using WinReg
end

import Base: convert, unsafe_convert, unsafe_string, Ptr

JULIA_COPY_STACKS = false

include("JNI.jl")
using .JNI
import .JNI.Threads
include("jvm.jl")
include("core.jl")
include("convert.jl")
include("reflect.jl")

Base.@deprecate_binding jnifunc JavaCall.JNI.jniref[]

function __init__()
    global JULIA_COPY_STACKS = get(ENV, "JULIA_COPY_STACKS", "") ∈ ("1", "yes")
    if ! Sys.iswindows()
        # On Windows, JULIA_COPY_STACKS is not needed and causes crash
        if VERSION ≥ v"1.1-" && VERSION < v"1.3-"
            @warn("JavaCall does not work correctly on Julia v$VERSION. \n" *
                    "Either use Julia v1.0.x, or v1.3.0 or higher.\n"*
                    "For 1.3 onwards, please also set the environment variable `JULIA_COPY_STACKS` to be `1` or `yes`")
        end
        if VERSION ≥ v"1.3-" && ! JULIA_COPY_STACKS
            @warn("JavaCall needs the environment variable `JULIA_COPY_STACKS` to be `1` or `yes`.\n"*
                  "Calling the JVM may result in undefined behavior.")
        end
    end
    Threads.resize_nthreads!(_jmc_cache)
end

end # module


text:
import .Grid

function neighbors_lambda(g::Grid.GridGraph)
    idx, pi, lam = Grid.neighboridx_lambda(g)
    return IncidenceIndex(idx, pi), lam
end

max_stars_repo_name: EQt/graphidx

max_issues_repo_name: EQt/graphidx

max_forks_repo_name: EQt/graphidx


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md

using Pkg, Documenter, Literate, JuliaFEM

# really?
juliafem_dir = abspath(joinpath(dirname(pathof(JuliaFEM)), ".."))

"""
    copy_docs(pkg_name)

Copy documentation of some package `pkg_name` to `docs/src/pkg_name`,
where `pkg_name` is the name of package. Return true if success, false
otherwise.

If package is undocumented, i.e. directory `docs/src` is missing,
but there still exists `README.md`, copy that file to
`docs/src/pkg_name/index.md`.
"""
function copy_docs(pkg_name)

    local pkg_path
    try
        pkg = getfield(JuliaFEM, Symbol(pkg_name))
        pkg_path = abspath(joinpath(pathof(pkg), "..", ".."))
    catch
        @info("Could not find package $pkg_name from JuliaFEM namespace, pkg_path unknown.")
        return false
    end

    src_dir = joinpath(pkg_path, "docs", "src")
    dst_dir = joinpath(juliafem_dir, "docs", "src", "packages", pkg_name)
    pkg_dir = joinpath(juliafem_dir, "docs", "src", "packages")
    isdir(pkg_dir) || mkpath(pkg_dir)

    # if can find pkg_name/docs/src =>
    # copy that to docs/src/packages/pkg_name
    if isdir(src_dir)
        isdir(dst_dir) || cp(src_dir, dst_dir)
        @info("Copied documentation of package $pkg_name from $src_dir succesfully.")
        return true
    end

    # if can find pkg_name/README.md =>
    # copy that to docs/src/packages/pkg_name/index.md
    readme_file = joinpath(pkg_path, "README.md")
    if isfile(readme_file)
        isdir(dst_dir) || mkpath(dst_dir)
        dst_file = joinpath(dst_dir, "index.md")
        isfile(dst_file) || cp(readme_file, dst_file)
        @info("Copied README.md of package $pkg_name from $readme_file succesfully.")
        return true
    end

    @warn("Cannot copy documentation of package $pkg_name from $src_dir: " *
          "No such directory exists. (Is the package in REQUIRE of JuliaFEM?)")

    return false

end

"""
    add_page!(part, page)

Add new page to documentation. Part can be USER_GUIDE, DEVELOPER_GUIDE or
PACKAGES. `page` is a string containing path to the page, or `Pair`
where key is the name of the page and value is the path to the page.

# Examples

If page starts with `#`, i.e. having title, page is automatically included
to documentation with that name:

```julia
add_page!(PACKAGES, "MyPackage/index.md")
```

Title can be changed using `Pair`, i.e.

```julia
add_page!(PACKAGES, "Theory" => "MyPackage/theory.md")
```

"""
function add_page!(dst, src)
    file = isa(src, Pair) ? src.second : src
    src_dir = joinpath(juliafem_dir, "docs", "src")
    if isfile(joinpath(src_dir, file))
        push!(dst, src)
        return true
    end
    if isfile(joinpath(src_dir, "packages", file))
        push!(dst, joinpath("packages", src))
        return true
    end
    @warn("Cannot add page $file: no such file")
    return false
end

"""
    generate_developers_guide()

Generate Developer's guide.
"""
function generate_developers_guide()
    DEVELOPER_GUIDE = []
    if copy_docs("FEMBase")
        add_page!(DEVELOPER_GUIDE, "FEMBase/mesh.md")
        add_page!(DEVELOPER_GUIDE, "FEMBase/fields.md")
        add_page!(DEVELOPER_GUIDE, "FEMBase/basis.md")
        add_page!(DEVELOPER_GUIDE, "FEMBase/integration.md")
        add_page!(DEVELOPER_GUIDE, "FEMBase/elements.md")
        add_page!(DEVELOPER_GUIDE, "FEMBase/problems.md")
        add_page!(DEVELOPER_GUIDE, "FEMBase/solvers.md")
        add_page!(DEVELOPER_GUIDE, "FEMBase/postprocessing.md")
        add_page!(DEVELOPER_GUIDE, "FEMBase/results.md")
        add_page!(DEVELOPER_GUIDE, "FEMBase/materials.md")
    end
    return DEVELOPER_GUIDE
end

"""
    generate_packages()

Generate single page description for packages by fetching the README.md, index.md
or similar from the package documentation. Returns a vector containing pages for
a documentation.
"""
function generate_packages()
    PACKAGES = []
    copy_docs("FEMBase")                && add_page!(PACKAGES, "FEMBase/index.md")
    copy_docs("FEMBasis")               && add_page!(PACKAGES, "FEMBasis/index.md")
    copy_docs("FEMQuad")                && add_page!(PACKAGES, "FEMQuad/index.md")
    copy_docs("FEMSparse")              && add_page!(PACKAGES, "FEMSparse/index.md")
    copy_docs("Materials")              && add_page!(PACKAGES, "Materials/index.md")
    copy_docs("AsterReader")            && add_page!(PACKAGES, "AsterReader/index.md")
    copy_docs("AbaqusReader")           && add_page!(PACKAGES, "AbaqusReader/index.md")
    copy_docs("LinearImplicitDynamics") && add_page!(PACKAGES, "LinearImplicitDynamics/index.md")
    copy_docs("HeatTransfer")           && add_page!(PACKAGES, "HeatTransfer/index.md")
    copy_docs("PlaneElasticity")        && add_page!(PACKAGES, "PlaneElasticity/index.md")
    copy_docs("FEMBeam")                && add_page!(PACKAGES, "FEMBeam/index.md")
    copy_docs("FEMCoupling")            && add_page!(PACKAGES, "FEMCoupling/index.md")
    copy_docs("FEMTruss")               && add_page!(PACKAGES, "FEMTruss/index.md")
    copy_docs("Mortar2D")               && add_page!(PACKAGES, "Mortar2D/index.md")
    copy_docs("Mortar3D")               && add_page!(PACKAGES, "Mortar3D/index.md")
    copy_docs("MortarContact2D")        && add_page!(PACKAGES, "MortarContact2D/index.md")
    copy_docs("MortarContact2DAD")      && add_page!(PACKAGES, "MortarContact2DAD/index.md")
    copy_docs("OptoMechanics")          && add_page!(PACKAGES, "OptoMechanics/index.md")
    copy_docs("Miniball")               && add_page!(PACKAGES, "Miniball/index.md")
    copy_docs("ModelReduction")         && add_page!(PACKAGES, "ModelReduction/index.md")
    copy_docs("NodeNumbering")          && add_page!(PACKAGES, "NodeNumbering/index.md")
    #copy_docs("Xdmf")                   && add_page!(PACKAGES, "Xdmf/index.md")
    copy_docs("UMAT")                   && add_page!(PACKAGES, "UMAT/index.md")
    return PACKAGES
end

"""
    license_stripper(s)

Strip license strings away from source file.
"""
function license_stripper(s)
    # return replace(s, r"# This file is a part of JuliaFEM.*?\n\n"sm, "")
    lines = split(s, '\n')
    function myfilt(line)
        occursin(line, "# This file is a part of JuliaFEM.") && return false
        occursin(line, "# License is MIT") && return false
        return true
    end
    new_s = join(filter(myfilt, lines), '\n')
    return new_s
end

"""
    generate_example(example_src, dst_dir)

Given excutable Julia script file, use Literate.jl to render file to Markdown
file. If there exists a directory with the same name than script, but without
extension, it will be copied to destination also. Function returns a path to
the generated example file

# Example

If `example_src` is `examples/my_example_analysis.jl`, and there exists a
directory `examples/my_example_analysis`, there will be a rendered file
`dst_dir/my_example_analysis.md` and additional directory
`dst_dir/my_example_analysis` containing complementary material like mesh,
result files and so on.
"""
function generate_example(example_src, dst_dir)
    dst_file = Literate.markdown(example_src, dst_dir;
                                 documenter=true,
                                 preprocess=license_stripper)

    # ex_dir = directory for complementary material
    complementary_dir = first(splitext(basename(example_src)))
    src_ex_dir = joinpath(juliafem_dir, "examples", complementary_dir)
    dst_ex_dir = joinpath(dst_dir, complementary_dir)
    if isdir(src_ex_dir) && !isdir(dst_ex_dir)
        cp(src_ex_dir, dst_ex_dir)
    end
    docs_src = joinpath(juliafem_dir, "docs", "src")
    return relpath(dst_file, docs_src)
end

"""
    get_example_files(ex_src)

Return all example files from directory `examples`.
"""
function get_example_files(ex_src)
    ex_flt(s) = endswith(s, ".jl") && !startswith(s, "test_")
    return filter(ex_flt, readdir(ex_src))
end

"""
    generate_examples()

Generate examples using Literate.jl. Returns a vector containing pages for
a documentation.
"""
function generate_examples()
    EXAMPLES = []
    ex_src = joinpath(juliafem_dir, "examples")
    ex_dst = joinpath(juliafem_dir, "docs", "src", "examples")
    for ex_file in get_example_files(ex_src)
        example_file = joinpath(ex_src, ex_file)
        markdown_file = generate_example(example_file, ex_dst)
        add_page!(EXAMPLES, markdown_file)
    end
    return EXAMPLES
end

# Collect all together
USER_GUIDE = []
EXAMPLES = generate_examples()
DEVELOPER = generate_developers_guide()
PACKAGES = generate_packages()
APIDOC = ["api.md"]

PAGES = [
         "Home" => "index.md",
         "User's guide" => USER_GUIDE,
         "Examples" => EXAMPLES,
         "Developer's guide" => DEVELOPER,
         "Description of packages" => PACKAGES,
         "API documentation" => APIDOC
        ]

@info("Pages in documentation", PAGES)

makedocs(modules=[JuliaFEM],
         format = Documenter.HTML(),
         checkdocs = :all,
         sitename = "JuliaFEM",
         analytics = "UA-83590644-1",
         pages = PAGES)


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

include("model_AliasPathType.jl")
include("model_AliasType.jl")
include("model_BasicDependency.jl")
include("model_DebugSetting.jl")
include("model_Dependency.jl")
include("model_Deployment.jl")
include("model_DeploymentExportResult.jl")
include("model_DeploymentExtended.jl")
include("model_DeploymentExtendedFilter.jl")
include("model_DeploymentListResult.jl")
include("model_DeploymentOperation.jl")
include("model_DeploymentOperationProperties.jl")
include("model_DeploymentOperationsListResult.jl")
include("model_DeploymentProperties.jl")
include("model_DeploymentPropertiesExtended.jl")
include("model_DeploymentValidateResult.jl")
include("model_ExportTemplateRequest.jl")
include("model_GenericResourceFilter.jl")
include("model_HttpMessage.jl")
include("model_Identity.jl")
include("model_ParametersLink.jl")
include("model_Plan.jl")
include("model_Provider.jl")
include("model_ProviderListResult.jl")
include("model_ProviderResourceType.jl")
include("model_Resource.jl")
include("model_ResourceGroup.jl")
include("model_ResourceGroupExportResult.jl")
include("model_ResourceGroupFilter.jl")
include("model_ResourceGroupListResult.jl")
include("model_ResourceGroupPatchable.jl")
include("model_ResourceGroupProperties.jl")
include("model_ResourceListResult.jl")
include("model_ResourceManagementErrorWithDetails.jl")
include("model_ResourceProviderOperationDisplayProperties.jl")
include("model_ResourcesMoveInfo.jl")
include("model_Sku.jl")
include("model_SubResource.jl")
include("model_TagCount.jl")
include("model_TagDetails.jl")
include("model_TagValue.jl")
include("model_TagsListResult.jl")
include("model_TargetResource.jl")
include("model_TemplateLink.jl")
include("model_GenericResource.jl")


text:
function blocks(A, rows)
    block_size = div(size(A,1),rows)
    blks = Array{SparseMatrixCSC, 2}(rows, rows)
    for x in 1:rows
        for y in 1:rows
            xr = (x-1)*block_size+1 : ((x == rows) ? size(A,1) : x * block_size)
            yr = (y-1)*block_size+1 : ((y == rows) ? size(A,1) : y * block_size)
            blks[x,y] = A[xr,yr]
        end
    end

    blks
end

max_stars_repo_name: ligmg/ligmg

max_issues_repo_name: ligmg/ligmg

max_forks_repo_name: ligmg/ligmg


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

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

# ╔═╡ a358b3ab-b1d4-471e-a9be-565e7ec339de
module includes
	using FromFile
	using Distances
	@from "../Oboe/Oboe.jl" using Oboe
	export Oboe
	########
	using HypertextLiteral
	include("./viz-ui.jl")
	include("./vega-specs.jl")
	include("./exporter.jl")
end

# ╔═╡ 94762b1c-4f03-43b9-b288-4f4b77894d3b
begin
	using VegaLite,VegaDatasets
	using DataFrames, CSV, Statistics, Revise, PlutoUI, HypertextLiteral

	Oboe = includes.Oboe
	UI = includes.VizUI
	Specs = includes.VegaSpecs
	Exp = includes.Exporter
	#====PREP==THE==DATA=================#

	fips = map(Oboe.normalizeOneDigitFIPS, ["6", "41", "53"])

	nsRaw = Oboe.censorFluteTractByFIPS(Oboe.rdWholeUS(), fips)

	APs = Oboe.getProcessedAPs(Oboe.rdBTS() |> Oboe.grpBTS, Oboe.rdAPs())

	ns = Oboe.assignDsgAPs(nsRaw,APs)

	#cheapest way of filtering for _relevant_ county IDs in plots
	ctys = Oboe.aggByCty(ns)
	insertcols!(ctys,1,:id =>
		map(r-> parse(Int,r.Ste*"000") + parse(Int,r.Cty),eachrow(ctys)) )

	dsgAPcodes = ns.IATA_Code |> unique

	dsgAPs = filter(r -> r.IATA_Code ∈ dsgAPcodes,eachrow(APs)) |> DataFrame
end

# ╔═╡ fe20d728-2096-11ec-2740-4368d1340f88
md"""
# Network Explorer

This is a collection of plots of the *transportation network*, i.e.
the set of airports, statistical divisions and the air and commuter connections
between them. The transportation network is created and processed by Oboe.
Edit the `fips` variable to modify the set of states that are included in the
processing.

---
Authors: Kara Ignatenko, 2021;
	     Yaroslav Salii, 2021

Changelog:
- 2021-10-09 v.1.0: First complete version
"""

# ╔═╡ 75bc7766-500b-45fc-b5dc-af144c17d947
begin
	picker = @bind selName Select([
		"bySte" => "Tracts and APs, with state boundaries",
		"byCty" => "Tracts and APs, with county and state boundaries",
		"pop"   => "Population by county"
	])
	@htl("""Select plot: $picker""")
end

# ╔═╡ 592f98f7-5cc1-448b-aa7d-0d7a0be98835
begin
	specs = Dict(
		"bySte" => Specs.pltTractsAndAPsBySte(ns, dsgAPs, fips),
		"byCty" => Specs.pltTractsAndAPsByCty(ns, dsgAPs, fips, ctys),
		"pop"   => Specs.pltPopCtyWithBorders(ctys, fips)
	)

	selectedPlot = get(specs, selName, @vlplot())

	UI.vegaEmbed(selectedPlot)
end

# ╔═╡ f2bfd623-223b-4249-a57a-924f359cbede
begin
	sel = @bind format Select(["svg", "pdf"])
	cbx = @bind exportingEnabled UI.booleanButton("Export")
	figDir = Exp.figDir
	@htl("""Select export format: $sel $cbx – Will be saved to $figDir""")
end

# ╔═╡ d84d3340-9fec-470e-9717-4261c5c0cce2
if exportingEnabled
	Exp.savePlot("$(selName)_", selectedPlot, format)
end

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
Distances = "b4f34e82-e78d-54a5-968a-f98e89d6e8f7"
FromFile = "ff7dd447-1dcb-4ce3-b8ac-22a812192de7"
HypertextLiteral = "ac1192a8-f4b3-4bfe-ba22-af5b92cd3ab2"
PlutoUI = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
Revise = "295af30f-e4ad-537b-8983-00126c2a3abe"
Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
VegaDatasets = "0ae4a718-28b7-58ec-9efb-cded64d6d5b4"
VegaLite = "112f6efa-9a02-5b7d-90c0-432ed331239a"

[compat]
CSV = "~0.9.4"
DataFrames = "~1.2.2"
Distances = "~0.10.4"
FromFile = "~0.1.1"
HypertextLiteral = "~0.9.0"
PlutoUI = "~0.7.12"
Revise = "~3.1.20"
VegaDatasets = "~2.1.1"
VegaLite = "~2.6.0"
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

[[ArgTools]]
uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"

[[Artifacts]]
uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"

[[Base64]]
uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"

[[CSV]]
deps = ["CodecZlib", "Dates", "FilePathsBase", "Mmap", "Parsers", "PooledArrays", "SentinelArrays", "Tables", "Unicode", "WeakRefStrings"]
git-tree-sha1 = "3a877c2fc5c9b88ed7259fd0bdb7691aad6b50dc"
uuid = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
version = "0.9.4"

[[ChainRulesCore]]
deps = ["Compat", "LinearAlgebra", "SparseArrays"]
git-tree-sha1 = "e8a30e8019a512e4b6c56ccebc065026624660e8"
uuid = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
version = "1.7.0"

[[CodeTracking]]
deps = ["InteractiveUtils", "UUIDs"]
git-tree-sha1 = "9aa8a5ebb6b5bf469a7e0e2b5202cf6f8c291104"
uuid = "da1fd8a2-8d9e-5ec2-8556-3022fb5608a2"
version = "1.0.6"

[[CodecZlib]]
deps = ["TranscodingStreams", "Zlib_jll"]
git-tree-sha1 = "ded953804d019afa9a3f98981d99b33e3db7b6da"
uuid = "944b1d66-785c-5afd-91f1-9de20f533193"
version = "0.7.0"

[[Compat]]
deps = ["Base64", "Dates", "DelimitedFiles", "Distributed", "InteractiveUtils", "LibGit2", "Libdl", "LinearAlgebra", "Markdown", "Mmap", "Pkg", "Printf", "REPL", "Random", "SHA", "Serialization", "SharedArrays", "Sockets", "SparseArrays", "Statistics", "Test", "UUIDs", "Unicode"]
git-tree-sha1 = "31d0151f5716b655421d9d75b7fa74cc4e744df2"
uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
version = "3.39.0"

[[CompilerSupportLibraries_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"

[[ConstructionBase]]
deps = ["LinearAlgebra"]
git-tree-sha1 = "f74e9d5388b8620b4cee35d4c5a618dd4dc547f4"
uuid = "187b0558-2788-49d3-abe0-74a17ed4e7c9"
version = "1.3.0"

[[Crayons]]
git-tree-sha1 = "3f71217b538d7aaee0b69ab47d9b7724ca8afa0d"
uuid = "a8cc5b0e-0ffa-5ad4-8c14-923d3ee1735f"
version = "4.0.4"

[[DataAPI]]
git-tree-sha1 = "cc70b17275652eb47bc9e5f81635981f13cea5c8"
uuid = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a"
version = "1.9.0"

[[DataFrames]]
deps = ["Compat", "DataAPI", "Future", "InvertedIndices", "IteratorInterfaceExtensions", "LinearAlgebra", "Markdown", "Missings", "PooledArrays", "PrettyTables", "Printf", "REPL", "Reexport", "SortingAlgorithms", "Statistics", "TableTraits", "Tables", "Unicode"]
git-tree-sha1 = "d785f42445b63fc86caa08bb9a9351008be9b765"
uuid = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
version = "1.2.2"

[[DataStructures]]
deps = ["Compat", "InteractiveUtils", "OrderedCollections"]
git-tree-sha1 = "7d9d316f04214f7efdbb6398d545446e246eff02"
uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
version = "0.18.10"

[[DataValueInterfaces]]
git-tree-sha1 = "bfc1187b79289637fa0ef6d4436ebdfe6905cbd6"
uuid = "e2d170a0-9d28-54be-80f0-106bbe20a464"
version = "1.0.0"

[[DataValues]]
deps = ["DataValueInterfaces", "Dates"]
git-tree-sha1 = "d88a19299eba280a6d062e135a43f00323ae70bf"
uuid = "e7dc6d0d-1eca-5fa6-8ad6-5aecde8b7ea5"
version = "0.4.13"

[[Dates]]
deps = ["Printf"]
uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"

[[DelimitedFiles]]
deps = ["Mmap"]
uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab"

[[Distances]]
deps = ["LinearAlgebra", "Statistics", "StatsAPI"]
git-tree-sha1 = "9f46deb4d4ee4494ffb5a40a27a2aced67bdd838"
uuid = "b4f34e82-e78d-54a5-968a-f98e89d6e8f7"
version = "0.10.4"

[[Distributed]]
deps = ["Random", "Serialization", "Sockets"]
uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"

[[DocStringExtensions]]
deps = ["LibGit2"]
git-tree-sha1 = "a32185f5428d3986f47c2ab78b1f216d5e6cc96f"
uuid = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
version = "0.8.5"

[[DoubleFloats]]
deps = ["GenericLinearAlgebra", "LinearAlgebra", "Polynomials", "Printf", "Quadmath", "Random", "Requires", "SpecialFunctions"]
git-tree-sha1 = "1c962cf7e75c09a5f1fbf504df7d6a06447a1129"
uuid = "497a8b3b-efae-58df-a0af-a86822472b78"
version = "1.1.23"

[[Downloads]]
deps = ["ArgTools", "LibCURL", "NetworkOptions"]
uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6"

[[ExprTools]]
git-tree-sha1 = "b7e3d17636b348f005f11040025ae8c6f645fe92"
uuid = "e2ba6199-217a-4e67-a87a-7c52f15ade04"
version = "0.1.6"

[[FileIO]]
deps = ["Pkg", "Requires", "UUIDs"]
git-tree-sha1 = "3c041d2ac0a52a12a27af2782b34900d9c3ee68c"
uuid = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549"
version = "1.11.1"

[[FilePaths]]
deps = ["FilePathsBase", "MacroTools", "Reexport", "Requires"]
git-tree-sha1 = "919d9412dbf53a2e6fe74af62a73ceed0bce0629"
uuid = "8fc22ac5-c921-52a6-82fd-178b2807b824"
version = "0.8.3"

[[FilePathsBase]]
deps = ["Dates", "Mmap", "Printf", "Test", "UUIDs"]
git-tree-sha1 = "6d4b609786127030d09e6b1ee0e2044ec20eb403"
uuid = "48062228-2e41-5def-b9a4-89aafe57970f"
version = "0.9.11"

[[FileWatching]]
uuid = "7b1f6079-737a-58dc-b8bc-7a2ca5c1b5ee"

[[Formatting]]
deps = ["Printf"]
git-tree-sha1 = "8339d61043228fdd3eb658d86c926cb282ae72a8"
uuid = "59287772-0a20-5a39-b81b-1366585eb4c0"
version = "0.4.2"

[[FromFile]]
git-tree-sha1 = "f4fce099f9e14cd12d91e29fd105ef891fb55297"
uuid = "ff7dd447-1dcb-4ce3-b8ac-22a812192de7"
version = "0.1.1"

[[Future]]
deps = ["Random"]
uuid = "9fa8497b-333b-5362-9e8d-4d0656e87820"

[[GenericLinearAlgebra]]
deps = ["LinearAlgebra", "Printf", "Random"]
git-tree-sha1 = "eddbb6ee8fe2c3244a2c973874a3179c3c4d3ac5"
uuid = "14197337-ba66-59df-a3e3-ca00e7dcff7a"
version = "0.2.6"

[[HTTP]]
deps = ["Base64", "Dates", "IniFile", "Logging", "MbedTLS", "NetworkOptions", "Sockets", "URIs"]
git-tree-sha1 = "14eece7a3308b4d8be910e265c724a6ba51a9798"
uuid = "cd3eb016-35fb-5094-929b-558a96fad6f3"
version = "0.9.16"

[[HypertextLiteral]]
git-tree-sha1 = "72053798e1be56026b81d4e2682dbe58922e5ec9"
uuid = "ac1192a8-f4b3-4bfe-ba22-af5b92cd3ab2"
version = "0.9.0"

[[IOCapture]]
deps = ["Logging", "Random"]
git-tree-sha1 = "f7be53659ab06ddc986428d3a9dcc95f6fa6705a"
uuid = "b5f81e59-6552-4d32-b1f0-c071b021bf89"
version = "0.2.2"

[[IniFile]]
deps = ["Test"]
git-tree-sha1 = "098e4d2c533924c921f9f9847274f2ad89e018b8"
uuid = "83e8ac13-25f8-5344-8a64-a9f2b223428f"
version = "0.5.0"

[[InteractiveUtils]]
deps = ["Markdown"]
uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"

[[Intervals]]
deps = ["Dates", "Printf", "RecipesBase", "Serialization", "TimeZones"]
git-tree-sha1 = "323a38ed1952d30586d0fe03412cde9399d3618b"
uuid = "d8418881-c3e1-53bb-8760-2df7ec849ed5"
version = "1.5.0"

[[InvertedIndices]]
git-tree-sha1 = "bee5f1ef5bf65df56bdd2e40447590b272a5471f"
uuid = "41ab1584-1d38-5bbf-9106-f11c6c58b48f"
version = "1.1.0"

[[IrrationalConstants]]
git-tree-sha1 = "f76424439413893a832026ca355fe273e93bce94"
uuid = "92d709cd-6900-40b7-9082-c6be49f344b6"
version = "0.1.0"

[[IterableTables]]
deps = ["DataValues", "IteratorInterfaceExtensions", "Requires", "TableTraits", "TableTraitsUtils"]
git-tree-sha1 = "70300b876b2cebde43ebc0df42bc8c94a144e1b4"
uuid = "1c8ee90f-4401-5389-894e-7a04a3dc0f4d"
version = "1.0.0"

[[IteratorInterfaceExtensions]]
git-tree-sha1 = "a3f24677c21f5bbe9d2a714f95dcd58337fb2856"
uuid = "82899510-4779-5014-852e-03e436cf321d"
version = "1.0.0"

[[JLLWrappers]]
deps = ["Preferences"]
git-tree-sha1 = "642a199af8b68253517b80bd3bfd17eb4e84df6e"
uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210"
version = "1.3.0"

[[JSON]]
deps = ["Dates", "Mmap", "Parsers", "Unicode"]
git-tree-sha1 = "8076680b162ada2a031f707ac7b4953e30667a37"
uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
version = "0.21.2"

[[JSONSchema]]
deps = ["HTTP", "JSON", "URIs"]
git-tree-sha1 = "2f49f7f86762a0fbbeef84912265a1ae61c4ef80"
uuid = "7d188eb4-7ad8-530c-ae41-71a32a6d4692"
version = "0.3.4"

[[JuliaInterpreter]]
deps = ["CodeTracking", "InteractiveUtils", "Random", "UUIDs"]
git-tree-sha1 = "e273807f38074f033d94207a201e6e827d8417db"
uuid = "aa1ae85d-cabe-5617-a682-6adf51b2e16a"
version = "0.8.21"

[[LazyArtifacts]]
deps = ["Artifacts", "Pkg"]
uuid = "4af54fe1-eca0-43a8-85a7-787d91b784e3"

[[LibCURL]]
deps = ["LibCURL_jll", "MozillaCACerts_jll"]
uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21"

[[LibCURL_jll]]
deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"]
uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0"

[[LibGit2]]
deps = ["Base64", "NetworkOptions", "Printf", "SHA"]
uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"

[[LibSSH2_jll]]
deps = ["Artifacts", "Libdl", "MbedTLS_jll"]
uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8"

[[Libdl]]
uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"

[[LinearAlgebra]]
deps = ["Libdl"]
uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

[[LogExpFunctions]]
deps = ["ChainRulesCore", "DocStringExtensions", "IrrationalConstants", "LinearAlgebra"]
git-tree-sha1 = "34dc30f868e368f8a17b728a1238f3fcda43931a"
uuid = "2ab3a3ac-af41-5b50-aa03-7779005ae688"
version = "0.3.3"

[[Logging]]
uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"

[[LoweredCodeUtils]]
deps = ["JuliaInterpreter"]
git-tree-sha1 = "491a883c4fef1103077a7f648961adbf9c8dd933"
uuid = "6f1432cf-f94c-5a45-995e-cdbf5db27b0b"
version = "2.1.2"

[[MacroTools]]
deps = ["Markdown", "Random"]
git-tree-sha1 = "5a5bc6bf062f0f95e62d0fe0a2d99699fed82dd9"
uuid = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
version = "0.5.8"

[[Markdown]]
deps = ["Base64"]
uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"

[[MbedTLS]]
deps = ["Dates", "MbedTLS_jll", "Random", "Sockets"]
git-tree-sha1 = "1c38e51c3d08ef2278062ebceade0e46cefc96fe"
uuid = "739be429-bea8-5141-9913-cc70e7f3736d"
version = "1.0.3"

[[MbedTLS_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"

[[Missings]]
deps = ["DataAPI"]
git-tree-sha1 = "bf210ce90b6c9eed32d25dbcae1ebc565df2687f"
uuid = "e1d29d7a-bbdc-5cf2-9ac0-f12de2c33e28"
version = "1.0.2"

[[Mmap]]
uuid = "a63ad114-7e13-5084-954f-fe012c677804"

[[Mocking]]
deps = ["ExprTools"]
git-tree-sha1 = "748f6e1e4de814b101911e64cc12d83a6af66782"
uuid = "78c3b35d-d492-501b-9361-3d52fe80e533"
version = "0.7.2"

[[MozillaCACerts_jll]]
uuid = "14a3606d-f60d-562e-9121-12d972cd8159"

[[MutableArithmetics]]
deps = ["LinearAlgebra", "SparseArrays", "Test"]
git-tree-sha1 = "3927848ccebcc165952dc0d9ac9aa274a87bfe01"
uuid = "d8a4904e-b15c-11e9-3269-09a3773c0cb0"
version = "0.2.20"

[[NetworkOptions]]
uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"

[[NodeJS]]
deps = ["Pkg"]
git-tree-sha1 = "905224bbdd4b555c69bb964514cfa387616f0d3a"
uuid = "2bd173c7-0d6d-553b-b6af-13a54713934c"
version = "1.3.0"

[[Nullables]]
git-tree-sha1 = "8f87854cc8f3685a60689d8edecaa29d2251979b"
uuid = "4d1e1d77-625e-5b40-9113-a560ec7a8ecd"
version = "1.0.0"

[[OpenLibm_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "05823500-19ac-5b8b-9628-191a04bc5112"

[[OpenSpecFun_jll]]
deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "13652491f6856acfd2db29360e1bbcd4565d04f1"
uuid = "efe28fd5-8261-553b-a9e1-b2916fc3738e"
version = "0.5.5+0"

[[OrderedCollections]]
git-tree-sha1 = "85f8e6578bf1f9ee0d11e7bb1b1456435479d47c"
uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
version = "1.4.1"

[[Parsers]]
deps = ["Dates"]
git-tree-sha1 = "9d8c00ef7a8d110787ff6f170579846f776133a9"
uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
version = "2.0.4"

[[Pkg]]
deps = ["Artifacts", "Dates", "Downloads", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"]
uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"

[[PlutoUI]]
deps = ["Base64", "Dates", "HypertextLiteral", "IOCapture", "InteractiveUtils", "JSON", "Logging", "Markdown", "Random", "Reexport", "UUIDs"]
git-tree-sha1 = "d1fb76655a95bf6ea4348d7197b22e889a4375f4"
uuid = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
version = "0.7.14"

[[Polynomials]]
deps = ["Intervals", "LinearAlgebra", "MutableArithmetics", "RecipesBase"]
git-tree-sha1 = "0bbfdcd8cda81b8144de4be8a67f5717e959a005"
uuid = "f27b6e38-b328-58d1-80ce-0feddd5e7a45"
version = "2.0.14"

[[PooledArrays]]
deps = ["DataAPI", "Future"]
git-tree-sha1 = "a193d6ad9c45ada72c14b731a318bedd3c2f00cf"
uuid = "2dfb63ee-cc39-5dd5-95bd-886bf059d720"
version = "1.3.0"

[[Preferences]]
deps = ["TOML"]
git-tree-sha1 = "00cfd92944ca9c760982747e9a1d0d5d86ab1e5a"
uuid = "21216c6a-2e73-6563-6e65-726566657250"
version = "1.2.2"

[[PrettyTables]]
deps = ["Crayons", "Formatting", "Markdown", "Reexport", "Tables"]
git-tree-sha1 = "0d1245a357cc61c8cd61934c07447aa569ff22e6"
uuid = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d"
version = "1.1.0"

[[Printf]]
deps = ["Unicode"]
uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"

[[Quadmath]]
deps = ["Printf", "Random", "Requires"]
git-tree-sha1 = "5a8f74af8eae654086a1d058b4ec94ff192e3de0"
uuid = "be4d8f0f-7fa4-5f49-b795-2f01399ab2dd"
version = "0.5.5"

[[REPL]]
deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"]
uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"

[[Random]]
deps = ["Serialization"]
uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

[[RecipesBase]]
git-tree-sha1 = "44a75aa7a527910ee3d1751d1f0e4148698add9e"
uuid = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
version = "1.1.2"

[[Reexport]]
git-tree-sha1 = "45e428421666073eab6f2da5c9d310d99bb12f9b"
uuid = "189a3867-3050-52da-a836-e630ba90ab69"
version = "1.2.2"

[[Requires]]
deps = ["UUIDs"]
git-tree-sha1 = "4036a3bd08ac7e968e27c203d45f5fff15020621"
uuid = "ae029012-a4dd-5104-9daa-d747884805df"
version = "1.1.3"

[[Revise]]
deps = ["CodeTracking", "Distributed", "FileWatching", "JuliaInterpreter", "LibGit2", "LoweredCodeUtils", "OrderedCollections", "Pkg", "REPL", "Requires", "UUIDs", "Unicode"]
git-tree-sha1 = "41deb3df28ecf75307b6e492a738821b031f8425"
uuid = "295af30f-e4ad-537b-8983-00126c2a3abe"
version = "3.1.20"

[[SHA]]
uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"

[[SentinelArrays]]
deps = ["Dates", "Random"]
git-tree-sha1 = "54f37736d8934a12a200edea2f9206b03bdf3159"
uuid = "91c51154-3ec4-41a3-a24f-3f23e20d615c"
version = "1.3.7"

[[Serialization]]
uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"

[[Setfield]]
deps = ["ConstructionBase", "Future", "MacroTools", "Requires"]
git-tree-sha1 = "fca29e68c5062722b5b4435594c3d1ba557072a3"
uuid = "efcf1570-3423-57d1-acb7-fd33fddbac46"
version = "0.7.1"

[[SharedArrays]]
deps = ["Distributed", "Mmap", "Random", "Serialization"]
uuid = "1a1011a3-84de-559e-8e89-a11a2f7dc383"

[[Sockets]]
uuid = "6462fe0b-24de-5631-8697-dd941f90decc"

[[SortingAlgorithms]]
deps = ["DataStructures"]
git-tree-sha1 = "b3363d7460f7d098ca0912c69b082f75625d7508"
uuid = "a2af1166-a08f-5f64-846c-94a0d3cef48c"
version = "1.0.1"

[[SparseArrays]]
deps = ["LinearAlgebra", "Random"]
uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"

[[SpecialFunctions]]
deps = ["ChainRulesCore", "IrrationalConstants", "LogExpFunctions", "OpenLibm_jll", "OpenSpecFun_jll"]
git-tree-sha1 = "793793f1df98e3d7d554b65a107e9c9a6399a6ed"
uuid = "276daf66-3868-5448-9aa4-cd146d93841b"
version = "1.7.0"

[[Statistics]]
deps = ["LinearAlgebra", "SparseArrays"]
uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"

[[StatsAPI]]
git-tree-sha1 = "1958272568dc176a1d881acb797beb909c785510"
uuid = "82ae8749-77ed-4fe6-ae5f-f523153014b0"
version = "1.0.0"

[[TOML]]
deps = ["Dates"]
uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"

[[TableShowUtils]]
deps = ["DataValues", "Dates", "JSON", "Markdown", "Test"]
git-tree-sha1 = "14c54e1e96431fb87f0d2f5983f090f1b9d06457"
uuid = "5e66a065-1f0a-5976-b372-e0b8c017ca10"
version = "0.2.5"

[[TableTraits]]
deps = ["IteratorInterfaceExtensions"]
git-tree-sha1 = "c06b2f539df1c6efa794486abfb6ed2022561a39"
uuid = "3783bdb8-4a98-5b6b-af9a-565f29a5fe9c"
version = "1.0.1"

[[TableTraitsUtils]]
deps = ["DataValues", "IteratorInterfaceExtensions", "Missings", "TableTraits"]
git-tree-sha1 = "78fecfe140d7abb480b53a44f3f85b6aa373c293"
uuid = "382cd787-c1b6-5bf2-a167-d5b971a19bda"
version = "1.0.2"

[[Tables]]
deps = ["DataAPI", "DataValueInterfaces", "IteratorInterfaceExtensions", "LinearAlgebra", "TableTraits", "Test"]
git-tree-sha1 = "1162ce4a6c4b7e31e0e6b14486a6986951c73be9"
uuid = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
version = "1.5.2"

[[Tar]]
deps = ["ArgTools", "SHA"]
uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e"

[[Test]]
deps = ["InteractiveUtils", "Logging", "Random", "Serialization"]
uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

[[TextParse]]
deps = ["CodecZlib", "DataStructures", "Dates", "DoubleFloats", "Mmap", "Nullables", "WeakRefStrings"]
git-tree-sha1 = "eb1f4fb185c8644faa2d18d14c72f2c24412415f"
uuid = "e0df1984-e451-5cb5-8b61-797a481e67e3"
version = "1.0.2"

[[TimeZones]]
deps = ["Dates", "Future", "LazyArtifacts", "Mocking", "Pkg", "Printf", "RecipesBase", "Serialization", "Unicode"]
git-tree-sha1 = "6c9040665b2da00d30143261aea22c7427aada1c"
uuid = "f269a46b-ccf7-5d73-abea-4c690281aa53"
version = "1.5.7"

[[TranscodingStreams]]
deps = ["Random", "Test"]
git-tree-sha1 = "216b95ea110b5972db65aa90f88d8d89dcb8851c"
uuid = "3bb67fe8-82b1-5028-8e26-92a6c54297fa"
version = "0.9.6"

[[URIParser]]
deps = ["Unicode"]
git-tree-sha1 = "53a9f49546b8d2dd2e688d216421d050c9a31d0d"
uuid = "30578b45-9adc-5946-b283-645ec420af67"
version = "0.4.1"

[[URIs]]
git-tree-sha1 = "97bbe755a53fe859669cd907f2d96aee8d2c1355"
uuid = "5c2747f8-b7ea-4ff2-ba2e-563bfd36b1d4"
version = "1.3.0"

[[UUIDs]]
deps = ["Random", "SHA"]
uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"

[[Unicode]]
uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"

[[Vega]]
deps = ["DataStructures", "DataValues", "Dates", "FileIO", "FilePaths", "IteratorInterfaceExtensions", "JSON", "JSONSchema", "MacroTools", "NodeJS", "Pkg", "REPL", "Random", "Setfield", "TableTraits", "TableTraitsUtils", "URIParser"]
git-tree-sha1 = "43f83d3119a868874d18da6bca0f4b5b6aae53f7"
uuid = "239c3e63-733f-47ad-beb7-a12fde22c578"
version = "2.3.0"

[[VegaDatasets]]
deps = ["DataStructures", "DataValues", "FilePaths", "IterableTables", "IteratorInterfaceExtensions", "JSON", "TableShowUtils", "TableTraits", "TableTraitsUtils", "TextParse"]
git-tree-sha1 = "c997c7217f37205c5795de8c797f8f8531890f1d"
uuid = "0ae4a718-28b7-58ec-9efb-cded64d6d5b4"
version = "2.1.1"

[[VegaLite]]
deps = ["Base64", "DataStructures", "DataValues", "Dates", "FileIO", "FilePaths", "IteratorInterfaceExtensions", "JSON", "MacroTools", "NodeJS", "Pkg", "REPL", "Random", "TableTraits", "TableTraitsUtils", "URIParser", "Vega"]
git-tree-sha1 = "3e23f28af36da21bfb4acef08b144f92ad205660"
uuid = "112f6efa-9a02-5b7d-90c0-432ed331239a"
version = "2.6.0"

[[WeakRefStrings]]
deps = ["DataAPI", "Parsers"]
git-tree-sha1 = "4a4cfb1ae5f26202db4f0320ac9344b3372136b0"
uuid = "ea10d353-3f73-51f8-a26c-33c1cb351aa5"
version = "1.3.0"

[[Zlib_jll]]
deps = ["Libdl"]
uuid = "83775a58-1f1d-513f-b197-d71354ab007a"

[[nghttp2_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d"

[[p7zip_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0"
"""

# ╔═╡ Cell order:
# ╠═fe20d728-2096-11ec-2740-4368d1340f88
# ╠═a358b3ab-b1d4-471e-a9be-565e7ec339de
# ╠═94762b1c-4f03-43b9-b288-4f4b77894d3b
# ╟─75bc7766-500b-45fc-b5dc-af144c17d947
# ╟─592f98f7-5cc1-448b-aa7d-0d7a0be98835
# ╟─f2bfd623-223b-4249-a57a-924f359cbede
# ╠═d84d3340-9fec-470e-9717-4261c5c0cce2
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


text:
module ProxyInterface
using Compat
export @forward

"""
    @forward T.x functions...
Define methods for `functions` on type `T`, which call the relevant function
on the field `x`.
# Example
```julia
struct Wrapper
    x
end
@forward Wrapper.x  Base.sqrt                                  # now sqrt(Wrapper(4.0)) == 2.0
@forward Wrapper.x  Base.length, Base.getindex, Base.iterate   # several forwarded functions are put in a tuple
@forward Wrapper.x (Base.length, Base.getindex, Base.iterate)  # equivalent to above
```
"""
macro forward(ex, fs)
  check = ex isa Expr && ex.head == :. && ex.args[1] isa Symbol && ex.args[2] isa QuoteNode
  check || error("Syntax: @forward T.x f, g, h")
  T = esc(ex.args[1])
  field = ex.args[2].value  # no escape needed for . syntax
  fs = (fs isa Expr && fs.head == :tuple) ? map(esc, fs.args) : [esc(fs)]
  :($([:($f(x::$T, args...; kwargs...) = (Base.@_inline_meta; $f(x.$field, args...; kwargs...)))
       for f in fs]...);
    nothing)
end

function iterator end
macro iterator(ProxyType)
  ProxyType = esc(ProxyType)
  quote
    Base.iterate(p::$ProxyType) = Base.iterate(ProxyInterface.iterator(p))
    Base.iterate(p::$ProxyType, state) = Base.iterate(ProxyInterface.iterator(p), state)
    Base.IteratorSize(A::Type{<:$ProxyType}) = Base.IteratorSize(ProxyInterface.iterator(A))
    Base.IteratorEltype(A::Type{<:$ProxyType}) = Base.IteratorEltype(ProxyInterface.iterator(A))
    Base.eltype(p::$ProxyType) = Base.eltype(ProxyInterface.iterator(p))
    Base.length(p::$ProxyType) = Base.length(ProxyInterface.iterator(p))
    Base.size(p::$ProxyType) = Base.size(ProxyInterface.iterator(p))
    Base.size(p::$ProxyType, d) = Base.size(ProxyInterface.iterator(p), d)
    Base.axes(p::$ProxyType) = Base.axes(ProxyInterface.iterator(p))  # analog to Base.Generator
    Base.ndims(p::$ProxyType) = Base.ndims(ProxyInterface.iterator(p))
    # we include foreach to iterator, as it directly corresponds to a for loop
    Base.foreach(f, p::$ProxyType) = Base.foreach(f, ProxyInterface.iterator(p))
    # we don't include Base.map, as it requires constructing the Proxy, which we don't know how to do it
  end
end

function indexing end
macro indexing(ProxyType)
  ProxyType = esc(ProxyType)
  quote
    Base.getindex(p::$ProxyType, args...)	= Base.getindex(ProxyInterface.indexing(p), args...)
    Base.firstindex(p::$ProxyType) = Base.firstindex(ProxyInterface.indexing(p))
    Base.lastindex(p::$ProxyType)	= Base.lastindex(ProxyInterface.indexing(p))
  end
end
macro indexing_mutable(ProxyType)
  ProxyType = esc(ProxyType)
  quote
    @indexing($ProxyType)
    Base.setindex!(p::$ProxyType, args...) = Base.setindex!(ProxyInterface.indexing(p), args...)
  end
end

function array end
macro array(ProxyType)
  ProxyType = esc(ProxyType)
  quote
    Base.size(p::$ProxyType) = Base.size(ProxyInterface.array(p))
    Base.getindex(p::$ProxyType, args...)	= Base.getindex(ProxyInterface.array(p), args...)
    Base.firstindex(p::$ProxyType) = Base.firstindex(ProxyInterface.array(p))
    Base.lastindex(p::$ProxyType)	= Base.lastindex(ProxyInterface.array(p))
    # optional
    Base.IndexStyle(A::Type{<:$ProxyType}) = Base.IndexStyle(ProxyInterface.array(A))
    Base.length(p::$ProxyType) = Base.length(ProxyInterface.array(p))
    Base.similar(p::$ProxyType, args...) = Base.similar(ProxyInterface.array(p), args...)
    # non-traditional
    Base.axes(p::$ProxyType) = Base.axes(ProxyInterface.array(p))
    # strides
    Base.strides(p::$ProxyType, args...) = Base.strides(ProxyInterface.array(p), args...)
    Base.unsafe_convert(p::$ProxyType, args...) = Base.strides(ProxyInterface.array(p), args...)
    Base.unsafe_convert(t::Type{Ptr{T}}, p::$ProxyType) where T = Base.unsafe_convert(t, ProxyInterface.array(p))
    # extra
    Base.in(a, p::$ProxyType) = Base.in(a, ProxyInterface.array(p))
  end
end

macro array_mutable(ProxyType)
  ProxyType = esc(ProxyType)
  quote
    @array($ProxyType)
    Base.setindex!(p::$ProxyType, args...) = Base.setindex!(ProxyInterface.array(p), args...)
    Base.empty!(p::$ProxyType) = Base.empty!(ProxyInterface.array(p))
  end
end

function dict end
macro dict(ProxyType)
  ProxyType = esc(ProxyType)
  quote
    Base.getindex(p::$ProxyType, args...)	= Base.getindex(ProxyInterface.dict(p), args...)
    Base.length(p::$ProxyType) = Base.length(ProxyInterface.dict(p))
    Base.in(a, p::$ProxyType) = Base.in(a, ProxyInterface.dict(p))
    Base.haskey(p::$ProxyType, key) = Base.haskey(ProxyInterface.dict(p), key)
    Base.keys(p::$ProxyType) = Base.keys(ProxyInterface.dict(p))
    Base.values(p::$ProxyType) = Base.values(ProxyInterface.dict(p))
  end
end

macro dict_mutable(ProxyType)
  ProxyType = esc(ProxyType)
  quote
    @dict($ProxyType)
    Base.setindex!(p::$ProxyType, args...) = Base.setindex!(ProxyInterface.dict(p), args...)
    Base.empty!(p::$ProxyType) = Base.empty!(ProxyInterface.dict(p))
  end
end
end # module


text:
for pkg in ["Utils"]
    ispath(Pkg.dir(pkg)) && (Pkg.checkout(pkg); continue)
    Pkg.clone("https://github.com/AStupidBear/$pkg.jl.git")
    Pkg.build(pkg)
end


text:
include("readData.jl")
include("readSparseData.jl")
include("saveSparseData.jl")
include("stochasticCVB0.jl")

#This file shows a demo of the SCVB0 algorithm on the NIPS data.

numWords = 13649
numTopics = 50
minibatchSize = 20
numDocumentIterations = 50000
time_limit = 120.0
burnIn = 1
alpha = 0.1
eta = 0.01
tau = 1000.0
kappa = 0.9
tau2= 1.0
kappa2 = 0.9
scale = 100.0
saveParams = true
saveCount = 100

dataFilename = "data/NIPS.txt";                 #one line per document, space-separated one-based dictionary indices for each consecutive word in the document.
sparseDataFilename = "data/NIPSsparse.txt";     #sparse file format, obtained from the dense data format via saveSparseData("data/NIPSsparse.txt", documents);
dictionaryFilename = "data/NIPSdict.txt";       #one line per dictionary word

readSparseDataset = true; #This script demonstrates reading either the non-sparse file format or the sparse file format, selected here

if readSparseDataset
    documents, dictionary = readSparseData(sparseDataFilename, dictionaryFilename);
else
    documents, dictionary = readData(dataFilename, dictionaryFilename)
    println("Converting to sparse format...")
    documents = wordVectorToSparseCounts(documents);
    println("Done.")
end

tic();
VBparams,saved_iters, saved_time, saved_topics, saved_alphas  = stochasticCVB0(documents,numWords,numTopics,
                              numDocumentIterations,time_limit,burnIn,minibatchSize,alpha,eta,
                              dictionary,tau,kappa,tau2,kappa2,scale,saveParams,saveCount);
t = toc();

phi = computeTopics(VBparams); #compute the normalized topic distributions


text:
# Wandb Artifacts
mutable struct WandbArtifact
    artifact::PyObject
end

WandbArtifact(args...; kwargs...) = WandbArtifact(wandb.Artifact(args...; kwargs...))

Base.show(io::IO, ::WandbArtifact) = Base.print(io, "WandbArtifact()")

for func in (:download, :get, :finalize, :wait)
    @eval begin
        Base.$(func)(wa::WandbArtifact, args...; kwargs...) = wa.artifact.$(func)(args...; kwargs...)
    end
end

for func in
    (:add, :add_file, :add_dir, :add_reference, :checkout, :delete, :get_added_local_path_name, :get_path, :logged_by,
     :new_file, :used_by, :verify, :save)
    @eval begin
        $(func)(wa::WandbArtifact, args...; kwargs...) = wa.artifact.$(func)(args...; kwargs...)
    end
end

Base.log(::WandbLogger, wa::WandbArtifact) = wandb.log_artifact(wa.artifact)


text:
##############################################################
"""
$(SIGNATURES)

Calculation of Bernoulli function via Horner scheme based on Taylor
coefficients around 0.
"""

function bernoulli_horner(x)
    y=x/47_900_160
    y=x*y
    y=x*(-1/1_209_600+y)
    y=x*y
    y=x*(1/30_240+y)
    y=x*y
    y=x*(-1/720+y)
    y=x*y
    y=x*(1/12+y)
    y=x*(-1/2+y)
    y=1+y
end

# Bernoulli thresholds optimized for Float64
const bernoulli_small_threshold=0.25
const bernoulli_large_threshold=40.0
##############################################################
"""
$(SIGNATURES)

Bernoulli function ``B(x)=\\frac{x}{e^x-1}`` for exponentially
fitted upwinding.

The name `fbernoulli` has been chosen to avoid confusion
with Bernoulli from JuliaStats/Distributions.jl

Returns a real number containing the result.
"""
function fbernoulli(x)
    if x<-bernoulli_large_threshold
        -x
    elseif x>bernoulli_large_threshold
	zero(x)
    else
        expx=exp(x)
        expxm1=expx-1.0
        if abs(expxm1)>bernoulli_small_threshold
            x/expxm1
        else
            bernoulli_horner(x)
        end
    end
end

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

"""
$(SIGNATURES)

Bernoulli function ``B(x)=\\frac{x}{e^x-1}`` for exponentially
fitted upwind, joint evaluation for positive and negative
argument

Usually, we need ``B(x), B(-x)`` togehter,
and it is cheaper to calculate them together.

Returns two real numbers containing the result for argument
`x` and argument `-x`.

The error in comparison with the evaluation of the original expression
with BigFloat is less than 1.0e-15
"""
function fbernoulli_pm(x)
    if x< -bernoulli_large_threshold
        return -x, zero(x)
    elseif x>bernoulli_large_threshold
	return zero(x),x
    else
        expx=exp(x)
        expxm1=expx-1.0
        if abs(expxm1)>bernoulli_small_threshold
            bp=x/expxm1
            bm=x/(1.0-1.0/expx)
            return bp,bm
        else
            y=bernoulli_horner(x)
            return y,x+y
        end
    end
end


text:
module BALUtils

using StaticArrays
using SparseArrays
using LinearAlgebra
using Rotations

export Camera, pose, axisangle, BA, readbal, visibility_graph, restrict, writebal, center, AngleAxis, num_cameras, num_points, num_observations

function rodrigues_rotate(angle_axis :: AbstractArray, point :: AbstractArray)
    theta = norm(angle_axis)
    axis = normalize(angle_axis)
    point * cos(theta) + cross(axis, point) * sin(theta) + axis * dot(axis, point) * (1 - cos(theta))
end

"""
Rotation represented as a 3 vector who's norm is the angle to rotate around the
normalized vector.
"""
struct AngleAxis{T <: Number}
    angleaxis :: SVector{3, T}
end

function AngleAxis(mat :: Array{<: Number, 2})
    AngleAxis(axisangle(mat))
end

function Base.:*(v :: AngleAxis, x :: AbstractArray)
    rodrigues_rotate(v.angleaxis, x)
end

# from https://math.stackexchange.com/questions/382760/composition-of-two-axis-angle-rotations
function Base.:*(x :: AngleAxis, y :: AngleAxis)
    alpha = norm(x.angleaxis)
    beta = norm(y.angleaxis)
    ax = normalize(x.angleaxis)
    ay = normalize(y.angleaxis)
    gamma = 2 * acos(cos(alpha/2) * cos(beta/2) - sin(alpha/2) * sin(beta/2) * dot(ax, ay))
    axis = (sin(alpha/2) * cos(beta/2) * ax + cos(alpha/2) * sin(beta/2) * ay - sin(alpha/2) * sin(beta/2)  * cross(ax, ay)) / sin(gamma/2)
    AngleAxis(axis * gamma)
end

Base.transpose(v :: AngleAxis) = AngleAxis(-v.angleaxis)
Base.adjoint(v :: AngleAxis) = transpose(v)

struct Camera{T}
    pose :: SVector{3,T}
    rotation :: SVector{3,T}
    intrinsics :: SVector{3,T}
end

"""
Construct a camera from a rotation, translation, instrinsics vector.
"""
function Camera(x :: AbstractArray)
    Camera(SVector{3,Float64}(x[4:6]), SVector{3,Float64}(x[1], x[2], x[3]), SVector{3,Float64}(x[7:9]))
end

rot(c :: Camera) = RodriguesVec(c.rotation...)
pose(c :: Camera) = c.pose
center(c :: Camera) = -(rot(c)' * pose(c))
Base.vec(c :: Camera) = vcat(c.rotation, c.pose, c.intrinsics)

function axisangle(x :: Array{Float64, 2})
    u = [x[3,2]-x[2,3], x[3,1]-x[1,3], x[2,1] - x[1,2]]
    angle = acos((trace(x)-1)/2)
    SVector{3,Float64}((u / norm(u) * 2 * sin(angle))...)
end

struct BA
    cameras :: Vector{Camera}
    observations :: Vector{Vector{Tuple{Int64,Float64,Float64}}}
    points :: Vector{SVector{3,Float64}}
end

function restrict(ba :: BA, inds; ignore_points=false, vis_thresh=2)
    cams = ba.cameras[inds]
    obs = ba.observations[inds]
    obs, points = if !ignore_points
        # get all points still visible
        ps = vcat(map(x -> map(y -> y[1], x), obs)...)
        # remove points visible by only one camera (or threshold)
        ps = filter(x -> count(i -> i == x, ps) >= vis_thresh, unique(ps))
        # remap point ids
        d = Dict(zip(ps, 1:length(ps)))
        points = ba.points[ps]
        obs = map(obs) do ob
            vs = []
            for (i, x, y) in ob
                # drop points
                if i in keys(d)
                    push!(vs, (d[i], x, y))
                end
            end
            vs
        end
        obs, points
    else
        obs, ba.points
    end
    BA(cams, obs, points)
end

num_cameras(ba :: BA) = length(ba.cameras)
num_points(ba :: BA) = length(ba.points)
num_observations(ba :: BA) = sum(length.(ba.observations))

function Base.show(io :: IO, ba :: BA)
    print(io, "Bundle adjustment problem with $(length(ba.cameras)) cameras, $(length(ba.points)) points, $(sum(length.(ba.observations))) observations")
end

function readbal(filename)
    isfile(filename) || SystemError("$filename is not a file")
    if splitext(filename)[2] == ".bbal"
        open(filename) do f
            num_cameras = read(f, UInt64) |> ntoh
            num_points = read(f, UInt64) |> ntoh
            num_observations = read(f, UInt64) |> ntoh

            obs = map(1:num_cameras) do i
                nobs = read(f, UInt64) |> ntoh
                obs = Vector{Tuple{Int64,Float64,Float64}}()
                sizehint!(obs, nobs)
                for j in 1:nobs
                    p = read(f, UInt64) |> ntoh
                    u = read(f, Float64) |> ntoh
                    v = read(f, Float64) |> ntoh
                    push!(obs, (p+1, u, v))
                end
                obs
            end

            cameras = map(1:num_cameras) do i
                ary = Array{Float64}(undef, 9)
                read!(f, ary)
                Camera(ntoh.(ary))
            end

            points = map(1:num_points) do i
                ary = Array{Float64}(undef, 3)
                read!(f, ary)
                SVector{3,Float64}(ntoh.(ary))
            end

            BA(cameras, obs, points)
        end
    else
        f = open(filename)
        line = readline(f)
        (num_cameras, num_points, num_observations) = map(x->parse(Int,x), split(line))

        obs = [Vector{Tuple{Int64,Float64,Float64}}() for _ in  1:num_cameras]

        for i in 1:num_observations
            (c_ind_, p_ind_, x_, y_) = split(readline(f))
            c_ind = parse(Int, c_ind_)+1
            p_ind = parse(Int, p_ind_)+1
            x = parse(Float64, x_)
            y = parse(Float64, y_)
            push!(obs[c_ind], (p_ind, x, y))
        end

        cameras = map(1:num_cameras) do i
            t = map(1:9) do j
                parse(Float64, readline(f))
            end
            Camera(t)
        end

        points = map(1:num_points) do i
            t = map(1:3) do j
                parse(Float64, readline(f))
            end
            SVector{3,Float64}(t...)
        end

        BA(cameras, obs, points)
    end
end

function visibility_graph(ba :: BA)
    is = []
    sizehint!(is, num_observations(ba))
    js = []
    sizehint!(js, num_observations(ba))
    for (i, obs) in enumerate(ba.observations)
        for (j, _, _) in obs
            push!(is, i)
            push!(js, j)
        end
    end
    m = sparse(is, js, ones(Int, length(js)))
    m * m'
end

function writebal(filepath :: AbstractString, ba :: BA)
    open(filepath, "w") do f
        if splitext(filepath)[2] == ".bbal"
            write(f, UInt64(num_cameras(ba)) |> hton)
            write(f, UInt64(num_points(ba)) |> hton)
            write(f, UInt64(num_observations(ba)) |> hton)

            for (cam, obs) in enumerate(ba.observations)
                write(f, UInt64(length(obs)) |> hton)
                for (point, x, y) in obs
                    write(f, UInt64(point - 1) |> hton)
                    write(f, Float64(x) |> hton)
                    write(f, Float64(y) |> hton)
                end
            end

            for c in ba.cameras
                write(f, hton.(Float64.(vec(c))))
            end

            for p in ba.points
                write(f, hton.(Float64.(p)))
            end
        else
            write(f, "$(num_cameras(ba)) $(num_points(ba)) $(num_observations(ba))\n")
            for (cam, obs) in enumerate(ba.observations)
                for (point, x, y) in obs
                    # zero indexed
                    write(f, "$(cam - 1) $(point - 1) $x $y\n")
                end
            end
            for cam in ba.cameras
                for v in vec(cam)
                    write(f, "$v ")
                end
                write(f, "\n")
            end
            for point in ba.points
                write(f, "$(point[1]) $(point[2]) $(point[3])\n")
            end
        end
    end
end

end # module