text:
using WhereTheWaterFlows
const WWF = WhereTheWaterFlows
using Test

# test examples
module Test_Examples # use a module to avoid name-space pollution
plotyes = false
for fl in readdir("../examples/")
    include(joinpath("../examples/", fl))
end
end

"An artificial DEM"
function dem1(x, y; withpit=false, randfac=0.0)
    out = - (x^2 - 1)^2 - (x^2*y - x - 1)^2 + 6 + 0.1*x + 3*y
    if withpit
        out -= 2*exp(-(x^2 + y^2)*50)
    end
    out += randfac*randn()
    return out<0 ? 0.0 : out
end

const _offset = [0.2, 0.5]

"DEM with a few maxs and mins"
function peaks(n=100)
    coords = range(-pi+_offset[1], pi-_offset[2], length=n)
    return coords, sin.(coords) .* cos.(coords')
end

"DEM with a few more maxs and mins"
function peaks2(n=100, randfac=0.0)
    coords = range(-pi+_offset[1], pi-_offset[2], length=n)
    return coords, sin.(coords) .* cos.(coords') .-
        0.7*(sin.(coords.+1) .* cos.(coords')).^8 .+
        0.01 * (coords.+coords') .+
        randfac*randn(n,n) # 0.02
end

"DEM with a few more maxs and mins.  NaN masked middle."
function peaks2_nan(n=100, randfac=0.0)
    coords = range(-pi+_offset[1], pi-_offset[2], length=n)
    dem = sin.(coords) .* cos.(coords') .-
        0.7*(sin.(coords.+1) .* cos.(coords')).^8 .+
        0.01 * coords.*coords' .+
        randfac*randn(n,n) # 0.02
    dem[(n÷2:n÷2+2).-8, n÷2+1:n÷2+4] .= NaN
    return coords, dem
end

"DEM with a few more maxs and mins.  NaN masked edge."
function peaks2_nan_edge(n=100, randfac=0.0)
    coords = range(-pi+_offset[1], pi-_offset[2], length=n)
    dem = sin.(coords) .* cos.(coords') .-
        0.7*(sin.(coords.+1) .* cos.(coords')).^8 .+
        0.01 * (coords.+coords') .+
        randfac*randn(n,n) # 0.02
    outflow = dem[n÷2, 1]
    dem[1, :] .= NaN
    dem[end, :] .= NaN
    dem[:, 1] .= NaN
    dem[:, end] .= NaN
    dem[n÷2, 1] = outflow
    return coords, dem
end

function dem_one_point()
    dem = zeros(4,5) .* NaN
    dem[2,2] = 1.0
    coords = range(0, 1.0, length=size(dem,1)), range(0, 1.0, length=size(dem,2))
    return coords, dem
end

function dem_two_points()
    coords, dem = dem_one_point()
    dem[3,4] = 1.0
    return coords, dem
end

function dem_patho1()
    dem = [NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN; NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN; NaN NaN NaN NaN NaN 224973.93798828125 NaN NaN NaN NaN NaN NaN NaN NaN NaN; NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN; NaN NaN NaN NaN NaN NaN 1.5264895434570312e6 780071.6101074219 611466.6186523438 68967.85522460938 NaN NaN NaN NaN NaN; NaN NaN NaN NaN NaN NaN NaN 340947.7197265625 1.7224892822265625e6 1.503103447265625e6 646871.5270996094 NaN NaN NaN NaN; NaN NaN NaN NaN NaN 493601.3024902344 NaN 953505.1794433594 1.33096515625e6 251280.46264648438 689901.5191650391 NaN NaN NaN NaN; NaN NaN NaN 1.210690463256836e6 2.2551471240234375e6 1.620159921875e6 1.9493891125488281e6 2.07443427734375e6 2.1659396948242188e6 10309.342041015625 NaN NaN NaN NaN NaN; NaN NaN NaN 1.4056021618652344e6 2.9576955078125e6 2.6524201977539062e6 2.3390813134765625e6 2.8080285302734375e6 2.6852572412109375e6 1.7304499780273438e6 NaN NaN NaN NaN NaN; NaN NaN NaN NaN 3.1962603833007812e6 3.308844228515625e6 3.3650431396484375e6 3.518876923828125e6 2.88937056640625e6 2.29438943359375e6 2.190683447265625e6 1.920630322265625e6 NaN NaN NaN; NaN NaN NaN NaN 2.521603028564453e6 3.1732824658203125e6 2.8601568969726562e6 3.4935072900390625e6 3.284232265625e6 2.7940991479492188e6 2.621620751953125e6 1.9656544458007812e6 NaN NaN NaN; NaN NaN NaN NaN 2.1466443212890625e6 2.5994008081054688e6 997498.37890625 3.0142672705078125e6 3.1849692797851562e6 2.6993770849609375e6 2.7337567944335938e6 2.1765316455078125e6 NaN NaN NaN; NaN NaN NaN NaN 1.4447062829589844e6 1.6352845727539062e6 297836.35498046875 2.5211809497070312e6 2.6398626635742188e6 2.2131605932617188e6 1.3856159045410156e6 NaN NaN NaN NaN; NaN NaN NaN NaN NaN NaN NaN NaN 932178.7219238281 NaN NaN NaN NaN NaN NaN; NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN]
    coords = range(0, 1.0, length=size(dem,1)), range(0, 1.0, length=size(dem,2))
    return coords, dem
end

function dem_patho2()
    dem = [2.4669351831054688e6 2.4912734814453125e6 2.4511637768554688e6 2.297378525390625e6 2.1467358154296875e6 2.0279176049804688e6 1.8322703247070312e6 1.6419496728515625e6 1.3741556689453125e6 985379.6069335938 170602.13500976562; 2.6378740991210938e6 2.621620751953125e6 2.4854557543945312e6 2.4079890673828125e6 2.2196343212890625e6 1.9963872729492188e6 1.8812639379882812e6 1.6405933081054688e6 1.3537354772949219e6 937561.1877441406 NaN; 2.7794645190429688e6 2.6392191967773438e6 2.6666443579101562e6 2.4175318090820312e6 2.2348992724609375e6 2.0938221411132812e6 1.7992342114257812e6 1.6129501953125e6 1.1581694873046875e6 650434.7094726562 NaN; 2.7840773608398438e6 2.764429228515625e6 2.6759121826171875e6 2.5908404150390625e6 2.393361572265625e6 2.2053372387695312e6 2.128048779296875e6 1.8252124182128906e6 1.2706359448242188e6 820004.8461914062 NaN; 2.818317265625e6 2.876910830078125e6 2.7989625903320312e6 2.6884156298828125e6 2.3920689697265625e6 2.2358143334960938e6 2.211969111328125e6 1.9571274853515625e6 1.3319587353515625e6 864203.271484375 NaN; 2.7831515209960938e6 2.7405954125976562e6 2.728760107421875e6 2.5505696166992188e6 2.4580276049804688e6 2.0936141577148438e6 2.1094039184570312e6 1.9852256274414062e6 1.5260443835449219e6 868576.2139892578 NaN; 2.6441077563476562e6 2.666107099609375e6 2.6324500732421875e6 2.3455210571289062e6 2.2824994604492188e6 2.2010346118164062e6 2.149140135498047e6 1.762273427734375e6 1.3110776196289062e6 NaN NaN; 2.4332154907226562e6 2.3887434301757812e6 2.3717401684570312e6 2.1986926538085938e6 1.96104419921875e6 1.7805842163085938e6 1.8096322509765625e6 1.395735947265625e6 682117.8509521484 NaN NaN; 2.1785780810546875e6 2.1734258032226562e6 2.1239198754882812e6 1.974646943359375e6 1.6406056872558594e6 1.2984880700683594e6 1.258363583984375e6 726651.7492675781 NaN NaN NaN; 1.8946451831054688e6 1.81941337890625e6 1.7843751513671875e6 1.7247107556152344e6 1.2050103857421875e6 93141.61499023438 535295.9533691406 NaN NaN NaN NaN; 1.5536723583984375e6 1.3856159045410156e6 1.3403697265625e6 1.2850621618652344e6 793587.7758789062 NaN NaN NaN NaN NaN NaN]
    coords = range(0, 1.0, length=size(dem,1)), range(0, 1.0, length=size(dem,2))
    return coords, dem
end

# Test low-level fns

# flowsinto
# """
# Tests whether a cell `J` with flowdir `dirJ` flows into cell `I`.
# """
# flowsinto(J::CartesianIndex, dirJ, I::CartesianIndex) = ind2dir(I-J) == dirJ
@testset "flowsinto" begin
    J = CartesianIndex(5,6)
    dirs = 0:9
    for i=[3,7], j=4:8, dir in dirs
        @test_throws BoundsError WWF.flowsinto(J, dir, CartesianIndex(i,j))
    end
    for i=3:7, j=[4,8], dir in dirs
        @test_throws BoundsError WWF.flowsinto(J, dir, CartesianIndex(i,j))
    end
    trueflow = Dict(1 => J+CartesianIndex(-1,-1),
                    2 => J+CartesianIndex(0,-1),
                    3 => J+CartesianIndex(1,-1),
                    4 => J+CartesianIndex(-1,0),
                    5 => J+CartesianIndex(0,0),
                    6 => J+CartesianIndex(1,0),
                    7 => J+CartesianIndex(-1,1),
                    8 => J+CartesianIndex(0,1),
                    9 => J+CartesianIndex(1,1)
                    )
    for (td, I) in trueflow
        for d in dirs
            if d==td
                @test WWF.flowsinto(J, d, I)
            else
                @test !WWF.flowsinto(J, d, I)
            end
        end
    end
end

@test WWF.dirnums==reverse([ 7 8 9
                             4 5 6
                             1 2 3]', dims=2)
tmp = convert(Matrix, WWF.dirnums)

for i=1:9
    dem = ones(3,3)
    dem[i] = 0
    dir, nout, nin, pits = WWF.d8dir_feature(dem, false)

    @test dir[2,2] == WWF.dirnums[i]
    @test dem[WWF.dir2ind(dir[2,2])+CartesianIndex(2,2)]==0

    @test WWF.d8dir_feature(dem, false)[1][2,2] == i
    @test maximum(nout)<=1
    @test minimum(nout)>=0
    @test maximum(nin)<=8
    @test minimum(nin)>=0
    @test sum(nout)==sum(nin)

    @test sum(dem[WWF.iterate_D9(CartesianIndex(2,2), dem)]) == 8

end

@testset "dir2ind" begin
    @test WWF.dir2ind(1) == CartesianIndex(-1,-1)
    @test WWF.dir2ind(2) == CartesianIndex(0,-1)
    @test WWF.dir2ind(3) == CartesianIndex(1,-1)
    @test WWF.dir2ind(4) == CartesianIndex(-1,0)
    @test WWF.dir2ind(5) == CartesianIndex(0,0)
    @test WWF.dir2ind(6) == CartesianIndex(1,0)
    @test WWF.dir2ind(7) == CartesianIndex(-1,1)
    @test WWF.dir2ind(8) == CartesianIndex(0,1)
    @test WWF.dir2ind(9) == CartesianIndex(1,1)
    for i=1:9
        @test i==WWF.ind2dir(WWF.dir2ind(i))
    end
end

@testset "diagonal gradient adjustment" begin
    dem = ones(4,4)
    dem[2,2] = 0.5
    dem[2,3] = 0.45
    area, slen, dir, nout, nin, pits, c, bnds = WWF.waterflows(dem);
    @test dir==[5  5  5  5
                5  8  1  5
                5  4  4  5 # without adjustment this row would be '5  7  4  5'
                5  5  5  5]
end

# Test high-level functions
@testset "DEM: dem1" begin
    dx = 0.9
    xs = -1.5:dx:1
    ys = -0.5:dx:3.0
    dem = dem1.(xs, ys')
    @test size(dem)==(length(xs), length(ys))

    area, slen, dir, nout, nin, pits, c, bnds = WWF.waterflows(dem, bnd_as_pits=false)
    @test area == [1.0 1.0 4.0 2.0; 1.0 1.0 1.0 1.0; 5.0 3.0 2.0 1.0]
    @test slen == [1 1 2 2; 1 1 1 1; 4 3 2 1]
    @test dir == Int8[5 8 5 5; 6 7 4 4; 5 2 2 2]
    @test nout == Bool[false true false false; true true true true; false true true true]
    @test nin == Int8[0 0 3 1; 0 0 0 0; 2 1 1 0]
    @test pits ==CartesianIndex{2}[CartesianIndex(1, 1), CartesianIndex(3, 1), CartesianIndex(1, 3), CartesianIndex(1, 4)]

    # non default cellarea
    area, _ = WWF.waterflows(dem, fill(10.0, size(dem)), bnd_as_pits=false)
    @test area == [10.0 10.0 40.0 20.0; 10.0 10.0 10.0 10.0; 50.0 30.0 20.0 10.0]
end

@testset "DEM: peaks2" begin
    xs, dem = peaks2()
    ys = xs
    @test size(dem)==(length(xs), length(ys))
    area, slen, dir, nout, nin, pits, c, bnds = WWF.waterflows(dem, drain_pits=false, bnd_as_pits=false);
    #plotarea_dem(xs, ys, dem, area, pits)
    @test length(pits) == 6
    @test maximum(slen)==86
    @test maximum(area)==4759
    @test length(unique(c))==6
    @test sum(c) == 33862
    @test sum(diff(c[:])) == 5
    @test sort(unique(c))[[1,end]] ==[1,6]
    @test all([c[pits[cc]]==cc  for cc=1:length(pits)]) # pit in catchment of same color

    area, slen, dir, nout, nin, pits, c, bnds = WWF.waterflows(dem, drain_pits=true, bnd_as_pits=false)
    # plotarea_dem(xs, ys, dem, area, pits)
    @test length(pits) == 4
    @test maximum(slen)==129
    @test maximum(area)==7702
    @test length(unique(c))==4
    @test sum(c) == 23855
    @test sum(diff(c[:])) == 3
    @test sort(unique(c))[[1,end]] ==[1,4]
    @test bnds isa Array{Array{CartesianIndex{2},1},1}
    @test all([c[pits[cc]]==cc  for cc=1:length(pits)]) # pit in catchment of same color
end

@testset "DEM: peaks2_nan" begin
    xs, dem = peaks2_nan()
    ys = xs
    nanlocs = findall(isnan.(dem))
    @test size(dem)==(length(xs), length(ys))

    area, slen, dir, nout, nin, pits, c, bnds = WWF.waterflows(dem, drain_pits=false, bnd_as_pits=false);
    # plotarea_dem(xs, ys, dem, area, pits)
    @test length(pits) == 7
    @test maximum(slen)==87
    @test maximum(area[.!isnan.(area)])==4725
    @test length(unique(c))==8
    @test sum(c) == 36447
    @test sum(diff(c[:])) == 6
    @test sort(unique(c))[[1,end]] ==[0,7]
    @test all(c[nanlocs].==0)
    @test all([c[pits[cc]]==cc  for cc=1:length(pits)]) # pit in catchment of same color

    area, slen, dir, nout, nin, pits, c, bnds = WWF.waterflows(dem, drain_pits=false, bnd_as_pits=true);
    # plotarea_dem(xs, ys, dem, area, pits)
    @test length(pits) == 408
    @test maximum(slen)==66
    @test maximum(area[.!isnan.(area)])==4457
    @test length(unique(c))==409
    @test sum(c) == 2079380
    @test sum(diff(c[:])) == 407
    @test sort(unique(c))[[1,end]] ==[0,408]
    @test all(c[ [n for n in nanlocs if !(n in pits)]].==0)
    @test all([c[pits[cc]]==cc  for cc=1:length(pits)]) # pit in catchment of same color

    # fill pits
    area, slen, dir, nout, nin, pits, c, bnds = WWF.waterflows(dem, drain_pits=true, bnd_as_pits=false);
    # plotarea_dem(xs, ys, dem, area, pits)
    @test length(pits) == 4
    @test maximum(slen)==142
    @test maximum(area[.!isnan.(area)])==7706
    @test length(unique(c))==5
    @test sum(c) == 23820
    @test sum(diff(c[:])) == 3
    @test sort(unique(c))[[1,end]] ==[0,4]
    @test bnds isa Array{Array{CartesianIndex{2},1},1}
    @test all(c[nanlocs].==0)
    @test all([c[pits[cc]]==cc  for cc=1:length(pits)]) # pit in catchment of same color

    area, slen, dir, nout, nin, pits, c, bnds = WWF.waterflows(dem, drain_pits=true, bnd_as_pits=true);
    # plotarea_dem(xs, ys, dem, area, pits)
    @test length(pits) == 406
    @test maximum(slen)==75
    @test maximum(area[.!isnan.(area)])==4458
    @test length(unique(c))==407
    @test sum(c) == 2067155
    @test sum(diff(c[:])) == 405
    @test sort(unique(c))[[1,end]] ==[0,406]
    @test bnds isa Array{Array{CartesianIndex{2},1},1}
    @test all(c[ [n for n in nanlocs if !(n in pits)]].==0)
    @test all([c[pits[cc]]==cc  for cc=1:length(pits)]) # pit in catchment of same color
end

@testset "DEM: peaks2_nan_edge" begin
    xs, dem = peaks2_nan_edge()
    ys = xs
    @test size(dem)==(length(xs), length(ys))

    area, slen, dir, nout, nin, pits, c, bnds = WWF.waterflows(dem, drain_pits=false, bnd_as_pits=false);
    # plotarea_dem(xs, ys, dem, area, pits)
    @test length(pits) == 6
    @test maximum(slen)==84
    @test maximum(area[.!isnan.(area)])==4658
    @test length(unique(c))==7
    @test sum(c) == 32486
    @test sum(diff(c[:])) == 0
    @test sort(unique(c))[[1,end]] ==[0,6]
    @test all([c[pits[cc]]==cc  for cc=1:length(pits)]) # pit in catchment of same color

    area, slen, dir, nout, nin, pits, c, bnds = WWF.waterflows(dem, drain_pits=false, bnd_as_pits=true);
    # plotarea_dem(xs, ys, dem, area, pits)
    @test length(pits) == 390
    @test maximum(slen)==64
    @test maximum(area[.!isnan.(area)])==4434
    @test length(unique(c))==391
    @test sum(c) == 2000894
    @test sum(diff(c[:])) == -1
    @test sort(unique(c))[[1,end]] ==[0,390]
    @test all([c[pits[cc]]==cc  for cc=1:length(pits)]) # pit in catchment of same color

    area, slen, dir, nout, nin, pits, c, bnds = WWF.waterflows(dem, drain_pits=true, bnd_as_pits=false);
    #plotarea_dem(xs, ys, dem, area, pits)
    @test length(pits) == 1
    @test maximum(slen)==195
    @test maximum(area[.!isnan.(area)])==9605
    @test length(unique(c))==2
    @test sum(c) == 9605
    @test sum(diff(c[:])) == 0
    @test sort(unique(c))[[1,end]] ==[0,1]
    @test bnds isa Array{Array{CartesianIndex{2},1},1}
    @test all([c[pits[cc]]==cc  for cc=1:length(pits)]) # pit in catchment of same color

    area, slen, dir, nout, nin, pits, c, bnds = WWF.waterflows(dem, drain_pits=true, bnd_as_pits=true);
    # plotarea_dem(xs, ys, dem, area, pits)
    @test length(pits) == 388
    @test maximum(slen)==81
    @test maximum(area[.!isnan.(area)])==4978
    @test length(unique(c))==389
    @test sum(c) == 1980928
    @test sum(diff(c[:])) == -1
    @test sort(unique(c))[[1,end]] ==[0,388]
    @test bnds isa Array{Array{CartesianIndex{2},1},1}
    @test all([c[pits[cc]]==cc  for cc=1:length(pits)]) # pit in catchment of same color
end

# bug with non-empty set
@testset "non_empty" begin
    (xs, ys), dem = dem_one_point()
    mask = .!isnan.(dem)
    #should not error
    area, slen, dir, nout, nin, pits, c, bnds = WWF.waterflows(dem, drain_pits=true, bnd_as_pits=false)
    @test mask[pits[1]]
    @test length(pits)==1
    @test all([c[pits[cc]]==cc  for cc=1:length(pits)]) # pit in catchment of same color
    #@test all(getindex.(Ref(mask), pits).==0) # tests that there are no interior pits left
    area, slen, dir, nout, nin, pits, c, bnds = WWF.waterflows(dem, drain_pits=true, bnd_as_pits=true)
    @test all(getindex.(Ref(mask), pits).==0)
    #area[isnan.(dem)] .= NaN; WWF.plotarea(xs, ys, area, pits)
    @test all([c[pits[cc]]==cc  for cc=1:length(pits)]) # pit in catchment of same color

    (xs, ys), dem = dem_two_points()
    mask = .!isnan.(dem)
    #should not error
    area, slen, dir, nout, nin, pits, c, bnds = WWF.waterflows(dem, drain_pits=true, bnd_as_pits=false)
    @test mask[pits[1]]
    @test mask[pits[2]]
    @test length(pits)==2
    @test all([c[pits[cc]]==cc  for cc=1:length(pits)]) # pit in catchment of same color
    area, slen, dir, nout, nin, pits, c, bnds = WWF.waterflows(dem, drain_pits=true, bnd_as_pits=true)
    @test all(getindex.(Ref(mask), pits).==0)
    #area[isnan.(dem)] .= NaN; WWF.plotarea(xs, ys, area, pits)
    @test all([c[pits[cc]]==cc  for cc=1:length(pits)]) # pit in catchment of same color

    (xs, ys), dem = dem_patho1()
    mask = .!isnan.(dem)
    #should not error
    area, slen, dir, nout, nin, pits, c, bnds = WWF.waterflows(dem, drain_pits=true, bnd_as_pits=false)
    @test length(pits)==2
    @test all([c[pits[cc]]==cc  for cc=1:length(pits)]) # pit in catchment of same color
    area, slen, dir, nout, nin, pits, c, bnds = WWF.waterflows(dem, drain_pits=true, bnd_as_pits=true)
    @test all(getindex.(Ref(mask), pits).==0)
    #area[isnan.(dem)] .= NaN; WWF.plotarea(xs, ys, area, pits)
    @test all([c[pits[cc]]==cc  for cc=1:length(pits)]) # pit in catchment of same color

    # should not have pit in the interior
    (xs, ys), dem = dem_patho2()
    mask = .!isnan.(dem)
    # also mask border points
    mask[1,:] .= false; mask[end,:] .= false; mask[:,1] .= false; mask[:,end] .= false
    # should not have a pit in the interior
    area, slen, dir, nout, nin, pits, c, bnds = WWF.waterflows(dem, drain_pits=true, bnd_as_pits=false)
    @test all(getindex.(Ref(mask), pits).==0)
    @test all([c[pits[cc]]==cc  for cc=1:length(pits)]) # pit in catchment of same color
    area, slen, dir, nout, nin, pits, c, bnds = WWF.waterflows(dem, drain_pits=true, bnd_as_pits=true)
    @test all(getindex.(Ref(mask), pits).==0)
    #area[isnan.(dem)] .= NaN; WWF.plotarea(xs, ys, area, pits)
    @test all([c[pits[cc]]==cc  for cc=1:length(pits)]) # pit in catchment of same color
end

@testset "NOFLOWer" begin
    xs, dem = peaks2_nan_edge()
    ys = xs

    area, slen, dir, nout, nin, pits, c, bnds = WWF.waterflows(dem, drain_pits=false, bnd_as_pits=true);
    @test sum(c.==0) == 8
    # add a NOFLOWer cell
    # (Note: this will not be consistent with the routing
    # (i.e. flow is still into NOFLOWer, i.e. NOFLOWer should be a pit) and thus
    # leads to bongous results.  Well, let's just test the bongous results anyway,
    # to make sure they stay the same bongous.)
    ind = CartesianIndex(43,30)
    dir[ind] = WWF.NOFLOWer
    # and drain the pits
    WWF.drainpits!(dir, nin, nout, pits, c, bnds, dem)
    area, slen, c = WWF.flowrouting_catchments(dir, pits, ones(size(dem)))
    @test sum(c.==0) == 42 # these are bongous...

    area_flow, dir_flow = WWF.waterflows(dem, drain_pits=true, bnd_as_pits=true)[[1,3]];
    @test sum(area.!=area_flow) == 79
    @test sum(dir.!=dir_flow) == 1
end

#################################
include("postproc.jl")

ext: jl


text:
using Pkg
Pkg.activate(normpath(joinpath(@__DIR__, ".."))) # Activate package from repo if it is not installed
using QuantExQASM

X = QuantExQASM.NCU.GateOps.default_gates["X"]

ctrl = collect(0:3)
aux = Array{Int64,1}()
tgt = length(ctrl)

cct = QuantExQASM.NCU.apply_ncu("q", ctrl, aux, tgt, X, 0)
println(QuantExQASM.Utils.format_string_nl(cct))


text:
include("measurement.jl")
include("comparison.jl")
include("projection.jl")
include("intersection.jl")
include("transformation.jl")


text:
#######################
## Timedata getindex ##
#######################

## Based on following DataFrame getindex behaviour:
## isa(df[1:2], DataFrame)
## isa(df[1:2, 3:4], DataFrame)
## isa(df[1, 3:4], DataFrame)
## isa(df[1], Array{Float64,1})        # loses column name
## isa(df[2:4, 1], DataArray{Float64,1})   # loses column name
## isa(df[2, 2], Float64)                  # loses column name

typealias ColumnIndex Union(Real, Symbol)
import Base.getindex

for t = (:Timedata, :Timenum, :Timematr)
    @eval begin

        ############################################
        ## already type preserving for DataFrames ##
        ############################################

        ## select multiple columns
        function getindex{T <: ColumnIndex}(td::$(t), col_inds::AbstractVector{T})
            valsDf = getindex(td.vals, col_inds) # dataframe
            return $(t)(valsDf, idx(td))
        end

        ## multiple rows, multiple columns
        function getindex{R <: Real, T <: ColumnIndex}(td::$(t),
                                                       row_inds::AbstractVector{R},
                                                       col_inds::AbstractVector{T})
            valsDf = getindex(td.vals, row_inds, col_inds) # dataframe
            dats = idx(td)[row_inds]
            return $(t)(valsDf, dats)
        end

        ## single row, multiple columns
        function getindex{T <: ColumnIndex}(td::$(t), row_ind::Real,
                                            col_inds::AbstractVector{T})
            valsDf = getindex(td.vals, row_ind, col_inds) # dataframe

            dat = [idx(td)[row_ind]]
            return $(t)(valsDf, dat)
        end

        #####################################################
        ## only type preserving for AbstractTimedata types ##
        #####################################################

        ## select single column
        function getindex(td::$(t), col_ind::ColumnIndex)
            vals = getindex(td.vals, col_ind) # dataarray or array

            ## manually get column name
            selected_column = td.vals.colindex[col_ind]
            name = names(td)[selected_column] # ASCIIString

            ## create respective dataframe
            if isa(vals, DataArray)
                valsDf = DataFrame(a = vals)
            else
                valsDf = convert(DataFrame, vals)
            end
            names!(valsDf, [name])           # names must be given as
            # array

            return $(t)(valsDf, idx(td))
        end

        ## multiple rows, single column
        function getindex{T <: Real}(td::$(t), row_inds::AbstractVector{T},
                                     col_ind::ColumnIndex)
            vals = getindex(td.vals, row_inds, col_ind) # dataarray or array

            ## manually get column name
            selected_column = td.vals.colindex[col_ind]
            name = names(td)[selected_column] # ASCIIString

            ## create respective dataframe
            if isa(vals, DataArray)
                valsDf = DataFrame(a = vals)
            else
                valsDf = convert(DataFrame, vals)
            end

            names!(valsDf, [name])           # names must be given as
            # array

            ## get idx
            dats = [idx(td)[row_inds]]

            return $(t)(valsDf, dats)
        end

        ## select single row and single column index
        function getindex(td::$(t), row_ind::Real, col_ind::ColumnIndex)
            val = getindex(td.vals, row_ind, col_ind) # single value / NA

            ## manually get column name
            selected_column = td.vals.colindex[col_ind]
            name = names(td)[selected_column] # ASCIIString

            ## create respective dataframe
            valDf = DataFrame(a = val)
            names!(valDf, [name])           # names must be given as
            # array

            ## single date needs to be transformed to DataArray
            dat = [idx(td)[row_ind]]

            return $(t)(valDf, dat)
        end

        #############################################################
        ## assumed as type preserving without further tests so far ##
        #############################################################

        ## getindex(td::$(t), ex::Expr) = getindex(td, with(td.vals, ex))
        ## getindex(td::$(t), ex::Expr, c::ColumnIndex) =
        ##     getindex(td, with(td.vals, ex), c)

        ## ## typealias ColumnIndex Union(Real, Symbol)
        ## getindex{T <: ColumnIndex}(td::$(t), ex::Expr, c::AbstractVector{T}) =
        ##     getindex(td, with(td.vals, ex), c)

        ## getindex(td::$(t), c::Real, ex::Expr) =
        ##     getindex(td, c, with(td.vals, ex))
        ## getindex{T <: Real}(td::$(t), c::AbstractVector{T}, ex::Expr) =
        ##     getindex(td, c, with(td.vals, ex))
        ## getindex(td::$(t), ex1::Expr, ex2::Expr) =
        ##     getindex(td, with(td.vals, ex1), with(td.vals, ex2))

        #########################
        ## indexing with dates ##
        #########################

        function getindex(td::$(t), date::Date)
            row_ind = idx(td) .== date
            return td[row_ind, :]
        end

        function getindex(td::$(t), date::Date, x::Any)
            row_ind = idx(td) .== date
            return getindex(td, row_ind, x)
        end

        function getindex(td::$(t), idxs::Array{Date,1})
            row_inds = findin(TimeData.idx(td), idxs)
            return td[row_inds, :]
        end

        function getindex(td::$(t), idxs::Array{Date,1}, x::Any)
            row_inds = findin(TimeData.idx(td), idxs)
            return td[row_inds, x]
        end

        function getindex(td::$(t), idx::StepRange{Date,Dates.Day})
            return getindex(td, [idx])
        end

        function getindex(td::$(t), idx::StepRange{Date,Dates.Day},
                          x::Any)
            return getindex(td, [idx], x)
        end

    end
end


text:
#=

  Lectures problem in Julia ConstraintSolver.js

  Biggs: Discrete Mathematics (2nd ed), page 187.
  """
  Suppose we wish to schedule six one-hour lectures, v1, v2, v3, v4, v5, v6.
  Among the the potential audience there are people who wish to hear both

   - v1 and v2
   - v1 and v4
   - v3 and v5
   - v2 and v6
   - v4 and v5
   - v5 and v6
   - v1 and v6

  How many hours are necessary in order that the lectures can be given
  without clashes?
  """

  Model created by Hakan Kjellerstrand, hakank@gmail.com
  See also my Julia page: http://www.hakank.org/julia/

=#

using ConstraintSolver, JuMP
using Cbc, GLPK, Ipopt
const CS = ConstraintSolver
include("constraints_utils.jl")

function lectures(problem,print_solutions=true,all_solutions=true)

    cbc_optimizer = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0)
    glpk_optimizer = optimizer_with_attributes(GLPK.Optimizer)
    ipopt_optimizer = optimizer_with_attributes(Ipopt.Optimizer)

    model = Model(optimizer_with_attributes(CS.Optimizer,   # "all_solutions"=> all_solutions,
                                                            "all_optimal_solutions"=>all_solutions,
                                                            "logging"=>[],

                                                            "traverse_strategy"=>:BFS,
                                                            # "traverse_strategy"=>:DFS,
                                                            # "traverse_strategy"=>:DBFS,

                                                            # "branch_split"=>:Smallest,
                                                            # "branch_split"=>:Biggest,
                                                            "branch_split"=>:InHalf,

                                                            # https://wikunia.github.io/ConstraintSolver.jl/stable/options/#branch_strategy-(:Auto)
                                                            "branch_strategy" => :IMPS, # default
                                                            # "branch_strategy" => :ABS, # Activity Based Search
                                                            # "activity.decay" => 0.999, # default 0.999
                                                            # "activity.max_probes" => 10, # default, 10
                                                            # "activity.max_confidence_deviation" => 20, # default 20

                                                            # "simplify"=>false,
                                                            # "simplify"=>true, # default

                                                            "time_limit"=>6,

                                                            # "backtrack" => false, # default true
                                                            # "backtrack_sorting" => false, # default true

                                                            # "lp_optimizer" => cbc_optimizer,
                                                            # "lp_optimizer" => glpk_optimizer,
                                                            # "lp_optimizer" => ipopt_optimizer,
                                        ))
    n, _ = size(problem)
    @variable(model, 1 <= x[1:n] <= n, Int)
    @variable(model, 1 <= max_t <= n, Int)

    for (l1,l2) in eachrow(problem)
        @constraint(model, x[l1] != x[l2])
    end

    # symmetry breaking
    for i in 1:n
        @constraint(model, x[i] <= i)
    end
    @constraint(model, x .<= max_t)
    @objective(model,Min,max_t)

    # Solve the problem
    optimize!(model)

    status = JuMP.termination_status(model)
    # println("status:$status")
    if status == MOI.OPTIMAL
        num_sols = MOI.get(model, MOI.ResultCount())
        println("num_sols:$num_sols\n")
        if print_solutions
            for sol in 1:num_sols
                println("solution #$sol")
                x_val = convert.(Integer,JuMP.value.(x; result=sol))
                max_t_val = convert.(Integer,JuMP.value.(max_t; result=sol))
                println("x:$x_val")
                println("max_t:$max_t_val")
                println()

            end
        end
    else
        println("status:$status")
    end

    return status
end

# The schedule requirements:
#   lecture a cannot be held at the same time as b
 lectures_problem = resize_matrix([[1, 2],
                                    [1, 4],
                                    [3, 5],
                                    [2, 6],
                                    [4, 5],
                                    [5, 6],
                                    [1, 6]])

@time lectures(lectures_problem)

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
"""
```
l = segment_image(Algorithm::GrowCut, feature_vector::AbstractArray, seeds::AbstractArray{Int, 3}, t1::Int=27, t2::Int=27; max_iter::Int=1000, converge_at::Int=1, rule::Tuple{Function, Function}=(colour_diff,find_maxC))
```
Same as [`segment_image`](@ref segment_image(::GrowCut, ::AbstractArray, ::AbstractArray{Int, 2}, ::Int, ::Int; ::Int, ::Int, ::Tuple{Function, Function}))
except runs on 3d images.

# Example

Create synthetic 3d image and segment:

```julia
using Images, PictureSegmentation

# Create image of cone with multiple colours
img=Array{RGB{Float64},3}(undef,100,100,100)
for i in CartesianIndices(img)
    if (i[1]-50)^2 + (i[2]-50)^2 < (i[3]/4)^2
        img[i]=RGB(i[3]/100,0,1-(i[3]/100))
    else
        img[i]=RGB(0,0,0)
    end
end

# Create seed array and select seed pixels
seeds=zeros(Int,axes(img))
set_seed_pixels(seeds, img)

# Segment image
l = segment_image(GrowCut(), img, seeds)

# Create new image from labels and display
img2 = convert_labels(l)
display_3d(img2)
```

See also: [`display_3d`](@ref display_3d(::AbstractArray{RGB{Float64}, 3}, ::Int, ::Bool))
and [`set_seed_pixels`](@ref set_seed_pixels(::Array{Int,3}, ::AbstractArray, ::Int; ::Bool).
"""
function segment_image(Algorithm::GrowCut, feature_vector::AbstractArray, seeds::AbstractArray{Int, 3}, t1::Int=27, t2::Int=27; max_iter::Int=1000, converge_at::Int=1, rule::Tuple{Function, Function}=(colour_differance,find_maximum))
    # initialize values
    calculate_adjusted_strength = rule[1]
    extra_function = rule[2]
    extras = extra_function(feature_vector)
    labels = copy(seeds)
    E = zeros(Int, axes(labels))
    θ = set_strength(labels)
    labelsₜ₊₁ = copy(seeds)
    Eₜ₊₁ = zeros(Int, axes(labels))
    θₜ₊₁ = set_strength(labels)
    regions = maximum(seeds)
    iter = 0
    count = length(labels)
    changed = Array{Tuple{UnitRange{Int},UnitRange{Int},UnitRange{Int}},1}(undef,1)
    changedₜ₊₁ = Array{Tuple{UnitRange{Int},UnitRange{Int},UnitRange{Int}},1}(undef,0)
    rₘ, cₘ, wₘ = size(labels)
    changed[1] = 1:rₘ, 1:cₘ, 1:wₘ

    # overflow warning
    if typeof(feature_vector) == Array{RGB, 3} && typeof(feature_vector) != Array{RGB{Float64}, 3}
        @warn "feature_vector is RGB and not of type RGB{Float64}. This may cause overflow errors."
    end

    # iterate till convergence or maximum iterations reached
    while count > converge_at && iter < max_iter
        iter += 1
        count = 0
        copyto!(θₜ₊₁, θ)
        copyto!(labelsₜ₊₁, labels)
        copyto!(Eₜ₊₁, E)

        # loop through active pixel neighbourhoods (all pixels on first iteration)
        for j in CartesianIndices(changed)
            @inbounds for i in CartesianIndices(changed[j])
                Nbound = max(i[1]-1,1):min(i[1]+1, rₘ), max(i[2]-1, 1):min(i[2]+1, cₘ), max(i[3]-1, 1):min(i[3]+1, wₘ)
                Eₜ₊₁[i] = count_enemy(labels, i, Nbound)
                if E[i] < t2
                    θₜ₊₁, labelsₜ₊₁, count, changedₜ₊₁ = update_pixel(θ, θₜ₊₁, labels, labelsₜ₊₁, feature_vector, E, i, extras, t1, count, Nbound, changedₜ₊₁, calculate_adjusted_strength)
                else
                    θₜ₊₁, labelsₜ₊₁, count, changedₜ₊₁ = occupy_pixel(θ, θₜ₊₁, labels, labelsₜ₊₁, feature_vector, i, extras, count, Nbound, changedₜ₊₁, calculate_adjusted_strength)
                end
            end
        end

        # update values
        # changed=changedₜ₊₁
        changedₜ₊₁ = similar(changedₜ₊₁, 0)
        E = Eₜ₊₁
        θ = θₜ₊₁
        labels = labelsₜ₊₁
    end
    return labels
end


text:
# Julia 0.6
function singleNumber(vector)
	reduce(xor, vector)
end

a = [1,1,2,3,3]
println(singleNumber(a))

max_stars_repo_name: cdluminate/MyNotes

max_issues_repo_name: cdluminate/MyNotes

max_forks_repo_name: cdluminate/MyNotes


text:
"""
    OutwardsPopulationDispersal <: SetNeighborhoodRule

    OutwardsPopulationDispersal(; kw...)
    OutwardsPopulationDispersal{R}(; kw...)
    OutwardsPopulationDispersal{R,W}(; kw...)

Implements deterministic dispersal from the current cell to populations in neighboring
cells.

This will make sense ecologically where cell populations are large,
otherwise a randomised kernel may be more suitable.

The result should be identical to those obtained substituting `OutwardsDispersal` for
[`InwardsDispersal`](@ref) but may be more efficient when a small number of cells are
occupied. Conversely, it will become less efficient when a large proportion of the grid
is occupied.

# Keywords

- `neighborhood`: Any DynamicGrids.jl `Neighborhood`, or an
    already constructed [`DispersalKernel`](@ref). Using this keyword means `radius` is
    ignored, and for a `DispersalKernel`, all other keywords are ignored.
- `neighborhood`: `Neighborhood` object specifying the range from the origin of the
    discretised dispersal kernal. Defaults to `Window(radius)`.
- `formulation`: kernel formulation object holding the exact form of the kernal.
    Default [`ExponentialKernel`](@ref).
- `cellsize`: the cell size of the discretised kernal (i.e. simulation grid size).
    Default is 1.0.
- `distancemethod`: [`DistanceMethod`](@ref) object for calculating distance between cells.
    The default is [`CentroidToCentroid`](@ref).

Pass grid name `Symbol`s to `R` and `W` type parameters to use specific grids.
"""
struct OutwardsDispersal{R,W,N<:AbstractKernelNeighborhood} <: SetNeighborhoodRule{R,W}
    neighborhood::N
end
function OutwardsDispersal{R,W}(; kw...) where {R,W}
    OutwardsDispersal{R,W}(DispersalKernel(; kw...))
end

@inline function applyrule!(data, rule::OutwardsDispersal{R,W}, N, I) where {R,W}
    N == zero(N) && return nothing
    sum = zero(N)
    for (offset, k) in zip(offsets(rule), kernel(rule))
        @inbounds propagules = N * k
        @inbounds add!(data[W], propagules, I .+ offset...)
        sum += propagules
    end
    @inbounds sub!(data[W], sum, I...)
    return nothing
end


text:
# AUTO GENERATED FILE - DO NOT EDIT

export ''_leafletcircle

"""
    ''_leafletcircle(;kwargs...)

A LeafletCircle component.

Keyword arguments:
- `id` (String; optional)
- `center` (required): . center has the following type: lists containing elements 'lng', 'lat'.
Those elements have the following types:
  - `lng` (Real; optional)
  - `lat` (Real; optional)
- `loading_state` (optional): . loading_state has the following type: lists containing elements 'is_loading', 'prop_name', 'component_name'.
Those elements have the following types:
  - `is_loading` (Bool; optional): Determines if the component is loading or not
  - `prop_name` (String; optional): Holds which property is loading
  - `component_name` (String; optional): Holds the name of the component that is loading
- `pathOptions` (optional)
- `radius` (Real; required)
"""
function ''_leafletcircle(; kwargs...)
        available_props = Symbol[:id, :center, :loading_state, :pathOptions, :radius]
        wild_props = Symbol[]
        return Component("''_leafletcircle", "LeafletCircle", "feffery_leaflet_components", available_props, wild_props; kwargs...)
end

max_stars_repo_name: CNFeffery/feffery-leaflet-components

max_issues_repo_name: CNFeffery/feffery-leaflet-components

max_forks_repo_name: CNFeffery/feffery-leaflet-components


text:
struct Vars
    t::Num
    θ::Vector{Num}
    ϕ::Vector{Num}
end

Vars(N) = Vars((@variables t, θ[1:N](t), ϕ[1:N](t))...)

struct RHS
    θ::Vector{Num}
    ϕ::Vector{Num}
end

Base.:+(x::RHS, y::RHS) = RHS(x.θ .+ y.θ, x.ϕ .+ y.ϕ)

function effectiveField(field::Array, vars::Vars, α)
    # field is Nx3 dimensional, given in cartesian coordinates
    N = length(vars.ϕ)
    ϕres = Vector{Num}(undef, N)
    θres = Vector{Num}(undef, N)
    for i in 1:N
        hx, hy, hz = field[i,:]
        θ = vars.θ[i]
        ϕ = vars.ϕ[i]
        ϕres[i] = hz + cos(ϕ) * (-(hx * cot(θ)) + hy * α * csc(θ)) - hy * cot(θ) * sin(ϕ) - hx * α * csc(θ) * sin(ϕ)
        θres[i] = (hy + hx * α * cos(θ)) * cos(ϕ) - hz * α * sin(θ) + (-hx + hy * α * cos(θ)) * sin(ϕ)
    end
    RHS(θres, ϕres)
end

function spinTorque(torque::Array, vars::Vars)
    # takes as input an array
    N = length(vars.ϕ)
    ϕres = Vector{Num}(undef, N)
    θres = Vector{Num}(undef, N)
    for i in 1:N
        ϕ = vars.ϕ[i]
        θ = vars.θ[i]
        ϕres[i] = csc(θ) * (torque[i,1] * sin(ϕ) - torque[i,2] * cos(ϕ))
        θres[i] = torque[i,3] * sin(θ) - cos(θ) * (torque[i,1] * cos(ϕ) + torque[i,2] * sin(ϕ))
    end
    return RHS(θres, ϕres)
end

function spinTorque(torquefunc::Function, vars::Vars)
    # takes as input an function with as parameter Particle and should return a 3D array
    mapping = mappingAtoms(unit_cells)
    N = length(mapping)
    res = zeros(N, 3)
    for (i, p) in mapping
        res[i, :] = torquefunc(p)
    end
    spintorque(res, vars)
end

function exchangeCoupling(coupling::Matrix, vars::Vars)
    # coupling is a NxN array which describes which sites couple to which
    # returns two arrays with θ and ϕ right hand sides of the EOM
    ϕres = Vector{Num}(undef, size(coupling)[1])
    θres = Vector{Num}(undef, size(coupling)[1])
    for i in 1:size(coupling)[1]
        tempθ = 0
        tempϕ = 0
        ϕi = vars.ϕ[i]
        θi = vars.θ[i]
        for j in 1:size(coupling)[2]
            ϕj = vars.ϕ[j]
            θj = vars.θ[j]
            tempθ += -coupling[i,j] * sin(θj) * sin(θi - θj)
            tempϕ += coupling[i,j] * (
                cos(θj) - cos(ϕi - ϕj) * cot(θi) * sin(θj)
            )
        end
        θres[i] = tempθ
        ϕres[i] = tempϕ
    end
    RHS(θres, ϕres)
end

function DMICoupling(coupling::Array, DMI::Array, vars::Vars)
    # coupling is a NxN array which describes which sites couple to which
    # DMI is a 3D vector that describes the DMI field
    # returns two arrays with θ and ϕ right hand sides of the EOM
    ϕres = Vector{Num}(undef, size(coupling)[1])
    θres = Vector{Num}(undef, size(coupling)[1])
    Dx, Dy, Dz = DMI
    for i in 1:size(coupling)[1]
        tempθ = 0
        tempϕ = 0
        ϕi = vars.ϕ[i]
        θi = vars.θ[i]
        for j in 1:size(coupling)[2]
            ϕj = vars.ϕ[j]
            θj = vars.θ[j]

            tempθ += coupling[i,j] * (
                Dz * cos(ϕi - ϕj) * sin(θj) - cos(θj) * (Dx * cos(ϕi) + Dy * sin(ϕi))
            )
            tempϕ += coupling[i,j] * (
                -(cos(θj) * cot(θi) * (Dy * cos(ϕi) - Dx * sin(ϕi))) +
                -  sin(θj) * (-(cos(ϕj) * (Dy + Dz * cot(θi) * sin(ϕi))) + (Dx + Dz * cos(ϕi) * cot(θi)) * sin(ϕj))
                )
        end
        θres[i] = tempθ
        ϕres[i] = tempϕ
    end
    RHS(θres, ϕres)
end

function setupEOM(vars::Vars, rhs::RHS)
    D = Differential(vars.t)
    [D.(vars.θ) .~ rhs.θ; D.(vars.ϕ) .~ rhs.ϕ]
end

function initBoltzmann(vars::Vars, λ)
    # boltzmann is proportional to θ (for now)
    # the angle ϕ is chosen randomly between 0 and 2pi
    N = length(vars.θ)
    dist = Exponential(λ / 2)
    init_θ = rand(dist, N)
    init_ϕ = rand(Float64, N) * 2π
    return [[vars.θ[i] => init_θ[i] for i in 1:N]; [vars.ϕ[i] => init_ϕ[i] for i in 1:N]]
end

max_stars_repo_name: pgunnink/LLGTools

max_issues_repo_name: pgunnink/LLGTools

max_forks_repo_name: pgunnink/LLGTools


text:
# Julia wrapper for header: /usr/local/include/libavutil/dict.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

export
    av_dict_get,
    av_dict_count,
    av_dict_set,
    av_dict_set_int,
    av_dict_parse_string,
    av_dict_copy,
    av_dict_free,
    av_dict_get_string

function av_dict_get(m,key,prev,flags::Integer)
    ccall((:av_dict_get,libavutil),Ptr{AVDictionaryEntry},(Ptr{AVDictionary},Cstring,Ptr{AVDictionaryEntry},Cint),m,key,prev,flags)
end

function av_dict_count(m)
    ccall((:av_dict_count,libavutil),Cint,(Ptr{AVDictionary},),m)
end

function av_dict_set(pm,key,value,flags::Integer)
    ccall((:av_dict_set,libavutil),Cint,(Ptr{Ptr{AVDictionary}},Cstring,Cstring,Cint),pm,key,value,flags)
end

function av_dict_set_int(pm,key,value::Int64,flags::Integer)
    ccall((:av_dict_set_int,libavutil),Cint,(Ptr{Ptr{AVDictionary}},Cstring,Int64,Cint),pm,key,value,flags)
end

function av_dict_parse_string(pm,str,key_val_sep,pairs_sep,flags::Integer)
    ccall((:av_dict_parse_string,libavutil),Cint,(Ptr{Ptr{AVDictionary}},Cstring,Cstring,Cstring,Cint),pm,str,key_val_sep,pairs_sep,flags)
end

function av_dict_copy(dst,src,flags::Integer)
    ccall((:av_dict_copy,libavutil),Cint,(Ptr{Ptr{AVDictionary}},Ptr{AVDictionary},Cint),dst,src,flags)
end

function av_dict_free(m)
    ccall((:av_dict_free,libavutil),Void,(Ptr{Ptr{AVDictionary}},),m)
end

function av_dict_get_string(m,buffer,key_val_sep::UInt8,pairs_sep::UInt8)
    ccall((:av_dict_get_string,libavutil),Cint,(Ptr{AVDictionary},Ptr{Cstring},UInt8,UInt8),m,buffer,key_val_sep,pairs_sep)
end

max_stars_repo_name: asafasban/AAVideoIO

max_issues_repo_name: asafasban/AAVideoIO

max_forks_repo_name: asafasban/AAVideoIO


text:
module DataAndPlotting

using HDF5
using Plots

#Wave parameters
const ki = 3
const kj = 3

#grid size
const ni = 3000
const nj = 2000

#i ranges from 1 to ni
#i - 1 ranges from 0 to ni - 1
#x ranges from 0 to 1
function index2coord(i, ni)
    @assert ni > 1
    @assert 1 <= i <= ni
    (i - 1) / (ni - 1)
end

#set up standing wave
export initialize
function initialize()
    arr = zeros(Float64, ni, nj)
    for j in 1:nj, i in 1:ni
        x = index2coord(i, ni)
        y = index2coord(j, nj)
        arr[i,j] = sin(2π*ki*x) * sin(2π*kj*y)
    end
    arr

end

#write data to file
export output
function output(arr::Array{Float64,2}, filename::String)
    h5write(filename, "data", arr)
end

#Read data from file
export input
function input(filename::String)
    arr = h5read(filename, "data") ## (1:ni, 1:nj)
end

#Create a plot
export makeplot
function makeplot(arr::Array{Float64,2}, filename::String)
    heatmap(arr, clim=(-1.0,+1.0), color=:viridis)
    savefig(filename)
end

end # module

max_stars_repo_name: brodriguezmilla/DataAndPlotting

max_issues_repo_name: brodriguezmilla/DataAndPlotting

max_forks_repo_name: brodriguezmilla/DataAndPlotting


text:
using Findpeaks
using Test

gaussian(x,μ,σ) = 1/sqrt(2*π)/σ*exp(-((x-μ)^2)/2/σ^2)
x = collect(range(1,stop=1000,length=1000))
μ = rand(x)
data = gaussian.(x,μ,rand(x)/3)

p = findpeaks(data)
@test length(p) == 1
@test p[1] == μ
p = findpeaks(data,x)
@test length(p) == 1
@test p[1] == μ


text:
module DecentralizedInternet

end # module


text:
using StaticStrings
using Test
using Random: shuffle, shuffle!

@testset "StaticStrings.jl" begin
    include("static.jl")

    include("padded.jl")
end


text:
if basename(pwd()) == "aoc"
    cd("2016/18")
end

function countsafe(str::AbstractString, n::Integer)
    row = [c == '^' for c in str]
    next = falses(size(row))
    total = 0
    for _ in 1:n
        total += count(!, row)
        next[1] = row[2]
        next[2:end-1] = row[1:end-2] .⊻ row[3:end]
        next[end] = row[end-1]
        row, next = next, row
    end
    total
end

@assert countsafe(".^^.^.^^^^", 10) == 38
@time countsafe(readline("input.txt"), 40)
@time countsafe(readline("input.txt"), 400000)

max_stars_repo_name: stellartux/advent-of-code

max_issues_repo_name: stellartux/advent-of-code

max_forks_repo_name: stellartux/advent-of-code


text:
module RadeonProRender

using GeometryBasics
using Colors
using CEnum
using RadeonProRender_jll
using Printf

assetpath(paths...) = normpath(joinpath(@__DIR__, "..", "assets", paths...))

include("LibRPR.jl")

using .RPR

include("highlevel-api.jl")
include("materials.jl")

export set!
export setradiantpower!
export setintensityscale!
export setportal!
export setalbedo!
export setturbidity!
export setscale!
export clear!
export transform!
export layeredshader
export TileIterator
export set_standard_tonemapping!
export setsubdivisions!
export setdisplacementscale!
export lookat!

export Tahoe, Northstar, Hybrid, HybridPro

end


text:
"""
    range(from::Date, rdate::RDate; inc_from=true, cal_mgr=nothing)
    range(from::Date, to::Date, rdate::RDate; inc_from=true, inc_to=true, cal_mgr=nothing)

The range provides a mechanism for iterating over a range of dates given a period. This can
provide a mechanism for getting an infinite range (from a given date) or appropriately
clipped.

```julia-repl
julia> collect(Iterators.take(range(Date(2017,1,25), rd"1d"), 3))
3-element Vector{Date}:
 2017-01-25
 2017-01-26
 2017-01-27
julia> collect(Iterators.take(range(Date(2017,1,25), rd"1d"; inc_from=false), 3))
3-element Vector{Date}:
 2017-01-26
 2017-01-27
 2017-01-28
julia> collect(range(Date(2019,4,17), Date(2019,4,22), rd"2d"))
3-element Vector{Date}:
 2019-04-17
 2019-04-19
 2019-04-21
```

Under the hoods, the range will `multiply` the period. Since non-periodic RDates will
always give back self when you multiply it allows us to set a reference point.
```julia-repl
julia> rd"1JAN2001+3m+3rd WED"
1JAN2001+3m[LDOM;PDOM]+3rd WED
julia> 3*rd"1JAN2001+3m+3rd WED"
1JAN2001+9m[LDOM;PDOM]+3rd WED
```

This provides the basic building blocks to come up with more complex functionality. For
example to get the next four [IMM dates](https://en.wikipedia.org/wiki/IMM_dates)

```julia-repl
julia> d = Date(2017,1,1)
julia> collect(Iterators.take(range(d, rd"1MAR+3m+3rd WED"), 4))
4-element Vector{Date}:
 2017-03-15
 2017-06-21
 2017-09-20
 2017-12-20
```
"""
struct RDateRange
    from::Dates.Date
    to::Union{Dates.Date, Nothing}
    period::RDate
    inc_from::Bool
    inc_to::Bool
    calendar_mgr::CalendarManager
end

function Base.iterate(iter::RDateRange, state=nothing)
    if state === nothing
        count = 0
        elem = apply(multiply(iter.period, count), iter.from, iter.calendar_mgr)
        while elem > iter.from
            count -= 1
            elem = apply(multiply(iter.period, count), iter.from, iter.calendar_mgr)
        end

        from_op = iter.inc_from ? Base.:< : Base.:<=
        while from_op(elem, iter.from)
            count += 1
            elem = apply(multiply(iter.period, count), iter.from, iter.calendar_mgr)
        end

        state = (elem, count)
    end
    elem, count = state
    op = iter.inc_to ? Base.:> : Base.:>=
    if iter.to !== nothing && op(elem,iter.to)
        return nothing
    end

    return (elem, (apply(multiply(iter.period, count+1), iter.from, iter.calendar_mgr), count+1))
end

Base.IteratorSize(::Type{RDateRange}) = Base.SizeUnknown()
Base.eltype(::Type{RDateRange}) = Dates.Date
function Base.range(from::Dates.Date, period::RDate; inc_from::Bool=true, cal_mgr::Union{CalendarManager,Nothing}=nothing)
    return RDateRange(from, nothing, period, inc_from, false, cal_mgr !== nothing ? cal_mgr : NullCalendarManager())
end

function Base.range(from::Dates.Date, to::Dates.Date, period::RDate; inc_from::Bool=true, inc_to::Bool=true, cal_mgr::Union{CalendarManager,Nothing}=nothing)
    return RDateRange(from, to, period, inc_from, inc_to, cal_mgr !== nothing ? cal_mgr : NullCalendarManager())
end


text:
using Stemmers

algs = stemmer_types()
@assert !isempty(algs)

for alg in algs
    stmr = Stemmer(alg)
    Stemmers.release(stmr)
end

test_cases = {
    "english" => {
        "working" => "work",
        "worker" => "worker"
    }
}

for (alg, test_words) in test_cases
    stmr = Stemmer(alg)
    for (n,v) in test_words
        @assert v == stem(stmr, n)
    end
end


text:
import Pkg
Pkg.activate(joinpath(@__DIR__, ".."))
import Juno
using Revise
using Beamforming, Seis
using Plots

Beamforming.verbose(false)

function f(; kwargs...)
    t = Seis.sample_data(:array)
    Beamforming.beamform(t, 1110, 1140, 6, 0.25; kwargs...)
end

Juno.@profiler f(; method=:phaseweight, n=3)
@time fk = f()

heatmap(fk.sx, fk.sy, fk.power', aspect_ratio=:equal,
    xlim=extrema(fk.sx), ylim=extrema(fk.sy),
    c=Plots.cgrad([:white,:blue,:green,:yellow,:orange,:red]))

max_stars_repo_name: anowacki/Beamforming

max_issues_repo_name: anowacki/Beamforming

max_forks_repo_name: anowacki/Beamforming


text:
function fill_refs!(refs::AbstractArray, X::AbstractArray,
                    breaks::AbstractVector, extend::Bool, nullok::Bool)
    n = length(breaks)
    lower = first(breaks)
    upper = last(breaks)

    @inbounds for i in eachindex(X)
        x = X[i]

        if extend && x == upper
            refs[i] = n-1
        elseif !extend && !(lower <= x < upper)
            throw(ArgumentError("value $x (at index $i) does not fall inside the breaks: adapt them manually, or pass extend=true"))
        else
            refs[i] = searchsortedlast(breaks, x)
        end
    end
end

function fill_refs!{T<:Nullable}(refs::AbstractArray, X::AbstractArray{T},
                                  breaks::AbstractVector, extend::Bool, nullok::Bool)
    n = length(breaks)
    lower = first(breaks)
    upper = last(breaks)

    @inbounds for i in eachindex(X)
        isnull(X[i]) && continue

        x = unsafe_get(X[i])

        if extend && x == upper
            refs[i] = n-1
        elseif !extend && !(lower <= x < upper)
            nullok || throw(ArgumentError("value $x (at index $i) does not fall inside the breaks: adapt them manually, or pass extend=true or nullok=true"))
            refs[i] = 0
        else
            refs[i] = searchsortedlast(breaks, x)
        end
    end
end

function fill_refs!(refs::AbstractArray, X::NullableArray,
                    breaks::AbstractVector, extend::Bool, nullok::Bool)
    n = length(breaks)
    lower = first(breaks)
    upper = last(breaks)

    @inbounds for i in eachindex(X)
        X.isnull[i] && continue

        x = X.values[i]

        if extend && x == upper
            refs[i] = n-1
        elseif !extend && !(lower <= x < upper)
            nullok || throw(ArgumentError("value $x (at index $i) does not fall inside the breaks: adapt them manually, or pass extend=true or nullok=true"))
            refs[i] = 0
        else
            refs[i] = searchsortedlast(breaks, x)
        end
    end
end

_extrema(X::Any) = extrema(X)
# NullableArrays provide a more efficient version with higher precedence
_extrema{T<:Nullable}(X::AbstractArray{T}) = (minimum(X), maximum(X))

unwrap(x::Any) = x
unwrap(x::Nullable) = x.value

"""
    cut(x::AbstractArray, breaks::AbstractVector;
        extend::Bool=false, labels::AbstractVector=[])

Cut a numeric array into intervals and return an ordered `CategoricalArray` indicating
the interval into which each entry falls. Intervals are of the form `[lower, upper)`,
i.e. the lower bound is included and the upper bound is excluded.

# Arguments
* `extend::Bool=false`: when `false`, an error is raised if some values in `x` fall
  outside of the breaks; when `true`, breaks are automatically added to include all
  values in `x`, and the upper bound is included in the last interval.
* `labels::AbstractVector=[]`: a vector of strings giving the names to use for the
  intervals; if empty, default labels are used.

    cut(x::AbstractArray{<:Nullable}, breaks::AbstractVector;
        extend::Bool=false, labels::AbstractVector=[]), nullok::Bool=false)

For nullable arrays, return a `NullableCategoricalArray`. If `nullok=true`, values outside
of breaks result in null values.
"""
function cut{T, N, U<:AbstractString}(x::AbstractArray{T, N}, breaks::AbstractVector;
                                      extend::Bool=false, labels::AbstractVector{U}=String[],
                                      nullok::Bool=false)
    if !issorted(breaks)
        breaks = sort(breaks)
    end

    if extend
        min_x, max_x = _extrema(x)
        if !isnull(min_x) && breaks[1] > unwrap(min_x)
            unshift!(breaks, unwrap(min_x))
        end
        if !isnull(max_x) && breaks[end] < unwrap(max_x)
            push!(breaks, unwrap(max_x))
        end
    end

    refs = Array{DefaultRefType, N}(size(x))
    try
        fill_refs!(refs, x, breaks, extend, nullok)
    catch err
        # So that the error appears to come from cut() itself,
        # since it refers to its keyword arguments
        if isa(err, ArgumentError)
            throw(err)
        else
            rethrow(err)
        end
    end

    n = length(breaks)
    if isempty(labels)
        from = map(x -> sprint(showcompact, x), breaks[1:n-1])
        to = map(x -> sprint(showcompact, x), breaks[2:n])
        levs = Vector{String}(n-1)
        for i in 1:n-2
            levs[i] = string("[", from[i], ", ", to[i], ")")
        end
        if extend
            levs[end] = string("[", from[end], ", ", to[end], "]")
        else
            levs[end] = string("[", from[end], ", ", to[end], ")")
        end
    else
        length(labels) == n-1 || throw(ArgumentError("labels must be of length $(n-1), but got length $(length(labels))"))
        levs = copy(labels)
    end

    pool = CategoricalPool(levs, true)
    if T <: Nullable
        NullableCategoricalArray{String, N, DefaultRefType}(refs, pool)
    else
        CategoricalArray{String, N, DefaultRefType}(refs, pool)
    end
end

"""
    cut(x::AbstractArray, ngroups::Integer;
        labels::AbstractVector=String[])

Cut a numeric array into `ngroups` quantiles, determined using
[`quantile`](@ref).
"""
cut{U<:AbstractString}(x::AbstractArray, ngroups::Integer;
                       labels::AbstractVector{U}=String[]) =
    cut(x, quantile(x, (1:ngroups-1)/ngroups); extend=true, labels=labels)


text:
factor_residual_variance(dfm::DynamicFactorModel) = sum(dfm.factor_residuals.^2)/apply(*, size(dfm.x))  # see page 201 of Bai Ng 2002
#factor_residual_variance(dfm::DynamicFactorModel) = sum(mapslices(x->x'x/length(x), dfm.factor_residuals, 1))/size(dfm.x, 2)  # the same as above
# and var(dfm.factor_residuals) is approximately the same as well but it probably corrects for small sample bias

function criterion_cumulative_variance(pca_result, threshold=0.95)   # simply use factors until a certain threshold of variance is reached
    pca_result.cumulative_variance .< threshold  # take threshold% of variance
end

# PCp criteria as defined on page 201 of Bai and Ng 2002
function criterion_PCp1(dfm::DynamicFactorModel)
    dfm_unrestricted = DynamicFactorModel(dfm.y, dfm.w, dfm.x)
    N_plus_T_by_NT = apply(+, size(dfm.x))/apply(*, size(dfm.x))
    factor_residual_variance(dfm) + dfm.number_of_factors*factor_residual_variance(dfm_unrestricted)*(N_plus_T_by_NT)*log(N_plus_T_by_NT^-1)
end
function criterion_PCp2(dfm::DynamicFactorModel)
    dfm_unrestricted = DynamicFactorModel(dfm.y, dfm.w, dfm.x)
    N_plus_T_by_NT = apply(+, size(dfm.x))/apply(*, size(dfm.x))
    factor_residual_variance(dfm) + sum(dfm.number_of_factors)*factor_residual_variance(dfm_unrestricted)*(N_plus_T_by_NT)*log(minimum(size(dfm.x)))
end
function criterion_PCp3(dfm::DynamicFactorModel)
    dfm_unrestricted = DynamicFactorModel(dfm.y, dfm.w, dfm.x)
    N_plus_T_by_NT = apply(+, size(dfm.x))/apply(*, size(dfm.x))
    factor_residual_variance(dfm) + sum(dfm.number_of_factors)*factor_residual_variance(dfm_unrestricted)*log(minimum(size(dfm.x)))/minimum(size(dfm.x))
end

# ICp criteria as defined on page 201 of Bai and Ng 2002
function criterion_ICp1(dfm::DynamicFactorModel)
    N_plus_T_by_NT = apply(+, size(dfm.x))/apply(*, size(dfm.x))
    log(factor_residual_variance(dfm)) + sum(dfm.number_of_factors)*(N_plus_T_by_NT)*log(N_plus_T_by_NT^-1)
end
function criterion_ICp2(dfm::DynamicFactorModel)
    N_plus_T_by_NT = apply(+, size(dfm.x))/apply(*, size(dfm.x))
    log(factor_residual_variance(dfm)) + sum(dfm.number_of_factors)*(N_plus_T_by_NT)*log(minimum(size(dfm.x)))
end
function criterion_ICp3(dfm::DynamicFactorModel)
    N_plus_T_by_NT = apply(+, size(dfm.x))/apply(*, size(dfm.x))
    log(factor_residual_variance(dfm)) + sum(dfm.number_of_factors)*log(minimum(size(dfm.x)))/minimum(size(dfm.x))
end

#  default criteria which are not quite consistent at least not when the number of factors is estimated TODO: why exactly is that?
function criterion_BIC(dfm::DynamicFactorModel)
    T, N = size(dfm.x)
    factor_residual_variance(dfm) + dfm.number_of_factors * log(T)/T
end


text:
"""
    MasterCalendar(startdate::Date,enddate::Date)

Produces a mastercalendar (Dict{Int64,Date}) of weekdays where clinic is open between the start and enddate
# Examples
```julia-repl
julia> MasterCalendar(Date("2019-01-01"),Date("2019-02-01"))
Dict{Int64,Date} with 24 entries:
  18 => 2019-01-18
  24 => 2019-01-24
  ⋮  => ⋮
```
"""
function MasterCalendar(startdate::Date,enddate::Date)
    Dict((date-startdate).value+1 => date for date = filter(td -> dayofweek(td) <= 5, startdate:Day(1):enddate))
end
"""
    MasterCalendar(mastercalendar::Dict{Int64,Date},startdate,enddate)

Produces a mastercalendar ('Dict{Int64,Date}') of weekdays where clinic is open.

# Examples
```julia-repl
julia> mc = MasterCalendar(Date("2019-01-01"),Date("2019-02-01"))
Dict{Int64,Date} with 24 entries:
  18 => 2019-01-18
  24 => 2019-01-24
  ⋮  => ⋮
  julia> MasterCalendar(mc,Date("2019-01-01"),Date("2019-01-10"))
  Dict{Int64,Date} with 6 entries:
    7 => 2019-01-07
    9 => 2019-01-09
    4 => 2019-01-04
    2 => 2019-01-02
    3 => 2019-01-03
    8 => 2019-01-08
```
"""
function MasterCalendar(mastercalendar::Dict{Int64,Date},startdate::Date,enddate::Date)
    filter(date -> startdate <date[2] < enddate,mastercalendar)
 end

@enum STATUS free = 1 booked = 2

@enum WEEK Odd = 1 Even = 2


text:
BUY = true
SELL = false


text:
using Flux, Optimisers, XLA
using Flux.Data.MNIST
using Flux: onehotbatch, crossentropy
using Flux.Data: DataLoader

images = reduce(hcat, [vec(Float64.(im)) for im in MNIST.images()])
images = reshape(images, 28, 28, 1, :)
labels = Float64.(onehotbatch(MNIST.labels(), 0:9))

opt = Optimisers.Descent(0.1)

function step(m, x, y)
  m̄, = gradient(m) do m
    @show crossentropy(m(x), y)
  end
  return opt(m, m̄)
end

m = Chain(
  Conv((3, 3), 1=>16, pad=(1,1), relu),
  MaxPool((2,2)),
  Conv((3, 3), 16=>32, pad=(1,1), relu),
  MaxPool((2,2)),
  Conv((3, 3), 32=>32, pad=(1,1), relu),
  MaxPool((2,2)),
  flatten,
  Dense(288, 10), softmax) |> f64

xstep = xla(step)

for (x, y) in DataLoader((images, labels), batchsize = 100)
  global m = xstep(m, x, y)
end


text:
# provide a structure to extract the underlying MDP of a RPOMDP

struct UnderlyingMDP{P <: RPOMDP, S, A} <: MDP{S, A}
    rpomdp::P
end

function UnderlyingMDP{S, A, O}(rpomdp::RPOMDP{S, A, O})
    P = typeof(rpomdp)
    return UnderlyingMDP{P, S, A}(rpomdp)
end

RPOMDPs.transition{P, S, A}(mdp::UnderlyingMDP{P, S, A}, s::S, a::A) = transition(mdp.rpomdp, s, a)
RPOMDPs.initial_state_distribution(mdp::UnderlyingMDP) = initial_state_distribution(mdp.rpomdp)
RPOMDPs.generate_s(mdp::UnderlyingMDP, s, a, rng::AbstractRNG) = generate_s(mdp.rpomdp, s, a, rng)
RPOMDPs.generate_sr(mdp::UnderlyingMDP, s, a, rng::AbstractRNG) = generate_sr(mdp.rpomdp, s, a, rng)
RPOMDPs.initial_state(mdp::UnderlyingMDP, rng::AbstractRNG) = initial_state(mdp.rpomdp, rng)
RPOMDPs.states(mdp::UnderlyingMDP) = states(mdp.rpomdp)
RPOMDPs.actions(mdp::UnderlyingMDP) = actions(mdp.rpomdp)
RPOMDPs.reward{P, S, A}(mdp::UnderlyingMDP{P, S, A}, s::S, a::A) = reward(mdp.rpomdp, s, a)
RPOMDPs.reward{P, S, A}(mdp::UnderlyingMDP{P, S, A}, s::S, a::A, sp::S) = reward(mdp.rpomdp, s, a, sp)
RPOMDPs.isterminal{P, S, A}(mdp ::UnderlyingMDP{P, S, A}, s::S) = isterminal(mdp.rpomdp, s)
RPOMDPs.discount(mdp::UnderlyingMDP) = discount(mdp.rpomdp)
RPOMDPs.n_actions(mdp::UnderlyingMDP) = n_actions(mdp.rpomdp)
RPOMDPs.n_states(mdp::UnderlyingMDP) = n_states(mdp.rpomdp)
RPOMDPs.state_index{P, S, A}(mdp::UnderlyingMDP{P, S, A}, s::S) = state_index(mdp.rpomdp, s)
RPOMDPs.state_index{P, A}(mdp::UnderlyingMDP{P, Int, A}, s::Int) = state_index(mdp.rpomdp, s) # fix ambiguity with src/convenience
RPOMDPs.state_index{P, A}(mdp::UnderlyingMDP{P, Bool, A}, s::Bool) = state_index(mdp.rpomdp, s)
RPOMDPs.action_index{P, S, A}(mdp::UnderlyingMDP{P, S, A}, a::A) = action_index(mdp.rpomdp, a)
RPOMDPs.action_index{P, S}(mdp::UnderlyingMDP{P,S, Int}, a::Int) = action_index(mdp.rpomdp, a)
RPOMDPs.action_index{P, S}(mdp::UnderlyingMDP{P,S, Bool}, a::Bool) = action_index(mdp.rpomdp, a)


text:
"""
NETALIGNMR
----------
    solve the network alignment problem with Klau's algorithm
"""

function netalignmr(S::SparseMatrixCSC{Int64,Int64},w::Vector{Float64},
                       a::Int64,b::Int64,li::Vector{Int64},lj::Vector{Int64},
                        gamma::Float64,stepm::Int64,rtype::Int64,maxiter::Int64,verbose::Bool)

  m = maximum(li)
  n = maximum(lj)
  Matching_Setup = bipartite_matching_setup(w,li,lj,m,n)
  tripi = Matching_Setup.tripi
  matm = Matching_Setup.m
  matn = Matching_Setup.n
  rp = Matching_Setup.rp
  ci = Matching_Setup.ci
  mperm = tripi[tripi.>0]

  U = spzeros(Float64,size(S,1),size(S,2))
  xbest = zeros(Float64,length(w))

  flower = 0.0                      #best lower bound on the solution
  fupper = Inf                      #best upper bound on the solution
  next_reduction_iteration = stepm  #reduction step

  hist = zeros(Float64,maxiter,7)
  iter = 1
  @inbounds for iter = 1:maxiter
    Q = (b/2)*S + U-U'
    Qt = Q'
    Qp = Qt.colptr
    Qr = Qt.rowval
    Qv = Qt.nzval
    M = size(Qt,1)
    N = size(Qt,2)
    nedges = length(li)

    all_matching = column_maxmatchsum(M,N,Qp,Qr,Qv,m,n,nedges,li,lj)
    q = all_matching.q
    qj = all_matching.mi
    qi = all_matching.mj
    medges = all_matching.medges

    SM = sparse(qi[1:medges],qj[1:medges],1,size(Q,1),size(Q,2))

    x = a*w + q
    ai = zeros(Float64,length(tripi))

    ai[tripi.>0] = x[mperm]

    M_output = MatrixNetworks.bipartite_matching_primal_dual(rp,ci,ai,matm,matn)
    mi = MatrixNetworks.edge_indicator(M_output,li,lj)
    val = M_output.weight

    # compute stats
    matchval = dot(mi,w)
    overlap = dot(mi,(S*mi)/2)
    card = M_output.cardinality
    f = a*matchval + b*overlap

    if val < fupper
        fupper = val
        next_reduction_iteration = iter+stepm
    end

    if f > flower
        flower = f
        itermark = "*"
        xbest = convert(Vector{Float64},mi)
    else
        itermark = " "
    end

    if rtype == 1
        # no work
    elseif rtype==2
        mw = S*x
        mw = a*w + b/2*mw

        ai = zeros(Float64,length(tripi))
        ai[tripi.>0] = mw[mperm]

        M_output2 = MatrixNetworks.bipartite_matching_primal_dual(rp,ci,ai,matm,matn)
        mx = MatrixNetworks.edge_indicator(M_output2,li,lj)
        card = M_output2.cardinality
        matchval = dot(w,mx)
        overlap = dot(mx,(S*mx)/2)

        f = a*matchval + b*overlap

        if f > flower
            flower = f
            itermark = "**"
            mi = mx
            xbest = mw
        end
    end

    # report on current iter
    hist[iter,1:end] = [norm(nonzeros(U),1), flower, fupper, f, matchval, card, overlap]

    # the below if statement causes type instability due to @printf being type instable
    if verbose
        @printf("%5s   %4i   %8.1e   %7.2f %7.2f %7.2f  %7.2f %7.2f %7i %7i\n",
            itermark, iter, norm(nonzeros(U),1),
            flower, fupper, val,
            f, matchval, card, overlap)
    end

    if iter == next_reduction_iteration
        gamma = gamma*0.5

        # the below if statement causes type instability due to @printf being type instable
        if verbose
            @printf("%5s   %4s   reducing step to %f\n", "", "", gamma);
        end
        if gamma < 1e-24
            break
        end
        next_reduction_iteration = iter+stepm
    end

    if (fupper-flower) < 1e-2
        break
    end

    Wtemp = sparse(1:length(mi),1:length(mi),gamma*mi)
    U = U - Wtemp*triu(SM) + tril(SM)'*Wtemp
    Utemp1 = spones(U)
    Utemp1 *= 0.5
    U = min(U,Utemp1)
    Utemp1 *= -1
    U = max(U,Utemp1)
  end
  status = zeros(Float64,2)

  st = (fupper-flower) < 1e-2
  status[1] = flower
  status[2] = fupper
  hist = hist[1:iter,:]
  return (xbest,st,status,hist)

end

############################
### Additional functions ###
############################

function netalignmr(S::SparseMatrixCSC{Int64,Int64},w::Vector{Float64},
                       a::Int64,b::Int64,li::Vector{Int64},lj::Vector{Int64},
                        gamma::Float64,stepm::Int64,rtype::Int64,maxiter::Int64)
    return netalignmr(S,w,a,b,li,lj,gamma,stepm,rtype,maxiter,false)
end

function netalignmr(S::SparseMatrixCSC{Int64,Int64},w::Vector{Float64},
                       a::Int64,b::Int64,li::Vector{Int64},lj::Vector{Int64},
                        gamma::Float64,stepm::Int64,rtype::Int64)
    return netalignmr(S,w,a,b,li,lj,gamma,stepm,rtype,1000,false)
end

function netalignmr(S::SparseMatrixCSC{Int64,Int64},w::Vector{Float64},
                       a::Int64,b::Int64,li::Vector{Int64},lj::Vector{Int64},
                        gamma::Float64,stepm::Int64)
    return netalignmr(S,w,a,b,li,lj,gamma,stepm,1,1000,false)
end

function netalignmr(S::SparseMatrixCSC{Int64,Int64},w::Vector{Float64},
                       a::Int64,b::Int64,li::Vector{Int64},lj::Vector{Int64},
                        gamma::Float64)
    return netalignmr(S,w,a,b,li,lj,gamma,25,1,1000,false)
end

function netalignmr(S::SparseMatrixCSC{Int64,Int64},w::Vector{Float64},
                       a::Int64,b::Int64,li::Vector{Int64},lj::Vector{Int64})
    return netalignmr(S,w,a,b,li,lj,0.4,25,1,1000,false)
end

max_stars_repo_name: nassarhuda/NetworkAlign

max_issues_repo_name: nassarhuda/NetworkAlign

max_forks_repo_name: nassarhuda/NetworkAlign


text:
function harmonic_oscillator(; X=Universe(2))
    A = [0.0 1.0; -1.0 0.0]
    prob = @ivp(x' = Ax, x(0) ∈ Singleton([1.0, 0.0]), x ∈ X)
    tspan = (0.0, 20.0)
    return prob, tspan
end


text:
struct LLGProblem{uType,tType,P,NP,K,ND} <: AbstractSDEProblem{uType,tType,true,ND}
    f::Function
    g::Function
    α::Number
    vars::VarsCart
    u0::uType
    tspan::tType
    kbT::Number
    p::P
    noise_rate_prototype::ND
    noise::NP
    kwargs::K
    seed::UInt64

    @add_kwonly function LLGProblem(A, α, vars, u0, tspan, p = NullParameters() ;
    kbT = 0., noise_rate_prototype = nothing,
        noise = nothing, seed = UInt64(0),
        kwargs...)
        g = LLG.noiseLLG(vars)
        f = convert(SciMLBase.SDEFunction{true}, A, g)
        _tspan = promote_tspan(tspan)
        new{typeof(u0),typeof(_tspan),typeof(p),typeof(noise),typeof(kwargs),typeof(nothing)}(f, g, α, vars, u0, tspan, kbT, p, noise_rate_prototype, noise, kwargs, seed)
    end
end

abstract type LLGAlgorithm <: StochasticDiffEqAlgorithm end
struct SIB <: LLGAlgorithm end

alg_order(alg::SIB) = 1
alg_compatible(prob::LLGProblem,alg::LLGAlgorithm) = true

@cache struct SIBCache{uType} <: StochasticDiffEqMutableCache
    up::uType
    Acache::uType
    N::Int64
    kbT::Float64
    α::Float64
    detM::Vector{Float64}
    detMx::Vector{Float64}
    detMy::Vector{Float64}
    detMz::Vector{Float64}
end

function alg_cache(alg::SIB, prob::LLGProblem, u, dW, dZ, p, rate_prototype, noise_rate_prototype, jump_prototype, uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits, uprev, f, t, dt, iip)
    Acache = similar(u)
    up = similar(u)
    N = size(u)[1]
    detM = zeros(Float64, N)
    detMx = zeros(Float64, N)
    detMy = zeros(Float64, N)
    detMz = zeros(Float64, N)
    SIBCache(up, Acache, N, prob.kbT, prob.α, detM, detMx, detMy, detMz)
end

function DiffEqBase.__solve(prob::LLGProblem,
    alg::LLGAlgorithm,
    timeseries = [],ts = [],ks = nothing, # needed for variable rate
recompile::Type{Val{recompile_flag}} = Val{true};
    kwargs...) where recompile_flag
    integrator = DiffEqBase.__init(prob, alg, timeseries, ts, recompile;kwargs...)
    solve!(integrator)
    integrator.sol
end

@fastmath @muladd function perform_step!(integrator, cache::SIBCache)
    @unpack t, dt, u, f, W = integrator
    # @unpack M, Mx, My, Mz, a, Acache, N, up  = cache
    @unpack Acache, N, up, α, kbT, detM, detMx, detMy, detMz = cache

    # do the predictor step:
    f(Acache, u, t)
    @. Acache = .5dt * (Acache + W.dW * √(2α * kbT))

    fdet(x11, x21, x31, x12, x22, x32, x13, x23, x33) = x11 * (x22 * x33 - x23 * x32) - x12 * (x21 * x33 - x23 * x31) + x13 * (x21 * x32 - x22 * x31)

    combined_cross(uloc, Ax, Ay, Az) = (uloc[1] + uloc[2] * Az - Ay * uloc[3], uloc[2] + uloc[3] * Ax - uloc[1] * Az, uloc[3] + uloc[1] * Ay - uloc[2] * Ax)

    @inbounds for i in 1:N
        uloc = @view u[i,:]
        Ax, Ay, Az = @view Acache[i,:]
        a1, a2, a3 = combined_cross(uloc, Ax, Ay, Az)
        detMx[i] =  fdet(a1, a2, a3, -Az, 1, Ax, Ay, -Ax, 1)
        detMy[i] = fdet(1, Az, -Ay, a1, a2, a3, Ay, -Ax, 1)
        detMz[i] = fdet(1, Az, -Ay, -Az, 1, Ax, a1, a2, a3)
        detM[i] = 1 / fdet(1, Az, -Ay, -Az, 1, Ax, Ay, -Ax, 1)
    end
    @. up[:,1] = detMx * detM
    @. up[:,2] = detMy * detM
    @. up[:,3] = detMz * detM

    # # and the final step:
    @. up = .5 * ( up + u)
    f(Acache,  up, t + .5dt)
    @. Acache = .5dt * (Acache + W.dW * √(2α * kbT))

    @inbounds for i in 1:N
        uloc = @view u[i,:]
        Ax, Ay, Az = @view Acache[i,:]
        a1, a2, a3 = combined_cross(uloc, Ax, Ay, Az)
        detMx[i] =  fdet(a1, a2, a3, -Az, 1, Ax, Ay, -Ax, 1)
        detMy[i] = fdet(1, Az, -Ay, a1, a2, a3, Ay, -Ax, 1)
        detMz[i] = fdet(1, Az, -Ay, -Az, 1, Ax, a1, a2, a3)
        detM[i] = 1 / fdet(1, Az, -Ay, -Az, 1, Ax, Ay, -Ax, 1)
    end
    @. u[:,1] = detMx * detM
    @. u[:,2] = detMy * detM
    @. u[:,3] = detMz * detM
    # @inbounds for i in 1:N
    #     # Ax, Ay, Az = @view Acache[:,i]
    #     # uloc = @view u[:,i]
    #     # # a = SVector{3}(uloc + cross(uloc, @view Acache[:,i]))
    #     # a .= uloc .+ cross(uloc, @view Acache[:,i])
    #     # # up[:,i] .= calc_up(Ax, Ay, Az, a)
    #     calc_up!(@view(up[:,i]), @view( Acache[:,i]), @view(u[:,i]))
# end
end
# @muladd function perform_step!(integrator, cache::SIBCache)
#     @unpack t, dt, u, f, W = integrator
#     # @unpack M, Mx, My, Mz, a, Acache, N, up  = cache
#     @unpack Acache, N, up, α, kbT = cache

#     # do the predictor step:
#     f(Acache, u, t)
#     Acache .+= W.dW * √(2α * kbT)
#     Acache .*= .5dt

#     fdet(x11, x21, x31, x12, x22, x32, x13, x23, x33) = x11 * (x22 * x33 - x23 * x32) - x12 * (x21 * x33 - x23 * x31) + x13 * (x21 * x32 - x22 * x31)

#     combined_cross(uloc, Ax, Ay, Az) = (uloc[1] + uloc[2] * Az - Az * uloc[2], uloc[2] + uloc[3] * Ay - uloc[1] * Az, uloc[3] + uloc[1] * Ay - uloc[2] * Ax)
#     function calc_up!(res, Acachei, uloc)
#         Ax, Ay, Az = Acachei
#         # a = uloc .+ cross(uloc, Acachei)
#         a1, a2, a3 = combined_cross(uloc, Ax, Ay, Az)
#         detMx =  fdet(a1, -Az, Ay, a2, 1, -Ax, a3, Ax, 1)

#         detMy = fdet(1, a1, Ay, Az,  a2, -Ax, -Ay, a3, 1)
#         detMz = fdet(1, -Az, a1, Az, 1, a2, -Ay, Ax, a3)
#         detM = fdet(1, -Az, Ay, Az, 1, -Ax, -Ay, Ax, 1)

#         # detMx = a[1] * (1 + Ax^2) + Az * (a[2] * Ax - a[3]) + Ay * (a[2] * Ax - a[3])
#         # detMy = det(@SMatrix  [  1  a[1] Ay
#         #         Az  a[2] -Ax
#         #         -Ay a[3] 1
#         #         ] )
#         # detMz = det(@SMatrix  [  1 -Az a[1]
#         #     Az 1 a[2]
#         #     -Ay Ax a[3]
#         #     ] )
#         # detM = det(@SMatrix [1 -Az Ay
#         # Az 1 -Ax
#         # -Ay Ax 1])

#         res .= [detMx / detM, detMy / detM, detMz / detM]
#     end

#     @inbounds for i in 1:N
#         calc_up!(@view(up[:,i]), @view( Acache[:,i]), @view(u[:,i]))
#     end

#     # and the final step:
#     f(Acache, .5 * (up + u), t + dt / 2)
#     Acache += W.dW * √(2α * kbT)
#     Acache *= dt / 2
#     @inbounds for i in 1:N
#         # Ax, Ay, Az = @view Acache[:,i]
#         # uloc = @view u[:,i]
#         # # a = SVector{3}(uloc + cross(uloc, @view Acache[:,i]))
#         # a .= uloc .+ cross(uloc, @view Acache[:,i])
#         # # up[:,i] .= calc_up(Ax, Ay, Az, a)
#         calc_up!(@view(u[:,i]), @view( Acache[:,i]), @view(u[:,i]))
#     end
# end
# @muladd function perform_step!(integrator, cache::SIBCache)
#     @unpack t, dt, u, f, W = integrator
#     # @unpack M, Mx, My, Mz, a, Acache, N, up  = cache
#     @unpack Acache, N, up, α, kbT  = cache

#     # do the predictor step:
#     f(Acache, u, t)
#     Acache += W.dW * √(2α * kbT)
#     Acache *= dt / 2

#     @inbounds for i in 1:N
#         Ax, Ay, Az = @view Acache[i,:]
#         uloc = @view u[i,:]
#         a = SVector{3}(uloc + cross(uloc, Acache[i,:]))
#         M = @SMatrix [1 -Az Ay
#             Az 1 -Ax
#         -Ay Ax 1]

#         Mx = @SMatrix [  a[1] -Az Ay
#         a[2] 1 -Ax
#         a[3] Ax 1
#         ]
#         My = @SMatrix  [  1  a[1] Ay
#                 Az  a[2] -Ax
#                 -Ay a[3] 1
#                 ]
#         Mz = @SMatrix  [  1 -Az a[1]
#                 Az 1 a[2]
#                 -Ay Ax a[3]
#                 ]

#         up[i,:] .= [det(Mx) / det(M), det(My) / det(M), det(Mz) / det(M)]
#         end

#         # # and the final step:
#     f(Acache, .5 * (up + u), t + dt / 2)
#     Acache += W.dW * √(2α * kbT)

#     Acache *= dt / 2

#     @inbounds for i in 1:N
#         Ax, Ay, Az = @view Acache[i,:]
#         uloc = @view u[i,:]
#         a = SVector{3}(uloc + cross(uloc, Acache[i,:]))
#         M = @SMatrix [1 -Az Ay
#             Az 1 -Ax
#         -Ay Ax 1]

#         Mx = @SMatrix [  a[1] -Az Ay
#         a[2] 1 -Ax
#         a[3] Ax 1
#         ]
#         My = @SMatrix  [  1  a[1] Ay
#                 Az  a[2] -Ax
#                 -Ay a[3] 1
#                 ]
#         Mz = @SMatrix  [  1 -Az a[1]
#                 Az 1 a[2]
#                 -Ay Ax a[3]
#                 ]

#         u[i,:] .= [det(Mx) / det(M), det(My) / det(M), det(Mz) / det(M)]
#     end
# end


text:
# A collection of polygon domains

export Rectangle,
    RegularPolygon,
    EquilateralTriangle,
    Square,
    Pentagon,
    Hexagon

# We first define an open line segment, this represents the side of a polygon
# given in counter clockwise order.

"An open line segment domain in `ℝ^2`."
struct OpenLineSegment{T} <: Domain{SVector{2,T}}
    A   ::  SVector{2,T}       # Two endpoints
    B   ::  SVector{2,T}

    OpenLineSegment{T}(A = SVector(-one(T), zero(T)),
            B = SVector(one(T), zero(T))) where T = new(A, B)
end

OpenLineSegment() = OpenLineSegment{Float64}()
OpenLineSegment(A::SVector{2,T}, B::SVector{2,T}) where T = OpenLineSegment{T}(A, B)
OpenLineSegment(x1::T, y1::T, x2::T, y2::T) where T = OpenLineSegment{T}(SVector(x1, y1), SVector(x2, y2))

indomain(x, ::OpenLineSegment) = error("Don't know how to evaluate indomain for an open line segment.")

"""
Parameterization of the open line segment domain, mapping from `[0,1]` to the open line segment.
"""
struct OpenLineSegmentMap{T} <: Parameterization{T}
    A   ::  SVector{2,T}       # Two endpoints
    B   ::  SVector{2,T}
end

hasparameterization(d::OpenLineSegment) = true
parameterization(d::OpenLineSegment) = OpenLineSegmentMap(d.A, d.B)

domain(m::OpenLineSegmentMap{T}) where T = zero(T)..one(T)
image(m::OpenLineSegmentMap) = OpenLineSegment(m.A, m.B)

_ols_applymap(A::SVector{2,T}, B::SVector{2,T}, t) where T =
    (1 - t)*A + t*B
applymap(m::OpenLineSegmentMap, t) = _ols_applymap(m.A, m.B, t)

function _ols_gradient(A::SVector{2,T}, B::SVector{2,T}, t) where T
    # As gradient we take the unit vector perpendicular on AB
    # facing to the left of AB
    t = B - A
    tu = norm(t) \ t
    SVector(tu[2], tu[1])
end
gradient(m::OpenLineSegmentMap, t) = _ols_gradient(A, B, t)

# Small utility function that probably is somewhere in the standard library
# Map a range to a unit range
function _map2ur(i, b, e)
    (i - b) / (e - b)
end

# A rectangle

"A rectangular domain in `ℝ^2`."
struct Rectangle{T} <: Domain{SVector{2,T}}
    ll :: SVector{2,T}   # Two endpoints
    ur :: SVector{2,T}

    function Rectangle{T}(A = SVector(-one(T), -one(T)), B = SVector(one(T), one(T))) where T
        new(SVector(min(A[1], B[1]), min(A[2], B[2])),
            SVector(max(A[1], B[1]), max(A[2], B[2])))
    end
end

Rectangle() = Rectangle{Float64}()
Rectangle(A::SVector{2,T}, B::SVector{2,T}) where T = Rectangle{T}(A, B)
Rectangle(x1::T, y1::T, x2::T, y2::T) where T = Rectangle{T}(SVector(x1, y1), SVector(x2, y2))

indomain(x::SVector{2,T}, d::Rectangle{T}) where T =
    d.ll[1] <= x[1] <= d.ur[1] && d.ll[2] <= x[2] <= d.ur[2]

"""
Parameterization of a rectangular domain, mapping from `[0,1)` to the rectangle.
"""
struct RectangleMap{T} <: Parameterization{T}
    ll :: SVector{2,T}
    ur :: SVector{2,T}
end

hasparameterization(d::Rectangle) = true
parameterization(d::Rectangle) = RectangleMap(d.ll, d.ur)

domain(m::RectangleMap{T}) where T = Interval{:closed,:open,T}(0, 1)
image(m::RectangleMap) = Rectangle(m.ll, m.ur)

function _getpart(m::RectangleMap, t)
    w, h = m.ur[1] - m.ll[1], m.ur[2] - m.ll[2]
    c = 2w + 2h

    te = w / c
    t < te && return m.ll, SVector(m.ur[1], m.ll[2]), 0, te

    ts, te = te, (w + h) / c
    t < te && return SVector(m.ur[1], m.ll[2]), m.ur, ts, te

    ts, te = te, (2w + h) / c
    t < te && return m.ur, SVector(m.ll[1], m.ur[2]), ts, te

    SVector(m.ll[1], m.ur[2]), m.ll, te, 1
end

function applymap(m::RectangleMap, t)
    s, e, st, et = _getpart(m, t)
    _ols_applymap(s, e, _map2ur(t, st, et))
end

function gradient(m::RectangleMap, t)
    s, e, st, et = _getpart(m, t)
    _ols_gradient(s, e, _map2ur(t, st, et))
end

# Regular polygons

"A regular polygon domain in `ℝ^2`."
struct RegularPolygon{T} <: Domain{SVector{2,T}}
    n      :: Int          # number of verticies
    center :: SVector{2,T} # center of polygon
    radius :: T            # the radius of the enclosing circle

    function RegularPolygon{T}(n, center=SVector(zero(T), zero(T)), radius=one(T)) where T
        if n < 3
            error("A regular polygon has 3 or more sides")
        end

        new(n, center, radius)
    end
end

RegularPolygon(n) = RegularPolygon{Float64}(n)
RegularPolygon(n::Int, center::SVector{2,T}, radius::T) where T = RegularPolygon{T}(n, center, radius)
RegularPolygon(n::Int, x::T, y::T, radius::T) where T = RegularPolygon{T}(n, SVector(x, y), radius)

function indomain(x::SVector{2,T}, d::RegularPolygon{T}) where T
    #=x -= d.center # Work relative to the origin
    Δθ = 2π / d.n
    θ0 = 3π/2 - Δθ/2=#

    # TODO
    error("Don't know how to evaluate indomain for regular polygon")
end

"""
Parameterization a regular polygon domain, mapping from `[0,1)` to the polygon.
"""
struct RegularPolygonMap{T} <: Parameterization{T}
    n      :: Int          # number of verticies
    center :: SVector{2,T} # center of the polygon
    radius :: T            # the radius of the enclosing circle
end

hasparameterization(d::RegularPolygon) = true
parameterization(d::RegularPolygon) = RegularPolygonMap(d.n, d.center, d.radius)

domain(m::RegularPolygonMap{T}) where T = Interval{:closed,:open,T}(0, 1)
image(m::RegularPolygonMap) = RegularPolygon(m.n, m.center, m.radius)

function _getpart(m::RegularPolygonMap, t)
    Δt, Δθ = 1 / m.n, 2π / m.n
    θ0 = 3π/2 - Δθ/2

    # The last segment is m.n - 1, min catches t = 1
    seg = min(floor(Int, t * m.n), m.n - 1)

    θ1 = θ0 + seg * Δθ
    y, x = sincos(θ1)
    s = m.radius * SVector(x, y) + m.center

    θ2 = θ1 + Δθ
    y, x = sincos(θ2)
    e = m.radius * SVector(x, y) + m.center

    s, e, seg * Δt, (seg + 1) * Δt
end

function applymap(m::RegularPolygonMap, t)
    s, e, st, et = _getpart(m, t)
    _ols_applymap(s, e, _map2ur(t, st, et))
end

function gradient(m::RegularPolygonMap, t)
    s, e, st, et = _getpart(m, t)
    _ols_gradient(s, e, _map2ur(t, st, et))
end

# A collection of regular polygons

EquilateralTriangle() = RegularPolygon(3)
EquilateralTriangle(center::SVector{2,T}, radius=one(T)) where T = RegularPolygon{T}(3, center, radius)
EquilateralTriangle(x::T, y::T, radius=one(T)) where T = RegularPolygon{T}(3, SVector(x, y), radius)

Square() = RegularPolygon(4)
Square(center::SVector{2,T}, radius=one(T)) where T = RegularPolygon{T}(4, center, radius)
Square(x::T, y::T, radius=one(T)) where T = RegularPolygon{T}(4, SVector(x, y), radius)

Pentagon() = RegularPolygon(5)
Pentagon(center::SVector{2,T}, radius=one(T)) where T = RegularPolygon{T}(5, center, radius)
Pentagon(x::T, y::T, radius=one(T)) where T = RegularPolygon{T}(5, SVector(x, y), radius)

Hexagon() = RegularPolygon(6)
Hexagon(center::SVector{2,T}, radius=one(T)) where T = RegularPolygon{T}(6, center, radius)
Hexagon(x::T, y::T, radius=one(T)) where T = RegularPolygon{T}(6, SVector(x, y), radius)


text:
using CLIMAParameters.Planet: cv_d, T_0

export BoundaryCondition, InitStateBC

export AtmosBC,
    Impenetrable,
    FreeSlip,
    NoSlip,
    DragLaw,
    Insulating,
    PrescribedTemperature,
    PrescribedEnergyFlux,
    Impermeable,
    ImpermeableTracer,
    PrescribedMoistureFlux,
    PrescribedTracerFlux

"""
    AtmosBC(momentum = Impenetrable(FreeSlip())
            energy   = Insulating()
            moisture = Impermeable()
            tracer  = ImpermeableTracer())

The standard boundary condition for [`AtmosModel`](@ref). The default options imply a "no flux" boundary condition.
"""
Base.@kwdef struct AtmosBC{M, E, Q, TR}
    momentum::M = Impenetrable(FreeSlip())
    energy::E = Insulating()
    moisture::Q = Impermeable()
    tracer::TR = ImpermeableTracer()
end

function boundary_state!(nf, atmos::AtmosModel, args...)
    atmos_boundary_state!(nf, atmos.boundarycondition, atmos, args...)
end

function boundary_state!(
    nf::Union{CentralNumericalFluxHigherOrder, CentralNumericalFluxDivergence},
    m::AtmosModel,
    x...,
)
    nothing
end

@generated function atmos_boundary_state!(
    nf,
    tup::Tuple,
    atmos,
    state⁺,
    aux⁺,
    n,
    state⁻,
    aux⁻,
    bctype,
    t,
    args...,
)
    N = fieldcount(tup)
    return quote
        Base.Cartesian.@nif(
            $(N + 1),
            i -> bctype == i, # conditionexpr
            i -> atmos_boundary_state!(
                nf,
                tup[i],
                atmos,
                state⁺,
                aux⁺,
                n,
                state⁻,
                aux⁻,
                bctype,
                t,
                args...,
            ), # expr
            i -> error("Invalid boundary tag")
        ) # elseexpr
        return nothing
    end
end

function atmos_boundary_state!(nf, bc::AtmosBC, atmos, args...)
    atmos_momentum_boundary_state!(nf, bc.momentum, atmos, args...)
    atmos_energy_boundary_state!(nf, bc.energy, atmos, args...)
    atmos_moisture_boundary_state!(nf, bc.moisture, atmos, args...)
    atmos_tracer_boundary_state!(nf, bc.tracer, atmos, args...)
end

function normal_boundary_flux_second_order!(
    nf,
    atmos::AtmosModel,
    fluxᵀn::Vars{S},
    n⁻,
    state⁻,
    diff⁻,
    hyperdiff⁻,
    aux⁻,
    state⁺,
    diff⁺,
    hyperdiff⁺,
    aux⁺,
    bctype::Integer,
    t,
    args...,
) where {S}
    atmos_normal_boundary_flux_second_order!(
        nf,
        atmos.boundarycondition,
        atmos,
        fluxᵀn,
        n⁻,
        state⁻,
        diff⁻,
        hyperdiff⁻,
        aux⁻,
        state⁺,
        diff⁺,
        hyperdiff⁺,
        aux⁺,
        bctype,
        t,
        args...,
    )
end
@generated function atmos_normal_boundary_flux_second_order!(
    nf,
    tup::Tuple,
    atmos::AtmosModel,
    fluxᵀn,
    n⁻,
    state⁻,
    diff⁻,
    hyperdiff⁻,
    aux⁻,
    state⁺,
    diff⁺,
    hyperdiff⁺,
    aux⁺,
    bctype,
    t,
    args...,
)
    N = fieldcount(tup)
    return quote
        Base.Cartesian.@nif(
            $(N + 1),
            i -> bctype == i, # conditionexpr
            i -> atmos_normal_boundary_flux_second_order!(
                nf,
                tup[i],
                atmos,
                fluxᵀn,
                n⁻,
                state⁻,
                diff⁻,
                hyperdiff⁻,
                aux⁻,
                state⁺,
                diff⁺,
                hyperdiff⁺,
                aux⁺,
                bctype,
                t,
                args...,
            ), #expr
            i -> error("Invalid boundary tag")
        ) # elseexpr
        return nothing
    end
end
function atmos_normal_boundary_flux_second_order!(
    nf,
    bc::AtmosBC,
    atmos::AtmosModel,
    args...,
)
    atmos_momentum_normal_boundary_flux_second_order!(
        nf,
        bc.momentum,
        atmos,
        args...,
    )
    atmos_energy_normal_boundary_flux_second_order!(
        nf,
        bc.energy,
        atmos,
        args...,
    )
    atmos_moisture_normal_boundary_flux_second_order!(
        nf,
        bc.moisture,
        atmos,
        args...,
    )
    atmos_tracer_normal_boundary_flux_second_order!(
        nf,
        bc.tracer,
        atmos,
        args...,
    )
end

include("bc_momentum.jl")
include("bc_energy.jl")
include("bc_moisture.jl")
include("bc_initstate.jl")
include("bc_tracer.jl")


text:
@testsuite "uniformscaling" (AT, eltypes)->begin
    eltypes = (ComplexF32, Float32)
    wrappers = (identity, UnitLowerTriangular, UnitUpperTriangular, LowerTriangular, UpperTriangular, Hermitian, Symmetric)

    @testset "UniformScaling $f $T1 $T2" for T1 in eltypes, T2 in eltypes, f in wrappers
        x = ones(T1, 5, 5)
        y = AT(x)

        xw = f(x)
        yw = f(y)

        J = one(T2) * I

        # TODO: remove @allowscalar as soon as there is a proper implementation for Symmetric / Hermitian
        @test @allowscalar collect(xw + J) ≈ collect(yw + J)

        # Missing methods in Base it seems... -(x - I) can be removed when Base supports (I - x)
        @test @allowscalar collect(-(xw - J)) ≈ collect(J - yw)
    end
end