text:
struct Sector{S<:AbstractCharge, N}
    charges :: NTuple{N, S}

    Sector{S,N}(cs) where {S,N} = new{S,N}(cs)
    Sector{S,N}(cs::S...) where {S<:AbstractCharge,N} = new{S,N}(cs)
    function Sector{S}(cs...) where {S}
        new{S,length(cs)}(Tuple(convert(S,c) for c in cs))
    end
    Sector(cs::NTuple{N, S}) where {S, N} = new{S,N}(cs)
    Sector(cs::S...) where {S<:AbstractCharge} = new{S, length(cs)}(cs)
end

@inline Base.length(::Sector{S,N}) where {S,N} = N
@inline function Base.sum(s::Sector{S, N}, dualinfo::NTuple{N, Bool}) where{S,N}
    r = zero(S)
    for i=1:N
        r += dualinfo[i] ? inv(s[i]) : s[i]
    end
    r
end
@inline Base.inv(s::Sector{S}) where{S} = Sector{S}([inv(c) for c in s]...)
Base.getindex(s::Sector, i) = getindex(s.charges, i)
Base.iterate(s::Sector) = iterate(s.charges)
Base.iterate(s::Sector, i) = iterate(s.charges, i)

@inline layout(space::NTuple{N, VectorSpace{S}},
               s::Sector{S, N}) where {S, N} =
                   Tuple(layout(space[i])[s[i]] for i=1:N)

"""
we always sort sectors based on charges The sorting is column
major. That means the leftmost (first) VSpace changes charge faster
(first)!
"""
@inline Base.isless(s1::Sector{S, N}, s2::Sector{S, N}) where {N, S} =
    reverse(s1.charges) < reverse(s2.charges)

# @inline function Base.isless(s1::Sector{S, N},
#                              s2::Sector{S, N}) where {N, S}
#     for n in N:-1:1
#         @inbounds a, b = s1[n], s2[n]
#         if a < b
#             return true
#         elseif a > b
#             return false
#         end
#     end
#     false
# end

function _check_sectordatadict(charge, space, data)
    sects, sizes = _allsectorsandsizes(charge, space)
    i = 1
    for (s,d) in data
        s == sects[i] || throw{SectorMismatch()}
        size(d) == sizes[i] || throw(SizeMismatch())
        i += 1
    end
end

"""
    _allsectorsandsizes(charge, Vs)

Generate all sectors for a tuple of VSpace that add up to a given
`charge` and are thus sectors are the corresponding SymTensor. This is
recursive function calling itself with Vs[1:N-1], etc. The output
sector is sorted by construction.
"""
function _allsectorsandsizes(charge::S, Vs::NTuple{N, VectorSpace{S}}) where{S, N}
    if N < 2
        _c = Vs[1].isdual ? inv(charge) : charge
        d = dim(Vs[1], _c)
        if d > 0
            return [Sector(_c)], [(d,)]
        end
        return Sector{S, 1}[], NTuple{1, Int}[]
    end

    sects = Sector{S, N}[]
    sizes = NTuple{N, Int}[]
    for (c,d) in Vs[N]
        _c = Vs[N].isdual ? inv(c) : c
        sect_head, size_head = _allsectorsandsizes(charge - _c, Vs[1:N-1])
        for i in eachindex(sect_head)
            push!(sects, Sector(sect_head[i].charges..., c))
            push!(sizes, (size_head[i]..., d))
        end
    end
    sects, sizes
end
@inline _allsectorsandsizes(c::Int, Vs::NTuple{N, VectorSpace{S}}) where{S,N} =
    _allsectorsandsizes(convert(S, c), Vs)

function Base.show(io::IO, s::Sector)
    if !get(io, :compact, false)
        show(io, typeof(s))
    end
    print(io, "(")
    separator = ""
    comma = ", "
    io2 = IOContext(io, :typeinfo => typeof(s))
    for c in s.charges
        print(io2, separator, c)
        separator=comma
    end
    print(io2, ")")
    return nothing
end

max_stars_repo_name: amiragha/mps

max_issues_repo_name: amiragha/mps

max_forks_repo_name: amiragha/mps


text:
function maptocolours(x,cmap)
    m=minimum(x)
    M=maximum(x)
    if m == M
        return [cmap[128] for i in 1:length(x)]
    end
    y = Int.(floor.(255*(x.-m)/(M-m)))
    return [cmap[i+1] for i in y]
end

"""
    Train model on data, against the loss. Plots the network (nodes/edges) with colours
    representing the relative connection strength (blue = negative, white = 0, red = positive)
    compared to all others (biases and weights done separately).

    Plots the current performance on data as well.

    Saves an mp4 of the process.
"""
function WatchNetwork(model::Flux.Chain,data, loss::Function, truth; Epochs = 100, opt=Descent(0.01),
        debug=false,minwidth=1.0,maxwidth=4.0,res=(1920,1080),fig=nothing, neuron_colour=RGBAf0(0,0,0,1),
        bias_colour=RGBAf0(0.5,0.5,0.5,1),fps=60,titlesize=30,filename="Makie.mp4"
)::Figure

    fig = Figure(resolution=res,fontsize=titlesize)

    main_title = Node("Network") # will update to track epochs
    axes = [Axis(fig[1,1],title=main_title)]

    title = Node("")

    if data != nothing
        push!(axes,Axis(fig[1,2],title=title))
    end

    time_ = Node(1)

    Neurons = Vector{Point2{Float64}}()
    Biases = Vector{Point2{Float64}}()
    Connections = Vector{Point2{Float64}}()
    ConnectionStrength = Vector{Float64}()
    nColours = Vector{PerceptualColourMaps.RGBA{Float64}}()
    bColours = Vector{PerceptualColourMaps.RGBA{Float64}}()

    cm = cmap("D1")

    for layer in model.layers
        if !(typeof(layer)<:Flux.Dense)
            @warn "Model must only have Dense layers"
	    @info "If your model has a function at the end\n like softmax, put it into the loss"
            @assert typeof(layer)<:Flux.Dense
        end
    end

    Ns = [(size(layer.weight,2),size(layer.weight,1)) for layer in model.layers]
    biases = true

    spacing_x = 8.0
    spacing_y = 4.0
    neuron_size = 1.0
    bias_size = 0.5
    bias_placement = [0*bias_size,+2*bias_size] # gets added to neuron position

    pos = [(1,1)]

    # maximum bias magnitude
    B = 0.0
    for l in model.layers
        for b in abs.(l.bias)
            if b > B
                B = b
            end
        end
    end

    # collect neurons and biases into arrays
    for l in 1:length(Ns)
        start = 0.0
        ns = Ns[l][1]
        if ns > 1
            if ns % 2 == 0
                start = -spacing_y * (ns / 2)+spacing_y/2.
            else
                start = -spacing_y * floor(ns / 2)
            end
        end
        for n in 1:Ns[l][1]
            # add a neuron
            x,y = spacing_x*(l-1),start+spacing_y*(n-1) # with padding
            neuron = Point2(x,y)
            push!(Neurons,neuron)
            push!(nColours,neuron_colour)
        end

        start = 0.0
        ns = Ns[l][2]
        if ns > 1
            if ns % 2 == 0
                start = -spacing_y * (ns / 2)+spacing_y/2.
            else
                start = -spacing_y * floor(ns / 2)
            end
        end
        if biases
            for n in 1:Ns[l][2]
                x,y = spacing_x*(l),start+spacing_y*(n-1)
                # add it's bias
                bias = Point2(x+bias_placement[1],y+bias_placement[2])
                push!(Biases,bias)
                push!(bColours,bias_colour)
                # add a connection for the bias
                push!(Connections,Point2(x,y))
                push!(Connections,Point2(x+bias_placement[1],y+bias_placement[2]))
                # relative magnitude*sign
                mag = sign(model.layers[l].bias[n])*maxwidth*(minwidth+abs.(model.layers[l].bias[n])/B)
                isnan(mag) ? mag = 0 : nothing
                push!(ConnectionStrength,mag)
                push!(ConnectionStrength,mag)
            end
        end

        # final layer
        if l == length(Ns)
            # draw output
            start = 0.0
            ns = Ns[l][2]
            if ns > 1
                if ns % 2 == 0
                    start = -spacing_y * (ns / 2)+spacing_y/2.
                else
                    start = -spacing_y * floor(ns / 2)
                end
            end
            for n in 1:Ns[l][2]
                # add a neuron
                x,y = spacing_x*(l),start+spacing_y*(n-1)

                neuron = Point2(x,y)
                push!(Neurons,neuron)
                push!(nColours,neuron_colour)
            end
        end

        # largest weight magnitude
        M = 0.0
        for l in model.layers
            for w in abs.(l.W)
                if w > M
                    M = w
                end
            end
        end

        # draw connections
        for n1 in 1:Ns[l][1]
            for n2 in 1: Ns[l][2]
                start1 = 0.0
                start2 = 0.0
                ns1 = Ns[l][1]
                if ns1 > 1
                    if ns1 % 2 == 0
                        start1 = -spacing_y * (ns1 / 2)+spacing_y/2.
                    else
                        start1 = -spacing_y * floor(ns1 / 2)
                    end
                end

                ns2 = Ns[l][2]
                if ns2 > 1
                    if ns2 % 2 == 0
                        start2 = -spacing_y * (ns2 / 2)+spacing_y/2.
                    else
                        start2 = -spacing_y * floor(ns2 / 2)
                    end
                end

                x1,y1 = spacing_x*(l-1),start1+spacing_y*(n1-1)
                x2,y2 = spacing_x*(l),start2+spacing_y*(n2-1)
                push!(Connections,Point2(x1,y1))
                push!(Connections,Point2(x2,y2))
                # relative magnitude*sign
                mag = sign(model.layers[l].weight[n2,n1]/M)*maxwidth*(minwidth+abs.(model.layers[l].weight[n2,n1])/M)
                isnan(mag) ? mag = 0 : nothing
                push!(ConnectionStrength,mag)
                push!(ConnectionStrength,mag)
            end
        end
    end

    # make the arrays into nodes so we can update them
    Neurons = Node(Neurons)
    Biases = Node(Biases)
    Connections = Node(Connections)
    ConnectionStrength = Node(abs.(ConnectionStrength))
    ConnectionColour = Node(maptocolours(ConnectionStrength[],cm))
    nColours = Node(nColours);
    bColours = Node(bColours);

    # plot neurons, biases, and weights
    scatter!(axes[1],Neurons,color=nColours,markersize=(neuron_size,neuron_size),markerspace=SceneSpace)
    scatter!(axes[1],Biases,color=bColours,markersize=(bias_size,bias_size),markerspace=SceneSpace)
    linesegments!(axes[1],Connections,linewidth=ConnectionStrength,color=ConnectionColour)

    ns = maximum([max(l[1],l[2]) for l in Ns])

    l = length(model.layers)

    # determine position of viewport
    # so the network is drawn centered

    mx = 0+bias_placement[1]-bias_size
    Mx = l*spacing_x + neuron_size

    my = 0
    My = 0

    if ns % 2 == 0
        my = -spacing_y * (0.5+ns / 2)-bias_placement[2] - bias_size
        My = spacing_y * (ns / 2)+ bias_placement[2] + bias_size
    else
        my = -spacing_y * floor(ns / 2) - bias_placement[2] - bias_size
        My = spacing_y * floor(ns / 2) + bias_placement[2] + bias_size
    end

    dx = Mx-mx
    dy = My-my

    if dx < dy
        Mx = Mx+(dy-dx)
    elseif dy < dx
        My = My +(dx-dy)
    end

    # re-center
    pad = floor(l/2.)*spacing_x
    if l % 2 != 0
        pad += spacing_x/2.0
    end

    c = ((Mx+mx)/2.0-pad,(My+my)/2.0)
    # set viewport
    axes[1].aspect = AxisAspect(1)
    limits!(axes[1],mx-c[1],Mx-c[1],my-c[2],My-c[2])
    if debug == false
        hidedecorations!(axes[1,1])
        hidespines!(axes[1,1])
    end
    if length(data[1][1])==1 && length(data[1][2])==1
        # 1-d problem plotting code
        p = [model(data[i][1])[1] for i in 1:length(data)]|>Node
        if data != nothing
            x = [data[i][1][1] for i in 1:length(data)]
            scatter!(axes[2],x,truth)
            lines!(axes[2],x,p,linewidth=2.0,color="red")
            e = sum([loss(model,data[i][1],data[i][2]) for i in 1:length(data)])/length(data)
            title[]="Error: $(round(e,digits=5))"
        end
    elseif length(data[1][1])==2 && length(data[1][2])==2
        # 2-d problem plotting code
        p=[model(data[i][1]) for i in 1:length(data)]|>Node
        if data != nothing
            axes[2].aspect = AxisAspect(1)
            hidedecorations!(axes[2])
            N = length(truth)|>floor|>Int
            n = sqrt(N)|>floor|>Int

            c = reshape([RGBAf0(x[1],x[2],0,1) for x in truth],n,n)
            pos = Vector{Point2}()
            col = Vector{RGBAf0}()
            for i in 1:size(c,1)
                for j in 1:size(c,2)
                    push!(pos,Point2(i-0.5,j-0.5))
                    push!(col,c[i,j])
                end
            end
            push!(axes,Axis(
                fig[1,3],title="Goal"
            ))
            axes[3].aspect = AxisAspect(1)
            hidedecorations!(axes[3])
            scatter!(axes[3],pos,marker=Rect2D(0,0,1,1), # like a heatmap
                    markersize=(1,1),markerspace=SceneSpace,
                    color = col
            )

            q = reshape([RGBAf0(x[1],x[2],0,1) for x in p[][1:N]],n,n)|>Node
            pos = Vector{Point2}()
            col = Vector{RGBAf0}()
            for i in 1:size(q[],1)
                for j in 1:size(q[],2)
                    push!(pos,Point2(i-0.5,j-0.5))
                    push!(col,q[][i,j])
                end
            end
            pos=Node(pos)
            col = Node(col)
            scatter!(axes[2],pos,marker=Rect2D(0,0,1,1),
                markersize=(1,1),markerspace=SceneSpace,
                color = col
            )
            e = sum([loss(model,data[i][1],data[i][2]) for i in 1:length(data)])/length(data)
            a = Flux.Losses.logitbinarycrossentropy([0,1],[1,0])
            e = (e-a) # normalise loss
            title[]="Error: $(round(e,digits=5))"
        end
    end

    prog = Progress(Epochs-1)

    # now begin training and recording

    record(fig, filename, collect(2:Epochs), framerate=fps) do k

        # train 1 epoch
        Flux.train!((x,y)->loss(model,x,y),Flux.params(model),data,opt)

        newConnectionStrength = []

        if data != nothing
            main_title[] = "Network (Epochs: $k)"
            if length(data[1][1])==1 && length(data[1][2])==1
                # update prediction node
                p[] = [model(data[i][1])[1] for i in 1:length(data)]
                e = sum([loss(model,data[i][1],data[i][2]) for i in 1:length(data)])/length(data)
                # title node with error
                title[] = "Error: $(round(e,digits=5))"
            elseif length(data[1][1])==2 && length(data[1][2])==2
                # update prediction node
                p[]=[model(data[i][1]) for i in 1:length(data)]
                N = length(truth)|>floor|>Int
                n = sqrt(N)|>floor|>Int
                q[] = reshape([RGBAf0(x[1],x[2],0,1) for x in p[][1:N]],n,n)
                newpos = []
                newcol = []
                for i in 1:size(q[],1)
                    for j in 1:size(q[],2)
                        push!(newpos,Point2(i-0.5,j-0.5))
                        push!(newcol,q[][i,j])
                    end
                end
                # update current performance node
                pos[] = newpos
                col[] = newcol
                e = sum([loss(model,data[i][1],data[i][2]) for i in 1:length(data)])/length(data)
                a = Flux.Losses.logitbinarycrossentropy([0,1],[0,1])
                e = (e-a)
                title[]="Error: $(round(e,digits=5))"
            end
        end

        # update colour nodes
        B = 0.0
        for l in model.layers
            for b in abs.(l.bias)
                if b > B
                    B = b
                end
            end
        end

        for l in 1:length(Ns)
            start = 0.0
            ns = Ns[l][1]

            if biases
                for n in 1:Ns[l][2]
                    mag = sign(model.layers[l].bias[n])*maxwidth*(minwidth+abs.(model.layers[l].bias[n])/B)
                    isnan(mag) ? mag = 0 : nothing
                    push!(newConnectionStrength,mag)
                    push!(newConnectionStrength,mag)
                end
            end

            M = 0.0
            for l in model.layers
                for w in abs.(l.W)
                    if w > M
                        M = w
                    end
                end
            end

            # draw connections
            for n1 in 1:Ns[l][1]
                for n2 in 1: Ns[l][2]
                    mag = sign(model.layers[l].weight[n2,n1]/M)*maxwidth*(minwidth+abs.(model.layers[l].weight[n2,n1])/M)
                    isnan(mag) ? mag = 0 : nothing
                    push!(newConnectionStrength,mag)
                    push!(newConnectionStrength,mag)
                end
            end
        end
        # propogate colour and magnitude updates
        ConnectionStrength[] = newConnectionStrength
        ConnectionColour[] = maptocolours(newConnectionStrength,cm)

        # display loss in progress
        e = sum([loss(model,data[i][1],data[i][2]) for i in 1:length(data)])/length(data)
        next!(prog,showvalues=[(:iter,k),(:loss,e)])
    end

    return fig
end

ext: jl

max_stars_repo_name: Jerboa-app/MLPVis

max_issues_repo_name: Jerboa-app/MLPVis

max_forks_repo_name: Jerboa-app/MLPVis


text:
# Butcher Tableau
const B = [[0.161],
        [-0.008480655492356992, 0.3354806554923570],
        [2.8971530571054944, -6.359448489975075, 4.362295432869581],
        [5.32586482843926, -11.74888356406283, 7.495539342889836, -0.09249506636175525],
        [5.86145544294642, -12.92096931784711, 8.159367898576159,  -0.07158497328140100,
        -0.02826905039406838],
        [0.09646076681806523, 0.01, 0.4798896504144996, 1.379008574103742,
        -3.290069515436081, 2.324710524099774]
        ]

# Step size fractions
const nodes = [0.161, 0.327, 0.9, 0.9800255409045097, 1, 1]

# Weights for 5th order method
const b5 = [0.001780011052226, 0.000816434459657, -0.007880878010262, 0.144711007173263,
            -0.582357165452555, 0.458082105929187, 1/66]
# Weights for 4th order method
const b4 = [0.09646076681806523, 0.01, 0.4798896504144996, 1.379008574103742,
                -3.290069515436081, 2.324710524099774, 0]
# Error estimate
const errest = b5 .- b4

#Interpolation coefficients
const interpC = hcat(
    [1.0, -2.763706197274826, 2.9132554618219126, -1.0530884977290216],
    [0, 0.13169999999999998, -0.2234, 0.1017],
    [0, 3.9302962368947516, -5.941033872131505, 2.490627285651253],
    [0, -12.411077166933676, 30.33818863028232, -16.548102889244902],
    [0, 37.50931341651104, -88.1789048947664, 47.37952196281928],
    [0, -27.896526289197286, 65.09189467479366, -34.87065786149661],
    [0, 1.5, -4, 2.5]
)


text:
@testset "Forcing" begin
    tolerances = Dict(Float16 => 1e-4, Float32 => 1e-5, Float64 => 1e-7, Real => 1e-7)
    @testset "Slippage one dimension with velocity" begin
        N = 8
        typelist = [Float16 Float32 Float64]
        for type in typelist
            delta_set = [1.0 0.5 0.0]
            sol_set = [1/110 2/115 1/20]
            # Test different delta values
            for delta in 1:3
                label = delta_set[delta]
                @testset "δ = $label" begin
                    moment = JuSwalbe.Macroquant(height=ones(type, N), velocity=fill(type(0.1), N), pressure=zeros(type, N), energy=zeros(type, N))
                    forces = JuSwalbe.Forces(slip=zeros(type, N), thermal=zeros(type, N), h∇p=zeros(type, N), bathymetry=zeros(type, N))
                    input = JuSwalbe.Inputconstants(δ=type(delta_set[delta]), τ=type(1.0))
                    slippage = computeslip(moment, forces, input)
                    @test isa(forces.slip, Array{type, 1})
                    @test length(slippage) == N

                    for i in 1:N
                        @test slippage[i] - type(sol_set[delta]) ≈ type(0.0) atol=tolerances[type]
                    end
                end
            end
        end
    end
    @testset "Slippage one dimension without velocity" begin
        N, M = 8, 9
        typelist = [Float16 Float32 Float64]
        for type in typelist
            # Test different delta values
            delta_set = [1.0 0.5 0.0]
            for delta in 1:3
                label = delta_set[delta]
                @testset "δ = $label" begin
                    moment = JuSwalbe.Macroquant(height=ones(type, N), velocity=zeros(type, N), pressure=zeros(type, N), energy=zeros(type, N))
                    forces = JuSwalbe.Forces(slip=zeros(type, N), thermal=zeros(type, N), h∇p=zeros(type, N), bathymetry=zeros(type, N))
                    input = JuSwalbe.Inputconstants(δ=type(delta_set[delta]), τ=type(1.0))
                    slippage = computeslip(moment, forces, input)
                    @test isa(forces.slip, Array{type, 1})
                    @test length(slippage) == N
                    for i in 1:N
                        @test slippage[i] ≈ type(0.0) atol=tolerances[type]
                    end
                end
            end
        end
    end
    @testset "Slippage two dimensions with velocity" begin
        N, M = 8, 9
        typelist = [Float16 Float32 Float64]
        for type in typelist
            # Test different delta values
            delta_set = [1.0 0.5 0.0]
            sol_set = [1/110 2/115 1/20]
            for delta in 1:3
                label = delta_set[delta]
                @testset "δ = $label" begin
                    vel = JuSwalbe.Twovector(x=fill(type(0.1), (N,M)), y=fill(type(-0.1), (N,M)))
                    zerovec = JuSwalbe.Twovector(x=zeros(type ,(N,M)), y=zeros(type ,(N,M)))
                    moment = JuSwalbe.Macroquant(height=ones(type, (N,M)), velocity=vel, pressure=zeros(type, (N,M)), energy=zeros(type, (N,M)))
                    forces = JuSwalbe.Forces(slip=zerovec, thermal=zerovec, h∇p=zerovec, bathymetry=zerovec)
                    input = JuSwalbe.Inputconstants(δ=type(delta_set[delta]), τ=type(1.0))
                    slippage = computeslip(moment, forces, input)
                    @test isa(zerovec, JuSwalbe.Twovector{Array{type, 2}})
                    @test isa(moment, JuSwalbe.Macroquant{Array{type, 2}, JuSwalbe.Twovector{Array{type, 2}}})
                    @test isa(forces, JuSwalbe.Forces{JuSwalbe.Twovector{Array{type, 2}}})
                    @test isa(input, JuSwalbe.Inputconstants)

                    for i in 1:N, j in 1:M
                        @test slippage.x[i,j] + type(sol_set[delta]) ≈ type(0.0) atol=tolerances[type]
                        @test slippage.y[i,j] - type(sol_set[delta]) ≈ type(0.0) atol=tolerances[type]
                    end
                end
            end
        end
    end
    @testset "Slippage two dimensions without velocity" begin
        N, M = 8, 9
        typelist = [Float16 Float32 Float64]
        for type in typelist
            # Test different delta values
            delta_set = [1.0 0.5 0.0]
            sol_set = [1/110 2/115 1/20]
            for delta in 1:3
                label = delta_set[delta]
                @testset "δ = $label" begin
                    vel = JuSwalbe.Twovector(x=fill(type(0.1), (N,M)), y=fill(type(-0.1), (N,M)))
                    zerovec = JuSwalbe.Twovector(x=zeros(type ,(N,M)), y=zeros(type ,(N,M)))
                    moment = JuSwalbe.Macroquant(height=ones(type, (N,M)), velocity=vel, pressure=zeros(type, (N,M)), energy=zeros(type, (N,M)))
                    forces = JuSwalbe.Forces(slip=zerovec, thermal=zerovec, h∇p=zerovec, bathymetry=zerovec)
                    input = JuSwalbe.Inputconstants(δ=type(delta_set[delta]), τ=type(1.0))
                    slippage = computeslip(moment, forces, input)

                    for i in 1:N, j in 1:M
                        @test slippage.x[i,j] + type(sol_set[delta]) ≈ type(0.0) atol=tolerances[type]
                        @test slippage.y[i,j] - type(sol_set[delta]) ≈ type(0.0) atol=tolerances[type]
                    end
                end
            end
        end
    end
    @testset "Thermo capillary forcing two dimensions" begin
        N, M = 8, 9
        typelist = [Float16 Float32 Float64]
        for type in typelist
            # Test different delta values
            mom = simplemoment2d(N, M; T=type)
            force = JuSwalbe.Forces(slip=JuSwalbe.Twovector(x=fill(type(0.1), (N,M)), y=fill(type(0.1), (N,M))), h∇p = simpleTwovector(N, M; T=type),
                                    thermal=JuSwalbe.Twovector(x=zeros(type, (N,M)), y=zeros(type, (N,M))), bathymetry=JuSwalbe.Twovector(x=zeros(type, (N,M)), y=zeros(type, (N,M))))
            @test isa(force, JuSwalbe.Forces{JuSwalbe.Twovector{Matrix{type}}})
            @test isa(force.thermal.x, Array{type, 2})
            @test isa(force.thermal.y, Array{type, 2})
            @test size(force.thermal.x) == (N,M)
            @test size(force.thermal.y) == (N,M)
        end
    end
    @testset "Thermo capillary forcing one dimension" begin
        N = 8
        typelist = [Float16 Float32 Float64]
        for type in typelist
            # Test different delta values
            mom = simplemoment1d(N; T=type)
            force = JuSwalbe.Forces(slip=fill(type(0.1), N) , h∇p = zeros(type, N),
                                    thermal=zeros(type, N), bathymetry=zeros(type, N))
            @test isa(force, JuSwalbe.Forces{Vector{type}})
            @test isa(force.thermal, Array{type, 1})
            @test length(force.thermal) == N
        end
    end
end

max_stars_repo_name: Zitzeronion/JuSwalbe

max_issues_repo_name: Zitzeronion/JuSwalbe


text:
module GLExtender

export @glfunc

using ModernGL
using Libdl

# ModernGL functions for Extension
gl_lib = C_NULL

isOpenLib() = gl_lib != C_NULL
function openLib(); if Sys.iswindows() && !isOpenLib(); gl_lib = Libdl.dlopen("opengl32"); end; gl_lib; end
function closeLib(); if Sys.iswindows() && isOpenLib(); Libdl.dlclose(gl_lib); gl_lib=C_NULL; end; end

const enable_opengl_debugging  = false

function debug_opengl_expr(func_name, args)
    if enable_opengl_debugging && func_name != :glGetError
        quote
            err = glGetError()
            if err != GL_NO_ERROR
                arguments = gl_represent.(tuple($(args...)))
                warn("OpenGL call to $($func_name), with arguments: $(arguments)
                Failed with error: $(GLENUM(err).name).")
            end
        end
    else
        :()
    end
end

# based on getCFun macro
macro glfunc(opengl_func)
    arguments = map(opengl_func.args[1].args[2:end]) do arg
        isa(arg, Symbol) && return Expr(:(::), arg)
        arg
    end
    # Get info out of arguments of `opengl_func`
    arg_names       = map(arg->arg.args[1], arguments)
    return_type     = opengl_func.args[2]
    input_types     = map(arg->arg.args[2], arguments)
    func_name       = opengl_func.args[1].args[1]
    func_name_sym   = Expr(:quote, func_name)
    func_name_str   = string(func_name)
    ptr_expr        = :(ModernGL.getprocaddress_e($func_name_str))
    #debug_opengl_expr = ModernGL.debug_opengl_expr

    wasOpenLib=true

    if Sys.iswindows() # windows has some function pointers statically available and some not, this is how we deal with it:
        wasOpenLib=isOpenLib()
        lib=openLib()
        ptr = Libdl.dlsym_e(lib, func_name)
        if (ptr != C_NULL)
            ptr_expr = :(($func_name_sym, "opengl32"))
            ret = quote
                function $func_name($(arg_names...))
                    result = ccall($ptr_expr, $return_type, ($(input_types...),), $(arg_names...))
                    $(debug_opengl_expr(func_name, arg_names))
                    result
                end
                $(Expr(:export, func_name))
            end
            return esc(ret)
        end
    end
    ptr_sym = gensym("$(func_name)_func_pointer")
    ret = quote
        $ptr_sym = ModernGL.GLFunc(C_NULL)
        function $func_name($(arg_names...))
            if $ptr_sym.p::Ptr{Cvoid} == C_NULL
                $ptr_sym.p::Ptr{Cvoid} = $ptr_expr
            end
            result = ccall($ptr_sym.p::Ptr{Cvoid}, $return_type, ($(input_types...),), $(arg_names...))
            $(debug_opengl_expr(func_name, arg_names))
            result
        end
        $(Expr(:export, func_name))
        end
    r=esc(ret)
    if !wasOpenLib closeLib() end
    r
end

"""
TODO
"""
function getExtensions()
  count = GLint[0]
  glGetIntegerv(GL_NUM_EXTENSIONS, count)
  count=count[]
  exts=[]
  for i=1:count
    name=unsafe_string(glGetStringi(GL_EXTENSIONS, i-1))
    push!(exts,name)
    #if strcmp(ccc, (const GLubyte *)"GL_ARB_debug_output") == 0
    #  # The extension is supported by our hardware and driver
    #  # Try to get the "glDebugMessageCallbackARB" function :
    #  glDebugMessageCallbackARB  = (PFNGLDEBUGMESSAGECALLBACKARBPROC) wglGetProcAddress("glDebugMessageCallbackARB");
    #end
  end
  exts
end
export getExtensions

# function setExtensionsCallback()
#   if strcmp(ccc, (const GLubyte *)"GL_ARB_debug_output") == 0
#    # The extension is supported by our hardware and driver
#    # Try to get the "glDebugMessageCallbackARB" function :
#    glDebugMessageCallbackARB  = (PFNGLDEBUGMESSAGECALLBACKARBPROC) wglGetProcAddress("glDebugMessageCallbackARB");
#   end
# end

end #GLExtender

max_stars_repo_name: Gilga/JuliaOpenGL

max_issues_repo_name: Gilga/JuliaOpenGL

max_forks_repo_name: Gilga/JuliaOpenGL


text:
module StickBreak
include("custom_grads.jl") # cumsum, cumprod

"""
x: real vector of dim K - 1
return: simplex of dim K
"""
function transform(x::T) where T
  size_x = size(x)
  ndim = ndims(x)
  dim_head = size_x[1:end-1]
  K = size_x[end] + 1
  ks = cumsum(one.(x), dims=ndim)
  z = sigmoid.(x .- log.(K .- ks))
  one_minus_z_cumprod = cumprod(1 .- z, dims=ndim)

  if ndim == 1
    ones_pad =  Tracker.collect([one.(layer(x, 1))])
    p = vcat(z, ones_pad) .* vcat(ones_pad, one_minus_z_cumprod)
  else
    ones_pad = one.(layer(x, 1))
    p = cat(z, ones_pad, dims=ndim) .* cat(ones_pad, one_minus_z_cumprod, dims=ndim)
  end

  return p
end

"""
x: real vector of dim K - 1
p: simplex of dim K
return: log abs value of determinant of x
"""
function logabsdetJ(x, p)
  ndim = ndims(x)
  K = size(p)[end]
  ks = cumsum(one.(x), dims=ndim)
  z = sigmoid.(x .- log.(K .- ks))
  detJ = sum(log1p.(-z) .+ log.(head(p)), dims=ndim)

  return detJ
end

end # StickBreak


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

using Markdown
using InteractiveUtils

# ╔═╡ 514ab1cc-fb98-11eb-2a1e-c7b264f7066c
begin
	let
		using Pkg
		Pkg.activate(mktempdir())
		Pkg.Registry.update()
		Pkg.add("PlutoUI")
		Pkg.add("Tar")
		Pkg.add("MLDataPattern")
		Pkg.add("Glob")
		Pkg.add("NIfTI")
		Pkg.add("DLPipelines")
		Pkg.add("DataAugmentation")
		Pkg.add("CairoMakie")
		Pkg.add("ImageCore")
		Pkg.add("DataLoaders")
		Pkg.add("Flux")
		Pkg.add("CUDA")
		Pkg.add(url="https://github.com/FluxML/FastAI.jl")
	end

	using PlutoUI
	using Tar
	using MLDataPattern
	using Glob
	using NIfTI
	using DLPipelines
	using DataAugmentation
	using DataAugmentation: OneHot, Image
	using CairoMakie
	using ImageCore
	using DataLoaders
	using Flux
	using CUDA
	using FastAI
end

# ╔═╡ 5b2ed761-eb57-44c2-9ed9-bc1e5240af16
TableOfContents()

# ╔═╡ 3cfef02c-00f4-48e9-90b1-4ad43472f397
md"""
## Load data
Part of the [Medical Decathlon Dataset](http://medicaldecathlon.com/)
"""

# ╔═╡ 92776990-345c-4abc-82ad-3cec0bee356f
data_dir = "/Users/daleblack/Google Drive/Datasets/Task02_Heart"

# ╔═╡ bf88491c-a431-4575-a88e-6cc80afce35e
function loadfn_label(p)
	a = NIfTI.niread(string(p)).raw
	convert_a = convert(Array{UInt8}, a)
	return convert_a
end

# ╔═╡ ff620c1f-e08f-4e8d-b538-31a48cea2e0b
function loadfn_image(p)
	a = NIfTI.niread(string(p)).raw
	convert_a = convert(Array{Float32}, a)
	return convert_a
end

# ╔═╡ ae8852a9-5b80-4bda-b58f-3ff3d927bfb4
begin

	niftidata_image(dir) = mapobs(loadfn_image, Glob.glob("*.nii*", dir))
	niftidata_label(dir) =  mapobs(loadfn_label, Glob.glob("*.nii*", dir))
	data = (
		niftidata_image(joinpath(data_dir, "imagesTr")),
		niftidata_label(joinpath(data_dir, "labelsTr")),
	)
end

# ╔═╡ 2fd15108-ced5-4dc2-abb4-ceb7f394e718
md"""
## Create `LearningMethod`
"""

# ╔═╡ d2b61982-7829-486e-a0d9-c02fb6d77e74
struct ImageSegmentationSimple <: DLPipelines.LearningMethod
    imagesize
end

# ╔═╡ 69278a55-81b9-4936-9146-f2510774ff22
image_size = (112, 112, 96)

# ╔═╡ 55f8c444-3486-4cd8-a1af-a3774ad55594
method = ImageSegmentationSimple(image_size)

# ╔═╡ c981defe-f337-47ac-8e7b-ea646b4cf482
md"""
### Set up `encode` pipelines
"""

# ╔═╡ 65570158-ee5a-437d-937f-e37f47c07823
begin
	function DLPipelines.encode(
			method::ImageSegmentationSimple,
			context::Training,
			(image, target)::Union{Tuple, NamedTuple}
			)

		tfm_proj = RandomResizeCrop(method.imagesize)
		tfm_im = DataAugmentation.compose(
			ImageToTensor(),
			NormalizeIntensity()
			)
		tfm_mask = OneHot()

		items = Image(ImageCore.colorview(Gray, image)), MaskMulti(target .+ 1)
		item_im, item_mask = apply(tfm_proj, (items))

		return apply(tfm_im, item_im), apply(tfm_mask, item_mask)
	end

	function DLPipelines.encode(
			method::ImageSegmentationSimple,
			context::Validation,
			(image, target)::Union{Tuple, NamedTuple}
			)

		tfm_proj = CenterResizeCrop(method.imagesize)
		tfm_im = DataAugmentation.compose(
			ImageToTensor(),
			NormalizeIntensity()
			)
		tfm_mask = OneHot()

		items = Image(ImageCore.colorview(Gray, image)), MaskMulti(target .+ 1)
		item_im, item_mask = apply(tfm_proj, (items))

		return apply(tfm_im, item_im), apply(tfm_mask, item_mask)
	end
end

# ╔═╡ c85b7bf9-ab6e-45bf-b3e1-a029d715f4d5
md"""
## Create `methoddataloaders`
"""

# ╔═╡ 7ea4e1fa-cbb4-42ee-9c6b-ee801c35f2f2
batchsize = 2

# ╔═╡ 8dc5c540-4545-40d5-9999-2a4d3a853f5c
traindl, validdl = methoddataloaders(data, method, batchsize)

# ╔═╡ f94431fa-89c3-4030-807e-cc61a9977e81
md"""
## Visualize
"""

# ╔═╡ 85de2def-01ee-4fff-bfac-73c8d30ff34d
traindl[1]

# ╔═╡ 92f78dea-2c53-4fe0-aca2-247dbded247a
md"""
## Create  model
"""

# ╔═╡ a458ad0a-b48d-4c31-9f1c-cbe8717cc901
begin
	# 3D layer utilities
	conv = (stride, in, out) -> Conv((3, 3, 3), in=>out, stride=stride, pad=SamePad())
	tran = (stride, in, out) -> ConvTranspose((3, 3, 3), in=>out, stride=stride, pad=SamePad())

	conv1 = (in, out) -> Chain(conv(1, in, out), BatchNorm(out), x -> leakyrelu.(x))
	conv2 = (in, out) -> Chain(conv(2, in, out), BatchNorm(out), x -> leakyrelu.(x))
	tran2 = (in, out) -> Chain(tran(2, in, out), BatchNorm(out), x -> leakyrelu.(x))
end

# ╔═╡ c082b7e4-71d9-4a27-a302-18e543492f3d
begin
	function unet3D(in_chs, lbl_chs)
		# Contracting layers
		l1 = Chain(conv1(in_chs, 4))
		l2 = Chain(l1, conv1(4, 4), conv2(4, 16))
		l3 = Chain(l2, conv1(16, 16), conv2(16, 32))
		l4 = Chain(l3, conv1(32, 32), conv2(32, 64))
		l5 = Chain(l4, conv1(64, 64), conv2(64, 128))

		# Expanding layers
		l6 = Chain(l5, tran2(128, 64), conv1(64, 64))
		l7 = Chain(Parallel(+, l6, l4), tran2(64, 32), conv1(32, 32))
		l8 = Chain(Parallel(+, l7, l3), tran2(32, 16), conv1(16, 16))
		l9 = Chain(Parallel(+, l8, l2), tran2(16, 4), conv1(4, 4))
		l10 = Chain(l9, conv1(4, lbl_chs))
	end
end

# ╔═╡ ca6a29dc-55b8-4238-9ea8-6da9e0372158
md"""
## Create `Learner`
"""

# ╔═╡ 63a851e5-33bd-4756-aa49-d4f872020afc
begin
	model = unet3D(1, 2)
	optimizer = Flux.ADAM(0.01)
	loss_function = Flux.Losses.dice_coeff_loss
end

# ╔═╡ 83c1f0db-0526-407e-a620-3d913061da66
learner = Learner(model, (traindl, validdl), optimizer, loss_function)

# ╔═╡ 103ceb8d-995c-483f-863c-dd4e1f397f3f
md"""
## Train
"""

# ╔═╡ d60a4eee-b796-4f4b-bef8-99df74db4d3a
# fitonecycle!(learner, 1)

# ╔═╡ Cell order:
# ╠═514ab1cc-fb98-11eb-2a1e-c7b264f7066c
# ╠═5b2ed761-eb57-44c2-9ed9-bc1e5240af16
# ╟─3cfef02c-00f4-48e9-90b1-4ad43472f397
# ╠═92776990-345c-4abc-82ad-3cec0bee356f
# ╠═bf88491c-a431-4575-a88e-6cc80afce35e
# ╠═ff620c1f-e08f-4e8d-b538-31a48cea2e0b
# ╠═ae8852a9-5b80-4bda-b58f-3ff3d927bfb4
# ╟─2fd15108-ced5-4dc2-abb4-ceb7f394e718
# ╠═d2b61982-7829-486e-a0d9-c02fb6d77e74
# ╠═69278a55-81b9-4936-9146-f2510774ff22
# ╠═55f8c444-3486-4cd8-a1af-a3774ad55594
# ╟─c981defe-f337-47ac-8e7b-ea646b4cf482
# ╠═65570158-ee5a-437d-937f-e37f47c07823
# ╟─c85b7bf9-ab6e-45bf-b3e1-a029d715f4d5
# ╠═7ea4e1fa-cbb4-42ee-9c6b-ee801c35f2f2
# ╠═8dc5c540-4545-40d5-9999-2a4d3a853f5c
# ╟─f94431fa-89c3-4030-807e-cc61a9977e81
# ╠═85de2def-01ee-4fff-bfac-73c8d30ff34d
# ╟─92f78dea-2c53-4fe0-aca2-247dbded247a
# ╠═a458ad0a-b48d-4c31-9f1c-cbe8717cc901
# ╠═c082b7e4-71d9-4a27-a302-18e543492f3d
# ╟─ca6a29dc-55b8-4238-9ea8-6da9e0372158
# ╠═63a851e5-33bd-4756-aa49-d4f872020afc
# ╠═83c1f0db-0526-407e-a620-3d913061da66
# ╟─103ceb8d-995c-483f-863c-dd4e1f397f3f
# ╠═d60a4eee-b796-4f4b-bef8-99df74db4d3a


text:
module _Simple

const guard_digits = 5
const eps_digits = 5

pi(digits::Integer) = begin
    o = big(10)^(digits + guard_digits + eps_digits)
    osqr = o^2

    # a = 1
    # b = 1 / sqrt(2)
    a = o
    b = osqr ÷ isqrt(osqr << 1)

    s = big(0)
    d = a - b

    n = 0
    ϵ = 10^eps_digits

    while d >= ϵ
        # s += 2^n * d^2
        s += (d^2) << n
        # a, b = (a + b) / 2, sqrt(a * b)
        a, b = (a + b) >> 1, isqrt(a * b)
        d = a - b
        n += 1
    end

    # 4*a^2 / (1 - s)
    (4*a*a*o) ÷ ((osqr - s) * big(10)^(guard_digits + eps_digits))
end

end


text:
# Julia wrapper for header: /Volumes/BOOTCAMP/programs/compa_libs/iup/include/iup_plot.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

@static Sys.iswindows() ? (const iup_plot = "iup_plot") : (const iup_plot = "libiup_plot")  # Name of IUP shared lib.

function IupPlotOpen()
  ccall( (:IupPlotOpen, iup_plot), Cvoid, (), )
end

function IupPlot()
    ccall((:IupPlot, iup_plot), Ptr{Ihandle}, (),)
end

function IupPlotBegin(ih::Ptr{Ihandle}, strXdata::Integer)
  ccall( (:IupPlotBegin, iup_plot), Cvoid, (Ptr{Ihandle}, Cint), ih, strXdata)
end

function IupPlotAdd(ih::Ptr{Ihandle}, x::Cdouble, y::Cdouble)
  ccall( (:IupPlotAdd, iup_plot), Cvoid, (Ptr{Ihandle}, Cdouble, Cdouble), ih, x, y)
end

function IupPlotAddStr(ih::Ptr{Ihandle}, x::String, y::Cdouble)
  ccall( (:IupPlotAddStr, iup_plot), Cvoid, (Ptr{Ihandle}, Ptr{UInt8}, Cdouble), ih, x, y)
end

function IupPlotEnd(ih::Ptr{Ihandle})
  ccall( (:IupPlotEnd, iup_plot), Cint, (Ptr{Ihandle},), ih)
end

function IupPlotLoadData(ih::Ptr{Ihandle}, filename::Ptr{UInt8}, strXdata::Integer)
  ccall( (:IupPlotLoadData, iup_plot), Cint, (Ptr{Ihandle}, Ptr{UInt8}, Cint), ih, filename, strXdata)
end

function IupPlotInsertStr(ih::Ptr{Ihandle}, ds_index::Integer, sample_index::Integer, x::Ptr{UInt8}, y::Cdouble)
  ccall( (:IupPlotInsertStr, iup_plot), Cvoid, (Ptr{Ihandle}, Cint, Cint, Ptr{UInt8}, Cdouble), ih, ds_index, sample_index, x, y)
end

function IupPlotInsert(ih::Ptr{Ihandle}, ds_index::Integer, sample_index::Integer, x::Cdouble, y::Cdouble)
  ccall( (:IupPlotInsert, iup_plot), Cvoid, (Ptr{Ihandle}, Cint, Cint, Cdouble, Cdouble), ih, ds_index, sample_index, x, y)
end

function IupPlotInsertStrPoints(ih::Ptr{Ihandle}, ds_index::Integer, sample_index::Integer, x::Ptr{Ptr{UInt8}}, y::Ptr{Cdouble}, count::Integer)
  ccall( (:IupPlotInsertStrPoints, iup_plot), Cvoid, (Ptr{Ihandle}, Cint, Cint, Ptr{Ptr{UInt8}}, Ptr{Cdouble}, Cint), ih, ds_index, sample_index, x, y, count)
end

function IupPlotInsertPoints(ih::Ptr{Ihandle}, ds_index::Integer, sample_index::Integer, x::Ptr{Cdouble}, y::Ptr{Cdouble}, count::Integer)
  ccall( (:IupPlotInsertPoints, iup_plot), Cvoid, (Ptr{Ihandle}, Cint, Cint, Ptr{Cdouble}, Ptr{Cdouble}, Cint), ih, ds_index, sample_index, x, y, count)
end

function IupPlotAddPoints(ih::Ptr{Ihandle}, ds_index::Integer, x::Ptr{Cdouble}, y::Ptr{Cdouble}, count::Integer)
  ccall( (:IupPlotAddPoints, iup_plot), Cvoid, (Ptr{Ihandle}, Cint, Ptr{Cdouble}, Ptr{Cdouble}, Cint), ih, ds_index, x, y, count)
end

function IupPlotAddStrPoints(ih::Ptr{Ihandle}, ds_index::Integer, x::Ptr{Ptr{UInt8}}, y::Ptr{Cdouble}, count::Integer)
  ccall( (:IupPlotAddStrPoints, iup_plot), Cvoid, (Ptr{Ihandle}, Cint, Ptr{Ptr{UInt8}}, Ptr{Cdouble}, Cint), ih, ds_index, x, y, count)
end

function IupPlotGetSample(ih::Ptr{Ihandle}, ds_index::Integer, sample_index::Integer, x::Ptr{Cdouble}, y::Ptr{Cdouble})
  ccall( (:IupPlotGetSample, iup_plot), Cvoid, (Ptr{Ihandle}, Cint, Cint, Ptr{Cdouble}, Ptr{Cdouble}), ih, ds_index, sample_index, x, y)
end

function IupPlotGetSampleStr(ih::Ptr{Ihandle}, ds_index::Integer, sample_index::Integer, x::Ptr{Ptr{UInt8}}, y::Ptr{Cdouble})
  ccall( (:IupPlotGetSampleStr, iup_plot), Cvoid, (Ptr{Ihandle}, Cint, Cint, Ptr{Ptr{UInt8}}, Ptr{Cdouble}), ih, ds_index, sample_index, x, y)
end

function IupPlotGetSampleSelection(ih::Ptr{Ihandle}, ds_index::Integer, sample_index::Integer)
  ccall( (:IupPlotGetSampleSelection, iup_plot), Cint, (Ptr{Ihandle}, Cint, Cint), ih, ds_index, sample_index)
end

function IupPlotSetSample(ih::Ptr{Ihandle}, ds_index::Integer, sample_index::Integer, x::Cdouble, y::Cdouble)
  ccall( (:IupPlotSetSample, iup_plot), Cvoid, (Ptr{Ihandle}, Cint, Cint, Cdouble, Cdouble), ih, ds_index, sample_index, x, y)
end

function IupPlotSetSampleStr(ih::Ptr{Ihandle}, ds_index::Integer, sample_index::Integer, x::Ptr{UInt8}, y::Cdouble)
  ccall( (:IupPlotSetSampleStr, iup_plot), Cvoid, (Ptr{Ihandle}, Cint, Cint, Ptr{UInt8}, Cdouble), ih, ds_index, sample_index, x, y)
end

function IupPlotSetSampleSelection(ih::Ptr{Ihandle}, ds_index::Integer, sample_index::Integer, selected::Integer)
  ccall( (:IupPlotSetSampleSelection, iup_plot), Cvoid, (Ptr{Ihandle}, Cint, Cint, Cint), ih, ds_index, sample_index, selected)
end

function IupPlotTransform(ih::Ptr{Ihandle}, x::Cdouble, y::Cdouble, cnv_x::Ptr{Cdouble}, cnv_y::Ptr{Cdouble})
    ccall( (:IupPlotTransform, iup_plot), Cvoid, (Ptr{Ihandle}, Cdouble, Cdouble, Ptr{Cdouble}, Ptr{Cdouble}), ih, x, y, cnv_x, cnv_y)
end

function IupPlotTransformTo(ih::Ptr{Ihandle}, cnv_x::Cdouble, cnv_y::Cdouble, x::Ptr{Cdouble}, y::Ptr{Cdouble})
    ccall( (:IupPlotTransformTo, iup_plot), Cvoid, (Ptr{Ihandle}, Cdouble, Cdouble, Ptr{Cdouble}, Ptr{Cdouble}), ih, cnv_x, cnv_y, x, y)
end
# COMMENTED because cdCanvas is not known here
#function IupPlotPaintTo(ih::Ptr{Ihandle}, cnv::Ptr{cdCanvas})
#  ccall( (:IupPlotPaintTo, iup_plot), Cvoid, (Ptr{Ihandle}, Ptr{cdCanvas}), ih, cnv)
#end


text:
function initialize!(integrator,cache::RichardsonEulerCache)
  integrator.kshortsize = 2
  @unpack k,fsalfirst = cache
  integrator.fsalfirst = fsalfirst
  integrator.fsallast = k
  resize!(integrator.k, integrator.kshortsize)
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
  integrator.f(integrator.fsalfirst,integrator.uprev,integrator.p,integrator.t) # For the interpolation, needs k at the updated point
end

function perform_step!(integrator,cache::RichardsonEulerCache,repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack k,fsalfirst,T = cache
  @unpack m = cache.tab

  @muladd @. u = uprev + dt*fsalfirst
  T[1,1] = copy(u)
  halfdt = dt/2
  for i in 2:m
    # Solve using Euler method
    @muladd @. u = uprev + halfdt*fsalfirst
    f(k, u, p, t+halfdt)
    for j in 2:2^(i-1)
      @muladd @. u = u + halfdt*k
      f(k, u, p, t+j*halfdt)
    end
    T[i,1] = copy(u)
    # Richardson Extrapolation
    for j in 2:i
      T[i,j] = ((2^(j-1))*T[i,j-1] - T[i-1,j-1])/((2^(j-1)) - 1)
    end
    halfdt = halfdt/2
  end

  # using extrapolated value of u
  @. u = T[m,m]

  f(k, u, p, t+dt)
end

function initialize!(integrator,cache::RichardsonEulerConstantCache)
  integrator.kshortsize = 2
  integrator.k = typeof(integrator.k)(undef, integrator.kshortsize)
  integrator.fsalfirst = integrator.f(integrator.uprev, integrator.p, integrator.t) # Pre-start fsal

  # Avoid undefined entries if k is an array of arrays
  integrator.fsallast = zero(integrator.fsalfirst)
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
end

function perform_step!(integrator,cache::RichardsonEulerConstantCache,repeat_step=false)
  @unpack t,dt,uprev,f,p = integrator
  @unpack m = cache
  @muladd u = @. uprev + dt*integrator.fsalfirst
  T = zeros(typeof(uprev), (m,m))
  T[1,1] = u

  halfdt = dt/2
  for i in 2:m
    # Solve using Euler method
    @muladd u = @. uprev + halfdt*integrator.fsalfirst
    k = f(u, p, t+halfdt)

    for j in 2:2^(i-1)
      @muladd u = @. u + halfdt*k
      k = f(u, p, t+j*halfdt)
    end
    T[i,1] = u

    # Richardson Extrapolation
    for j in 2:i
      T[i,j] = ((2^(j-1))*T[i,j-1] - T[i-1,j-1])/((2^(j-1)) - 1)
    end
    halfdt = halfdt/2
  end

  k = f(T[m,m], p, t+dt)

  integrator.fsallast = k
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
  # Use extrapolated value of u
  integrator.u = T[m,m]
end


text:
include("../src/core.jl")
include("../src/integration.jl")

using Test

new = createIntegral

@testset "Validation" begin
    # Everything ok
    @test new(:(x-1), :x) ≠ nothing
    # Differential with no vars
    @test_throws ErrorException new(:(x-1), :(10 + 4*3))
    # Differential with multiple vars
    @test_throws ErrorException new(:(x^2), :(x*y))
end

@testset "Indefinite" begin
    calc(args...) = eval(new(args...))
    @test calc(:(sin(x)), :x) == :(-cos(x))
end

@testset "Definite" begin

end


text:
ia> [RpcApi.SimpleOrganism(0x000000000005fd85,:(function (o::Creature.Organism,)
     function func_92(var_90::Int8=-3,var_91::Bool=false)
         var_90 = var_90 + var_90
         return var_90
     end
     function func_78(var_76::Int8=114,var_77::Int8=15)
         local var_93::Int16 = -12977
         return var_76
     end
     function func_52(var_51::Int8=-36)
         local var_87::Int64 = 5892437190061028745
         Creature.stepDown(o)
         Creature.eatDown(o,Int(var_51))
         for i = 1:var_51
         end
         Creature.stepLeft(o)
         Creature.stepUp(o)
         Creature.eatUp(o,Int(var_51))
         var_87 = var_87 / var_87
         Creature.eatUp(o,Int(var_51))
         Creature.stepUp(o)
         Creature.eatUp(o,Int(var_51))
         Creature.eatUp(o,Int(var_51))
         Creature.stepRight(o)
         Creature.eatDown(o,Int(var_51))
         return var_51
     end
     local var_1::Int8 = 121
     var_1 = func_52(var_1)
     for i = 1:var_1
         local var_68::String = "1kR3kL1t"
         local var_15::Bool = true
         local var_8::Int16 = 23832
         local var_6::Int8 = 66
         var_15 = var_15 & var_15
         Creature.eatRight(o,Int(var_6))
         Creature.eatRight(o,Int(var_6))
         o.mem[var_8] = var_8
         var_6 = func_52(var_6)
         var_6 = func_52(var_6)
         var_8 = if haskey(o.mem,var_8)
                 o.mem[var_8]
             else
                 var_8
             end
         Creature.eatRight(o,Int(var_6))
         var_15 = var_15 | var_15
         Creature.stepUp(o)
         Creature.eatLeft(o,Int(var_6))
         Creature.stepLeft(o)
         var_8 = if haskey(o.mem,var_8)
                 o.mem[var_8]
             else
                 var_8
             end
         Creature.stepUp(o)
         Creature.stepDown(o)
         var_6 = func_52(var_6)
         Creature.eatDown(o,Int(var_6))
         var_15 = Bool(abs(var_15 - var_15))
     end
     Creature.eatDown(o,Int(var_1))
     Creature.stepLeft(o)
     Creature.stepRight(o)
     var_1 = var_1 / var_1
     Creature.eatUp(o,Int(var_1))
     Creature.eatDown(o,Int(var_1))
     Creature.stepDown(o)
     Creature.stepDown(o)
     var_1 = var_1 + var_1
     Creature.stepDown(o)
     var_1 = var_1 + var_1
     for i = 1:var_1
         local var_63::Bool = false
         Creature.stepDown(o)
         Creature.stepRight(o)
         Creature.stepDown(o)
         Creature.stepDown(o)
         Creature.stepRight(o)
     end
     return true
 end),[100,300,95,0,1,1,1],7,4654,8,4673824,0x0161fac0,Dict(-7040=>-7040,10405=>10405,-17872=>-17872,2=>2,23832=>23832,1=>1),[23,34],62),RpcApi.SimpleOrganism(0x000000000005fda3,:(function (o::Creature.Organism,)
     function func_78(var_76::Int8=114,var_77::Int8=15)
         return var_76
     end
     function func_52(var_51::Int8=-36)
         local var_90::Int16 = -3304
         local var_87::Int64 = 5892437190061028745
         Creature.stepDown(o)
         for i = 1:var_51
         end
         Creature.stepLeft(o)
         Creature.eatUp(o,Int(var_51))
         Creature.stepUp(o)
         Creature.eatUp(o,Int(var_51))
         Creature.eatUp(o,Int(var_51))
         Creature.eatUp(o,Int(var_51))
         Creature.eatUp(o,Int(var_51))
         Creature.eatUp(o,Int(var_51))
         Creature.eatDown(o,Int(var_51))
         Creature.stepRight(o)
         if var_87 === var_87
         end
         return var_51
     end
     local var_1::Int8 = 121
     for i = 1:var_1
     end
     for i = 1:var_1
         local var_68::String = "1kR3kL1t"
         local var_15::Bool = true
         local var_8::Int16 = 23832
         local var_6::Int8 = 66
         var_6 = func_52(var_6)
         var_6 = func_78(var_6,var_6)
         var_8 = if haskey(o.mem,var_8)
                 o.mem[var_8]
             else
                 var_8
             end
         Creature.eatUp(o,Int(var_6))
         var_6 = func_52(var_6)
         var_6 = func_52(var_6)
         Creature.eatLeft(o,Int(var_6))
         var_15 = var_15 | var_15
         Creature.stepUp(o)
         Creature.eatLeft(o,Int(var_6))
         Creature.stepLeft(o)
         var_68 = var_68[1:if length(var_68) > length(var_68)
                         0
                     else
                         length(var_68) - length(var_68)
                     end]
         var_6 = func_52(var_6)
         Creature.stepUp(o)
         Creature.stepDown(o)
         var_6 = func_52(var_6)
         var_8 = if haskey(o.mem,var_8)
                 o.mem[var_8]
             else
                 var_8
             end
         Creature.eatUp(o,Int(var_6))
         var_15 = Bool(abs(var_15 - var_15))
     end
     if var_1 < var_1
     end
     Creature.stepRight(o)
     Creature.stepRight(o)
     var_1 = var_1 / var_1
     Creature.eatUp(o,Int(var_1))
     Creature.stepDown(o)
     Creature.stepDown(o)
     var_1 = var_1 + var_1
     var_1 = var_1 + var_1
     for i = 1:var_1
         local var_63::Bool = false
         Creature.stepDown(o)
         Creature.stepRight(o)
         Creature.stepRight(o)
         Creature.stepDown(o)
         Creature.stepRight(o)
         Creature.stepDown(o)
         Creature.stepRight(o)
     end
     return true
 end),[100,300,95,0,1,1,1],7,6203,1,4208108,0x0163f128,Dict(1=>1,10405=>10405,-7040=>-7040,-17872=>-17872,2=>2,23832=>23832,-3304=>-3304),[20,15],61),RpcApi.SimpleOrganism(0x000000000005fd9b,:(function (o::Creature.Organism,)
     function func_52(var_51::Int8=-36)
         local var_87::Int64 = 5892437190061028745
         Creature.stepDown(o)
         Creature.eatRight(o,Int(var_51))
         if var_51 > var_51
         end
         Creature.stepLeft(o)
         Creature.stepUp(o)
         for i = 1:var_51
         end
         Creature.eatUp(o,Int(var_51))
         Creature.eatUp(o,Int(var_51))
         Creature.eatUp(o,Int(var_51))
         Creature.eatDown(o,Int(var_51))
         for i = 1:var_51
         end
         Creature.eatUp(o,Int(var_51))
         Creature.stepRight(o)
         return var_51
     end
     local var_1::Int8 = 121
     var_1 = func_52(var_1)
     for i = 1:var_1
     end
     for i = 1:var_1
         local var_68::String = "1kR3kL1t"
         local var_15::Bool = true
         local var_8::Int16 = 23832
         local var_6::Int8 = 66
         var_15 = var_15 & var_15
         Creature.eatRight(o,Int(var_6))
         Creature.eatRight(o,Int(var_6))
         Creature.eatUp(o,Int(var_6))
         var_6 = func_52(var_6)
         var_15 = Bool(abs(var_15 - var_15))
         var_8 = if haskey(o.mem,var_8)
                 o.mem[var_8]
             else
                 var_8
             end
         Creature.eatRight(o,Int(var_6))
         var_15 = var_15 | var_15
         Creature.stepUp(o)
         Creature.eatLeft(o,Int(var_6))
         Creature.stepLeft(o)
         Creature.eatRight(o,Int(var_6))
         Creature.stepUp(o)
         Creature.stepDown(o)
         var_8 = if haskey(o.mem,var_8)
                 o.mem[var_8]
             else
                 var_8
             end
         var_6 = func_52(var_6)
         Creature.eatDown(o,Int(var_6))
         var_15 = Bool(abs(var_15 - var_15))
     end
     Creature.eatDown(o,Int(var_1))
     var_1 = var_1 / var_1
     var_1 = var_1 / var_1
     Creature.eatUp(o,Int(var_1))
     Creature.eatDown(o,Int(var_1))
     Creature.stepRight(o)
     Creature.eatUp(o,Int(var_1))
     var_1 = var_1 + var_1
     var_1 = var_1 + var_1
     for i = 1:var_1
         local var_63::Bool = false
         Creature.stepDown(o)
         Creature.stepRight(o)
         Creature.stepRight(o)
         Creature.stepDown(o)
         Creature.stepDown(o)
         Creature.stepDown(o)
         Creature.stepRight(o)
     end
     return true
 end),[100,300,95,0,1,1,1],7,4859,1,4045603,0x0160f4a0,Dict(-7040=>-7040,10405=>10405,-17872=>-17872,2=>2,23832=>23832,1=>1),[15,18],60),RpcApi.SimpleOrganism(0x000000000005fcc0,:(function (o::Creature.Organism,)
     local var_87::String = "o7IXGFtJ"
     local var_85::String = "uYWnLeWF"
     function func_50(var_48::Int64=8246173392054715175,var_49::Bool=false)
         return var_48
     end
     function func_19(var_17::String="A06aRuJT",var_18::Int64=-4494005246402814635)
         var_17 = var_17[1:if length(var_17) > length(var_17) > 0
                         0
                     else
                         div(length(var_17),length(var_17))
                     end]
         return var_17
     end
     local var_1::Int8 = 121
     Creature.eatDown(o,Int(var_1))
     Creature.stepUp(o)
     Creature.stepDown(o)
     Creature.eatUp(o,Int(var_1))
     for i = 1:var_1
         local var_79::Int64 = 8722268545643689674
         local var_52::String = "3Zf0ESw4"
         local var_37::Int8 = 66
         local var_9::Bool = true
         local var_6::Int8 = 66
         Creature.eatUp(o,Int(var_37))
         Creature.eatUp(o,Int(var_6))
         Creature.stepLeft(o)
         Creature.stepDown(o)
         Creature.eatRight(o,Int(var_6))
         var_52 = var_52[if length(var_52) > length(var_52)
                         1
                     else
                         length(var_52)
                     end:if length(var_52) > length(var_52) > 0
                         0
                     else
                         end
                     end]
         Creature.eatRight(o,Int(var_6))
         Creature.eatLeft(o,Int(var_6))
         Creature.stepLeft(o)
         Creature.eatLeft(o,Int(var_6))
         Creature.eatDown(o,Int(var_6))
         Creature.stepLeft(o)
         Creature.eatLeft(o,Int(var_37))
         Creature.eatRight(o,Int(var_37))
         Creature.stepDown(o)
         Creature.eatDown(o,Int(var_37))
         Creature.stepDown(o)
         Creature.eatDown(o,Int(var_6))
         var_79 = var_79 - var_79
         Creature.stepDown(o)
         Creature.stepLeft(o)
         Creature.eatUp(o,Int(var_6))
     end
     for i = 1:var_1
         local var_88::Int16 = 31438
         local var_86::Int8 = -69
         Creature.stepRight(o)
         Creature.stepRight(o)
         Creature.stepRight(o)
         Creature.stepUp(o)
         Creature.stepLeft(o)
         var_86 = var_86 / var_86
         var_88 = var_88 * var_88
         Creature.stepRight(o)
         Creature.eatUp(o,Int(var_86))
     end
     Creature.stepDown(o)
     for i = 1:var_1
         Creature.stepLeft(o)
     end
     Creature.eatLeft(o,Int(var_1))
     for i = 1:var_1
         local var_75::Int8 = 63
         Creature.eatLeft(o,Int(var_75))
         Creature.eatLeft(o,Int(var_75))
     end
     Creature.stepDown(o)
     Creature.stepLeft(o)
     return true
 end),[100,300,95,0,1,1,1],13,5410,13,4038340,0x01678698,Dict(7511=>7511,23832=>23832,-24995=>-24995,1=>1),[249,216],60),RpcApi.SimpleOrganism(0x000000000005fcea,:(function (o::Creature.Organism,)
     local var_85::String = "uYWnLeWF"
     function func_50(var_48::Int64=8246173392054715175,var_49::Bool=false)
         return var_48
     end
     function func_19(var_17::String="A06aRuJT",var_18::Int64=-4494005246402814635)
         return var_17
     end
     local var_1::Int8 = 121
     Creature.eatRight(o,Int(var_1))
     Creature.stepDown(o)
     Creature.eatUp(o,Int(var_1))
     var_85 = var_85[1:if length(var_85) > length(var_85)
                     0
                 else
                     length(var_85) - length(var_85)
                 end]
     for i = 1:var_1
         local var_79::Int64 = 8722268545643689674
         local var_52::String = "3Zf0ESw4"
         local var_37::Int8 = 66
         local var_9::Bool = true
         local var_8::Int16 = 23832
         local var_6::Int8 = 66
         Creature.eatUp(o,Int(var_37))
         Creature.eatUp(o,Int(var_6))
         Creature.stepLeft(o)
         Creature.stepDown(o)
         Creature.eatRight(o,Int(var_6))
         Creature.eatLeft(o,Int(var_6))
         Creature.eatRight(o,Int(var_6))
         Creature.eatLeft(o,Int(var_6))
         Creature.stepLeft(o)
         Creature.eatLeft(o,Int(var_6))
         Creature.eatDown(o,Int(var_6))
         Creature.stepLeft(o)
         Creature.eatLeft(o,Int(var_37))
         var_8 = if haskey(o.mem,var_8)
                 o.mem[var_8]
             else
                 var_8
             end
         Creature.eatRight(o,Int(var_37))
         Creature.stepDown(o)
         Creature.eatDown(o,Int(var_37))
         Creature.stepDown(o)
         Creature.eatRight(o,Int(var_6))
         var_37 = var_6 / var_37
         var_8 = if haskey(o.mem,var_8)
                 o.mem[var_8]
             else
                 var_8
             end
         Creature.stepLeft(o)
         Creature.eatRight(o,Int(var_37))
         var_8 = var_8 + var_8
         Creature.eatUp(o,Int(var_6))
     end
     for i = 1:var_1
         local var_86::Int8 = -69
         Creature.stepRight(o)
         Creature.stepRight(o)
         Creature.stepUp(o)
         Creature.stepRight(o)
         Creature.stepRight(o)
     end
     for i = 1:var_1
         local var_90::String = "fElvFwMz"
         Creature.stepLeft(o)
     end
     Creature.eatLeft(o,Int(var_1))
     for i = 1:var_1
         local var_75::Int8 = 63
         Creature.eatLeft(o,Int(var_75))
         Creature.eatRight(o,Int(var_75))
         Creature.eatLeft(o,Int(var_75))
     end
     Creature.stepDown(o)
     Creature.eatUp(o,Int(var_1))
     return true
 end),[100,300,95,0,1,1,1],13,5410,13,3824216,0x0167b248,Dict(7511=>7511,23832=>23832,-24995=>-24995,1=>1),[430,216],58),RpcApi.SimpleOrganism(0x000000000005fc82,:(function (o::Creature.Organism,)
     function func_78(var_76::Int8=114,var_77::Int8=15)
         return var_76
     end
     function func_52(var_51::Int8=-36)
         local var_90::Int16 = -3304
         local var_87::Int64 = 5892437190061028745
         Creature.stepDown(o)
         for i = 1:var_51
             Creature.stepRight(o)
         end
         Creature.stepLeft(o)
         Creature.eatUp(o,Int(var_51))
         Creature.stepUp(o)
         Creature.eatUp(o,Int(var_51))
         Creature.eatUp(o,Int(var_51))
         Creature.eatUp(o,Int(var_51))
         o.mem[var_90] = var_90
         Creature.eatDown(o,Int(var_51))
         Creature.stepRight(o)
         if var_87 === var_87
         end
         return var_51
     end
     local var_1::Int8 = 121
     for i = 1:var_1
     end
     for i = 1:var_1
         local var_68::String = "1kR3kL1t"
         local var_15::Bool = true
         local var_8::Int16 = 23832
         local var_6::Int8 = 66
         var_6 = func_52(var_6)
         var_6 = func_78(var_6,var_6)
         var_8 = if haskey(o.mem,var_8)
                 o.mem[var_8]
             else
                 var_8
             end
         Creature.eatUp(o,Int(var_6))
         var_6 = func_52(var_6)
         var_6 = func_52(var_6)
         Creature.eatLeft(o,Int(var_6))
         var_15 = var_15 | var_15
         var_8 = var_8 - var_8
         Creature.eatLeft(o,Int(var_6))
         Creature.stepLeft(o)
         var_68 = var_68[1:if length(var_68) > length(var_68)
                         0
                     else
                         length(var_68) - length(var_68)
                     end]
         var_6 = func_52(var_6)
         Creature.stepUp(o)
         Creature.stepDown(o)
         var_6 = func_52(var_6)
         var_8 = if haskey(o.mem,var_8)
                 o.mem[var_8]
             else
                 var_8
             end
         Creature.eatUp(o,Int(var_6))
         var_15 = Bool(abs(var_15 - var_15))
     end
     if var_1 < var_1
     end
     Creature.stepRight(o)
     Creature.stepRight(o)
     var_1 = var_1 / var_1
     Creature.eatUp(o,Int(var_1))
     Creature.stepDown(o)
     Creature.stepDown(o)
     var_1 = var_1 + var_1
     var_1 = var_1 + var_1
     for i = 1:var_1
         local var_63::Bool = false
         Creature.stepDown(o)
         Creature.stepRight(o)
         Creature.stepRight(o)
         Creature.stepDown(o)
         Creature.stepRight(o)
         Creature.stepDown(o)
         Creature.stepRight(o)
     end
     return true
 end),[100,300,95,0,1,1,1],7,6203,1,3624857,0x01640700,Dict(1=>1,10405=>10405,-7040=>-7040,-17872=>-17872,2=>2,23832=>23832,-3304=>-3304),[638,335],61),RpcApi.SimpleOrganism(0x000000000005fd2e,:(function (o::Creature.Organism,)
     local var_85::String = "uYWnLeWF"
     function func_50(var_48::Int64=8246173392054715175,var_49::Bool=false)
         local var_89::Int64 = -4049598015957677781
         return var_48
     end
     function func_19(var_17::String="A06aRuJT",var_18::Int64=-4494005246402814635)
         if var_18 != var_18
         end
         return var_17
     end
     local var_1::Int8 = 121
     Creature.eatRight(o,Int(var_1))
     Creature.stepDown(o)
     Creature.eatUp(o,Int(var_1))
     var_85 = var_85[1:if length(var_85) > length(var_85)
                     0
                 else
                     length(var_85) - length(var_85)
                 end]
     for i = 1:var_1
         local var_52::String = "3Zf0ESw4"
         local var_37::Int8 = 66
         local var_9::Bool = true
         local var_8::Int16 = 23832
         local var_6::Int8 = 66
         Creature.eatUp(o,Int(var_37))
         Creature.eatUp(o,Int(var_6))
         Creature.stepLeft(o)
         Creature.stepDown(o)
         Creature.eatRight(o,Int(var_6))
         Creature.eatLeft(o,Int(var_6))
         Creature.eatRight(o,Int(var_6))
         Creature.eatLeft(o,Int(var_6))
         Creature.stepLeft(o)
         Creature.eatLeft(o,Int(var_6))
         Creature.eatDown(o,Int(var_6))
         Creature.stepLeft(o)
         Creature.eatLeft(o,Int(var_37))
         var_8 = if haskey(o.mem,var_8)
                 o.mem[var_8]
             else
                 var_8
             end
         Creature.eatRight(o,Int(var_37))
         Creature.stepDown(o)
         Creature.eatDown(o,Int(var_37))
         Creature.stepDown(o)
         Creature.eatRight(o,Int(var_6))
         var_37 = var_6 / var_37
         var_8 = if haskey(o.mem,var_8)
                 o.mem[var_8]
             else
                 var_8
             end
         Creature.eatRight(o,Int(var_6))
         Creature.eatRight(o,Int(var_37))
         var_8 = var_8 + var_8
         Creature.eatUp(o,Int(var_6))
     end
     for i = 1:var_1
         local var_86::Int8 = -69
         Creature.stepRight(o)
         Creature.stepRight(o)
         Creature.stepUp(o)
         Creature.stepRight(o)
         Creature.stepRight(o)
     end
     for i = 1:var_1
         Creature.stepLeft(o)
     end
     Creature.eatLeft(o,Int(var_1))
     for i = 1:var_1
         local var_75::Int8 = 63
         Creature.eatLeft(o,Int(var_75))
         var_75 = var_75 + var_75
         Creature.eatLeft(o,Int(var_75))
     end
     Creature.stepDown(o)
     Creature.eatUp(o,Int(var_1))
     return true
 end),[100,300,95,0,1,1,1],13,5410,13,3217757,0x01681f80,Dict(7511=>7511,23832=>23832,-24995=>-24995,1=>1),[509,215],58)]

max_stars_repo_name: tmptrash/jevo

max_issues_repo_name: tmptrash/jevo

max_forks_repo_name: tmptrash/jevo


text:
include("./shared.jl")
include("../scripts-part2/2D_SWE/shallow_water_2D_LF_xpu_mpi.jl")
include("../scripts-part2/2D_SWE/shallow_water_2D_LF_xpu.jl")
include("../scripts-part2/2D_SWE/shallow_water_2D.jl")

@testset "Height H between baseline and 2D xpu/mpi implementation of dam break in y-direction" begin
      H_mpi, xc = shallow_water_2D_xpu_mpi(; nx = n, ny = n, dam1D_x = false)
	@reset_parallel_stencil()
	H_xpu = shallow_water_2D_xpu(; nx = n, ny = n, dam1D_x = false)
	@reset_parallel_stencil()
      H_cpu = shallow_water_2D(; nx = n, ny = n, dam_x = false)
      @test isapprox(H_xpu, H_cpu)
	@test isapprox(H_mpi, H_cpu)
end;


text:
using Distances

using Test
using LinearAlgebra
using Random
using Statistics

include("F64.jl")
include("test_dists.jl")


text:
@doc raw"""
Holds the tableau of a Runge-Kutta method
```math
\begin{aligned}
Q_{n,i} &= q_{n} + h \sum \limits_{j=1}^{s} a_{ij} \, v(t_{n} + c_j \Delta t, Q_{n,j}) , &
q_{n+1} &= q_{n} + h \sum \limits_{i=1}^{s} b_{i}  \, v(t_{n} + c_j \Delta t, Q_{n,i}) , \\
\end{aligned}
```

Parameters:
 * `T`: datatype of coefficient arrays

Fields:
 * `name`: symbolic name of the tableau
 * `o`: order of the method
 * `s`: number of stages
 * `a`: coefficients $a_{ij}$ with $ 1 \le i,j \le s$
 * `b`: weights $b_{i}$  with $ 1 \le i \le s$
 * `c`: nodes $c_{i}$  with $ 1 \le i \le s$

Constructors:
```julia
Tableau{T}(name, o, s, a, b, c)
Tableau{T}(name, o, a, b, c)
Tableau(name::Symbol, o::Int, s::Int, a::AbstractMatrix, b::AbstractVector, c::AbstractVector)
Tableau(name::Symbol, o::Int, a::AbstractMatrix, b::AbstractVector, c::AbstractVector)
Tableau(name::Symbol, o::Int, t::AbstractMatrix)
```

The last constructor accepts an $(s+1) \times (s+1)$ array that holds the whole tableau in the form
of a Butcher tableau, i.e.,

 c | a
---|---
   | b

"""
struct Tableau{T}
    name::Symbol
    o::Int
    s::Int

    a::Matrix{T}
    b::Vector{T}
    c::Vector{T}

    â::Matrix{T}
    b̂::Vector{T}
    ĉ::Vector{T}

    function Tableau{T}(name,o,s,a,b,c) where {T}
        @assert s > 0 "Number of stages must be > 0"
        @assert s == size(a,1) == size(a,2) == length(b) == length(c)

        ã = convert(Matrix{T}, a)
        b̃ = convert(Vector{T}, b)
        c̃ = convert(Vector{T}, c)

        â = similar(ã)
        b̂ = similar(b̃)
        ĉ = similar(c̃)

        â .= a .- ã
        b̂ .= b .- b̃
        ĉ .= c .- c̃

        new(name,o,s,ã,b̃,c̃,â,b̂,ĉ)
    end

    function Tableau{T}(name,o,a,b,c) where {T}
        Tableau{T}(name,o,length(c),a,b,c)
    end
end

Tableau(name::Symbol, o::Int, s::Int, a::AbstractMatrix{AT}, b::AbstractVector{BT}, c::AbstractVector{CT}) where {AT,BT,CT} = Tableau{promote_type(AT,BT,CT)}(name,o,s,a,b,c)
Tableau(name::Symbol, o::Int, a::AbstractMatrix, b::AbstractVector, c::AbstractVector) = Tableau(name,o,length(c),a,b,c)

function Tableau(name::Symbol, o::Int, t::AbstractMatrix{T}) where {T}
    @assert size(t,1) == size(t,2)

    local s = size(t,1)-1
    local a = copy(t[1:s, 2:s+1])
    local b = copy(t[s+1, 2:s+1])
    local c = copy(t[1:s, 1    ])

    Tableau{T}(name, o, s, a, b, c)
end

Base.hash(tab::Tableau, h::UInt) = hash(tab.o, hash(tab.s, hash(tab.a, hash(tab.b, hash(tab.c, hash(tab.â, hash(tab.b̂, hash(tab.ĉ, hash(:Tableau, h)))))))))

Base.:(==)(tab1::Tableau, tab2::Tableau) = (tab1.o == tab2.o
                                         && tab1.s == tab2.s
                                         && tab1.a == tab2.a
                                         && tab1.b == tab2.b
                                         && tab1.c == tab2.c
                                         && tab1.â == tab2.â
                                         && tab1.b̂ == tab2.b̂
                                         && tab1.ĉ == tab2.ĉ)

Base.isapprox(tab1::Tableau, tab2::Tableau; kwargs...) = (
                                            tab1.o == tab2.o
                                         && tab1.s == tab2.s
                                         && isapprox(tab1.a, tab2.a; kwargs...)
                                         && isapprox(tab1.b, tab2.b; kwargs...)
                                         && isapprox(tab1.c, tab2.c; kwargs...))

Base.isequal(tab1::Tableau{T1}, tab2::Tableau{T2}) where {T1,T2} = (tab1 == tab2 && T1 == T2 && tab1.name == tab2.name)

Base.eltype(::Tableau{T}) where {T} = T

function to_array(tab::Tableau{T}) where {T}
    local s = tab.s
    local arr = zeros(T, s+1, s+1)
    arr[1:s, 1    ] .= tab.c
    arr[s+1, 2:s+1] .= tab.b
    arr[1:s, 2:s+1] .= tab.a
    return arr
end

Base.convert(::Type{Tableau}, t::AbstractMatrix; name::Symbol=:nonamespecified, o::Int=0) = Tableau(name, o, t)
Base.convert(::Type{Matrix{T}}, tab::Tableau) where {T} = convert(Matrix{T}, to_array(tab))
Base.convert(::Type{Matrix}, tab::Tableau{T}) where {T} = convert(Matrix{T}, tab)
Base.convert(::Type{Array{T}}, tab::Tableau) where {T} = convert(Matrix{T}, tab)
Base.convert(::Type{Array}, tab::Tableau{T}) where {T} = convert(Matrix{T}, tab)

"""
```julia
from_file(dir::AbstractString, name::AbstractString)
```

Read Runge-Kutta tableau from the file `<name>.tsv` in the directory `dir`.
"""
function from_file(dir::AbstractString, name::AbstractString)
    file = string(dir, "/", name, ".tsv")

    # Reads and parse Tableau metadata from file
    f = open(file, "r")
    header = readline(f)
    close(f)

    if header[1] == '#'
        header = split(header[2:end])
    else
        header = ()
    end

    if length(header) ≥ 1
        o = Base.parse(Int, header[1])
    else
        o = 0
    end

    if length(header) ≥ 2
        s = Base.parse(Int, header[2])
    else
        s = 0
    end

    if length(header) ≥ 3
        T = Core.eval(Main, Meta.parse(header[3]))
    else
        T = Float64
    end

    # TODO Read data in original format (e.g., Rational).
    #      For this we need to save tableaus as jld or hdf5.
    # tab_array = readdlm(file, T)
    tab_array = readdlm(file, comments=true)

    if s == 0
        s = size(tab_array, 1) - 1
    end

    @assert s == size(tab_array, 1) - 1 == size(tab_array, 2) - 1

    @info("Reading Runge-Kutta tableau $(name) with $(s) stages and order $(o) from file\n$(file)")

    Tableau(Symbol(name), o, tab_array)
end

"""
```julia
to_file(dir::AbstractString, tab::Tableau)
```

Write Runge-Kutta tableau to the file `<tab.name>.tsv` in the directory `dir`.
"""
function to_file(dir::AbstractString, tab::Tableau{T}) where {T}
    # tab_array = zeros(T, S+1, S+1)
    # tab_array[1:S, 2:S+1] = tab.a
    # tab_array[S+1, 2:S+1] = tab.b
    # tab_array[1:S, 1] = tab.c
    # tab_array[S+1, 1] = tab.order

    tab_array = convert(Matrix, tab)
    header = string("# ", tab.o, " ", tab.s, " ", T, "\n")
    file   = string(dir, "/", tab.name, ".tsv")

    @info("Writing Runge-Kutta tableau $(tab.name) with $(tab.s) stages and order $(tab.o) to file\n$(file)")

    f = open(file, "w")
    write(f, header)
    writedlm(f, float(tab_array))
    close(f)

    # TODO Write data in original format (e.g., Rational).
end

name(tab::Tableau) = tab.name
order(tab::Tableau) = tab.o
nstages(tab::Tableau) = tab.s
eachstage(tab::Tableau) = 1:tab.s
coefficients(tab::Tableau) = tab.a
weights(tab::Tableau) = tab.b
nodes(tab::Tableau) = tab.c

isexplicit(tab::Tableau) = istrilstrict(tab.a) && tab.c[1] == 0
isimplicit(tab::Tableau) = !isexplicit(tab)
isdiagnonallyimplicit(tab::Tableau) = tab.s != 1 && !istrilstrict(tab.a) && istril(tab.a)
isfullyimplicit(tab::Tableau) = (tab.s == 1 && tab.a[1,1] != 0) || (!istrilstrict(tab.a) && !istril(tab.a))

const tf_butcher_tableau = TextFormat(
    up_right_corner     = ' ',
    up_left_corner      = ' ',
    bottom_left_corner  = ' ',
    bottom_right_corner = ' ',
    up_intersection     = ' ',
    left_intersection   = ' ',
    right_intersection  = ' ',
    middle_intersection = ' ',
    bottom_intersection = ' ',
    column              = '│',
    row                 = ' '
)

"""
```julia
Base.show(io::IO, tab::Tableau)
```

Pretty-print Runge-Kutta tableau.
"""
function Base.show(io::IO, tab::Tableau)
    print(io, "\nRunge-Kutta Tableau $(tab.name) with $(tab.s) stages and order $(tab.o):\n")
    arr = convert(Matrix{Any}, tab)
    arr[tab.s+1,1] = ""
    pretty_table(io, arr,
                    tf = tf_butcher_tableau,
                    vlines = [1],
                    body_hlines = [tab.s],
                    body_hlines_format = ('─','┼','─','─'),
                    equal_columns_width = true,
                    noheader = true)
end

"""
```julia
Base.show(io::IO, ::MIME"text/markdown", tab::Tableau)
```

Generate and print a nice markdown table for the Runge-Kutta tableau.
"""
function Base.show(io::IO, ::MIME"text/markdown", tab::Tableau)
    show(io, "text/markdown", Markdown.parse("Runge-Kutta Tableau $(tab.name) with $(tab.s) stages and order $(tab.o):"))

    tab_arr = convert(Matrix{Any}, tab)
    tab_arr[tab.s+1,1] = ""

    strio = IOBuffer()
    pretty_table(strio, tab_arr,
                    backend = :latex,
                    vlines = [1],
                    hlines = [tab.s],
                    noheader = true)
    tab_latex = String(take!(strio))

    tab_markdown = replace(tab_latex, "tabular" => "array")
    tab_markdown = replace(tab_markdown, "\\begin{table}" => "```math")
    tab_markdown = replace(tab_markdown, "\\end{table}" => "```")

    print(io, tab_markdown)
end


text:
# Bio.Intervals
# =============
#
# Module for genomic intervals.
#
# This file is a part of BioJulia.
# License is MIT: https://github.com/BioJulia/Bio.jl/blob/master/LICENSE.md

module Intervals

export
    Strand,
    Interval,
    leftposition,
    hasleftposition,
    rightposition,
    hasrightposition,
    IntervalCollection,
    IntervalStream,
    STRAND_NA,
    STRAND_POS,
    STRAND_NEG,
    STRAND_BOTH,
    coverage,
    isoverlapping,
    seqname,
    hasseqname,
    metadata,
    strand,
    BEDReader,
    BEDWriter,
    BEDMetadata,
    BEDInterval,
    BigBedReader,
    BigBedWriter,
    GFF3,
    GFF3Reader,
    GFF3Metadata,
    GFF3Interval,
    hasfasta,
    getfasta,
    directives,
    Tabix,
    overlapchunks,
    tryread!,
    isfilled,
    MissingFieldException

import Base.Collections: heappush!, heappop!
import Bio.Exceptions: MissingFieldException
import Bio.Ragel: Ragel, tryread!
import Bio.StringFields: StringField
import BufferedStreams: BufferedStreams, BufferedInputStream, upanchor!
import ColorTypes: RGB
import IntervalTrees
import Libz
importall Bio

include("strand.jl")
include("interval.jl")
include("stream_buffer.jl")
include("intervalcollection.jl")
include("intervalstream.jl")
include("index/index.jl")

# Parsing file types
include("bed/bed.jl")
include("bbi/bbi.jl")
include("gff3/gff3.jl")

include("deprecated.jl")

end # module Intervals


text:
mutable struct Note
    interval::Interval
    beats::Int
end

function Note(; interval = Interval(), beats = 1)
    # convert strings and rationals to intervals first
    Note(interval, beats)
end

const NOTE_REGEX = r"(?<interval>.*) for (?<beats>.*)"

function parse(::Type{Note}, text::AbstractString; line_number = line_number)
    a_match = match(NOTE_REGEX, text)
    if a_match === nothing
        throw_parse_error(text, "note", line_number)
    else
        beats_string = a_match["beats"]
        beats = tryparse(Int, beats_string)
        if beats === nothing
            throw_parse_error(beats_string, "beats", line_number)
        else
            Note(parse(Interval, a_match["interval"]; line_number = line_number), beats)
        end
    end
end

function print(io::IO, note::Note)
    print(io, note.interval)
    print(io, " for ")
    print(io, note.beats)
end

# TODO: propertynames?

@inline function getproperty(note::Note, property_name::Symbol)
    if property_name === :numerator ||
       property_name === :denominator ||
       property_name === :octave
        getproperty(note.interval, property_name)
    else
        getfield(note, property_name)
    end
end

@inline function setproperty!(note::Note, property_name::Symbol, value)
    if property_name === :numerator ||
       property_name === :denominator ||
       property_name === :octave
        setproperty!(note.interval, property_name, value)
    else
        setfield!(
            note,
            property_name,
            convert(fieldtype(typeof(note), property_name), value),
        )
    end
end

max_stars_repo_name: bramtayl/Justly

max_issues_repo_name: bramtayl/Justly

max_forks_repo_name: bramtayl/Justly


text:
"""
KLDSampleSize(k::Int, ζ, η::Float64 = 0.05)

Return the minimum sample size in order to achieve an error at most ζ with a 1-η level of confidence according to KLD-Sampling.
"""
function KLDSampleSize(k::Int, ζ, η::Float64 = 0.05)
    a = (max(2.0, k) - 1.0) / 2.0
    b = 1.0/(a*9.0)
    return (1.0-b+sqrt(b)*quantile(Normal(), 1-η))^3.0*a/ζ
end


text:
module ServiceTests
include("test_daemon.jl")
include("test_filesystemcache.jl")
include("test_bucketcachedatasource.jl")
include("test_clibucket.jl")
end # module ServiceTests

max_forks_repo_name: seung-lab/SimpleTasks


text:
@testset "CRS" begin

    epsg_nl_code = 28992
    epsg_nl_str = "EPSG:28992"
    epsg_nl_wkt = "PROJCS[\"Amersfoort / RD New\",GEOGCS[\"Amersfoort\",DATUM[\"Amersfoort\",SPHEROID[\"Bessel 1841\",6377397.155,299.1528128,AUTHORITY[\"EPSG\",\"7004\"]],TOWGS84[565.2369,50.0087,465.658,-0.406857,0.350733,-1.87035,4.0812],AUTHORITY[\"EPSG\",\"6289\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"4289\"]],PROJECTION[\"Oblique_Stereographic\"],PARAMETER[\"latitude_of_origin\",52.15616055555555],PARAMETER[\"central_meridian\",5.38763888888889],PARAMETER[\"scale_factor\",0.9999079],PARAMETER[\"false_easting\",155000],PARAMETER[\"false_northing\",463000],UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]],AXIS[\"X\",EAST],AXIS[\"Y\",NORTH],AUTHORITY[\"EPSG\",\"28992\"]]"

    merc_proj_string = "+proj=merc +lat_ts=56.5 +ellps=GRS80"
    merc_wkt = "PROJCS[\"unnamed\",GEOGCS[\"GRS 1980(IUGG, 1980)\",DATUM[\"unknown\",SPHEROID[\"GRS80\",6378137,298.257222101]],PRIMEM[\"Greenwich\",0],UNIT[\"degree\",0.0174532925199433]],PROJECTION[\"Mercator_2SP\"],PARAMETER[\"standard_parallel_1\",56.5],PARAMETER[\"central_meridian\",0],PARAMETER[\"false_easting\",0],PARAMETER[\"false_northing\",0]]"

    @testset "Set CRS on GeoArray" begin
        ga = GeoArray(rand(10,10))
        epsg!(ga, 28992)
        @test ga.crs == epsg_nl_wkt

        epsg!(ga, "EPSG:28992")
        @test ga.crs == epsg_nl_wkt
    end

    @testset "Proj string" begin
        result = GeoArrays.str2wkt(merc_proj_string)
        @test result == merc_wkt
    end

    @testset "EPSG code" begin
        result = GeoArrays.epsg2wkt(epsg_nl_code)
        @test result == epsg_nl_wkt
    end

    @testset "Projection string checking" begin
        # EPSG string handling
        result = GeoArrays.str2wkt(epsg_nl_str)
        @test result == epsg_nl_wkt

        # Proj string handling
        result = GeoArrays.str2wkt(merc_proj_string)
        @test result == merc_wkt

        # WKT string handling
        result = GeoArrays.str2wkt(merc_wkt)
        @test result == merc_wkt

        # Invalid string
        @test_throws ArgumentError GeoArrays.str2wkt("INVALID")

    end

    @testset "Check projection of file" begin
        ga = GeoArrays.read("gdalworkshop/world.tif")
        @test ga.crs == "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0],UNIT[\"degree\",0.0174532925199433],AUTHORITY[\"EPSG\",\"4326\"]]"
    end

end


text:
function simulation(model::Dict{String,Any}, trainingdata::DataFrame;
    speciesindex::Int64 = 9)::Tuple{Array{Float64,2},Array{Float64,2}}

    # define contents -
    SF = 1e9 # all concentrations are in nmol/L

    # get size of the training data set -
    (ℛ, 𝒞) = size(trainingdata)

    # setup -
    species_archive = Array{Array{Float64,1},1}(undef, ℛ)
    T_archive = Array{Array{Float64,1},1}(undef, ℛ)

    # main loop -
    for i ∈ 1:ℛ

        # build new model -
        dd = deepcopy(model)

        # setup static -
        sfa = dd["static_factors_array"]
        sfa[1] = (2.5e-9) * SF        # 1 TFPI
        sfa[2] = (3.4e-6) * SF        # 2 AT
        sfa[3] = (5e-12) * SF         # 3 TF
        sfa[6] = 0.005                # 6 TRAUMA

        # grab the multiplier from the data -
        ℳ = dd["number_of_dynamic_states"]
        xₒ = zeros(ℳ)
        xₒ[1] = (1.4e-6) * (trainingdata[i, :II] / 100.0)     # 1 FII
        xₒ[2] = (1e-8) * (trainingdata[i, :VII] / 100.0)      # 2 FVII
        xₒ[3] = (2e-8) * (trainingdata[i, :V] / 100.0)        # 3 FV
        xₒ[4] = (1.6e-7) * (trainingdata[i, :X] / 100.0)      # 4 FX
        xₒ[5] = (7e-10) * (trainingdata[i, :VIII] / 100.0)    # 5 FVIII
        xₒ[6] = (9e-8) * (trainingdata[i, :IX] / 100.0)       # 6 FIX
        xₒ[7] = (1e-8) * (trainingdata[i, :XI] / 100.0)       # 7 FXI
        xₒ[8] = (1e-8) * (trainingdata[i, :XII] / 100.0)      # 8 FXII
        xₒ[9] = (1e-13)  # 9 FIIa
        xₒ = SF * xₒ # convert to nmol

        # setup the solver -
        tspan = (0.0, 20.0)
        prob = ODEProblem(balances, xₒ, tspan, dd; saveat = 0.01)
        soln = solve(prob)

        # smuch -
        T = soln.t
        X = hcat(soln.u...)

        # capture -
        species_archive[i] = X[speciesindex, :]
        T_archive[i] = T
    end

    # return -
    return (T_archive, species_archive)
end

function objective(model::Dict{String,Any}, trainingdata::DataFrame)

    # run the simulation -
    (T,X) = simulation(model, trainingdata)

    # compute the error vector -
    # think about this ...
    return nothing
end

max_stars_repo_name: varnerlab/UVM-TopDown-LegacyModel

max_issues_repo_name: varnerlab/UVM-TopDown-LegacyModel

max_forks_repo_name: varnerlab/UVM-TopDown-LegacyModel


text:
using LinearAlgebra, FFTW

function create_testfuncs(g::OneDGrid)
    if nx <= 8; error("nx must be > 8"); end

    m = 5
    k0 = g.k[2] # Fundamental wavenumber

    # Test function
    φ = π/3
    f1 = cos.(m*k0*g.x .+ φ)
    f1h = fft(f1)
    f1hr = rfft(f1)

    f1hr_mul = zeros(Complex{eltype(g.x)}, g.nkr)
    mul!(f1hr_mul, g.rfftplan, f1)

    # Theoretical values of the fft and rfft
    f1h_th = zeros(Complex{Float64}, size(f1h))
    f1hr_th = zeros(Complex{Float64}, size(f1hr))

    for i in 1:g.nk
      if abs(real(g.k[i])) == m*k0
        f1h_th[i] = -exp(sign(real(g.k[i]))*im*φ)*g.nx/2
      end
    end

    for i in 1:g.nkr
      if abs(real(g.k[i]))==m*k0
        f1hr_th[i] = -exp(sign(real(g.kr[i]))*im*φ)*g.nx/2
      end
    end

    f1, f1h, f1hr, f1hr_mul, f1h_th, f1hr_th
end

function create_testfuncs(g::TwoDGrid)
    if nx <= 8; error("nx must be > 8"); end

    m, n = 5, 2
    k0 = g.k[2] # fundamental wavenumber
    l0 = g.l[2] # fundamental wavenumber

    # some test functions
    f1 = @. cos(m*k0*g.X) * cos(n*l0*g.Y)
    f2 = @. sin(m*k0*g.X + n*l0*g.Y)

    f1h = fft(f1)
    f2h = fft(f2)
    f1hr = rfft(f1)
    f2hr = rfft(f2)

    f1hr_mul = zeros(Complex{eltype(g.x)}, (g.nkr, g.nl))
    f2hr_mul = zeros(Complex{eltype(g.x)}, (g.nkr, g.nl))
    mul!(f1hr_mul, g.rfftplan, f1)
    mul!(f2hr_mul, g.rfftplan, f2)

    # Theoretical results
    f1h_th = zeros(Complex{eltype(g.x)}, size(f1h))
    f2h_th = zeros(Complex{eltype(g.x)}, size(f2h))
    f1hr_th = zeros(Complex{eltype(g.x)}, size(f1hr))
    f2hr_th = zeros(Complex{eltype(g.x)}, size(f2hr))

    for j in 1:g.nl, i in 1:g.nk
      if ( abs(real(g.K[i, j])) == m*k0 && abs(real(g.L[i, j])) == n*l0 )
        f1h_th[i, j] = - g.nx*g.ny/4
      end
      if ( real(g.K[i, j]) == m*k0 && real(g.L[i, j]) == n*l0 )
        f2h_th[i, j] = -g.nx*g.ny/2
      elseif ( real(g.K[i, j]) == -m*k0 && real(g.L[i, j]) == -n*l0 )
        f2h_th[i, j] = g.nx*g.ny/2
      end
    end
    f2h_th = -im*f2h_th;

    for j in 1:g.nl, i in 1:g.nkr
      if ( abs(g.Kr[i, j])==m*k0 && abs(g.L[i, j])==n*l0 )
        f1hr_th[i, j] = - g.nx*g.ny/4
      end
      if ( real(g.Kr[i, j]) == m*k0 && real(g.L[i, j]) == n*l0 )
        f2hr_th[i, j] = -g.nx*g.ny/2
      elseif ( real(g.Kr[i, j]) == -m*k0 && real(g.L[i, j]) == -n*l0 )
        f2hr_th[i, j] = g.nx*g.ny/2
      end
    end
    f2hr_th = -im*f2hr_th;

    f1, f2, f1h, f2h, f1hr, f2hr, f1hr_mul, f2hr_mul,
      f1h_th, f1hr_th, f2h_th, f2hr_th
end

# -----------------------------------------------------------------------------
# FFT's TEST FUNCTIONS

rtol = 1e-12

function test_fft_cosmx(g::OneDGrid)
    f1, f1h, f1hr, f1hr_mul, f1h_th, f1hr_th = create_testfuncs(g)
    isapprox(f1h, f1h_th, rtol=rtol)
end

function test_rfft_cosmx(g::OneDGrid)
    f1, f1h, f1hr, f1hr_mul, f1h_th, f1hr_th = create_testfuncs(g)
    isapprox(f1hr, f1hr_th, rtol=rtol)
end

function test_rfft_mul_cosmx(g::OneDGrid)
    f1, f1h, f1hr, f1hr_mul, f1h_th, f1hr_th = create_testfuncs(g)
    isapprox(f1hr_mul, f1hr_th, rtol=rtol)
end

function test_fft_cosmxcosny(g::TwoDGrid)
    f1, f2, f1h, f2h, f1hr, f2hr, f1hr_mul, f2hr_mul,
            f1h_th, f1hr_th, f2h_th, f2hr_th = create_testfuncs(g)
    isapprox(f1h, f1h_th, rtol=rtol)
end

function test_rfft_cosmxcosny(g::TwoDGrid)
    f1, f2, f1h, f2h, f1hr, f2hr, f1hr_mul, f2hr_mul,
            f1h_th, f1hr_th, f2h_th, f2hr_th = create_testfuncs(g)
    isapprox(f1hr, f1hr_th, rtol=rtol)
end

function test_rfft_mul_cosmxcosny(g::TwoDGrid)
    f1, f2, f1h, f2h, f1hr, f2hr, f1hr_mul, f2hr_mul,
            f1h_th, f1hr_th, f2h_th, f2hr_th = create_testfuncs(g)
    isapprox(f1hr_mul, f1hr_th, rtol=rtol)
end

function test_fft_sinmxny(g::TwoDGrid)
    f1, f2, f1h, f2h, f1hr, f2hr, f1hr_mulPk, f2hr_mul,
            f1h_th, f1hr_th, f2h_th, f2hr_th = create_testfuncs(g)
    isapprox(f2h, f2h_th, rtol=rtol)
end

function test_rfft_sinmxny(g::TwoDGrid)
    f1, f2, f1h, f2h, f1hr, f2hr, f1hr_mul, f2hr_mul,
            f1h_th, f1hr_th, f2h_th, f2hr_th = create_testfuncs(g)
    isapprox(f2hr, f2hr_th, rtol=rtol)
end

function test_rfft_mul_sinmxny(g::TwoDGrid)
    f1, f2, f1h, f2h, f1hr, f2hr, f1hr_mul, f2hr_mul,
            f1h_th, f1hr_th, f2h_th, f2hr_th = create_testfuncs(g)
    isapprox(f2hr_mul, f2hr_th, rtol=rtol)
end

# Test 1D grid
nx = 32             # number of points
Lx = 2π             # Domain width
g1 = OneDGrid(nx, Lx)

# Test 2D rectangular grid
nx, ny = 32, 64     # number of points
Lx, Ly = 2π, 3π     # Domain width
g2 = TwoDGrid(nx, Lx, ny, Ly)

@test test_fft_cosmx(g1)
@test test_rfft_cosmx(g1)
@test test_rfft_mul_cosmx(g1)

@test test_fft_cosmxcosny(g2)
@test test_rfft_cosmxcosny(g2)
@test test_rfft_mul_cosmxcosny(g2)
@test test_fft_sinmxny(g2)
@test test_rfft_sinmxny(g2)
@test test_rfft_mul_sinmxny(g2)