text:
#Delete clusters
function b1_cb_template(widgetptr::Ptr,user_data::Tuple{Single_Channel})

    sc, = user_data
    clus=sc.buf.selected_clus

    if (clus<1) #do nothing if zeroth cluster selected
    else
        delete_cluster(sc.temp,clus)
        deleteat!(sc.sort_list,sc.total_clus+1)
        sc.total_clus -= 1
        sc.buf.selected_clus = 0
        selmodel=Gtk.GAccessor.selection(sc.sort_tv)
        Gtk.select!(selmodel, Gtk.iter_from_index(sc.sort_list,1))
        sc.buf.c_changed=true
    end
    nothing
end

function delete_cluster(c::ClusterTemplate,n)

    for i=1:size(c.templates,1)
        c.templates[i,n] = 0.0
        c.sig_max[i,n] = 0.0
        c.sig_min[i,n] = 0.0
    end

    c.tol[n]=0.0

    if n == c.num
        c.num -= 1
    else
        for i=n:(c.num-1)
            for j=1:size(c.templates,1)
                c.templates[j,i]=c.templates[j,i+1]
                c.sig_max[j,i]=c.sig_max[j,i+1]
                c.sig_min[j,i]=c.sig_min[j,i+1]
            end
            c.tol[i]=c.tol[i+1]
        end
        c.num -= 1
    end

    nothing
end

#Add Unit
function b2_cb_template(widget::Ptr,user_data::Tuple{Single_Channel})

    sc, = user_data

    #Add total number of units and go to that unit
    sc.total_clus += 1
    sc.temp.num += 1

    sc.buf.selected_clus = sc.total_clus
    push!(sc.sort_list,(sc.total_clus,))

    selmodel=Gtk.GAccessor.selection(sc.sort_tv)
    Gtk.select!(selmodel, Gtk.iter_from_index(sc.sort_list, sc.total_clus+1))

    nothing
end

function b3_cb_template(widget::Ptr,user_data::Tuple{Single_Channel,SortView})

    sc, sortview_widgets = user_data

    mybutton = convert(Button, widget)

    if sc.pause
        Gtk.GAccessor.active(sc.pause_button,false)
    end

    if !sc.hold
        clear_c2(sc.c2,sc.spike)
        sc.ctx2=Gtk.getgc(sc.c2)
        sc.ctx2s=copy(sc.ctx2)
        sc.buf.ind=1
        sc.buf.count=1
        sc.hold=true
        sc.pause=false
        change_button_label(mybutton,"Stop Collection")
    else

        Gtk.GAccessor.active(sc.pause_button,true)
        change_button_label(mybutton,"Collect Templates")

        if sc.buf.count==size(sc.buf.spikes,2)
            sc.buf.ind=sc.buf.count
        end

        if visible(sortview_widgets.win)
            recalc_features(sortview_widgets)
            replot_sort(sortview_widgets)
        end
    end

    nothing
end

function b4_cb_template(widgetptr::Ptr,user_data::Tuple{Single_Channel})

    sc, = user_data
    widget = convert(ToggleButton, widgetptr)

    #Untoggle
    if !getproperty(widget,:active,Bool)
        sc.buf.replot=true
    else
        #Toggling, so draw template
        if sc.pause
            draw_template_bounds(sc)
        end
    end

    nothing
end

function draw_template_bounds(sc::Single_Channel)

    ctx = sc.ctx2
    clus=sc.buf.selected_clus

     if clus>0
            s=sc.s
            o=sc.o

            Cairo.translate(ctx,0.0,sc.h2/2)
            Gtk.scale(ctx,sc.w2/sc.wave_points,s)

            move_to(ctx,1.0,sc.temp.templates[1,clus]+(sc.temp.sig_max[1,clus]*sc.temp.tol[clus])-o)

            for i=2:size(sc.temp.templates,1)
                y=sc.temp.templates[i,clus]+(sc.temp.sig_max[i,clus]*sc.temp.tol[clus])-o
                line_to(ctx,i,y)
            end

            y=sc.temp.templates[end,clus]-(sc.temp.sig_min[end,clus]*sc.temp.tol[clus])-o

            line_to(ctx,size(sc.temp.templates,1),y)

            for i=(size(sc.temp.templates,1)-1):-1:1
                y=sc.temp.templates[i,clus]-(sc.temp.sig_min[i,clus]*sc.temp.tol[clus])-o
                line_to(ctx,i,y)
            end

            close_path(ctx)

            select_color(ctx,clus+1)
            set_line_width(ctx,3.0)
            stroke_preserve(ctx)

            select_color(ctx,clus+1,.5)
            fill(ctx)
        end

    identity_matrix(ctx)
end

function check_cb_template(widget::Ptr,user_data::Tuple{Single_Channel})

    sc, = user_data

    mycheck=convert(CheckButton,widget)

    if getproperty(mycheck,:active,Bool)
        sc.sort_cb=true
    else
        sc.sort_cb=false
    end

    if sc.sort_cb
        draw_templates(sc)
    end

    nothing
end

function draw_templates(sc::Single_Channel)

    ctx = sc.ctx2s
    mywidth=width(ctx)
    myheight=height(ctx)

    s=sc.s
    o=sc.o

    Cairo.translate(ctx,0.0,myheight/2)
    Gtk.scale(ctx,mywidth/sc.wave_points,s)

    for clus=1:sc.total_clus

        move_to(ctx,1.0,(sc.temp.templates[1,clus])-o)

        for i=2:size(sc.temp.sig_max,1)
            y=sc.temp.templates[i,clus]-o
            line_to(ctx,i,y)
        end

        select_color(ctx,clus+1)
        set_line_width(ctx,3.0)
        stroke(ctx)
    end
    identity_matrix(ctx)
    nothing
end

function get_template_dims(sc::Single_Channel,clus)

    ctx=Gtk.getgc(sc.c3)
    mywidth=width(ctx)

    total_clus = max(sc.total_clus+1,5)

    xbounds=linspace(0.0,mywidth,total_clus+1)

    (xbounds[clus],xbounds[clus+1],0.0,130.0)
end

#=
Calculates the mean and bounds for a collection of spikes that
1) are not hidden by the mask vector
2) meet some condition

examples of this condition could be whether the spikes are contained within a rubberband, or if they are part of a certain cluster already
=#
function make_cluster(input::Array{T,2},mask,count,condition) where T

    hits=0
    mymean=zeros(Float64,size(input,1)-1)
    mysum=zeros(Int64,size(input,1)-1)
    mybounds=zeros(Float64,size(input,1)-1,2)

    for i=1:count
        if (condition[i])&&(mask[i])
            hits+=1
            for ii=1:length(mymean)
                mysum[ii] += input[ii,i]
                if hits==1
                    mybounds[ii,1]=input[ii,i]
                    mybounds[ii,2]=input[ii,i]
                else
                    if input[ii,i]>mybounds[ii,1]
                        mybounds[ii,1]=input[ii,i]
                    end
                    if input[ii,i]<mybounds[ii,2]
                        mybounds[ii,2]=input[ii,i]
                    end
                end
            end
        end
    end

    if hits==0
        hits=1
    end

    for i=1:length(mymean)
        mymean[i] = mysum[i]/hits
        mybounds[i,2] = abs(mymean[i]-mybounds[i,2])
        mybounds[i,1] = abs(mybounds[i,1] - mymean[i])
    end

    (mymean,mybounds)
end

function change_cluster(c::ClusterTemplate,mymean::Array{Float64,1},mystd::Array{Float64,2},n)

    for i=1:length(mymean)
        c.templates[i,n] = mymean[i]
        c.sig_max[i,n] = mystd[i,1]
        c.sig_min[i,n] = mystd[i,2]
    end

    nothing
end

function template_cluster(sc::Single_Channel,clus,mymean::Array{Float64,1},mymin::Array{Float64,1},mymax::Array{Float64,1},tol::Float64)

    @inbounds for i=1:sc.buf.ind

        mymisses=0
        for j=1:length(mymean)
            if (sc.buf.spikes[j,i]<(mymean[j]-(mymin[j]*tol)))|(sc.buf.spikes[j,i]>(mymean[j]+(mymax[j]*tol)))
                mymisses+=1
                if mymisses>1
                    break
                end
            end
        end
        if mymisses<1 #If passes template matching, set as unit
            sc.buf.clus[i]=clus
        elseif sc.buf.clus[i]==clus #If did not pass, but was previously, set to noise cluster
            sc.buf.clus[i]=0
        end
    end

    nothing
end

function send_clus(s::Array{T,1},sc::Single_Channel) where T<:Sorting
    s[sc.spike].c=deepcopy(sc.temp)
    nothing
end

function send_clus(s::DArray{T,1,Array{T,1}},sc::Single_Channel) where T<:Sorting
    (nn,mycore)=get_thres_id(s,sc.spike)
    remotecall_wait(((x,tt,num)->localpart(x)[num].c=tt),mycore,s,sc.temp,nn)
    nothing
end


text:
using Test
using Libxc
using DFTK: xc_fallback!

@testset "Fallback LDA" begin
    for func_name in (:lda_x, :lda_c_vwn)
        func = Libxc.Functional(func_name)

        # Create reference
        ρ = abs.(randn(100))
        Eref  = similar(ρ)
        Vref  = similar(ρ)
        V2ref = similar(ρ)
        Libxc.evaluate!(func, rho=ρ, zk=Eref, vrho=Vref, v2rho2=V2ref)

        # Compute in DFTK in elevated precision
        ρbig = Array{BigFloat}(ρ)
        E    = similar(ρbig)
        V    = similar(ρbig)
        V2   = similar(ρbig)
        xc_fallback!(func, Val(:lda), ρbig; zk=E, vrho=V, v2rho2=V2)

        @test E  ≈ Eref  atol=5e-15
        @test V  ≈ Vref  atol=5e-15
        @test V2 ≈ V2ref atol=5e-15
    end
end

ext: jl


text:
tobinary = Dict("0" => "0000", "1" => "0001", "2" => "0010", "3" => "0011", "4" => "0100", "5" => "0101", "6" => "0110", "7" => "0111", "8" => "1000", "9" => "1001", "A" => "1010", "B" => "1011", "C" => "1100", "D" => "1101", "E" => "1110", "F" => "1111")
version_sum = 0

function parseInput(filename)
    line = readlines(filename*".txt")
    binline = ""
    for c in line[1]
        b = tobinary[string(c)]
        binline = binline * b
    end
    return binline
end

function interpretHeader(binline, startidx)
    version = parse(Int, binline[startidx:startidx+2], base = 2)
    type = parse(Int, binline[startidx+3:startidx+5], base = 2)
    return version, type
end

function interpretLiteral(binline, startidx)
    #pkglength = 6
    fivebits = binline[startidx:startidx+4]
    number = ""
    stopdesignator = false
    while stopdesignator != true
        number = number * fivebits[2:5]
        #pkglength += 5
        if fivebits[1] == '0'
            stopdesignator = true
        end
        startidx += 5
        if startidx+4 < length(binline)
            fivebits = binline[startidx:startidx+4]
        end
    end
    number = parse(Int, number, base = 2)

    return number, startidx
end

function readInstructions(binline, idx)
    global version_sum
    packet = Dict()
    v, t = interpretHeader(binline, idx)
    version_sum += v
    packet["version"] = v
    packet["type"] = t
    idx = idx + 6
    if packet["type"] == 4
        packet["val"], idx = interpretLiteral(binline, idx)
    else
        packet["i"] = parse(Int, binline[idx], base = 2)
        idx += 1
        if packet["i"] == 0
            packet["length"] = parse(Int, binline[idx:idx+14], base = 2)
            idx = idx + 15
            p_end = idx + packet["length"]
            if haskey(packet, "subpackages")
                sp = sp = [packet["subpackages"]]
            else
                sp = []
            end
            while idx < p_end
                subpack, idx = readInstructions(binline, idx)
                sp = push!(sp, subpack)
                packet["subpackages"] = sp
            end
        else
            packet["sub_count"] = parse(Int, binline[idx:idx+10], base = 2)
            idx = idx + 11
            if haskey(packet, "subpackages")
                sp = sp = [packet["subpackages"]]
            else
                sp = []
            end
            for _ in 1:packet["sub_count"]
                subpack, idx = readInstructions(binline, idx)
                sp = push!(sp, subpack)
                packet["subpackages"] = sp
            end
        end
        subp_vals = [p["val"] for p in packet["subpackages"]]

        if packet["type"] == 0
            #sum
            packet["val"] = sum(subp_vals)
        elseif packet["type"] == 1
            #product
            packet["val"] = prod(subp_vals)
        elseif packet["type"] == 2
            #minimum
            packet["val"] = minimum(subp_vals)
        elseif packet["type"] == 3
            #maximum
            packet["val"] = maximum(subp_vals)
        elseif packet["type"] == 5
            #greater than
            if subp_vals[1] > subp_vals[2]
                packet["val"] = 1
            else
                packet["val"] = 0
            end
        elseif packet["type"] == 6
            #less than
            if subp_vals[1] < subp_vals[2]
                packet["val"] = 1
            else
                packet["val"] = 0
            end
        elseif packet["type"] == 7
            #equal
            if subp_vals[1] == subp_vals[2]
                packet["val"] = 1
            else
                packet["val"] = 0
            end
        end
    end

    return packet, idx
end

function solve(filename)
    bin = parseInput(filename)
    packets, _ = readInstructions(bin, 1)
    println(version_sum)
    #println(packets)
    println(packets["val"])
end

solve("literal")
version_sum = 0
 solve("op1")
version_sum = 0
solve("op2")
version_sum = 0
solve("input")


text:
using Roots, Plots

f(x) = x^3 - x - 1
@show fzeros(f, -10, 10)
println()

y = -10:0.1:10
plot(y, f.(y))

@show fzero(f, [-2,2])
println()

@show fzero(f, 1.0)
println()

f(x) = exp(x) - x^4
fp(x) = exp(x) - 4x^3
fpp(x) = exp(x) - 12x^2

Roots.newton(f, fp, 8)        # 8.613169456441398

halley(f, fp, fpp, 8)

secant_method(f, [8, 8.5])


text:
"""
    JDF.metadata(indir)

Load the metadata associated with the JDF in `indir`
"""
jdfmetadata(indir) = begin
    open(joinpath(indir, "metadata.jls")) do io
        deserialize(io)
    end
end

"""
    JDF.metadata(indir)

Load the metadata associated with the JDF in `indir`
"""
metadata(jdf::JDFFile) = jdfmetadata(path(jdf))

"""
    JDF.names(indir)

Load the column names associated with the JDF in `indir`

# Examples
```julia
using JDF, DataFrames
JDF.save(DataFrame(a = 1:3, b = 1:3), "plsdel.jdf")
JDF.names("plsdel.jdf")
```
"""
names(jdf::JDFFile) = metadata(jdf).names
# names(indir) = jdfmetadata(indir).names

"""
    JDF.size(indir)

Returns the JDF's `size`
"""
# JDF.size(indir) = begin
#     m = metadata(indir)
#     (m.rows, length(m.names))
# end
JDF.size(jdf::JDFFile) = begin
    m = metadata(jdf)
    (m.rows, length(m.names))
end

# JDF.size(indir, v) = JDF.size(indir, Val(v))
JDF.size(jdf::JDFFile, val) = JDF.size(jdf)[val]

max_stars_repo_name: xiaodaigh/JDF

max_issues_repo_name: xiaodaigh/JDF

max_forks_repo_name: xiaodaigh/JDF


text:
# This was formerly a part of Julia. License is MIT: http://julialang.org/license

import Base.Order.Reverse

@testset "PriorityQueue" begin

    # Test dequeing in sorted order.
    function test_issorted!(pq::PriorityQueue, priorities, rev=false)
        last = dequeue!(pq)
        while !isempty(pq)
            value = dequeue!(pq)
            if !rev
                @test priorities[last] <= priorities[value]
            else
                @test priorities[value] <= priorities[last]
            end
            value = last
        end
    end

    function test_isrequested!(pq::PriorityQueue, keys)
        i = 0
        while !isempty(pq)
            krqst =  keys[i+=1]
            krcvd = dequeue!(pq, krqst)
            @test krcvd == krqst
        end
    end

    @testset "Constructors" begin
        pmax = 1000
        n = 10000
        r = rand(1:pmax, n)
        priorities = Dict(zip(1:n, r))

        ks, vs = 1:n+1, rand(1:pmax, n)

        @testset "building from Dict" begin
            pq1 = PriorityQueue(priorities)
            test_issorted!(pq1, priorities)

            pq2 = PriorityQueue(priorities)
            test_isrequested!(pq2, 1:n)

            priorities = Dict(zip(ks, vs))
            pq4 = PriorityQueue(priorities)
            test_issorted!(pq4, priorities)
        end

        @testset "parameterized constructor tests pq5" begin
            pq5 = PriorityQueue{Int,Int}()
            for (k,v) in priorities
                pq5[k] = v
            end
            test_issorted!(pq5, priorities)
        end

        @testset "parameterized constructor tests pq6" begin
            pq6 = PriorityQueue{Int,Int}(priorities)
            test_issorted!(pq6, priorities)
        end

        @testset "parameterized constructor tests pq7" begin
            pq7 = PriorityQueue{Int,Int}(Reverse, priorities)
            test_issorted!(pq7, priorities, true)
        end

        @testset "parameterized constructor tests pq8" begin
            pq8 = PriorityQueue{Int,Int}(Reverse)
            for (k,v) in priorities
                pq8[k] = v
            end
            test_issorted!(pq8, priorities, true)
        end

        @testset "construction from pairs" begin
            @testset "pq9" begin
                pq9 = PriorityQueue('a'=>1, 'b'=>2)
                @test peek(pq9) == ('a'=>1)
            end

            @testset "pq10" begin
                pq10 = PriorityQueue(Reverse, 'a'=>1, 'b'=>2)
                @test peek(pq10) == ('b'=>2)
            end

            @testset "pq11" begin
                pq11 = PriorityQueue(Pair{Char}['a'=>1,'b'=>2])
                @test peek(pq11) == ('a'=>1)
            end
        end

        @testset "duplicate key => ArgumentError" begin
            @test_throws ArgumentError PriorityQueue('a'=>1, 'a'=>2)
        end

        @testset "Not a pair/tuple => ArguentError" begin
            @test_throws ArgumentError PriorityQueue(['a'])
            @test_throws ArgumentError PriorityQueue(Reverse, ['a'])
            @test_throws ArgumentError PriorityQueue{Char,Int}(Base.Order.Reverse, ['a'])
        end

        @testset "Silly test" begin
            @test_throws ArgumentError PriorityQueue(Reverse, Reverse)
        end

    end

    @testset "PriorityQueueMethods" begin
        pq1 = PriorityQueue('a'=>1, 'b'=>2)

        @testset "peek/get/dequeue!/get!" begin
            @test peek(pq1) == ('a'=>1)
            @test get(pq1, 'a', 0) == 1
            @test get(pq1, 'c', 0) == 0
            @test get!(pq1, 'b', 20) == 2
            @test dequeue!(pq1) == 'a'
            @test dequeue!(pq1) == 'b'
            @test get!(pq1, 'c', 0) == 0
            @test peek(pq1) == ('c'=>0)
            @test get!(pq1, 'c', 3) == 0
        end

        pmax = 1000
        n = 10000
        ks, vs = 1:n, rand(1:pmax, n)
        priorities = Dict(zip(ks, vs))

        @testset "peek" begin
            pq1 = PriorityQueue(priorities)
            lowpri = findmin(vs)
            @test peek(pq1)[2] == pq1[ks[lowpri[2]]]
        end

        @testset "enqueue error throw" begin
            ks, vs = 1:n, rand(1:pmax, n)
            pq = PriorityQueue(zip(ks, vs))
            @test_throws ArgumentError enqueue!(pq, 1, 10)
        end

        @testset "Iteration" begin
            pq = PriorityQueue(priorities)
            pq2 = PriorityQueue()
            for kv in pq
                enqueue!(pq2, kv)
            end
            @test pq == pq2
        end

        @testset "enqueing pairs via enqueue!" begin
            pq = PriorityQueue()
            for kv in priorities
                enqueue!(pq, kv)
            end
            test_issorted!(pq, priorities)
        end

        @testset "enqueing values via enqueue!" begin
            pq = PriorityQueue()
            for (k, v) in priorities
                enqueue!(pq, k, v)
            end
            test_issorted!(pq, priorities)
        end

        @testset "enqueing via assign" begin
            pq = PriorityQueue()
            for (k, v) in priorities
                pq[k] = v
            end
            test_issorted!(pq, priorities)
        end

        @testset "changing priorities" begin
            pq = PriorityQueue()
            for (k, v) in priorities
                pq[k] = v
            end

            for _ in 1:n
                k = rand(1:n)
                v = rand(1:pmax)
                pq[k] = v
                priorities[k] = v
            end

            test_issorted!(pq, priorities)
        end

        @testset "dequeuing" begin
            pq = PriorityQueue(priorities)
            @test_throws KeyError dequeue!(pq, 0)

            @test 10 == dequeue!(pq, 10)
            while !isempty(pq)
                @test 10 != dequeue!(pq)
            end
        end

        @testset "dequeuing pair" begin
            priorities2 = Dict(zip('a':'e', 5:-1:1))
            pq = PriorityQueue(priorities2)
            @test_throws KeyError dequeue_pair!(pq, 'g')
            @test dequeue_pair!(pq) == ('e'=> 1)
            @test dequeue_pair!(pq, 'b') == ('b'=>4)
            @test length(pq) == 3
        end

        @testset "delete!" begin
            pq = PriorityQueue(Base.Order.Forward, "a"=>2, "b"=>3, "c"=>1)
            pq_out = delete!(pq, "b")
            @test pq === pq_out
            @test Set(collect(pq)) == Set(["a"=>2, "c"=>1])
        end

        @testset "empty!" begin
            pq = PriorityQueue(Base.Order.Forward, "a"=>2, "b"=>3, "c"=>1)
            @test !isempty(pq)
            empty!(pq)
            @test isempty(pq)
            enqueue!(pq, "a"=>2)
            @test length(pq) == 1
        end
    end

    @testset "Iteration" begin
        pq = PriorityQueue(["a" => 10, "b" => 5, "c" => 15])
        @test collect(pq)         == ["b" => 5, "a" => 10, "c" => 15]
        @test collect(keys(pq))   == ["b", "a", "c"]
        @test collect(values(pq)) == [5, 10, 15]
    end

    @testset "UnorderedIteration" begin
        pq = PriorityQueue(["a" => 10, "b" => 5, "c" => 15])
        res = Pair{String, Int}[]
        next = iterate(pq, false)
        while next !== nothing
            p, i = next
            push!(res, p)
            next = iterate(pq, i)
        end
        @test Set(res) == Set(["a" => 10, "b" => 5, "c" => 15])
        empty!(res)
        next = iterate(pq, true)
        while next !== nothing
            p, i = next
            push!(res, p)
            next = iterate(pq, i)
        end
        @test res == ["b" => 5, "a" => 10, "c" => 15]
    end

    # Copy and merge operations in PriorityQueue utilize unordered
    # iteration
    @testset "Copy and Merge" begin
        v = [i for i=10000:-1:1]
        pq1 = PriorityQueue(zip(v, 10 .* v))
        @test collect(merge!(Dict(), pq1)) != collect(pq1)
        pq = PriorityQueue(["a" => 10, "b" => 5, "c" => 15])
        sd = SortedDict(["d"=>6, "e"=>4])
        merge!(sd, pq)
        @test collect(sd) == ["a" => 10, "b" => 5, "c" => 15, "d" => 6, "e" => 4]
        d = merge!(Dict("d"=> 6), pq)
        @test Set(collect(d)) == Set(["c" => 15, "b" => 5, "a" => 10, "d" => 6])
        sd = SortedDict(["c"=>6, "e"=>4])
        merge!(+, sd, pq)
        @test collect(sd) == ["a" => 10, "b" => 5, "c" => 21, "e" => 4]
        pq1 = empty(pq)
        @test length(pq1) == 0
        pq1 = copy(pq)
        @test pq1 isa PriorityQueue
    end

    @testset "LowLevelHeapOperations" begin
        pmax = 1000
        n = 10000
        r = rand(1:pmax, n)
        priorities = Dict(zip(1:n, r))

        @testset "low level heap operations" begin
            xs = heapify!([v for v in values(priorities)])
            @test issorted([heappop!(xs) for _ in length(priorities)])
        end

        @testset "heapify/issorted" begin
            xs = heapify(10:-1:1)
            @test issorted([heappop!(xs) for _ in 1:10])
        end

        @testset "heappush!" begin
            xs = Vector{Int}()
            for priority in values(priorities)
                heappush!(xs, priority)
            end
            @test issorted([heappop!(xs) for _ in length(priorities)])
        end

        @testset "isheap" begin
            @test isheap([1, 2, 3], Base.Order.Forward)
            @test !isheap([1, 2, 3], Base.Order.Reverse)
        end
    end

end # @testset "PriorityQueue"

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
include("../src/HMat.jl");

function laplacian2d(nx,ny)
    M = speye(nx*ny);
    for x=1:nx
        for y=1:ny
            s = x+(y-1)*nx;
            M[s,s] = 4;
            if x > 1
                M[s,s-1] = -1;
            end
            if x < nx
                M[s,s+1] = -1;
            end
            if y > 1
                M[s,s-nx] = -1;
            end
            if y < ny
                M[s,s+nx] = -1;
            end
        end
    end
    return M;
end

n = 32;
A = full(laplacian2d(n,n));

EPS = 1e-6;
MaxRank = 8;
minn = 4;

Z2C = Z2Cmapper([n,n],minn,reshape(1:n^2,n,n));
A = full(A[Z2C,Z2C]);

@time AH = HMat2dd2h(A, [n,n], [n,n], STANDARD, [0,0], [0,0], 1, EPS, MaxRank, minn);

v = rand(n^2,5);
uext  = A*v;
@time uhmat = hmatvec(AH,v);
@time uhmat = hmatvec(AH,v);
@printf("hmatvec error: %.3e\n",norm(uext-uhmat)/norm(uext));

uext  = A'*v;
@time uhmat = hmatTvec(AH,v);
@time uhmat = hmatTvec(AH,v);
@printf("hmatTvec error: %.3e\n",norm(uext-uhmat)/norm(uext));

uext  = v'*A;
@time uhmat = hvecmat(v',AH);
@time uhmat = hvecmat(v',AH);
@printf("hvecmat error: %.3e\n",norm(uext-uhmat)/norm(uext));

uext  = v'*A';
@time uhmat = hvecmatT(v',AH);
@time uhmat = hvecmatT(v',AH);
@printf("hvecmatT error: %.3e\n",norm(uext-uhmat)/norm(uext));

@time AD = hh2d(AH);
@printf("hh2d error: %.3e\n",norm(AD-A)/norm(A));

BH = hcopy(AH);
@time CH = hadd(AH,BH);
@time CH = hadd(AH,BH);
CD = hh2d(CH);
@printf("hadd error: %.3e\n",norm(CD-A-A)/norm(2*A));

UMat = rand(n^2,MaxRank);
VMat = rand(n^2,MaxRank);
@time CH = hadd(AH,UMat,VMat);
CD = hh2d(CH);
@printf("hadd error: %.3e\n",norm(CD-A-UMat*VMat')/norm(A+UMat*VMat'));

alpha = 3.1;
@time CH = hscale(alpha,AH);
CD = hh2d(CH);
@printf("hscale error: %.3e\n",norm(CD-alpha*A)/norm(alpha*A));

BH = hcopy(AH);
@time CH = hmul(AH,BH);
CD = hh2d(CH);
AA = A*A;
@printf("hmul error: %.3e\n",norm(CD-A*A)/norm(A*A));
Err = CD-A*A;

@time Anormfro = hnorm(AH,"fro");
@printf("hnorm fro error: %.3e\n",norm(Anormfro-vecnorm(A))/vecnorm(A));

max_stars_repo_name: YingzhouLi/HMat

max_issues_repo_name: YingzhouLi/HMat


text:
using Discord
include("base.jl")

function adjust(client::Client,m::Message,hp::Int)
    global management_dict
    channel_id_str = string(m.channel_id)

    if !haskey(management_dict,channel_id_str)
        return reply(client,m,"初期設定が完了しておりません．`.set`コマンドにて初期設定を完了させてください")
    end

    if !haskey(management_dict[channel_id_str],"remain_attack_channel")
        return reply(client,m,"残凸把握板の登録が完了しておりません．`.set`コマンドにて初期設定を完了させてください")
    end

    if  !haskey(management_dict[channel_id_str],"output_channel")
        return reply(client,m,"bot出力板の登録が完了しておりません．`.set`コマンドにて初期設定を完了させてください")
    end

    if hp <= 0
        return reply(client,m,"ボスを変更する場合には`.start`コマンドを使用して下さい")
    end
    boss_number = string(management_dict[channel_id_str]["boss_number"])

    management_dict[channel_id_str]["reserve"][boss_number]["remain_hp"] = hp
    management_dict[channel_id_str]["reserve"][boss_number]["plan_remain_hp"] = calc_plan_remain_hp(management_dict[channel_id_str]["reserve"][boss_number])

    output_message_content = create_output_message(management_dict[channel_id_str])
    edit_message(client,management_dict[channel_id_str]["output_channel"],management_dict[channel_id_str]["output_message"],content=output_message_content)

    save_management_dict(management_dict)
    create(client, Reaction, m, '👍')

end

max_stars_repo_name: hibibol/clan-battle

max_issues_repo_name: hibibol/clan-battle

max_forks_repo_name: hibibol/clan-battle


text:
module MLLibrary

# package code goes here
  files = ["DataParser","Statistics","LearningAlgorithms","DataSampling","DataManipulation","Utils","LinearAlgebra"]
  for file in files
    include(joinpath(file,"$file.jl"))
  end

end # module

max_stars_repo_name: ajermaky/CustomMLLibrary

max_issues_repo_name: ajermaky/CustomMLLibrary

max_forks_repo_name: ajermaky/CustomMLLibrary


text:
function setup_etc(backend::GPUBackend, layer::ChannelPoolingLayer, inputs, blobs)
  if isa(layer.pooling, Pooling.Max)
    masks = Array{CuPtr}(length(inputs))
    for i = 1:length(inputs)
      masks[i] = CUDA.cualloc(Csize_t, length(blobs[i]))
    end
    etc = masks
  elseif isa(layer.pooling, Pooling.Mean)
    integrals = Array{CuPtr}(length(inputs))
    for i = 1:length(inputs)
      integrals[i] = CUDA.cualloc(eltype(inputs[i]), prod(size(inputs[i])[1:end-1]))
    end
    etc = integrals
  else
    etc = nothing
  end
  return etc
end
function shutdown_etc(backend::GPUBackend, state::ChannelPoolingLayerState)
  if isa(state.layer.pooling, Pooling.Max)
    map(CUDA.free, state.etc)
  elseif isa(state.layer.pooling, Pooling.Mean)
    map(CUDA.free, state.etc)
  else
    error("Unknown pooling $(state.layer.pooling)")
  end
end

function forward(backend::GPUBackend, state::ChannelPoolingLayerState, inputs::Vector{Blob})
  forward(backend, state.layer.pooling, state, inputs)
end
function forward(backend::GPUBackend, pool::StdPoolingFunction,
    state::ChannelPoolingLayerState, inputs::Vector{Blob})

  for i = 1:length(inputs)
    input = inputs[i]
    output = state.blobs[i]

    if isa(pool, Pooling.Max)
      cuda_max_channel_pooling_forward(backend, input, output, state.etc[i], state.layer, state.op_dims[i])
    elseif isa(pool, Pooling.Mean)
      cuda_mean_channel_pooling_forward(backend, input, output, state.etc[i], state.layer, state.op_dims[i])
    else
      error("Pooling for $pool not implemented yet")
    end
  end
end

function backward(backend::GPUBackend, state::ChannelPoolingLayerState, inputs::Vector{Blob}, diffs::Vector{Blob})
  backward(backend, state.layer.pooling, state, inputs, diffs)
end

function backward(backend::GPUBackend, pool::StdPoolingFunction, state::ChannelPoolingLayerState,
    inputs::Vector{Blob}, diffs::Vector{Blob})

  for i = 1:length(inputs)
    diff = diffs[i]
    if !isa(diff, NullBlob)
      if isa(pool, Pooling.Max)
        cuda_max_channel_pooling_backward(backend, diff, state.blobs_diff[i], state.etc[i], state.layer, state.op_dims[i])
      elseif isa(pool, Pooling.Mean)
        cuda_mean_channel_pooling_backward(backend, diff, state.blobs_diff[i], state.layer, state.op_dims[i])
      else
        error("Pooling for $pool not implemented yet")
      end
    else
      continue # nothing to do if not propagating back
    end
  end
end

function cuda_mean_channel_pooling_forward{T}(backend::GPUBackend, input::CuTensorBlob{T},
    output::CuTensorBlob{T}, integral::CuPtr, layer, op_dim)

  spatial_dim_T, channels, num = split_dims(input, op_dim)
  pooled_chann = size(output, op_dim)
  one = convert(T, 1)
  neg_one = convert(T, -1)
  scale = convert(T, 1.0/layer.kernel)

  spatial_dim = spatial_dim_T * sizeof(T)
  fea_dim = spatial_dim * channels
  output_fea_dim = spatial_dim * pooled_chann

  for n = 1:num
    input_ptr = convert(Ptr{T}, input.ptr.p + fea_dim*(n-1))
    output_ptr = convert(Ptr{T}, output.ptr.p + output_fea_dim*(n-1))
    integral_ptr = convert(Ptr{T}, integral.p)

    # compute integral image
    CuBLAS.copy(backend.cublas_ctx, T, spatial_dim_T, input_ptr, 1, integral_ptr, 1)
    for c = 2:channels
      CuBLAS.copy(backend.cublas_ctx, T, spatial_dim_T, input_ptr + (c-1)*spatial_dim, 1,
          integral_ptr + (c-1)*spatial_dim, 1)
      CuBLAS.axpy(backend.cublas_ctx, spatial_dim_T, one, integral_ptr + (c-2)*spatial_dim, 1,
          integral_ptr + (c-1)*spatial_dim, 1)
    end

    for pc = 1:pooled_chann
      cstart = (pc-1)*layer.stride - layer.pad[1] + 1
      cend   = min(cstart + layer.kernel - 1, channels)
      cstart = max(1, cstart)

      output_ptr_pc = output_ptr + (pc-1)*spatial_dim

      CuBLAS.copy(backend.cublas_ctx, T, spatial_dim_T, integral_ptr + (cend-1)*spatial_dim, 1,
          output_ptr_pc, 1)
      if cstart > 1
        CuBLAS.axpy(backend.cublas_ctx, spatial_dim_T, neg_one,
            integral_ptr + (cstart-2)*spatial_dim, 1, output_ptr_pc, 1)
      end
      CuBLAS.scal(backend.cublas_ctx, spatial_dim_T, scale, output_ptr_pc, 1)
    end
  end
end

function cuda_mean_channel_pooling_backward{T}(backend::GPUBackend, input::CuTensorBlob{T},
    output::CuTensorBlob{T}, layer, op_dim)

  spatial_dim_T, channels, num = split_dims(input, op_dim)
  pooled_chann = size(output, op_dim)
  scale = 1/convert(T, layer.kernel)

  fill!(input, 0)

  spatial_dim = spatial_dim_T * sizeof(T)
  fea_dim = spatial_dim * channels
  output_fea_dim = spatial_dim * pooled_chann

  for n = 1:num
    input_ptr = convert(Ptr{T}, input.ptr.p) + fea_dim*(n-1)
    output_ptr = convert(Ptr{T}, output.ptr.p) + output_fea_dim*(n-1)

    for pc = 1:pooled_chann
      cstart = (pc-1)*layer.stride - layer.pad[1] + 1
      cend   = min(cstart + layer.kernel - 1, channels)
      cstart = max(1, cstart)
      output_ptr_pc = output_ptr + (pc-1)*spatial_dim

      for c = cstart:cend
        CuBLAS.axpy(backend.cublas_ctx, spatial_dim_T, scale, output_ptr_pc, 1,
            input_ptr + (c-1)*spatial_dim, 1)
      end
    end
  end
end

function cuda_geometry_max_chann_pool(sp_dim::Int, num::Int)
  x_block = round(Int, ceil(convert(Float64, num)/CUDA.THREADS_PER_BLOCK_X));
  y_block = 1;
  z_block = round(Int, ceil(convert(Float64, sp_dim)/CUDA.THREADS_PER_BLOCK_Z));
  return ((x_block,y_block,z_block),
          (CUDA.THREADS_PER_BLOCK_X,1,CUDA.THREADS_PER_BLOCK_Z))

end
function cuda_max_channel_pooling_forward{T}(backend::GPUBackend, input::CuTensorBlob{T},
    output::CuTensorBlob{T}, mask::CuPtr, layer, op_dim)

  sp_dim, channels, num = split_dims(input, op_dim)
  pooled_chann = size(output, op_dim)

  cuda_dim = cuda_geometry_max_chann_pool(sp_dim, num);
  if T == Float32
    kernel = backend.mocha.max_channel_pooling_forward_float
  elseif T == Float64
    kernel = backend.mocha.max_channel_pooling_forward_double
  else
    error("Unsupported data type for channel pooling: $T")
  end

  CUDA.launch(kernel, cuda_dim..., (input.ptr.p, output.ptr.p, mask.p, sp_dim, channels, num,
      pooled_chann, layer.kernel, layer.stride, layer.pad[1]))
end

function cuda_max_channel_pooling_backward{T}(backend::GPUBackend, input::CuTensorBlob{T},
    output::CuTensorBlob{T}, mask::CuPtr, layer, op_dim)

  sp_dim, channels, num = split_dims(input, op_dim)
  pooled_chann = size(output, op_dim)

  cuda_dim = cuda_geometry_max_chann_pool(sp_dim, num);
  if T == Float32
    kernel = backend.mocha.max_channel_pooling_backward_float
  elseif T == Float64
    kernel = backend.mocha.max_channel_pooling_backward_double
  else
    error("Unsupported data type for channel pooling: $T")
  end
  erase!(input)

  CUDA.launch(kernel, cuda_dim..., (input.ptr.p, output.ptr.p, mask.p, sp_dim, channels, num,
      pooled_chann, layer.kernel, layer.stride, layer.pad[1]))
end


text:
# # Heat Equation + Slab Tutorial

#  - load external packages:
import LinearAlgebra, UnPack
import ClimaCore: Fields, Domains, Topologies, Meshes, DataLayouts, Operators, Geometry, Spaces

using Base: show_supertypes
using OrdinaryDiffEq: ODEProblem, solve, SSPRK33

using Logging: global_logger
using TerminalLoggers: TerminalLogger

using RecursiveArrayTools

using OrdinaryDiffEq

using Statistics

#src ## Setup Logging Information
global_logger(TerminalLogger())
const CI = !isnothing(get(ENV, "CI", nothing))

# ## Define Parameters
#  - Global Constants
const FT = Float64;

#  - Experiment-specific Parameters
parameters = (
    ## atmos parameters
    zmin_atm = FT(0.0), # height of atm stack bottom [m]
    zmax_atm = FT(1.0), # height of atm stack top [m]
    n = 15,  # number of elements in atm stack
    μ = FT(0.0001), # diffusion coefficient [m^2 / s]
    T_top = FT(280.0), # fixed temperature at the top of the domain_atm [K]
    ## slab ice parameters
    ρc_ml = FT(4e6), # density times heat transfer coefficient for mixed layer [J / m2 / K ]
    F0_base = FT(120), # ice base transfer coefficient [W / m2 / K]
    T_base = FT(273.16), # ice base temperature [K]
    L_ice = FT(3e8), # latent heat coefficient for ice [J / m3]
    h_ml = FT(1), # mixed layer depth [m]
    T_freeze = FT(273.16), # temperature at freezing point [K]
    k_ice = FT(2), # thermal conductivity of ice [W / m / K]
    ## coupling parameters
    λ = FT(0.001), # coupling coefficient (J / K / s / m2)
)

# ## Define Model Functions

# - Model 1 (atm) Equations
"""
    ∑tendencies_atm!(du, u, (parameters, T_sfc), t)

Heat diffusion equation
    dT/dt =  ∇ μ ∇ T
    where
        (- μ ∇ T) = 0           at z = zmax_atm
        (- μ ∇ T) = F_sfc       at z = zmin_atm

We also use this model to calculate and accumulate the downward surface fluxes, F_sfc:
    F_sfc = - λ * (T_sfc - T1)
    d(F_integrated)/dt  = F_sfc
    where
        F_integrated is reset to 0 at the beginning of each coupling cycle
        T1 = atm temperature near the surface (here assumed equal to the first model level)
"""
function ∑tendencies_atm!(du, u, p, t)

    FT = eltype(u)
    T = u.T # u.x = vector of prognostic variables from DifferentialEquations
    F_sfc = calculate_flux(p.Ya.T_sfc[1], parent(T)[1], p.p)

    ## set BCs

    bcs_bottom = Operators.SetValue(Geometry.WVector(F_sfc ./ p.p.λ))
    bcs_top = Operators.SetValue(Geometry.WVector(FT(0)))

    gradc2f = Operators.GradientC2F()
    gradf2c = Operators.DivergenceF2C(bottom = bcs_bottom, top = bcs_top) # # Neumann BC (face-to-center)
    ## tendency calculations
    @. du.T = gradf2c(p.p.μ * gradc2f(T)) # dT/dt
    du.F .= -F_sfc[1] # d(F_integrated)/dt

end

get_∂F_atm∂T_sfc(p) = -p.λ # TODO: in GCM this will be: 4 σ T_sfc^3 + ||uh|| ρ cp (C_H + L(q_sat1 - q_sat) * dT^(-1))

# - Model 2 (ice) Equations
"""
    solve_ice!(dT_sfc, T_sfc, (parameters, F_accumulated), t)

slab RHS with an implicit solve ice and explicit (forward Euler) solve for ocean

"""
function solve_ice!(integ, Δt)

    Y = integ.u
    Ya = integ.p.Ya
    p = integ.p.p

    FT = eltype(Y)
    ocean_qflux = FT(0)

    # prognostic
    T_ml = Y.T_ml
    T_sfc = Y.T_sfc
    h_ice = Y.h_ice

    # auxiliary
    F_atm = Ya.F_atm

    ∂F_atm∂T_sfc = get_∂F_atm∂T_sfc(p) # this will be passed from atmos/SF.jl

    ΔT_ml = similar(F_atm)
    Δh_ice = similar(F_atm)
    F_ice = similar(F_atm)
    ΔT_sfc = similar(F_atm)

    # ice thickness and mixed layer temperature changes due to atmosphereic and ocean fluxes
    if h_ice[1] > 0 # ice-covered
        @. ΔT_ml = -(p.F0_base * (T_ml - p.T_base) + ocean_qflux) * Δt / (p.h_ml * p.ρc_ml)
        @. Δh_ice = (F_atm - p.F0_base * (T_ml - p.T_base)) * Δt / p.L_ice
    else # ice-free
        @. ΔT_ml = -(F_atm + ocean_qflux) * Δt / (p.h_ml * p.ρc_ml)
        @. Δh_ice = 0
    end

    # adjust if transition to ice-covered
    if (T_ml[1] + ΔT_ml[1] < p.T_freeze)
        @. Δh_ice = Δh_ice - (T_ml + ΔT_ml - p.T_freeze) * (p.h_ml * p.ρc_ml) / p.L_ice
        @. ΔT_ml = p.T_freeze - T_ml
    end

    # adjust if transition to ice-free
    if ((h_ice[1] > 0) & (h_ice[1] + Δh_ice[1] <= 0))
        @. ΔT_ml = ΔT_ml - (h_ice + Δh_ice) * p.L_ice / (p.h_ml * p.ρc_ml)
        @. Δh_ice = -h_ice
    end

    # solve for T_sfc
    if (h_ice[1] + Δh_ice[1] > 0) #  surface is ice-covered
        # if ice covered, solve implicity (for now one Newton iteration: ΔT_s = - F(T_s) / dF(T_s)/dT_s )
        @. F_ice = p.k_ice / (h_ice + Δh_ice) * (p.T_base - T_sfc)
        @. ΔT_sfc = (-F_atm + F_ice) / (p.k_ice / (h_ice + Δh_ice) + ∂F_atm∂T_sfc)
        if (T_sfc[1] + ΔT_sfc[1] > p.T_freeze)
            @. ΔT_sfc = p.T_freeze - T_sfc
        end
        # surface is ice-covered, so update T_sfc as ice surface temperature
        @. Y.T_sfc += ΔT_sfc
    else # ice-free, so update T_sfc as mixed layer temperature
        @. Y.T_sfc = T_ml + ΔT_ml
    end

    # update state
    @. Y.T_ml += ΔT_ml
    @. Y.h_ice += Δh_ice

    @. Ya.ice_mask = h_ice[1] > 0 ? h_ice[1] * FT(1) : h_ice[1] * FT(0)

    return
end

function ∑tendencies_ice_stub(du, u, p, t)
    dY = du
    FT = eltype(dY)

    solve_ice!((; u = u, p = p), p.Δt) # timestepping outside of DeffEq (but DeffEq still used here for saving vars in `integ.sol`)

    @. dY.T_ml = FT(0)
    @. dY.h_ice = FT(0)
    @. dY.T_sfc = FT(0)

end

# - Surface Flux Calculation (coarse bulk formula)
calculate_flux(T_sfc, T1, parameters) = -parameters.λ * (T_sfc - T1);

# - Coupler Communication Functions
# These functions export / import / transform variables
# These functions are now just place holders for coupler transformations (e.g. regridding, masking, etc)
coupler_get(x) = copy(x);
coupler_put(x) = copy(x);

# ## Model Initialization
# - initialize atm model domain and grid
domain_atm = Domains.IntervalDomain(
    Geometry.ZPoint{FT}(parameters.zmin_atm),
    Geometry.ZPoint{FT}(parameters.zmax_atm);
    boundary_tags = (:bottom, :top),
);
mesh_atm = Meshes.IntervalMesh(domain_atm, nelems = parameters.n); # struct, allocates face boundaries to 5,6: atmos
center_space_atm = Spaces.CenterFiniteDifferenceSpace(mesh_atm); # collection of the above, discretises space into FD and provides coords

# - initialize prognostic variables, either as ClimaCore's Field objects or as Arrays
T_atm_0 = Fields.ones(FT, center_space_atm) .* FT(265);
Y_atm = Fields.FieldVector(T = T_atm_0, F = [FT(0)])
Y_ice = Fields.FieldVector(T_sfc = [parameters.T_freeze], h_ice = [FT(0)], T_ml = [parameters.T_freeze])

# - initialize auxiliary variables
ics_aux = (; atm = (; T_sfc = copy(Y_ice.T_sfc)), ice = (; ice_mask = [FT(0)], F_atm = [FT(0)]))

# - specify timestepping information
stepping = (;
    Δt_min = 0.01,
    timerange = (0.0, 10.0),
    Δt_coupler = 0.01,
    odesolver = Euler(), #SSPRK33(),
    nsteps_atm = 1, # number of timesteps of atm per coupling cycle
    nsteps_ice = 1, # number of timesteps of ice per coupling cycle
);

# ## Define the sequential coupling loop
function coupler_solve!(stepping, ics_aux, parameters)
    t = 0.0
    Δt_min = stepping.Δt_min
    Δt_coupler = stepping.Δt_coupler
    t_start = stepping.timerange[1]
    t_end = stepping.timerange[2]

    ## SETUP ATMOS
    ## put all prognostic variable arrays into a vector and ensure that solve can partition them
    prob_atm = ODEProblem(∑tendencies_atm!, Y_atm, (t_start, t_end), (; p = parameters, Ya = ics_aux.atm))
    integ_atm = init(prob_atm, stepping.odesolver, dt = Δt_min, saveat = 1 * Δt_min)

    ## SETUP ICE
    prob_ice =
        ODEProblem(∑tendencies_ice_stub, Y_ice, (t_start, t_end), (; p = parameters, Ya = ics_aux.ice, Δt = Δt_coupler))
    integ_ice = init(prob_ice, Euler(), dt = Δt_coupler, saveat = 1 * Δt_coupler)

    ## coupler stepping
    for t in (t_start:Δt_coupler:t_end)

        ## STEP ATMOS
        ## pre_atmos
        integ_atm.u.F .= [0.0] # surface flux to be accumulated
        integ_atm.p.Ya.T_sfc .= coupler_get(integ_ice.u.T_sfc) # integ_atm.p is the parameter vector of an ODEProblem from DifferentialEquations

        ## run atmos
        ## NOTE: use (t - integ_atm.t) here instead of Δt_coupler to avoid accumulating roundoff error in our timestepping.
        OrdinaryDiffEq.step!(integ_atm, Δt_coupler, true)

        # ## no post_atmos

        ## STEP ICE
        ## pre_ice

        integ_ice.p.Ya.F_atm .= coupler_get(integ_atm.u.F) / Δt_coupler
        Δt_coupler_ice = Δt_coupler

        ## run ice
        OrdinaryDiffEq.step!(integ_ice, Δt_coupler_ice, true)

        # ## no post ice

    end

    return integ_atm, integ_ice
end;

# ## Run the Coupler Model Simulation
integ_atm, integ_ice = coupler_solve!(stepping, ics_aux, parameters);
sol_atm, sol_ice = integ_atm.sol, integ_ice.sol;

# ## Postprocessing and Visualization

# Each integrator output (`sol_atm`, `sol_ice`), contains the DifferentialEquations variable `.u` (the name is hard coded).

ENV["GKSwstype"] = "nul"
import Plots
Plots.GRBackend()

show_plots = isdefined(Main, :SHOWPLOTS) ? SHOWPLOTS : true

path = joinpath(dirname(@__FILE__), "images/")
mkpath(path);

# - Vertical profile at start and end
t0_ = parent(sol_atm.u[1].T)[:, 1];
tend_ = parent(sol_atm.u[end].T)[:, 1];
z_centers = parent(Fields.coordinate_field(center_space_atm))[:, 1];
show_plots ?
Plots.png(
    Plots.plot(
        [t0_ tend_],
        z_centers,
        title = "model 1: atm",
        labels = ["t=0" "t=end"],
        xlabel = "T (K)",
        ylabel = "z (m)",
    ),
    joinpath(path, "tc1_f1.png"),
) : nothing
# ![](images/tc1_f1.png)

# - Conservation: absolute "energy" of both models with time
# convert to the same units (analogous to energy conservation, assuming that is both domains density=1 and thermal capacity=1)
ice_sfc_h_t = [sum(parent(u.h_ice)[:]) for u in sol_ice.u] .* parameters.h_ml * parameters.L_ice;
ice_sfc_u_t = [sum(parent(u.T_sfc)[:]) for u in sol_ice.u] .* parameters.h_ml * parameters.ρc_ml .- ice_sfc_h_t;
atm_sum_u_t =
    [sum(parent(u.T)[:]) for u in sol_atm.u] .* (parameters.zmax_atm .- parameters.zmin_atm) ./ parameters.n .* parameters.λ;
v1 = ice_sfc_u_t .- ice_sfc_u_t[1];
v2 = atm_sum_u_t .- atm_sum_u_t[1];
show_plots ?
Plots.png(
    Plots.plot(
        sol_ice.t,
        [v1 v2 v1 + v2],
        labels = ["ice" "atm" "tot"],
        xlabel = "time (s)",
        ylabel = "pseudo-energy (J / m2)",
    ),
    joinpath(path, "tc1_f2.png"),
) : nothing
# ![](images/tc1_f2.png)

# - Conservation: relative error with time
total = atm_sum_u_t + ice_sfc_u_t;
rel_error = (total .- total[1]) / mean(total);
show_plots ?
Plots.png(
    Plots.plot(sol_ice.t, rel_error, labels = ["tot"], xlabel = "time (s)", ylabel = "relative error"),
    joinpath(path, "tc1_f3.png"),
) : nothing
# ![](images/tc1_f3.png)

#src # - Animation
#src anim = Plots.@animate for u in sol_atm.u
#src     Plots.plot(u.x[1], xlim=(220,280))
#src end
#src Plots.mp4(anim, joinpath(path, "heat.mp4"), fps = 10)
#src
#src function linkfig(figpath, alt = "")
#src     # buildkite-agent upload figpath
#src     # link figure in logs if we are running on CI
#src     if get(ENV, "BUILDKITE", "") == "true"
#src         artifact_url = "artifact://$figpath"
#src         print("\033]1338;url='$(artifact_url)';alt='$(alt)'\a\n")
#src     end
#src end
#src
#src dirname = "heat"
#src linkfig("output/$(dirname)/heat_end.png", "Heat End Simulation")

max_stars_repo_name: CliMA/ClimaCoupler


text:
module EMLsum
	using DanaTypes
	using DotPlusInheritance
	using Reexport
	@reexport using ...types.EMLtypes
	import EMLtypes.length
	include("sum/Sum.jl")
end


text:
export getHelmholtzOperator, getHelmholtzABL

"""
	H = getHelmholtzOperator(m,gamma,omega,Mesh)

	build Helmholtz operator
		H = Laplace - omega^2 * MassMatrix(m,gamma)

	Input:
		m                  - model
		gamma              - attenuation
		omega::Number      - frequency
		Mesh::RegularMesh  - mesh

	Output
		H::SparseMatrixCSC
"""
function getHelmholtzOperator(m,gamma,omega::Number,Mesh::RegularMesh)

Lap   = getNodalLaplacianMatrix(Mesh)
M     = getMassMatrix(vec(m),vec(gamma),Mesh)
# Get the Helmholtz operator (note the sign)
H = Lap .- omega^2 * M
return H
end

function getHelmholtzABL(Msh::RegularMesh,NeumannAtFirstDim::Bool,ABLpad::Array{Int64},ABLamp::Float64)
  h = Msh.h;
  n = Msh.n;
  pad = ABLpad;
  ntup = tuple(n...);

  if Msh.dim==2
	gamma = zeros(ntup);
	b_bwd1 = ((pad[1]:-1:1).^2)./pad[1]^2;
	b_bwd2 = ((pad[2]:-1:1).^2)./pad[2]^2;

	b_fwd1 = ((1:pad[1]).^2)./pad[1]^2;
	b_fwd2 = ((1:pad[2]).^2)./pad[2]^2;
	I1 = (n[1] - pad[1] + 1):n[1];
	I2 = (n[2] - pad[2] + 1):n[2];

	if NeumannAtFirstDim==false
		gamma[:,1:pad[2]] += ones(n[1],1)*b_bwd2';
		gamma[1:pad[1],1:pad[2]] -= b_bwd1*b_bwd2';
		gamma[I1,1:pad[2]] -= b_fwd1*b_bwd2';
	end

	gamma[:,I2] +=  ones(n[1],1)*b_fwd2';
	gamma[1:pad[1],:] += b_bwd1*ones(1,n[2]);
	gamma[I1,:] += b_fwd1*ones(1,n[2]);
	gamma[1:pad[1],I2] -= b_bwd1*b_fwd2';
	gamma[I1,I2] -= b_fwd1*b_fwd2';
	gamma *= ABLamp;
  else
	x1 = linspace(-1,1,n[1]);
	x2 = linspace(-1,1,n[2]);
	x3 = linspace( 0,1,n[3]);
	X1,X2,X3 = ndgrid(x1,x2,x3);
	padx1 = ABLpad[1];
	padx2 = ABLpad[2];
	padx3 = ABLpad[3];
	gammaL = (X1 - x1[padx1]).^2;
	gammaL[padx1+1:end,:,:] = 0.0
	gammaR = (X1 - x1[end-padx1+1]).^2
	gammaR[1:end-padx1,:,:] = 0.0

	gammat = gammaL + gammaR;
	gammat = gammat/maximum(gammat);
	gamma = copy(gammat);
	gammat[:] = 0.0;

	gammaL = (X2 - x2[padx2]).^2;
	gammaL[:,padx2+1:end,:] = 0.0
	gammaR = (X2 - x2[end-padx2+1]).^2
	gammaR[:,1:end-padx2,:] = 0.0

	gammat = gammaL + gammaR
	gammat = gammat/maximum(gammat);
	gamma += gammat;

	gammat = (X3 - x3[end-padx3+1]).^2
	gammat[:,:,1:end-padx3] = 0.0
	gammat = gammat/maximum(gammat);
	gamma += gammat;
	gamma *= ABLamp;
	gamma[gamma.>=ABLamp] = ABLamp;
  end
  return gamma;
end


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    Neighborhood

A neighborhood is a geometry that is not attached to any specific
point in space, and is free to slide over a domain of interest.
"""
abstract type Neighborhood end

# ----------------
# IMPLEMENTATIONS
# ----------------

include("neighborhoods/metricball.jl")


text:
"""
# Description

Walk down an expression recursively while rewriting interpolation expressions
into the interpolated symbol.

# Arguments

1. `e::Any`: An expression.

# Return Values

1. `_::Any`: A new expression in which all splicing has been removed.

# Examples

```
julia> unescape_locals(Expr(:call, :+, :a, Expr(:\$, :b)))
:(a + $(Expr(:escape, :b)))
```
"""
function unescape_locals(@nospecialize(e::Any))
    if isa(e, Expr) && e.head == :$
        # TODO: Reject if e.args[1] is not a symbol?
        # If doing this, do it in validate.
        e.args[1]
    elseif isa(e, Expr)
        Expr(e.head, unescape_locals.(e.args)...)
    else
        e
    end
end


text:
"""
AXEstimator
An abstract type for estimators.
"""
abstract type AXEstimator end

# Methods
"""
coef(fit::AXEstimator)
A method to get the coefficient from an AXEstimator object.
"""
function coef(fit::AXEstimator)
    return fit.β
end

"""
predict(fit::AXEstimator)
A method to predict data based on an AXEstimator object.
"""
function predict(fit::AXEstimator, newdata = nothing)
    isnothing(newdata) ? prediction = fit.X * fit.β : prediction = newdata * fit.β
    return(prediction)
end


text:
# This package provides colours from the popular nord project, which can be found at:
#     <https://github.com/arcticicestudio/nord>

module Nord

using Colors, ColorTypes, FixedPointNumbers

# Dictionaries of colours
export NordColours, PolarNight, SnowStorm, Frost, Aurora

const nord0 = colorant"#2E3440"
const nord1 = colorant"#3B4252"
const nord2 = colorant"#434C5E"
const nord3 = colorant"#4C566A"
const nord4 = colorant"#D8DEE9"
const nord5 = colorant"#E5E9F0"
const nord6 = colorant"#ECEFF4"
const nord7 = colorant"#8FBCBB"
const nord8 = colorant"#88C0D0"
const nord9 = colorant"#81A1C1"
const nord10 = colorant"#5E81AC"
const nord11 = colorant"#BF616A"
const nord12 = colorant"#D08770"
const nord13 = colorant"#EBCB8B"
const nord14 = colorant"#A3BE8C"
const nord15 = colorant"#B48EAD"

const black = nord0
const grey = nord3
const white = nord6
const blue = nord9
const red = nord11
const orange = nord12
const yellow = nord13
const green = nord14

const info = white
const debug = grey
const pass = green
const warn = yellow
const fail = red
const error = red

const NordColours = Dict("nord0" => nord0,
                         "nord1" => nord1,
                         "nord2" => nord2,
                         "nord3" => nord3,
                         "nord4" => nord4,
                         "nord5" => nord5,
                         "nord6" => nord6,
                         "nord7" => nord7,
                         "nord8" => nord8,
                         "nord9" => nord9,
                         "nord10" => nord10,
                         "nord11" => nord11,
                         "nord12" => nord12,
                         "nord13" => nord13,
                         "nord14" => nord14,
                         "nord15" => nord15,
                         )

const PolarNight = Dict("nord0" => nord0,
                        "nord1" => nord1,
                        "nord2" => nord2,
                        "nord3" => nord3
                        )

const SnowStorm = Dict("nord4" => nord4,
                       "nord5" => nord5,
                       "nord6" => nord6
                       )

const Frost = Dict("nord7" => nord7,
                   "nord8" => nord8,
                   "nord9" => nord9,
                   "nord10" => nord10,
                  )

const Aurora = Dict("nord11" => nord11,
                    "nord12" => nord12,
                    "nord13" => nord13,
                    "nord14" => nord14,
                    "nord15" => nord15,
                   )

end # module


text:
import CUDArt
using Base.Test

##########################
# Low-level memory tests #
##########################
devlist = CUDArt.devices(dev->true)
for dev in devlist
    CUDArt.device(dev)
    p = CUDArt.malloc(Uint8, 1)
    p2 = CUDArt.malloc(Uint16, 100)
    CUDArt.free(p)
    CUDArt.free(p2)
    a = rand(5,3)
    g = CUDArt.CudaArray(a)
    gp = CUDArt.CudaPitchedArray(a)
    CUDArt.free(g)
    CUDArt.free(gp)
    # Also test finalizer calls
    g = CUDArt.CudaArray(a)
    g = CUDArt.CudaPitchedArray(a)
    CUDArt.device_reset(dev)
end

gc()  # check for finalizer errors

# Test free when the allocated device is not the current device
if length(devlist) > 1
    CUDArt.device(devlist[1])
    p1a = CUDArt.malloc(Uint8, 1)
    p1 = CUDArt.malloc(Uint16, 100)
    CUDArt.device(devlist[2])
    p2 = CUDArt.malloc(Uint16, 100)
    CUDArt.free(p1)
    CUDArt.free(p2)
    CUDArt.free(p1a)
    CUDArt.device_reset(devlist[1])
    CUDArt.device_reset(devlist[2])
end

gc()

# Check whether gc/finalizers lead to bad behavior across
# device resets (the pointer address is reproducible)
CUDArt.device(0)
a = rand(5,3)
g = CUDArt.CudaArray(a)
CUDArt.device_reset(0)
@test isempty(CUDArt.cuda_ptrs)
CUDArt.device(0)
g = CUDArt.CudaArray(a)
# In case the next gc() test yields an error, the next two lines let
# us do some archaeology
dictcopy = deepcopy(CUDArt.cuda_ptrs)
gptrcopy = copy(pointer(g))
gc()   # Check that this doesn't delete the new g
@test !isempty(CUDArt.cuda_ptrs)
h_g = CUDArt.to_host(g)
CUDArt.free(g)
@test isempty(CUDArt.cuda_ptrs)
@test h_g == a
CUDArt.device_reset(0)

gc()


text:
using Test
using DRAGO
using Random

Random.seed!(1776)

@testset "Test" begin
    @test true
end


text:
##
using Random
using Plots
using LinearAlgebra
theme(:dark)
##
Random.seed!(7)
X = randn(2,100)
A = randn(2,2)
Y = A*X
scatter(Y[1,:], Y[2,:], label="", size=(300,300),color=3)
xlabel!("x")
ylabel!("y")
##
savefig("pca-2d-fig.pdf")

##
using Interact
using Statistics
using Printf
settheme!(:nativehtml)
mp = @manipulate for θ in 0.0:2π/100:2π
    v = [cos(θ),sin(θ)]
    z = vec(v'*Y)
    p = plot(
      scatter(Y[1,:], Y[2,:], label="",color=3, aspect_ratio=:equal),
      histogram(z,label=""),
      size=(900,450),tickfont=16,titlefont=20
      )
    plot!(p[1],2*[0,v[1]],2*[0,v[2]],label="",linewidth=2)
    scatter!(p[1], v[1].*z, v[2].*z, color=2, label="",
      markerstrokewidth=0,alpha=0.3, markersize=4)
    xlims!(p[2],(-5,5))
    ylims!(p[2],(0,30))
    xlims!(p[1],(-4,4))
    ylims!(p[1],(-4.5,4.5))
    title!(p[1],@sprintf("Angle = %.2f", θ))
    title!(p[2],@sprintf("Variance = %.2f", var(z)))
end
## If the figure doesn't show, run
using Interact
ui = button()
display(ui)
##
θs = 0.0:π/1000:π
plot(collect(θs), map(x -> var(vec([cos(x),sin(x)]'*Y)), θs),
  label="", size=(500,250))
xlabel!("Angle")
ylabel!("Variance")
##
savefig("pca-2d-angle-variance.pdf")
## PCA is a way to get this angle automatically!
function mypca(X::AbstractMatrix, donormalize::Bool = false)
  means = vec(mean(X, dims=2)) # take the average of each coordinate
  C = X .- means # center each coordinate
  if donormalize
    for i=1:size(C,1)
      normalize!(@view C[i,:])
    end
  end
  U,S,V = svd(C)
end
U,S,V = mypca(Y, false)
U[:,1]
##
θpca = atan(-U[2,1],-U[1,1])
##
plot(collect(θs), map(x -> var(vec([cos(x),sin(x)]'*Y)), θs),
label="", size=(500,250))
xlabel!("Angle")
ylabel!("Variance")
annotate!(θpca, var(vec([cos(θpca),sin(θpca)]'*Y)), text("*", 24, :red))
##
savefig("pca-2d-angle-variance-opt.pdf")
##
scatter(Y[1,:], Y[2,:], label="", size=(300,300),color=3)
plot!(2*[0,-U[1,1]],2*[0,-U[2,1]],label="",linewidth=2, color=5)
##
savefig("pca-2d-fig-PCA.pdf")
##
S
## Save a set of points
pyplot()
theme(:default)
scatter(Y[1,:], Y[2,:], label="", size=(300,300),color=3,
  framestyle=:none, background=:transparent)
savefig("pca-2d-data.pdf")


text:
# Write precompiles for userimg.jl
function write(io::IO, pc::Vector)
    for ln in pc
        println(io, ln)
    end
end

function write(filename::AbstractString, pc::Vector)
    path, fn = splitdir(filename)
    if !isdir(path)
        mkpath(path)
    end
    open(filename, "w") do io
        write(io, pc)
    end
    nothing
end

"""
    write(prefix::AbstractString, pc::Dict; always::Bool = false)

Write each modules' precompiles to a separate file.  If `always` is
true, the generated function will always run the precompile statements
when called, otherwise the statements will only be called during
package precompilation.
"""
function write(prefix::AbstractString, pc::Dict; always::Bool = false)
    if !isdir(prefix)
        mkpath(prefix)
    end
    for (k, v) in pc
        open(joinpath(prefix, "precompile_$k.jl"), "w") do io
            if any(str->occursin("__lookup", str), v)
                println(io, lookup_kwbody_str)
            end
            println(io, "function _precompile_()")
            !always && println(io, "    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing")
            for ln in sort(v)
                println(io, "    ", ln)
            end
            println(io, "end")
        end
    end
end


text:
struct PlaneWaveHH3DTD{T,P,F} <: TDFunctional{T}
    direction::P
    speed_of_light::T
    signature::F
    amplitude::T
end

function planewave(direction, speedoflight::Number, signature, amplitude=one(speedoflight)::Number)
    PlaneWaveHH3DTD(direction, speedoflight, signature, amplitude)
end

*(a, f::PlaneWaveHH3DTD) = PlaneWaveHH3DTD(f.direction, f.speed_of_light, f.signature, a*f.amplitude)

function(f::PlaneWaveHH3DTD)(r,t)
    k = f.direction
    u = sum(k[i]*r[i] for i in 1:length(k))
    #u = dot(f.direction, r)
    a = f.amplitude
    c = f.speed_of_light
    h = f.signature
    a * h(c*cartesian(t)[1] - u)
end

function gradient(f::PlaneWaveHH3DTD)
    @assert f.amplitude ≈ 1
    PlaneWaveMWTD(
        f.direction,
        -f.direction,
        f.speed_of_light,
        derive(f.signature)
    )
end

struct DotTraceHH{T,F} <: TDFunctional{T}
    field::F
end

dot(::NormalVector, field::TDFunctional) = DotTraceHH{scalartype(field), typeof(field)}(field)

function (f::DotTraceHH)(p,t)
    n = normal(p)
    # x = cartesian(p)
    return dot(n, f.field(p,t))
end


text:
module VortexHelperPurpleBooster
using ..Ahorn, Maple

@mapdef Entity "VortexHelper/PurpleBooster" PurpleBooster(x::Integer, y::Integer, lavender::Bool = false)

const placements = Ahorn.PlacementDict(
    "Booster (Purple) (Vortex Helper)" => Ahorn.EntityPlacement(
        PurpleBooster,
        "point"
    ),
    "Booster (Lavender) (Vortex Helper)" => Ahorn.EntityPlacement(
        PurpleBooster,
        "point",
        Dict{String, Any}(
            "lavender" => true
        )
    )
)

boosterSprite = "objects/VortexHelper/slingBooster/slingBooster00"
lavenderBoosterSprite = "objects/VortexHelper/lavenderBooster/boosterLavender00"

function Ahorn.selection(entity::PurpleBooster)
    x, y = Ahorn.position(entity)

    return Ahorn.Rectangle(x - 9, y - 9, 18, 18)
end

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::PurpleBooster, room::Maple.Room)
    Ahorn.drawSprite(ctx, Bool(get(entity.data, "lavender", false)) ? lavenderBoosterSprite : boosterSprite, 0, 0)
end

end

max_stars_repo_name: catapillie/VortexHelper

max_issues_repo_name: catapillie/VortexHelper

max_forks_repo_name: catapillie/VortexHelper


text:
module Experiments

include("MiniOmega.jl")

end # module

max_stars_repo_name: zenna/expect

max_issues_repo_name: zenna/expect

max_forks_repo_name: zenna/expect


text:
using PowerFlowData
using Documenter

DocMeta.setdocmeta!(PowerFlowData, :DocTestSetup, :(using PowerFlowData); recursive=true)

makedocs(;
    modules=[PowerFlowData],
    authors="Nick Robinson <npr251@gmail.com> and contributors",
    repo="https://github.com/nickrobinson251/PowerFlowData.jl/blob/{commit}{path}#{line}",
    sitename="PowerFlowData.jl",
    format=Documenter.HTML(canonical="https://nickrobinson251.github.io/PowerFlowData.jl"),
    pages=["Home" => "index.md" ],
    strict=true,
    checkdocs=:exports,
)

deploydocs(;
    repo="github.com/nickrobinson251/PowerFlowData.jl",
    devbranch="main",
    push_preview=true,
)


text:
"""
4. Initialization scheme from Mei et al
"""
struct IterativeInitializationMeiEtAl{T <: Real} <: InitializationStrategy
    τ::T
    ϵ::T
end
IterativeInitializationMeiEtAl() = IterativeInitialization()
IterativeInitialization() = IterativeInitializationMeiEtAl(1.0, 1E-7)

function initialize(
    strategy::IterativeInitializationMeiEtAl,
    q::Quadrature,
    problem::FluidFlowProblem,
    cm = SRT;
    process_method = ProcessingMethod(problem, true, 1000),
)
    # Initialize f such that it has a constant density for all nodes, We don't
    # care about the velocity here since this will be set in the collision method
    f = [
        q.weights[f_idx]
        for x_idx in 1:(problem.NX), y_idx in 1:(problem.NY), f_idx in 1:length(q.weights)
    ]

    # Simultae the lattice botlzmann method with a special collision operator that uses
    # u_0(x) as the velocity of each lattice node
    # The method stops if the density has converged
    model = LatticeBoltzmannModel(
        copy(f),
        copy(f),
        q,
        IterativeInitializationCollisionModel(q, strategy.τ, problem),
        boundary_conditions(problem),
        ProcessIterativeInitialization(strategy.ϵ, problem, process_method),
    )

    simulate(model, 1:10000)

    return model.f_stream
end


text:
using ITensors
using Test

@testset "Empty ITensor storage operations" begin
  i, j, k = Index.(2, ("i", "j", "k"))

  A = ITensor(i, j)
  B = ITensor(j, k)
  C = A * B
  @test hassameinds(C, (i, k))
  @test storage(C) isa
        ITensors.EmptyStorage{ITensors.EmptyNumber,<:ITensors.Dense{ITensors.EmptyNumber}}

  A = ITensor(Float64, i, j)
  B = ITensor(j, k)
  C = A * B
  @test hassameinds(C, (i, k))
  @test storage(C) isa ITensors.EmptyStorage{Float64,<:ITensors.Dense{Float64}}

  A = ITensor(i, j)
  B = ITensor(ComplexF64, j, k)
  C = A * B
  @test hassameinds(C, (i, k))
  @test storage(C) isa ITensors.EmptyStorage{ComplexF64,<:ITensors.Dense{ComplexF64}}

  A = ITensor(Float64, i, j)
  B = ITensor(ComplexF64, j, k)
  C = A * B
  @test hassameinds(C, (i, k))
  @test storage(C) isa ITensors.EmptyStorage{ComplexF64,<:ITensors.Dense{ComplexF64}}
end

@testset "Empty ITensor storage addition" begin
  i, j = Index.((2, 3))

  A = ITensor(i, j)
  B = randomITensor(j, i)

  C = A + B
  @test inds(C) == (i, j)
  @test C ≈ B

  C = B + A
  @test inds(C) == (j, i)
  @test C ≈ B
end

@testset "Empty QN ITensor storage operations" begin
  i = Index([QN(0) => 1, QN(1) => 1])
  A = ITensor(i', dag(i))

  @test storage(A) isa ITensors.EmptyStorage{
    ITensors.EmptyNumber,<:ITensors.BlockSparse{ITensors.EmptyNumber}
  }

  C = A' * A

  @test hassameinds(C, (i'', i))
  @test storage(C) isa ITensors.EmptyStorage{
    ITensors.EmptyNumber,<:ITensors.BlockSparse{ITensors.EmptyNumber}
  }

  B = randomITensor(dag(i), i')

  C = A' * B

  @test hassameinds(C, (i'', i))
  @test storage(C) isa ITensors.EmptyStorage{Float64,<:ITensors.BlockSparse{Float64}}

  C = B' * A

  @test hassameinds(C, (i'', i))
  @test storage(C) isa ITensors.EmptyStorage{Float64,<:ITensors.BlockSparse{Float64}}

  C = B + A
  @test inds(C) == inds(B)
  @test C ≈ B

  @test_broken A + B
end

nothing


text:
using LibGit2

function mergefiles(mainfile, depfiles)
    pkg_folder = normpath(@__DIR__, "../..")
    src_folder = joinpath(pkg_folder, "src")

    repo = LibGit2.GitRepo(pkg_folder)
    hash = LibGit2.GitHash(LibGit2.peel(LibGit2.GitCommit, LibGit2.head(repo))) |> string

    main = """
    #=
    AUTO-GENERATED FILE - DO NOT EDIT

    This file is derived from the following fork of the NormalHermiteSplines.jl package:

        https://github.com/jondeuce/NormalHermiteSplines.jl#$(hash)

    As it is not possible to depend on a package fork, the above module is included here verbatim.

    The `LICENSE.md` file contents from the original repository follows:

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

    MIT License

    Copyright (c) 2021 Igor Kohanovsky

    Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:

    The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.

    THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
    =#

    $(readchomp(joinpath(src_folder, mainfile)))
    """
    for (i, depfile) in enumerate(depfiles)
        maybenewline = i == length(depfiles) ? "" : "\n"
        main = replace(
            main,
            "include(\"./$(depfile)\")" => chomp(
                """
                ####
                #### $(depfile)
                ####

                $(readchomp(joinpath(src_folder, depfile)))
                """
            ) * maybenewline
        )
    end
    return main
end

mainfile = "NormalHermiteSplines.jl"
depfiles = [
    "ReproducingKernels.jl",
    "GramMatrix.jl",
    "Splines.jl",
    "Utils.jl",
    "Interpolate.jl",
]

outfile = isempty(ARGS) ? joinpath(pwd(), mainfile) : ARGS[1]
open(outfile; write = true) do io
    # Make sure we write 64bit integer in little-endian byte order
    write(io, mergefiles(mainfile, depfiles))
end


text:
using SafeTestsets

# Use all tests or just the active one
use_all_tests = false
active_test = "Default OptMap test"

@time begin
    # Run all of the tests
    if use_all_tests
        println("Running all tests ...")

    # Only run the active test
    else
        println("Only running active test ($active_test) ...")
        @time @safetestset "Default OptMap test" begin
            include("test_nloptmap.jl")
        end
    end

end

max_stars_repo_name: nmaterise/optmap

max_issues_repo_name: nmaterise/optmap

max_forks_repo_name: nmaterise/optmap


text:
# for things that are not actually related to the problem itself

"""
    try_get(collection, inds...)

Like `getindex`, except returns `nothing` if not found or
if the collection itsself is `nothing`.
"""
try_get(::Nothing, inds...) = nothing
try_get(x, inds...) = get(x, inds, nothing)


text:
"""
GenX: An Configurable Capacity Expansion Model
Copyright (C) 2021,  Massachusetts Institute of Technology
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.
A complete copy of the GNU General Public License v2 (GPLv2) is available
in LICENSE.txt.  Users uncompressing this from an archive may not have
received this license file.  If not, see <http://www.gnu.org/licenses/>.
"""

@doc raw"""
	discharge(EP::Model, inputs::Dict, EnergyShareRequirement::Int)
This module defines the power decision variable $\Theta_{y,t} \forall y \in \mathcal{G}, t \in \mathcal{T}$, representing energy injected into the grid by resource $y$ by at time period $t$.
This module additionally defines contributions to the objective function from variable costs of generation (variable O&M plus fuel cost) from all resources $y \in \mathcal{G}$ over all time periods $t \in \mathcal{T}$:
```math
\begin{aligned}
	Obj_{Var\_gen} =
	\sum_{y \in \mathcal{G} } \sum_{t \in \mathcal{T}}\omega_{t}\times(\pi^{VOM}_{y} + \pi^{FUEL}_{y})\times \Theta_{y,t}
\end{aligned}
```
"""
function discharge(EP::Model, inputs::Dict, EnergyShareRequirement::Int)

	println("Discharge Module")

	dfGen = inputs["dfGen"]

	G = inputs["G"]     # Number of resources (generators, storage, DR, and DERs)
	T = inputs["T"]     # Number of time steps
	Z = inputs["Z"]     # Number of zones
	### Variables ###

	# Energy injected into the grid by resource "y" at hour "t"
	@variable(EP, vP[y=1:G,t=1:T] >=0);

	### Expressions ###

	## Objective Function Expressions ##

	# Variable costs of "generation" for resource "y" during hour "t" = variable O&M plus fuel cost
	@expression(EP, eCVar_out[y=1:G,t=1:T], (inputs["omega"][t]*(dfGen[!,:Var_OM_Cost_per_MWh][y]+inputs["C_Fuel_per_MWh"][y,t])*vP[y,t]))
	#@expression(EP, eCVar_out[y=1:G,t=1:T], (round(inputs["omega"][t]*(dfGen[!,:Var_OM_Cost_per_MWh][y]+inputs["C_Fuel_per_MWh"][y,t]), digits=RD)*vP[y,t]))
	# Sum individual resource contributions to variable discharging costs to get total variable discharging costs
	@expression(EP, eTotalCVarOutT[t=1:T], sum(eCVar_out[y,t] for y in 1:G))
	@expression(EP, eTotalCVarOut, sum(eTotalCVarOutT[t] for t in 1:T))

	# Add total variable discharging cost contribution to the objective function
	EP[:eObj] += eTotalCVarOut

	# ESR Policy
	if EnergyShareRequirement >= 1

		@expression(EP, eESRDischarge[ESR=1:inputs["nESR"]], sum(inputs["omega"][t]*dfGen[!,Symbol("ESR_$ESR")][y]*EP[:vP][y,t] for y=dfGen[findall(x->x>0,dfGen[!,Symbol("ESR_$ESR")]),:R_ID], t=1:T)
						- sum(inputs["dfESR"][:,ESR][z]*inputs["omega"][t]*inputs["pD"][t,z] for t=1:T, z=findall(x->x>0,inputs["dfESR"][:,ESR])))

		EP[:eESR] += eESRDischarge
	end

	return EP

end

max_stars_repo_name: Betristor/LTESOM

max_issues_repo_name: Betristor/LTESOM

max_forks_repo_name: Betristor/LTESOM


text:
"""
---
FBMC waveform structure
# --- Syntax
- nFFT	  : FFT size
- K		  : Overlapping factor
- allocatedSubcarriers	: Vector of allocated subcarriers [Array{Int}]
# ---
# v 1.0 - Robin Gerzaguet.
"""
struct StrucFBMC<: Waveform
	nFFT::Int;
	K::Int;
	allocatedSubcarriers::Array{Int};
end

"""
---
Initiate FBMC structure
# --- Syntax
  fbmc = initFBMC(nFFT,K,allocatedSubcarriers)
# --- Input parameters
- nFFT	  : FFT size
- K		  : Overlapping factor
- allocatedSubcarriers	: Vector of allocated subcarriers [Array{Int}]
# --- Output parameters
- fbmc	  : FBMC structure [StrucFBMC]
# ---
# v 1.0 - Robin Gerzaguet.
"""
function initFBMC(nFFT,K,allocatedSubcarriers)
	# ---Checking FFT size
	if maximum(allocatedSubcarriers) > nFFT
		error("Subcarrier allocation is impossible");
	end
	# --- Create the FBMC structure
	return StrucFBMC(nFFT,K,allocatedSubcarriers)
end

"""
---
Apply OQAM pre-processing to incoming matrix qamMat of size nbSubcarriers x nbSymb
# --- Syntax
oqamMat = oqamConversion(qamMat)
# --- Input parameters
- qamMat : Input complex qam Matrix [Array{Complex{Float64,nbSubcarriers,nbSymb}}]
# --- Output parameters
- oqamMat : OQAM matrix (pure real and pure imag. alterns) [Array{Complex{Float64,nbSubcarriers,nbSymb*2}}]
# ---
# v 1.0 - Robin Gerzaguet.
"""
function oqamMapping(qamMat)
	# --- Size of incoming matrix :w
	(nbSubcarriers,nbSymb)	  = size(qamMat);
	# Real and imag split:
	oqamMat	  = zeros(Complex{Float64},nbSubcarriers,2*nbSymb);
	# --- Init matrix
	for iS = 1 : 1 : nbSubcarriers
		# Mask to keep real or imag part
		pattern		= (-1im)^( mod( mod(iS,2) + 1,2));
		pattern2	= (-1im)^( mod( mod(iS,2) + 2,2));
		# 1 // 1i shift
		scaleReIm	= (1im)^( mod( mod(iS,2) + 1,2));
		scaleReIm2	= (1im)^( mod( mod(iS,2) + 2,2));
		for iN = 1 : 1 : nbSymb
			# OQAM processing
			oqamMat[iS,2*(iN-1)+1]	= scaleReIm*real(pattern * qamMat[iS,iN]);
			oqamMat[iS,2*(iN-1)+2]	= scaleReIm2*real(pattern2* qamMat[iS,iN]);
		end
	end
	return oqamMat;
end

"""
---
Returns the PHYDIAS time domain impulse response of desired FBMC filter parametrized by its overlapping factor and FFT size
# --- Syntax
p = getFBMCFilter(K,nFFT,type);
# --- Input parameters
- K	  : Overlapping factor
- nFFT: FFT size
# --- Output parameters
- p	  : FBMC time impulse response
# ---
# v 1.0 - Robin Gerzaguet.
"""
function getFBMCFilter(K,nFFT)
	if K == 2
		p1			= sqrt(2)/2;
		P			= p1;
	elseif K== 3
		p1			= 0.91143783;
		P			= [p1;sqrt(1-p1^2)];
	elseif K== 4
		p1			= 0.97195983;
		p2			= 1/sqrt(2);
		p3          = sqrt(1-p1^2);
		P			= [p1;p2;p3];
	elseif K== 5
		# Parameters values are arbitrarily chosen
		p0          = 1;
		P1          = 0.998;
		P2          = 0.844;
		P           = [P1;P2;sqrt(1-P2^2);sqrt(1-P1^2)];
	elseif K== 6
		p0          = 1;
		p1          = 0.99722723;
		p2          = 0.94136732;
		P           = [ p1;p2;sqrt(2)/2;sqrt(1-p2^2);sqrt(1-p1^2)];
	elseif K== 8
		p0          = 1;
		p1          = 0.99988389;
		p2          = 0.99315513;
		p3          = 0.92708081;
		P          = [ p1;p2;p3;sqrt(2)/2;sqrt(1-p3^2);sqrt(1-p2^2);sqrt(1-p1^2)];
	else
		error("Overlapping factor should be 2, 3, 4, 6 or 8");
	end
	mSize = K*nFFT;
	t	  = collect(1:mSize);
	p	  = zeros(Float64,mSize);
	kV	  = 1:K-1;
	for n = t
		p[n] = 1 + 2 * sum( (-1).^(kV) .* P .* cos.( 2*pi*kV ./ (K*nFFT) * (n-1) )   );
	end
	return p./K;
end

"""
---
Generate a FBMC-OQAM signal in time domain, based on input complex matrix (before OQAM processing) and FBMC parameters.

Transmitter is based on Polyphase Network implementation, with PHYDIAS filter (parametrized by overlapping factor).
# --- Syntax
fbmcSigGen(qamMat,nFFT,K,allocatedSubcarriers)
# --- Input parameters
- qamMat	: Complex QAM Time-Frequency matrix [Array{Complex{Float64}}]
- nFFT		: FFT size  [Int]
- K			: Overlapping factor [Int]
- allocatedSubcarriers	: Vector of allocated subcarriers [Array{Int}]
# --- Output parameters
- sigId		: Complex baseband signal in time domain [Array{Complex{Float64}},nbEch]; nbEch = (2*nbSymb-1)*nFFT*/2+nFFT*K)
# ---
# v 1.0 - Robin Gerzaguet.
"""
function fbmcSigGen(qamMat,nFFT,K,allocatedSubcarriers)
	# Set 4 core for FFT computation
	#FFTW.set_num_threads(4)
	# --- Getting parameters
	nbSymb		= size(qamMat,2); 		# --- Applied on x symbols
	# --- Mapping to nFFT elements
	oqamMat		= oqamMapping(qamMat);
	oqamCurr	= zeros(Complex{Float64},nFFT,2*nbSymb);
	oqamCurr[allocatedSubcarriers,:] = oqamMat;
	# --- Switch to time domain
	sTmp 	= ifft(oqamCurr,1);
	# --- Init FBMC filter
	p = getFBMCFilter(K,nFFT);
	δ = 1/2;					# Space between 2 symbols
	y = zeros(Complex{Float64},Int((2*nbSymb-1)*nFFT*δ +nFFT*K));
	for iN = 1 : 1 : 2*nbSymb
		# Repeat IFFT and Phydias apodisation
		sigApod	  = p .* repeat(sTmp[:,iN],K);
		# Overlap and add
		y[ (iN-1)*Int(nFFT*δ)  .+ collect(1:( nFFT*K)) ] = y[ (iN-1)*Int(nFFT*δ)  .+ collect(1:( nFFT*K)) ] .+ sigApod;
	end
	return y
end

# --- MD is waveform structure is given
function fbmcSigGen(qamMat,fbmc::StrucFBMC)
	return fbmcSigGen(qamMat,fbmc.nFFT,fbmc.K,fbmc.allocatedSubcarriers);
end


text:
function simpleregressiongp(x, y; iterations = iterations, seed = 1, numberofrestarts = 1, kernel = kernel)

    @assert(length(x) == length(y))

    rg = MersenneTwister(seed)

    N = length(x)

    JITTER = 1e-8

    function unpack(param)

        @assert(length(param) == 3)

        exp.(param)

    end

    function objective(param)

        local ρ, α, σ² = unpack(param)

        local K = calculatekernelmatrix(kernel, ρ, α, x) + σ²*I

        -logpdf(MvNormal(zeros(N), K), y)

    end

    safeobj = safewrapper(objective)

    opt = Optim.Options(iterations = iterations, show_trace = false, show_every = 10)

    solve() = optimize(safeobj, randn(rg, 3)*3, NelderMead(), opt)

    solutions = [solve() for _ in 1:numberofrestarts]

    bestindex = argmin([s.minimum for s in solutions])

    ρopt, αopt, σ²opt = unpack(solutions[bestindex].minimizer)

    K = calculatekernelmatrix(kernel, ρopt, αopt, x) + σ²opt*I

    function predicttest(xtest)

        # dimensions: N × Ntest
        kB✴ = calculatekernelmatrix(kernel, ρopt, αopt, x, xtest)

        # Ntest × Ntest
        cB = calculatekernelmatrix(kernel, ρopt, αopt, xtest)

        # full predictive covariance
        Σpred = cB - kB✴' * (K \ kB✴)+ JITTER*I

        makematrixsymmetric!(Σpred)

        # predictive mean

        μpred = kB✴' * (K \ y)

        return μpred, Σpred

    end

    return solutions[bestindex].minimum, predicttest, (ρopt, αopt, σ²opt)

end


text:
# Calculate structure factor
function calc_strfact( Xpos::Array{Float64,2}, Nspecies::Int,
    atm2species::Array{Int}, G::Array{Float64,2} )

    Ng = size(G)[2]
    Na = size(Xpos)[2]
    println("Natoms = $(Na)")
    Sf = zeros(Complex128,Ng,Nspecies)
    for ia = 1:Na
        isp = atm2species[ia]
        for ig = 1:Ng
            GX = Xpos[1,ia]*G[1,ig] +
                 Xpos[2,ia]*G[2,ig] +
                Xpos[3,ia]*G[3,ig]
            Sf[ig,isp] = Sf[ig,isp] + cos(GX) - im*sin(GX)
        end
    end
    return Sf
end


text:
include("../flux_components.jl");

# Multigrid Helmholtz Shifted Laplacian Methods

function get_helmholtz_matrices!(kappa, omega, gamma; alpha=0.5)
    shifted_laplacian_matrix = kappa .* kappa .* omega .* (omega .- (im .* gamma) .- (im .* omega .* alpha))
    helmholtz_matrix = kappa .* kappa .* omega .* (omega .- (im .* gamma))

    return shifted_laplacian_matrix, helmholtz_matrix
end

function jacobi_helmholtz_method!(n, m, h, x, b, matrix; max_iter=1, w=0.8, use_gmres_alpha=0)
    for i in 1:max_iter
        y = helmholtz_chain!(real(reshape(x, n-1, m-1, 1, 1)), matrix; h=h) + im*helmholtz_chain!(imag(reshape(x, n-1, m-1, 1, 1)), matrix; h=h)
        residual = b - y[:,:,1,1]
        d = r_type(4.0 / h^2) .- matrix
        alpha = r_type(w) ./ d
        x = x + alpha .* residual
    end
    return x
end

function jacobi_helmholtz_method_channels!(n, m, h, x, b, matrix, matrixch; max_iter=1, w=0.8, use_gmres_alpha=0)
    for i in 1:max_iter
        y = helmholtz_chain_channels!(x, matrix; h=h)
        residual = b - y
        d = r_type(4.0 / h^2) .- matrixch
        alpha = r_type(w) ./ d
        x = x + alpha .* residual
    end
    return x
end

function v_cycle_helmholtz!(n, m, h, x, b, kappa, omega, gamma; u = 1, v1_iter = 1, v2_iter = 10, use_gmres_alpha = 0, alpha= 0.5, log = 0, level = nothing)

    shifted_laplacian_matrix, helmholtz_matrix = get_helmholtz_matrices!(kappa, omega, gamma; alpha=r_type(alpha))

    # Relax on Ax = b v1_iter times with initial guess x
    x = jacobi_helmholtz_method!(n, m, h, x, b, shifted_laplacian_matrix; max_iter=v1_iter, use_gmres_alpha=use_gmres_alpha)

    if( n % 2 == 0 && n > 4 && m % 2 == 0 && m > 4 && (level == nothing || level > 1))

        # Compute residual on fine grid
        x_matrix = reshape(x, n-1, m-1, 1, 1)
        residual_fine = b - (helmholtz_chain!(real(x_matrix), helmholtz_matrix; h=h) + im*helmholtz_chain!(imag(x_matrix), helmholtz_matrix; h=h))[:,:,1,1] #) helmholtz_chain!(x_matrix, helmholtz_matrix; h=h)[:,:,1,1]

        # Compute residual, kappa and gamma on coarse grid
        residual_coarse = down(reshape(real(residual_fine), n-1, m-1, 1, 1)|>pu)[:,:,1,1] + im * down(reshape(imag(residual_fine), n-1, m-1, 1, 1)|>pu)[:,:,1,1]
        kappa_coarse = down(reshape(kappa, n-1, m-1, 1, 1)|>pu)[:,:,1,1]
        gamma_coarse = down(reshape(gamma, n-1, m-1, 1, 1)|>pu)[:,:,1,1]

        # Recursive operation of the method on the coarse grid
        n_coarse = size(residual_coarse,1)+1
        m_coarse = size(residual_coarse,2)+1
        x_coarse = zeros(c_type,n_coarse-1, m_coarse-1)|>pu

        for i = 1:u
            x_coarse, helmholtz_matrix_coarse = v_cycle_helmholtz!(n_coarse, m_coarse, h*2, x_coarse, residual_coarse, kappa_coarse, omega, gamma_coarse; use_gmres_alpha = use_gmres_alpha,
                                                                    u=u, v1_iter=v1_iter, v2_iter=v2_iter, log=log, level = (level == nothing ? nothing : (level-1)))
        end
        x_coarse_matrix = reshape(x_coarse, n_coarse-1, m_coarse-1, 1, 1)

        # Correct
        fine_error = up(real(x_coarse_matrix)|>pu)[:,:,1,1] + im * up(imag(x_coarse_matrix)|>pu)[:,:,1,1]
        x = x + fine_error

        if log == 1
            r1 = residual_fine
            r2 = b - reshape(helmholtz_chain!(reshape(x, n-1, m-1, 1, 1), helmholtz_matrix; h=h), n-1, m-1)
            println("n = $(n), norm of x = $(norm(x)), norm of fine_error = $(norm(fine_error)), residual before vcycle =$(norm(r1)/norm(b)), residual after vcycle =$(norm(r2)/norm(b)), level =$(level)")
        end
    else
        # Coarsest grid
        A_Coarsest(v::a_type) = vec(helmholtz_chain!(reshape(real(v), n-1, m-1, 1, 1), shifted_laplacian_matrix; h=h) + im*helmholtz_chain!(reshape(imag(v), n-1, m-1, 1, 1), shifted_laplacian_matrix; h=h)) # vec(helmholtz_chain!(reshape(v, n-1, n-1, 1, 1), shifted_laplacian_matrix; h=h))
        M_Coarsest(v::a_type) = M_Jacobi(n, m, h, x, shifted_laplacian_matrix, 1, v; use_gmres_alpha=use_gmres_alpha)
        x,flag,err,iter,resvec = fgmres_func(A_Coarsest, vec(b), v2_iter, tol=1e-15, maxIter=1,
                                                    M=M_Coarsest, x=vec(x), out=-1, flexible=true)
        x = reshape(x, n-1, m-1)
    end

    # Relax on Ax = b v1_iter times with initial guess x
    x = jacobi_helmholtz_method!(n, m, h, x, b, shifted_laplacian_matrix; max_iter=v1_iter, use_gmres_alpha=use_gmres_alpha)

    return x, helmholtz_matrix
end

function v_cycle_helmholtz!(n, m, h, x, b, h_matrix_level1, sl_matrix_level1, h_matrix_level2, sl_matrix_level2, h_matrix_level3, sl_matrix_level3; u = 1, v1_iter = 1, v2_iter = 10, use_gmres_alpha = 0, alpha= 0.5, log = 0, level = nothing)

    if level == 3
        h_matrix = h_matrix_level3
        sl_matrix = sl_matrix_level3
    elseif level == 2
        h_matrix = h_matrix_level2
        sl_matrix = sl_matrix_level2
    else
        h_matrix = h_matrix_level1
        sl_matrix = sl_matrix_level1
    end

    # Relax on Ax = b v1_iter times with initial guess x
    x = jacobi_helmholtz_method!(n, m, h, x, b, sl_matrix; max_iter=v1_iter, use_gmres_alpha=use_gmres_alpha)

    if( n % 2 == 0 && n > 4 && m % 2 == 0 && m > 4 && (level == nothing || level > 1))

        # Compute residual on fine grid
        x_matrix = reshape(x, n-1, m-1, 1, 1)
        residual_fine = b - (helmholtz_chain!(real(x_matrix), h_matrix; h=h) + im*helmholtz_chain!(imag(x_matrix), h_matrix; h=h))[:,:,1,1] #) helmholtz_chain!(x_matrix, helmholtz_matrix; h=h)[:,:,1,1]

        # Compute residual, kappa and gamma on coarse grid
        residual_coarse = down(reshape(real(residual_fine), n-1, m-1, 1, 1)|>pu)[:,:,1,1] + im * down(reshape(imag(residual_fine), n-1, m-1, 1, 1)|>pu)[:,:,1,1]

        # Recursive operation of the method on the coarse grid
        n_coarse = size(residual_coarse,1)+1
        m_coarse = size(residual_coarse,2)+1
        x_coarse = zeros(c_type,n_coarse-1, m_coarse-1)|>pu

        for i = 1:u
            x_coarse, _ = v_cycle_helmholtz!(n_coarse, m_coarse, h*2, x_coarse, residual_coarse, h_matrix_level1, sl_matrix_level1, h_matrix_level2, sl_matrix_level2, h_matrix_level3, sl_matrix_level3; use_gmres_alpha = use_gmres_alpha,
                                                                    u=u, v1_iter=v1_iter, v2_iter=v2_iter, log=log, level = (level == nothing ? nothing : (level-1)))
        end
        x_coarse_matrix = reshape(x_coarse, n_coarse-1, m_coarse-1, 1, 1)

        # Correct
        fine_error = up(real(x_coarse_matrix)|>pu)[:,:,1,1] + im * up(imag(x_coarse_matrix)|>pu)[:,:,1,1]
        x = x + fine_error

        if log == 1
            r1 = residual_fine
            r2 = b - reshape(helmholtz_chain!(reshape(x, n-1, n-1, 1, 1), h_matrix; h=h), n-1, n-1)
            println("n = $(n), norm of x = $(norm(x)), norm of fine_error = $(norm(fine_error)), residual before vcycle =$(norm(r1)/norm(b)), residual after vcycle =$(norm(r2)/norm(b)), level =$(level)")
        end
    else
        # Coarsest grid
        A_Coarsest(v::a_type) = vec(helmholtz_chain!(real(reshape(v, n-1, m-1, 1, 1)), sl_matrix; h=h) + im*helmholtz_chain!(imag(reshape(v, n-1, m-1, 1, 1)), sl_matrix; h=h)) # vec(helmholtz_chain!(reshape(v, n-1, n-1, 1, 1), shifted_laplacian_matrix; h=h))
        M_Coarsest(v::a_type) = vec(jacobi_helmholtz_method!(n, m, h, x, reshape(v, n-1, m-1), sl_matrix)) #         M_Jacobi(n, h, x, sl_matrix, 1, v; use_gmres_alpha=use_gmres_alpha)
        x,flag,err,iter,resvec = fgmres_func(A_Coarsest, vec(b), v2_iter, tol=1e-15, maxIter=1,
                                                    M=M_Coarsest, x=vec(x), out=-1, flexible=true)
        x = reshape(x, n-1, m-1)
    end

    # Relax on Ax = b v1_iter times with initial guess x
    x = jacobi_helmholtz_method!(n, m, h, x, b, sl_matrix; max_iter=v1_iter, use_gmres_alpha=use_gmres_alpha)

    return x, h_matrix
end

function M_Jacobi(n, m, h, x, matrix, iterations, v;use_gmres_alpha=0)
    b = reshape(v, n-1, m-1)
    x = jacobi_helmholtz_method!(n, m, h, x, b, matrix; max_iter=iterations, use_gmres_alpha=use_gmres_alpha)
    return vec(x)
end

# Eran Code
function absorbing_layer!(gamma::Array,pad,ABLamp;NeumannAtFirstDim=false)

    n=size(gamma)

    #FROM ERAN ABL:

    b_bwd1 = ((pad[1]:-1:1).^2)./pad[1]^2;
	b_bwd2 = ((pad[2]:-1:1).^2)./pad[2]^2;

	b_fwd1 = ((1:pad[1]).^2)./pad[1]^2;
	b_fwd2 = ((1:pad[2]).^2)./pad[2]^2;
	I1 = (n[1] - pad[1] + 1):n[1];
	I2 = (n[2] - pad[2] + 1):n[2];

	if NeumannAtFirstDim==false
		gamma[:,1:pad[2]] += ones(n[1],1)*b_bwd2'.*ABLamp;
		gamma[1:pad[1],1:pad[2]] -= b_bwd1*b_bwd2'.*ABLamp;
		gamma[I1,1:pad[2]] -= b_fwd1*b_bwd2'.*ABLamp;
	end

	gamma[:,I2] +=  (ones(n[1],1)*b_fwd2').*ABLamp;
	gamma[1:pad[1],:] += (b_bwd1*ones(1,n[2])).*ABLamp;
	gamma[I1,:] += (b_fwd1*ones(1,n[2])).*ABLamp;
	gamma[1:pad[1],I2] -= (b_bwd1*b_fwd2').*ABLamp;
	gamma[I1,I2] -= (b_fwd1*b_fwd2').*ABLamp;

    return gamma
end

function fgmres_v_cycle_helmholtz!(n, m, h, b, kappa, omega, gamma; restrt=30, maxIter=10)
    shifted_laplacian_matrix, helmholtz_matrix = get_helmholtz_matrices!(kappa, omega, gamma; alpha=0.5)
    A(v) = vec(helmholtz_chain!(reshape(real(v), n-1, m-1, 1, 1), helmholtz_matrix; h=h) + im*helmholtz_chain!(reshape(imag(v), n-1, m-1, 1, 1), helmholtz_matrix; h=h)) # vec(helmholtz_chain!(reshape(v, n-1, n-1, 1, 1), helmholtz_matrix; h=h))
    function M(v)
        v = reshape(v, n-1, m-1)
        x = zeros(gmres_type,n-1,m-1)|>pu
        x, = v_cycle_helmholtz!(n, m, h, x, v, kappa, omega, gamma; u=1,
                    v1_iter = 1, v2_iter = 20, alpha=0.5, log = 0, level = 3)
        return vec(x)
    end
    x = zeros(gmres_type,n-1,m-1)|>pu
    x,flag,err,iter,resvec = fgmres_func(A, vec(b), restrt, tol=1e-30, maxIter=maxIter,
                                                    M=M, x=vec(x), out=-1, flexible=true)

    return reshape(x, n-1, m-1)
end

# not in use:
# function v_cycle_helmholtz_unet!(model, n, h, x, b, kappa, omega, gamma; u = 1, v1_iter = 1, v2_iter = 10, use_gmres_alpha = 0, alpha= 0.5, log = 0, level = nothing)
#
#     shifted_laplacian_matrix, helmholtz_matrix = get_helmholtz_matrices!(kappa, omega, gamma; alpha=alpha)
#
#     # Relax on Ax = b v1_iter times with initial guess x
#     x = jacobi_helmholtz_method!(n, h, x, b, shifted_laplacian_matrix; max_iter=v1_iter, use_gmres_alpha=use_gmres_alpha)
#
#     if( n % 2 == 0 && n > 4 && (level == nothing || level > 1))
#
#         # Compute residual on fine grid
#         x_matrix = reshape(x, n-1, n-1, 1, 1)
#         residual_fine = b - helmholtz_chain!(x_matrix, helmholtz_matrix; h=h)[:,:,1,1]
#
#         # Compute residual, kappa and gamma on coarse grid
#         residual_coarse = down(reshape(residual_fine, n-1, n-1, 1, 1))[:,:,1,1]
#         kappa_coarse = down(reshape(kappa, n-1, n-1, 1, 1))[:,:,1,1]
#         gamma_coarse = down(reshape(gamma, n-1, n-1, 1, 1))[:,:,1,1]
#
#         # Recursive operation of the method on the coarse grid
#         n_coarse = size(residual_coarse,1)+1
#         x_coarse = zeros(c_type,n_coarse-1, n_coarse-1)
#
#         for i = 1:u
#             x_coarse, helmholtz_matrix_coarse = v_cycle_helmholtz!(n_coarse, h*2, x_coarse, residual_coarse, kappa_coarse, omega, gamma_coarse; use_gmres_alpha = use_gmres_alpha,
#                                                                     u=u, v1_iter=v1_iter, v2_iter=v2_iter, log=log, level = (level == nothing ? nothing : (level-1)))
#         end
#         x_coarse_matrix = reshape(x_coarse, n_coarse-1, n_coarse-1, 1, 1)
#
#         # Correct
#         fine_error = up(x_coarse_matrix)[:,:,1,1]
#         x = x + fine_error
#
#         if log == 1
#             r1 = residual_fine
#             r2 = b - reshape(helmholtz_chain!(reshape(x, n-1, n-1, 1, 1), helmholtz_matrix; h=h), n-1, n-1)
#             println("n = $(n), norm of x = $(norm(x)), norm of fine_error = $(norm(fine_error)), residual before vcycle =$(norm(r1)/norm(b)), residual after vcycle =$(norm(r2)/norm(b)), level =$(level)")
#         end
#     else
#         # Coarsest grid
#         A_Coarsest(v) = vec(helmholtz_chain!(reshape(v, n-1, n-1, 1, 1), helmholtz_matrix; h=h))
#
#         function M_Unet(r)
#             r = reshape(r, n-1, n-1)
#             rj = reshape(r, n-1, n-1)
#             e = zeros(c_type, n-1, n-1)
#
#             # Relax
#             ej = jacobi_helmholtz_method!(n, h, e, r, shifted_laplacian_matrix)
#             rj = r - reshape(A(ej), n-1, n-1)
#
#             input = complex_grid_to_channels!(reshape(rj , n-1, n-1, 1, 1) , n)
#             input = cat(input, reshape(kappa, n-1, n-1, 1, 1), reshape(gamma, n-1, n-1, 1, 1), dims=3)
#             e_unet = (model(u_type.(input)|>cgpu)|>cpu)
#             e_unet = reshape((e_unet[:,:,1,1] + im*e_unet[:,:,2,1]), n-1, n-1, 1, 1) .* (h^2)
#             e = reshape(e_unet, n-1, n-1)
#
#             # Correction
#             e = ej + e
#             # Relax
#             e = jacobi_helmholtz_method!(n, h, e, r, shifted_laplacian_matrix)
#
#             return vec(e)
#         end
#
#         x,flag,err,iter,resvec = fgmres_func(A_Coarsest, vec(b), v2_iter, tol=1e-15, maxIter=1,
#                                                     M=M_Unet, x=vec(x), out=-1, flexible=true)
#         x = reshape(x, n-1, n-1)
#     end
#
#     # Relax on Ax = b v1_iter times with initial guess x
#     x = jacobi_helmholtz_method!(n, h, x, b, shifted_laplacian_matrix; max_iter=v1_iter, use_gmres_alpha=use_gmres_alpha)
#
#     return x, helmholtz_matrix
# end
# function compute_gmres_alpha!(residual, helmholtz_matrix, n, h)
#     helmholtz_residual = helmholtz_chain!(reshape(residual, n-1, n-1, 1, 1), helmholtz_matrix; h=h)[:,:,1,1]
#     residual_vector = reshape(residual,length(residual),1)
#     helmholtz_residual_vector = reshape(helmholtz_residual,length(helmholtz_residual),1)
#     alpha = (residual_vector' * helmholtz_residual_vector) ./ (helmholtz_residual_vector' * helmholtz_residual_vector)
#     return alpha
# end


text:
const USE_GPU = false
using ParallelStencil
using ParallelStencil.FiniteDifferences2D
@static if USE_GPU
    @init_parallel_stencil(CUDA, Float64, 2)
else
    @init_parallel_stencil(Threads, Float64, 2)
end
using Plots, Printf

macro H_dx(ix, iy)
    esc(:(H[$ix+1, $iy] - H[$ix, $iy]))
end
macro H_dy(ix, iy)
    esc(:(H[$ix, $iy+1] - H[$ix, $iy]))
end
macro Hv_dy(ix, iy)
    esc(:(H[$ix, $iy+1] * v[$ix, $iy+1] - H[$ix, $iy] * v[$ix, $iy]))
end
macro Hv_py(ix, iy)
    esc(:(H[$ix, $iy+1] * v[$ix, $iy+1] + H[$ix, $iy] * v[$ix, $iy]))
end
macro Huv_px(ix, iy)
    esc(:(H[$ix+1, $iy] * u[$ix+1, $iy] * v[$ix+1, $iy] + H[$ix, $iy] * u[$ix, $iy] * v[$ix, $iy]))
end
macro Huv_dx(ix, iy)
    esc(:(H[$ix, $iy+1] * v[$ix, $iy+1] * u[$ix, $iy+1] - H[$ix, $iy] * v[$ix, $iy] * u[$ix, $iy]))
end

"""
	compute_H!(H_new, H, u, v, dt_dx, dt_dy, size_H1_2, λx_2, λy_2)

ParallelStencil kernel that computes a time step of the shallow water 2D
process.
"""
@parallel_indices (ix, iy) function compute_H!(
    H_new,
    H,
    u,
    v,
    dt_dx,
    dt_dy,
    size_H1_2,
    λx_2,
    λy_2,
)
    if (ix <= size_H1_2 && iy <= size(H, 2))
        temp = H[ix+1, iy] -
               dt_dx * (
            1 / 2 * (H[ix+2, iy] * u[ix+2, iy] - H[ix, iy] * u[ix, iy]) -
            λx_2 * (@H_dx(ix + 1, iy) - @H_dx(ix, iy))
        )
        if (iy == size(H, 2))
            H_new[ix+1, iy] = (
                temp + dt_dy * (1 / 2 * @Hv_py(ix + 1, iy - 1) - λy_2 * @H_dy(ix + 1, iy - 1))
            )
        elseif (iy == 1)
            H_new[ix+1, iy] = (
                temp - dt_dy * (1 / 2 * @Hv_py(ix + 1, iy) - λy_2 * @H_dy(ix + 1, iy))
            )
        else
            H_new[ix+1, iy] = (
                temp - dt_dy * (1 / 2 * (@Hv_dy(ix + 1, iy) + @Hv_dy(ix + 1, iy - 1)) - λy_2 * (@H_dy(ix + 1, iy) - @H_dy(ix + 1, iy - 1)))
            )
        end
    end
    return
end

"""
    compute_u!(u_new, u, H_new, H, dt_dx, dt_dy, size_H1_2, g)

ParallelStencil kernel that computes the Lax-Friedrichs for the momentum in x-direction of the shallow water
equations.
"""
@parallel_indices (ix, iy) function compute_u!(
    u_new,
    u,
    H_new,
    H,
    dt_dx,
    dt_dy,
    size_H1_2,
    g,
)
    if (ix <= size_H1_2 && iy <= size(H, 2))
        u_new[ix+1, iy] =
            (
                1 / 2 * (H[ix, iy] * u[ix, iy] + H[ix+2, iy] * u[ix+2, iy]) -
                dt_dx * 1 / 2 * (
                    H[ix+2, iy] * u[ix+2, iy]^2 + 0.5 * g * H[ix+2, iy]^2 -
                    H[ix, iy] * u[ix, iy]^2 - 0.5 * g * H[ix, iy]^2
                )
            ) / H_new[ix+1, iy]
    end
    return
end

"""
    compute_v!(v_new, v, H_new, H, u, dt_dx, dt_dy, size_H2_2, g)

ParallelStencil kernel that computes the Lax-Friedrichs for the momentum in y-direction of the shallow water
equations.
"""
@parallel_indices (ix, iy) function compute_v!(
    v_new,
    v,
    H_new,
    H,
    u,
    dt_dx,
    dt_dy,
    size_H2_2,
    g,
)
    if (ix <= size(H, 1) && iy <= size_H2_2)
        temp = 1 / 2 * (H[ix, iy] * v[ix, iy] + H[ix, iy+2] * v[ix, iy+2]) - dt_dy * 1 / 2 * (
            H[ix, iy+2] * v[ix, iy+2]^2 - H[ix, iy] * v[ix, iy]^2 +
            0.5 * g * H[ix, iy+2]^2 - 0.5 * g * H[ix, iy]^2
        )
        if (ix == size(H, 1))
            v_new[ix, iy+1] =
                (
                    temp + dt_dx * 1 / 2 * @Huv_px(ix - 1, iy + 1)
                ) / H_new[ix, iy+1]
        elseif (ix == 1)
            v_new[ix, iy+1] =
                (
                    temp - dt_dx * 1 / 2 * @Huv_px(ix, iy + 1)
                ) / H_new[ix, iy+1]
        else
            v_new[ix, iy+1] =
                (
                    temp - dt_dx * 1 / 2 * (@Huv_dx(ix + 1, iy + 1) + @Huv_dx(ix - 1, iy + 1))
                ) / H_new[ix, iy+1]

        end
    end
    return
end

"""
	bc_x!(A)

Copy boundary conditions in x dimension of array `A`.
"""
@parallel_indices iy function bc_x!(A::Data.Array)
    A[1, iy] = A[2, iy]
    A[end, iy] = A[end-1, iy]
    return
end

"""
	shallow_water_2D_xpu(; nx, ny, dam1D_x=true, dam2D=false, do_visu=false)

2D shallow water equations solver for an instantaneous dam break.
The Lax-Friedrichs Method was applied to the continuity equation.
Geometry (length of 40 meters, width of 20 meters) and initial conditions
(half of domain has initial water level of 5 meters, other half is dry)
matches BASEMENT version 2.8 Test Case H_1 "Dam break in a closed channel."
Solution of momentum equations requires division by H. Therefore a minimum
H must be set throughout the domain, to avoid numerical instability.

# Arguments
    - `nx`: number of discretised cells for x dimension.
    - `ny`: number of discretised cells for y dimension.
    - `dam1D_x`: if true, 1D dam break in x-direction, else 1D dam break in y-direction. dam2D needs to be false to activate 1D dam break.
    - `dam2D`: if true, 2D dam break
    - `do_visu`: if true, each physical time step will be ploted.

# Return values
    - `H`: The solution arrray (Water surface height in m).
"""
@views function shallow_water_2D_xpu(;
    # Numerics
    nx = 512,
    ny = 256,
    # Modeling
    dam1D_x = true,
    dam2D = false,
    # Visualisation
    do_visu = false
)
    # Physics
    Lx, Ly = 40.0, 20.0
    g = 9.81
    u_max = 5 #from review of results with very small timesteps
    H_init = 5.0
    # Numerics
    nout = 100
    H_min = 1.0
    # Derived numerics
    dx, dy = Lx / nx, Ly / ny
    nt = 10000
    xc, yc = LinRange(dx / 2, Lx - dx / 2, nx), LinRange(dy / 2, Ly - dy / 2, ny)

    # Array initialisation
    H = @zeros(nx, ny) .+ H_min
    if (dam1D_x && !dam2D)
        H[1:round(Int64(nx / 2)), :] = @zeros(round(Int64(nx / 2)), ny) .+ H_init #1D dam break in x-direction
    elseif (!dam1D_x && !dam2D)
        H[:, 1:round(Int64(ny / 2))] = zeros(Float64, nx, round(Int64(ny / 2))) .+ H_init #1D dam break in y direction
    else
        H[1:round(Int64(nx/2)), 1:round(Int64(ny/2))] = zeros(Float64, round(Int64(nx/2)), round(Int64(ny/2))) .+ H_init #2D dam break
    end
    H_new = copy(H)
    u = @zeros(nx, ny)
    u_new = copy(u)
    v = @zeros(nx, ny)
    v_new = copy(v)
    size_H1_2, size_H2_2 = size(H, 1) - 2, size(H, 2) - 2
    t_tic = 0.0
    niter = 0
    time = 0
    # Time loop
    for it = 0:nt
        if (it == 10)
            t_tic = Base.time()
            niter = 0
        end
        dt = dx / (max(maximum(abs.(u)), maximum(abs.(v)), u_max)) / 4.1
        time = time + dt
        λx = dx / 2 / dt #x-direction signal speed for Lax-Friedrichs scheme
        λy = dy / 2 / dt #y-direction signal speed for Lax-Friedrichs scheme
        dt_dx, dt_dy = dt / dx, dt / dy
        λx_2, λy_2 = λx / 2, λy / 2

        #Lax-Friedrichs for Continuity
        @parallel compute_H!(H_new, H, u, v, dt_dx, dt_dy, size_H1_2, λx_2, λy_2)

        #Lax-Friedrichs for Momentum in x-direction
        @parallel compute_u!(u_new, u, H_new, H, dt_dx, dt_dy, size_H1_2, g)

        #Lax-Friedrichs for Momentum in y-direction
        @parallel compute_v!(v_new, v, H_new, H, u, dt_dx, dt_dy, size_H2_2, g)

        #No-flux boundary conditions
        @parallel 1:size(H, 2) bc_x!(H_new)
        u_new[1, :] .= 0
        u_new[end, :] .= 0
        v_new[:, 1] .= 0
        v_new[:, end] .= 0

        H .= H_new
        u .= u_new
        v .= v_new

        niter += 1
        if do_visu && (it % nout == 0)
            z_max = maximum(H)
            z_min = minimum(H)
            plot(xc, yc, Array(H)', st = :surface, xlims = (xc[1], xc[end]), ylims = (yc[1], yc[end]), zlims = (z_min, z_max), c = :blues, colorbar = false, xlabel = "Lx", ylabel = "Ly", zlabel = "water surface elevation (m)", title = "time = $(round(time, sigdigits=3)) s")
        end
    end
    t_toc = Base.time() - t_tic
    A_eff = (3 * 2) / 1e9 * nx * ny * sizeof(Float64)  # Effective main memory access per iteration [GB]
    t_it = t_toc / niter                      # Execution time per iteration [s]
    T_eff = A_eff / t_it                       # Effective memory throughput [GB/s]
    @printf(
        "Time = %1.3f sec, T_eff = %1.2f GB/s (niter = %d)\n",
        t_toc,
        round(T_eff, sigdigits = 3),
        niter
    )
    return Array(H)
end

# shallow_water_2D_xpu(;nx=128, ny=64, do_visu=true)