text:
module ReconAcquisitions
using FFTW
using Gadgetron
using Gadgetron.MRD
using Gadgetron.Stream
using Gadgetron.Default
using PartialFunctions
using LinearAlgebra
using Setfield

fft_scaling(xsize,dims=[1]) = sqrt(reduce((*),map(d->xsize[d],dims)))

ortho_fft(x,dims=[1]) = fft(x,dims) ./ Float32(fft_scaling(size(x),dims))

ortho_ifft(x,dims=[1]) = bfft(x,dims) ./ Float32(fft_scaling(size(x),dims))

cfft(x, dims=[1]) = ifftshift(ortho_fft(fftshift(x,dims),dims),dims)

cifft(x, dims=[1]) = fftshift(ortho_ifft(ifftshift(x,dims),dims),dims)

function noise_adjustment(header::MRDHeader, noise_acq::Acquisition)

    noise_bandwidth = header.acquisitionSystemInformation.relativeReceiverNoiseBandwidth |>   default(0.793f0)

    function whitening_transform(noise::Array{ComplexF32,2})
        mat = Float32(1.0 / (size(noise)[1] - 1)) * Hermitian(adjoint(noise)*noise)
        cholesky(mat).U |> inv
    end

    whitening_matrix = whitening_transform(noise_acq.data)
    noise_dwell_time = noise_acq.header.sample_time_us

    scaling_factor(acq::MRD.Acquisition) =
        sqrt(2 * acq.header.sample_time_us * noise_bandwidth / noise_dwell_time)

    apply_whitening(acq::MRD.Acquisition) = MRD.Acquisition(
        acq.header,
        scaling_factor(acq) * acq.data * whitening_matrix,
        acq.trajectory,
    )

    return Map(apply_whitening; spawn=true, buffer_size=Inf)
end

function remove_oversampling(header::MRDHeader)
    encoding_space = header.encoding[1].encodedSpace.matrixSize
    recon_space = header.encoding[1].reconSpace.matrixSize

    if encoding_space == recon_space
        return Map(identity)
    end

    x0 = (encoding_space.x - recon_space.x) ÷ 2 + 1
    x1 = x0 + recon_space.x -1

    function crop_acquisition(acq::MRD.Acquisition)
        img_space = cifft(acq.data)
        img_space = img_space[x0:x1, :]
        acq = @set acq.header.center_sample = recon_space.x ÷ 2
        acq.data = cfft(img_space)
        return acq
    end

    return Map(crop_acquisition; spawn=true, buffer_size=Inf)

end

function make_kspace(header::MRDHeader)

	matrix_size = header.encoding[1].encodedSpace.matrixSize

	function inner(acqs)
		number_of_samples, number_of_channels = size(acqs[1].data)
		buffer = zeros(ComplexF32,number_of_samples,matrix_size.y,matrix_size.z,number_of_channels)
		for acq in acqs
			buffer[:,acq.header.idx.kspace_encode_step_1+1,acq.header.idx.kspace_encode_step_2+1,:] = acq.data
		end
		return acqs[1].header,buffer

	end

	return Map(inner)
end

function reconstruct_image(header::MRDHeader)

    field_of_view = header.encoding[1].reconSpace.fieldOfView_mm
    fov = (field_of_view.x,field_of_view.y,field_of_view.z)

	recon(x) = cifft(x,[1,2,3])
    coil_combine(x) = sqrt.(sum(abs2,x,dims=4))

	inner(index, reference::AcquisitionHeader, buffer) =   Image(ImageHeader(reference; image_index=index, image_type = MRD.ImageType.magnitude, field_of_view=fov),  coil_combine(recon(buffer)))
    inner_splat(args ) = inner(args[1],args[2][1],args[2][2])
	return Map(inner_splat)
end

function reconstruct_acquisitions(connection)
    start = time()
    header = connection.header
    noise_data = connection |>
       TakeWhile(acq -> :is_noise_measurement ∈ acq.header.flags) |> collect |> last
    noise_whitener = noise_adjustment(header, noise_data)

    connection |>  noise_whitener |>
    remove_oversampling(header) |>
    SplitBy(acq -> :last_in_slice ∈ acq.header.flags, keepend = true) |>
    make_kspace(header) |>
    Enumerate() |>
    reconstruct_image(header) .|>
    push! $ connection

    tot_time = time()-start
    @info "Total recon time $tot_time"
end

Base.precompile(Tuple{typeof(reconstruct_acquisitions),Gadgetron.MRDChannel})   # time: 0.5008968

end


text:
module SchwarzChristoffel

export SchwarzChristoffelDerivate, SchwarzChristoffelMap, segment, integral, finv

using Bernstein
using OrdinaryDiffEq
using FastGaussQuadrature
using FixedPointNumbers
using Parameters: @unpack
using Roots
using StaticArrays

const NQUAD = 4

include("fargextremum.jl")

######################
# Low-level functions

"""
    schwarzchristoffelderivative(x,β,xv; skip = ())

Evluate the function `prod((xv.-x).^.-β)`.

Indices listed in skip are omitted from the product.
"""
function schwarzchristoffelderivative(x,β,xv; skip=())
    f = one(ComplexF64)
    for (i,(xk,βk)) in enumerate(zip(x,β))
        i in skip && continue
        f *= complex(xv-xk)^-βk
    end
    return f
end

function assemble_quadcache(β)
    qjac = map(βi->(SVector{NQUAD}.(gaussjacobi(NQUAD,0,-βi)), SVector{2NQUAD}.(gaussjacobi(2NQUAD,0,-βi))), β)
    qleg = (SVector{NQUAD}.(gausslegendre(NQUAD)), SVector{2NQUAD}.(gausslegendre(2NQUAD)))
    qinf = sum(β) > 1 ?
        (SVector{NQUAD}.(gaussjacobi(NQUAD,0,sum(β)-2)), SVector{2NQUAD}.(gaussjacobi(2NQUAD,0,sum(β)-2))) :
        nothing
    return (
        jacobi_quadrules=qjac,
        legendre_quadrule=qleg,
        inf_quadrule=qinf
    )
end

function adaptive_semijacobi_integral(f,x,β,Δx,tol,(qjac,qleg))
    Δx == 0 && return zero(ComplexF64)

    t = zero(Q0f63)
    I = zero(ComplexF64)

    Δt,ΔI = shrink_integral(
        Δt -> begin
            Iq = ntuple(i->
                semijacobi_integral(
                    f, x,β,Δx*-float(Δt),
                    qjac[i]
                ),
                Val(2)
            )
            return Iq[2],abs(Iq[1]-Iq[2])
        end,
        Q0f63(-1), tol
    )
    t += Δt; I += ΔI
    while t > -1
        Δt,ΔI = shrink_integral(
            Δt -> begin
                Iq = ntuple(i->
                    legendre_integral(
                        xv->f(xv)*complex(xv-x)^-β,
                        x+Δx*-float(t), Δx*-float(Δt),
                        qleg[i]
                    ),
                    Val(2)
                )
                return Iq[2],abs(Iq[1]-Iq[2])
            end,
            max(-1-t, Δt+Δt), tol
        )
        t += Δt; I += ΔI
    end
    return I::ComplexF64
end

legendre_integral(f,x,Δx, (xq,wq)) =
    Δx/2 * mapreduce(
        ((xqk,wqk),) -> f(x+Δx*(xqk+1)/2)*wqk,
        +, zip(xq,wq)
    )
semijacobi_integral(f,x,β,Δx, (xq,wq)) =
    complex(Δx/2)^(-β) * legendre_integral(f,x,Δx, (xq,wq))

function shrink_integral(IEfun,Δt, tol)
    I,E = IEfun(Δt)
    while E > sqrt(tol)*abs(I)
        Δt *= Q0f63(0.5)
        Δt == 0 && error("Could not converge Schwarz-Christoffel integral to desired accuracy")
        I,E = IEfun(Δt)
    end
    return Δt,I
end

struct InfSegmentMap
    x0::Float64
    b::Float64
end

function InfSegmentMap(xv,x1,x2)
    s = xv < x1 ? -1 : 1
    InfSegmentMap(xv,s*sqrt(xv*(xv-x1-x2) + x1*x2))
end

(csm::InfSegmentMap)(y) = csm.b/y-csm.b+csm.x0
grad(csm::InfSegmentMap) = y->-csm.b/y^2
finv(csm::InfSegmentMap) = x->csm.b/(x - csm.x0 + csm.b)

#############################
# SchwarzChristoffelDerivate

struct SchwarzChristoffelDerivate{X,B,QC}
    x::X
    β::B
    tol::Float64
    quadcache::QC
end

function SchwarzChristoffelDerivate(x,β; tol = eps())
    idx = sortperm(x)
    x = x[idx]
    β = β[idx]
    SchwarzChristoffelDerivate(x,β,tol,assemble_quadcache(β))
end

hasfiniteintegral(f::SchwarzChristoffelDerivate) = sum(f.β) > 1
nsegments(f::SchwarzChristoffelDerivate) = length(f.x)+1

"""
    (f::SchwarzChristoffelDerivate)(xv; skip = ())

Evluate the function `prod((xv.-x).^.-β)`.

Indices listed in skip are omitted from the product.
"""
(f::SchwarzChristoffelDerivate)(xv; skip=()) = schwarzchristoffelderivative(f.x,f.β,xv; skip=skip)

"""
    segment(f::SchwarzChristoffelDerivate, k)

Integral of `f` from `x[k]` to `x[k+1]`.
"""
segment(f::SchwarzChristoffelDerivate,k) = segment(f,one,k)
function segment(f::SchwarzChristoffelDerivate,w,k)
    @unpack x = f

    if k in 1:length(x)-1
        xm = (x[k]+x[k+1])/2
    elseif k == 0
        xm = InfSegmentMap(x[1],x[2],x[end])(0.5)
    elseif k == length(x)
        xm = InfSegmentMap(x[end],x[1],x[end-1])(0.5)
    else
        throw(ArgumentError("attempt to compute segment $k of $(length(x)+1)-segment SchwarzChristoffelDerivate"))
    end
    return integral(f,w,k,xm) - integral(f,w,k+1,xm)
end

"""
    integral(f::SchwarzChristoffelDerivate, k, xv)

Integral of `f` from `x[k]` to `xv`.
"""
integral(f::SchwarzChristoffelDerivate,k,xv) = integral(f,one,k,xv)
function integral(f::SchwarzChristoffelDerivate,w,k,xv)
    @unpack x,β,tol,quadcache = f
    @assert k in 0:length(x)+1
    qleg = quadcache.legendre_quadrule

    if k in 1:length(x)
        qjac = quadcache.jacobi_quadrules[k]
        return adaptive_semijacobi_integral(xx->w(xx)*f(xx,skip=(k,)),x[k],β[k],xv-x[k],tol,(qjac,qleg))
    else
        @assert isreal(xv)
        qinf = quadcache.inf_quadrule
        @assert !isnothing(qinf)
        m = InfSegmentMap(xv,x[1],x[end])
        return adaptive_semijacobi_integral(
            y->w(m(y))*f(m(y))*grad(m)(y),
            0, 2-sum(β), 1,
            tol,(qinf,qleg)
        )
    end
end

########################
# SchwarzChristoffelMap

struct SchwarzChristoffelMap{F<:SchwarzChristoffelDerivate,W,Z,Zinf}
    f::F
    w::W
    z::Z
    zinf::Zinf
end

function SchwarzChristoffelMap(f::SchwarzChristoffelDerivate, w = one)
    @unpack x = f
    z = similar(x,ComplexF64,)
    z[end] = 0
    for k = length(x)-1:-1:1
        z[k] = z[k+1] - segment(f,w,k)
    end
    zinf = hasfiniteintegral(f) ? segment(f,w,length(x)) : Inf
    return SchwarzChristoffelMap(f,w,z,zinf)
end

SchwarzChristoffelMap(x,β, args...; kwargs...) = SchwarzChristoffelMap(SchwarzChristoffelDerivate(x,β; kwargs...),args...)

function segment(F::SchwarzChristoffelMap,k)
    if k in 1:nsegments(F)-2
        return F.z[k+1]-F.z[k]
    elseif k == 0
        return F.z[1] - F.zinf
    elseif k == nsegments(F)-1
        return F.zinf - F.z[end]
    else
        throw(ArgumentError("attempt to compute segment $k of $(length(x)+1)-segment SchwarzChristoffelMap"))
    end
end
nsegments(F::SchwarzChristoffelMap) = nsegments(F.f)
grad(F::SchwarzChristoffelMap) = x->F.w(x)*F.f(x)

(F::SchwarzChristoffelMap)(xv::Real) = F(complex(xv))
function (F::SchwarzChristoffelMap)(xv::Complex)
    @unpack f,w,z,zinf = F
    @unpack x = f

    abs(xv) == Inf && return zinf

    k = fargmax(1:length(x)) do k
        xv == x[k] && return Inf

        # Map from x[k]..xv to -1..1
        ϕ = xx->(xx-x[k])/(xv-x[k]) - (xx-xv)/(x[k]-xv)

        k == 1 && return Bernstein.radius(ϕ(x[2]))
        k == length(x) && return Bernstein.radius(ϕ(x[end-1]))
        return min(
            Bernstein.radius(ϕ(x[k-1])),
            Bernstein.radius(ϕ(x[k+1]))
        )
    end
    return z[k] + integral(f,w,k,xv)
end

###############################
# InverseSchwarzChristoffelMap

struct InverseSchwarzChristoffelMap{SCM,X,D,Z}
    F::SCM
    x::X
    d::D
    z::Z
end

finv(F::SchwarzChristoffelMap, x = default_starting_points(F)) =
    InverseSchwarzChristoffelMap(F,x,conj.(sign.(F.f.(x))),F.(x))

const NSAMPLES = 3

function default_starting_points(F::SchwarzChristoffelMap)
    @unpack x = F.f
    n = length(x)
    midpoints = LinRange(0,1, 2NSAMPLES+1)[2:2:end]
    xs = Vector{Float64}(undef, NSAMPLES*(n+1))
    for k = 1:n-1
        xs[NSAMPLES*(k-1) .+ (1:NSAMPLES)] .= x[k] .+ (x[k+1]-x[k]).*midpoints
    end
    xs[NSAMPLES*(n-1) .+ (1:NSAMPLES)] .= InfSegmentMap(x[1],x[2],x[end]).(midpoints)
    xs[NSAMPLES*n .+ (1:NSAMPLES)] .= InfSegmentMap(x[end],x[1],x[end-1]).(midpoints)
    return xs
end

function find_starting_point(x,d,z,ze)
    @assert length(x) == length(d) == length(z)
    k = fargmin(1:length(x)) do k
        imag(d[k]*(ze-z[k])) < 0 && return Inf
        return abs(ze - z[k])
    end
    return x[k],z[k]
end

function (iF::InverseSchwarzChristoffelMap)(ze)
    @unpack F,x,d,z = iF
    xs,zs = find_starting_point(x,d,z,ze)
    x̃e = finv_ode(F,xs,zs,ze)
    return finv_newton(F,x̃e,ze)
end

function finv_ode(F::SchwarzChristoffelMap,xs,zs,ze)
    prob = ODEProblem(
        (x,_,z)->begin
            dx = (ze-zs)/grad(F)(x)
            # Make sure rounding errors do not push us into the lower half plane
            return imag(x) == 0 ? complex(real(dx),max(0,imag(dx))) : dx
        end,
        ComplexF64(xs), (0.0,1.0)
    )
    sol = solve(
        prob, Tsit5(),
        reltol = 1e-3,
        save_everystep = false
    )
    return sol(1)
end

function finv_newton(F::SchwarzChristoffelMap,x̃e,ze)
    fd = (x->F(x) - ze), (x->grad(F)(x))
    return find_zero(
        fd,
        complex(x̃e),
        Roots.Newton(),
        rtol=F.f.tol
    )
end

end

ext: jl


text:
using Clang, Clang.cindex
# We'll replace some of the wrap functions, and we need to call the following:
import Clang.wrap_c: repr_jl, rep_type, rep_args, name_safe

const HEADER = "/usr/local/cuda-8.0/targets/x86_64-linux/include/cudnn.h"

context=wrap_c.init(output_file="gen_libcudnn.jl",
                    common_file="gen_libcudnn_h.jl",
                    header_library=x->"libcufft",
                    clang_includes=["/usr/include"],
                    header_wrapped=(x,y)->(contains(x,"cuda")))

# context.options = wrap_c.InternalOptions(true)  # wrap structs, too
path = "/usr/local/cuda-8.0/targets/x86_64-linux/include"
headers = [joinpath(path,"cudnn.h")]

# Customize the wrap function for functions. This was copied
# from Clang/src/wrap_c.jl, with the following customizations:
#   - error-check each function that returns a cudaError_t
#   - omit types from function prototypes
function wrap_c.wrap(buf::IO, funcdecl::FunctionDecl, libname::String)
    function print_args(buf::IO, cursors, types)
        i = 1
#        for (c,t) in zip(cursors,types)
#            print(buf, name_safe(c), "::", t)
#            (i < length(cursors)) && print(buf, ", ")
#            i += 1
#        end
        for c in cursors
            print(buf, name_safe(c))
            (i < length(cursors)) && print(buf, ", ")
            i += 1
        end
    end

    cu_spelling = spelling(funcdecl)

    funcname = spelling(funcdecl)
    arg_types = cindex.function_args(funcdecl)
    args = [x for x in search(funcdecl, ParmDecl)]
    arg_list = tuple( [repr_jl(x) for x in arg_types]... )
    ret_type = repr_jl(return_type(funcdecl))

    print(buf, "function ")
    print(buf, spelling(funcdecl))
    print(buf, "(")
    print_args(buf, args, [myrepr_jl(x) for x in arg_types])
    println(buf, ")")
    print(buf, "  ")
    ret_type == "cufftResult" && print(buf, "checkerror(")
    print(buf, "ccall( (:", funcname, ", ", libname, "), ")
    print(buf, rep_type(ret_type))
    print(buf, ", ")
    print(buf, rep_args(arg_list), ", ")
    for (i,arg) in enumerate(args)
        print(buf, name_safe(arg))
        (i < length(args)) && print(buf, ", ")
    end
    ret_type == "cufftResult" && print(buf, ")")
    println(buf, ")")
    println(buf, "end")
end

function myrepr_jl(x)
    str = repr_jl(x)
    return (str == "Ptr{Cint}") ? "Array{Cint}" : str
end

wrap_c.wrap_c_headers(context, headers)

# run(wrap_c.init(headers = [HEADER],
#                 # index = None,
#                 # common_file="types.jl",
#                 output_file="libcudnn.jl",
#                 # output_dir = "",
#                 # clang_args = String[],
#                 # clang_includes = ["/usr/lib/clang/3.4.2/include","/usr/local/cuda/include","/usr/include"],
#                 clang_includes = ["/usr/usc/clang/default/include", "/usr/usc/cuda/default/include"],
#                 # clang_diagnostics = true,
#                 header_wrapped=(header,cursorname)->(contains(cursorname,"cudnn")),
#                 header_library=x->"libcudnn",
#                 # header_outputfile = None,
#                 cursor_wrapped=(cursorname,cursor)->!isempty(cursorname),
#                 # options = InternalOptions(),
#                 # rewriter = x->x,
#                 ))


text:
using Documenter
using Cosmology

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

makedocs(;
    modules = [Cosmology],
    authors = "Julia Astro",
    repo = "https://github.com/JuliaAstro/Cosmology.jl/blob/{commit}{path}#L{line}",
    sitename = "Cosmology.jl",
    format = Documenter.HTML(;
        prettyurls = get(ENV, "CI", "false") == "true",
        canonical = "https://juliaastro.github.io/Cosmology.jl",
        assets = String[],),
    pages = [
        "Home" => "index.md",
        "API/Reference" => "api.md"
    ],)

deploydocs(;
    repo = "github.com/JuliaAstro/Cosmology.jl",)


text:
abstract type ButcherTableau end

"""
    RungeKuttaMethod(α, β, γ, order)

Simple RungeKuttaMethod which can be used to solve a set of differential equations.
It's required that the number of rows and columns of γ match the dimensions of α and β respectively.
"""
struct RungeKuttaMethod <: ButcherTableau
    α::Vector{<:Real}
    β::Vector{<:Real}
    γ::Matrix{<:Real}
    order::Int
    RungeKuttaMethod(α, β, γ, order) = (length(α) == size(γ, 1)) & (length(β) == size(γ,2)) ?
        new(α, β, γ, order) : error("Not a valid Butcher Tableau")
end

"""
    AdaptiveRungeKuttaMethod(α, β, γ, order)

Extended version of RungeKuttaMethod where methods of two different orders are used to approximate the
error of each calculation. Allows for each step to be within a specified tolerance.
"""
struct AdaptiveRungeKuttaMethod <: ButcherTableau
    α::Vector{<:Real}
    β::Matrix{<:Real}
    γ::Matrix{<:Real}
    order::Int
    AdaptiveRungeKuttaMethod(α, β, γ, order) = length(α) == size(γ, 1) & size(β, 2) == size(γ,2) & size(β, 1) == 2 ?
        new(α, β, γ, order) : error("Not a valid Butcher Tableau")
end

RK4 = RungeKuttaMethod(
    [1/6, 1/3, 1/3, 1/6], # α
    [0, 1/2, 1/2, 1], # β
    [[0, 1/2, 0, 0] [0, 0, 1/2, 0] [0, 0, 0, 1] [0, 0, 0, 0]], # γ
    4 # order
)


text:
using WindowManager
using ModernGL
using TimeManager
using MatrixMath
using Quaternions
using StaticArrays

using JLGEngine.LibGL
using JLGEngine.ModelManager
using JLGEngine.ShaderManager
using JLGEngine.TransformManager
using JLGEngine.CameraManager
using JLGEngine.RenderManager

include("objects/light.jl")

SUN = Light()
SPOTLIGHT = Light()

SCRIPTTICK=0 #0.01
FRAMES = 0
TICKS = 0

fokusChange=false

_RELEASE = 0
_PRESS = 1
_REPEAT = 2
WINDOW = nothing
lastKey = 0
projection = IMat4x4f
camera = IMat4x4f
model = IMat4x4f
randomTick=false
isFocus=true
lastTime=0
currentFrames = 0
max_fps = 0
max_fmps = 0
time = 0
renderReady=false
cursorpos = Vec2f(0,0)
old_cursorpos = Vec2f(0,0)
windowsize = Vec2f(800,600)
fullscreen = false
keyPressed = false
mouseKeyPressed = false
modelID = 0
planeModel = nothing
cubeModel = nothing
getVertexBuffer = ()->nothing
getUVBuffer = ()->nothing

angle = Vec2f(0,0)
old_Angle = Vec2f(0,0)
r360 = pi * 2;

right = 0
up = 0
forward = 0
zoom = 0

plane = nothing
cube = nothing
fpsCam = nothing

#@printf("OnKey2: %i, %i\n", 'c','C')

function clearScreen()
	glClearColor(0.0, 0.0, 0.0, 0.0)
	glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT | GL_STENCIL_BUFFER_BIT)
	WindowManager.update(WINDOW)
	glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT | GL_STENCIL_BUFFER_BIT)
end

function OnWindowFocus(focused::Bool)
	if !fokusChange return end
	global isFocus = focused
	t=focused ? "" : " - PAUSE"
	WindowManager.title(WINDOW,"$(WINDOW.name)$t")
	setProgramTick(isFocus ? SCRIPTTICK : 0.1)
	setRenderUpdate(isFocus)
	clearScreen()
	println("Fokus: ", focused)
end

function OnWindowResize(size::Array{Int32,1})
	global windowsize=Vec2f(size[1],size[2])
	println("Window Size $(WINDOW.size)")
	CameraManager.setRatio(fpsCam, Float32(windowsize[1]/windowsize[2]))
end

keyFB = 0
keyLR = 0
keyUD = 0

wireframe = false

function OnKey(key::Number, scancode::Number, action::Number, mods::Number, unicode::Number)
  #println("OnKey: '$key', '",Char(unicode),"'($unicode) '$action'")

	if key == 290 && action == 1 # F1 = 290
		global fullscreen=!fullscreen
		WindowManager.fullscreen(WINDOW,fullscreen)

	elseif key == 291 && action == 1 # F2
		global wireframe
		wireframe = !wireframe
		println("Wireframe: ",wireframe)

	elseif key == 293 && action == 1 # F4 = 293
		#this.callList[:reload]()
		#this.callList[:shaderRefresh]()
		#whos()
		println(current_module())

	elseif unicode == UInt32('w')
		global keyFB = (action > 0)?1:0
	elseif unicode == UInt32('s')
		global keyFB = (action > 0)?-1:0
	elseif unicode == UInt32('d')
		global keyLR = (action > 0)?1:0
	elseif unicode == UInt32('a')
		global keyLR = (action > 0)?-1:0
	elseif unicode == UInt32(32)
		global keyUD = (action > 0)?1:0
	elseif unicode == UInt32('c') || key == 341
		global keyUD = (action > 0)?-1:0

	elseif unicode == UInt32('m') && action > 0
		global modelID += 1
		if modelID > 1 modelID = 0 end
		if modelID == 0 cube.model=planeModel end
		if modelID == 1 cube.model=cubeModel end
  end

	global keyPressed = action != 0

	lastKey = key
	#this.storage[:lastKey] = key

end

function OnMouseKey(key::Number, action::Number, mods::Number)
	#println("Mouse: ", key, ", ", action)
	global mouseKeyPressed = action != 0
  global mouseKeyReleased = action == 0

	#global rightMouseKeyPressed = (key == 1 && action == 1)
	#rightMouseKeyReleased = (key == 1 && action == 0)

	if mouseKeyPressed #rightMouseKeyPressed
		WindowManager.cursor(WINDOW, :CURSOR_DISABLED)
		global old_cursorpos = cursorpos
	end

	if mouseKeyReleased #rightMouseKeyReleased
		WindowManager.cursor(WINDOW, :CURSOR_NORMAL)
	end

end

function OnMousePos(pos::AbstractVector)
	global cursorpos = Vec2f(pos[1]/WINDOW.size[1],pos[2]/WINDOW.size[2])

	#println("Mouse: ", cursorpos)

	if mouseKeyPressed OnCamRotate() end #rightMouseKeyPressed
end

function OnMouseMove(pos::AbstractVector)
  #println("Mouse: ", pos)
end

function OnCamRotate()
	mx = cursorpos[1] - old_cursorpos[1]
	my = cursorpos[2] - old_cursorpos[2]
	global old_cursorpos = cursorpos

	fpsCam.transform.rotation+=Vec3f(-mx*2,my*2,0f0) #Vec3f((-mx+0.5f0),(-0.5f0+my),0f0)
	CameraManager.updateRotation(fpsCam)
end

function OnMove(key::Symbol, m::Number)
		forward = MatrixMath.forwardVector4(fpsCam.rotationMat)
		right = MatrixMath.rightVector4(fpsCam.rotationMat)
		up = MatrixMath.upVector4(fpsCam.rotationMat)

		if key == :FORWARD
			fpsCam.transform.position += forward*(m*0.05f0)
		elseif key == :RIGHT
			fpsCam.transform.position -= right*(m*0.05f0) #+Vec3f(-right*0.02f0,-up*0.02f0,forward*0.02f0)
		elseif key == :UP
			fpsCam.transform.position -= up*(m*0.05f0)
		end
		CameraManager.updateTranslation(fpsCam)
end

OnShaderPropertyUpdate = (
	"iDummy" => () -> 0
	,"iResolution" => () -> windowsize
	,"cursorpos" => () -> cursorpos
	,"key" => () -> keyPressed || mouseKeyPressed
	,"iTessLevel" => () -> UInt32(65)
	,"iGlobalTime" => () -> Float32(programTime())
	,"programTime" => () -> Float32(scriptTime())
	,"programTick" => () -> Float32(SCRIPTTICK)
	,"frames" => () -> Float32(currentFrames)
	#,"iProjection" => () -> fpsCam.projectionMat
	#,"iView" => () -> fpsCam.viewMat
	#,"iModel" => () -> IMat4x4f
	#,"iUV" => () -> RenderManager.is(:plane) ? RenderManager.getData(:VERTEX) : RenderManager.getData(:UV)
	#,"iNORMAL" => () -> RenderManager.is(:plane) ? RenderManager.getData(:VERTEX) : RenderManager.getData(:NORMAL)
	,"iMVP" => () -> RenderManager.getMVP()
	,"iModel" => () -> fpsCam.modelMat
	,"iVP" => () -> fpsCam.projectionMat*fpsCam.viewMat
	,"iCameraPos" => () -> fpsCam.transform.position
	,"iPos" => () -> cube.transform.position
	,"iLightPos" => () -> SPOTLIGHT.position #getLightData(SPOTLIGHT)
	,"iLightColor" => () -> SPOTLIGHT.color
	#,"tex" => () -> Int32(0)
)

function OnShaderReload(path)
	println("Reloaded Shader ",path,".")

	#chomp(readline())
end

function OnPreRender()
	global rcount
	global renderReady = true

	glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT | GL_STENCIL_BUFFER_BIT)
	glViewport(0, 0, windowsize[1], windowsize[2])

	# Configure global settings
	#glDisable(GL_CULL_FACE)
	glEnable(GL_CULL_FACE)
	glFrontFace(GL_CCW)
	glCullFace(GL_BACK)

	glEnable(GL_DEPTH_TEST)
	glDepthFunc(GL_LESS)
	glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA)

	glDepthMask(GL_TRUE)
	glColorMask(GL_TRUE, GL_TRUE, GL_TRUE, GL_TRUE)
	glDepthRange(0.1f0,	1f0) #glDepthRange(0.0001f0,	1000f0)
	glClearColor(0.75, 0.75, 0.75, 0.1)

	glPatchParameteri(GL_PATCH_VERTICES, 4)

	#glBlendEquation(GL_FUNC_ADD)
	#glBlendFunc(GL_SRC_ALPHA,GL_ONE_MINUS_SRC_ALPHA)
	#glBlendFunc(GL_ONE_MINUS_SRC_ALPHA,GL_SRC_ALPHA)
	#glBlendFunc(GL_SRC_ALPHA, GL_SRC_COLOR)

	#glDisable(GL_BLEND)
	glEnable(GL_BLEND)
	#glBlendFunc(GL_ONE, GL_ONE)
	glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA)

	showFrames()
end

function OnPostRender()
	#println("------------------------------------")
end

function OnRender()
	RenderManager.renderAll()
end

function OnPostDraw()
	r=RenderManager.getSelected()
	if r.id == :FG ||  r.id == :BG
		#glBlendFuncSeparate(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA, GL_ONE, GL_ONE_MINUS_SRC_ALPHA)
		glEnable(GL_DEPTH_TEST)
	elseif r.id == :iplane
		#glDepthMask(GL_TRUE)
	elseif r.id == :cube
		#glEnable(GL_DEPTH_TEST)
		#glDepthFunc(GL_LESS)
		#glDisable(GL_CULL_FACE)
	end
	#println(r.id)
end

function OnPreDraw()
	r=RenderManager.getSelected()
	#println(r.id) #cube,iplane,plane

	if r.id == :FG || r.id == :BG
		glDisable(GL_DEPTH_TEST)
		#glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA)
		#glBlendFuncSeparate(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA, GL_ONE, GL_ONE_MINUS_SRC_COLOR);
	elseif r.id == :iplane
		#glDepthMask(GL_FALSE)
	elseif r.id == :plane
	elseif r.id == :cube
		#glDepthMask(GL_TRUE)
		#glColorMask(GL_FALSE, GL_FALSE, GL_FALSE, GL_TRUE)
		#glDepthFunc(GL_LESS)
		#glEnable(GL_CULL_FACE)
		#glCullFace(GL_FRONT)
		#glCullFace(GL_BACK)
		#glFrontFace(GL_CCW)

		#glDisable(GL_DEPTH_TEST)
		#glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA)
		#glBlendFuncSeparate(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA, GL_ONE, GL_ONE_MINUS_SRC_COLOR);
	end

	if wireframe glPolygonMode(GL_FRONT_AND_BACK, GL_LINE)
	else glPolygonMode(GL_FRONT_AND_BACK, GL_FILL)
	end

	#if !isFocus return end

	#if haskey(this.storage, :lastKey)
		#DRIVER.setShaderPropberty("key", this.storage[:lastKey] == UInt32('k') )
		#this.storage[:lastKey]=0
	#end
	#lastKey = 0

	#glActiveTexture(GL_TEXTURE0)
	#glBindTexture(GL_TEXTURE_2D, texture)

	#if isFocus RENDER()
	#else OnRenderPause()
	#end

	#model = mgl32.HomogRotate3D(float32(angle), mgl32.Vec3{0, 1, 0})

	#pointer([convert(Ptr{GLchar}, pointer(source))]))

	# Load the texture
	# texture = GRAPHICSDRIVER().CreateTexture(string(Environment.PATHS[:ROOT],"square.png"))
	# texture, err := GoniGL.NewTexture("square.png")
	# if err != nil {
	# 	panic(err)
	# }

	# draw
end

function OnRenderPause()

end

prevTime1 = Ref(0.0)
prevTime2 = Ref(0.0)

function showFrames()
	global TICKS
	global FRAMES
	global max_fps
	global max_fmps
	global WINDOW

	fps = FRAMES
	fpms = fps > 0 ? (1000.0 / fps) : 0

	# MAX
	if (max_fps < fps)
		max_fps = fps
		max_fmps = fpms
	end

	lastFrames = FRAMES
	FRAMES += 1

	#global currentFrames = fps/max_fps

	if OnTime(1.0,prevTime1)
		#global SCRIPTTICK
		#est = (1/SCRIPTTICK) * 0.6
		WindowManager.title(WINDOW,"$(WINDOW.name) - FPS $(round(fps, 2)) [$(round(max_fps, 2))] | FMPS $(round(fpms, 2)) [$(round(max_fmps, 2))] - Ticks $(TICKS)")
		global currentFrames = fps/max_fps
		FRAMES = 0
		TICKS = 0
	end

	if randomTick && OnTime(5.0,prevTime2)
		r = [0.1,0.05,0.01,0.005,0.001,0.0005,0.0001,0.00005,0.00001]
		global SCRIPTTICK=r[rand(1:length(r))] #0.00016*rand(1:100)
		setProgramTick(SCRIPTTICK)
		#if isFocus this.callList[:setProgramTick](SCRIPTTICK) end
		max_fps=0
		max_fmps=0
		FRAMES = 0
	end
end

prevTime3 = Ref(0.0)
function OnUpdate()
	global TICKS
	TICKS += 1
	#cube.transform.position = Vec3f(sin(FRAMES*0.01f0),0,0)

	if isFocus
		if OnTime(0.01, prevTime3)
			cube.transform.rotation = Vec3f(0.01f0,0.02f0,0.03f0) * programTime() * 10
			if keyFB != 0 OnMove(:FORWARD,keyFB) end
			if keyLR != 0 OnMove(:RIGHT,keyLR) end
			if keyUD != 0 OnMove(:UP,keyUD) end
		end
	end
end

function OnStart()
	println("Reload Entities...")

	robj_bg=RenderManager.create(:BG)
	robj_fg=RenderManager.create(:FG)

	robj_plane=RenderManager.create(:plane)
	robj_iplane=RenderManager.create(:iplane)

	robj_cube=RenderManager.create(:cube)
	robj_cube_sphere=RenderManager.create(:cube_sphere)

	robj_icosphere=RenderManager.create(:icosphere)

	robj_loaded_cube=RenderManager.create(:loaded_cube)
	robj_loaded_monkey=RenderManager.create(:loaded_monkey)
	robj_loaded_bannana=RenderManager.create(:loaded_bannana)

	global plane=robj_plane
	global cube=robj_cube_sphere

	println("Reload Shaders...")
	listenOnShaderReload = (p)->OnShaderReload(ShaderManager.getSource(p).path)
	ShaderManager.setListenerOnShaderPropertyUpdate(:OnShaderPropertyUpdate, OnShaderPropertyUpdate)

	shader_tess = ShaderManager.create(:Prog1, "shaders/tess.glsl", listenOnShaderReload)
	shader_mandel = ShaderManager.create(:Prog2, "shaders/mandelbrot.glsl", listenOnShaderReload) #mandelbrot, uv
	shader_color = ShaderManager.create(:Prog3, "shaders/color.glsl", listenOnShaderReload)
	shader_uv = ShaderManager.create(:Prog4, "shaders/uv.glsl", listenOnShaderReload)
	shader_texture = ShaderManager.create(:Prog5, "shaders/texture.glsl", listenOnShaderReload)
	shader_sphere_texture = ShaderManager.create(:Prog6, "shaders/sphere_texture.glsl", listenOnShaderReload)
	shader_fg = ShaderManager.create(:FG, "shaders/fg.glsl", listenOnShaderReload)
	shader_bg = ShaderManager.create(:BG, "shaders/bg.glsl", listenOnShaderReload)

	#ShaderManager.API.reload(shader_sphere_texture)

	robj_loaded_cube.shaderProgram=shader_texture
	robj_loaded_monkey.shaderProgram=shader_mandel
	robj_loaded_bannana.shaderProgram=shader_mandel

	robj_plane.shaderProgram=shader_tess
	robj_iplane.shaderProgram=shader_color
	robj_cube.shaderProgram=shader_mandel
	robj_icosphere.shaderProgram=shader_mandel
	robj_cube_sphere.shaderProgram=shader_sphere_texture

	robj_fg.shaderProgram=shader_fg
	robj_bg.shaderProgram=shader_bg

	#RenderManager.setMode(robj_iplane, :POLYGONMODE, [:FRONT_AND_BACK,:LINE])
	#RenderManager.setMode(robj_iplane, :POLYGONMODE, [:FRONT_AND_BACK,:FILL])

	println("Reload Models...")
	model_loaded_cube = ModelManager.create(:cube)
	model_loaded_monkey = ModelManager.create(:monkey)
	model_loaded_bannana = ModelManager.create(:bannana)

	model_iplane = ModelManager.create(:iplane)
	model_plane = ModelManager.create(:plane)
	model_cube = ModelManager.create(:cube)
	model_icosphere = ModelManager.create(:icosphere)

	global planeModel = model_plane
	global cubeModel = model_cube

	robj_loaded_cube.model = model_loaded_cube
	robj_loaded_monkey.model = model_loaded_monkey
	robj_loaded_bannana.model = model_loaded_bannana

	robj_iplane.model = model_iplane
	robj_plane.model = model_plane
	robj_cube.model = model_cube
	robj_cube_sphere.model = model_cube
	robj_icosphere.model = model_icosphere

	robj_fg.model=model_iplane
	robj_bg.model=model_iplane

	model_loaded_cube.source.path = "models/box.obj"
	model_loaded_monkey.source.path = "models/suzanne.obj"
	model_loaded_bannana.source.path = "models/Bananne_Mario.obj"

	print("MODEL LOADED CUBE: "); @time ModelManager.addMesh(model_loaded_cube, :MODEL)
	print("MODEL LOADED MONKEY: "); @time ModelManager.addMesh(model_loaded_monkey, :MODEL)
	print("MODEL LOADED BANNANA: "); @time ModelManager.addMesh(model_loaded_bannana, :MODEL)

	print("MODEL IPLANE: "); @time ModelManager.addMesh(model_iplane, :IPLANE)
	print("MODEL PLANE: "); @time ModelManager.addMesh(model_plane, :PLANE)
	print("MODEL CUBE: "); @time ModelManager.addMesh(model_cube, :CUBE2)
	print("MODEL ICOSPHERE: "); @time ModelManager.addMesh(model_icosphere, :ICOSPHERE)

	println("Reload Textures...")
	RenderManager.addTexture(robj_loaded_cube,"textures/square.png")
	textures_cube = robj_loaded_cube.textures

	println("Reload Cameras...")
	camera_fps = CameraManager.create(:FPS)
	camera_bg = CameraManager.create(:BG)

	global fpsCam = camera_fps

	robj_loaded_cube.camera = camera_fps
	robj_loaded_monkey.camera = camera_fps
	robj_loaded_bannana.camera = camera_fps

	robj_iplane.camera = camera_fps
	robj_plane.camera = camera_fps
	robj_cube.camera = camera_fps
	robj_cube_sphere.camera = camera_fps
	robj_icosphere.camera = camera_fps

	robj_fg.camera=camera_bg
	robj_bg.camera=camera_bg

	println("Reload Transformations...")
	robj_fg.transform.scaling = Vec3f(0.25,0.25,1)
	robj_fg.transform.position = Vec3f(0.75,0.75,-1)
	robj_bg.transform.position = Vec3f(0,0,-1)

	robj_cube.transform.position = Vec3f(-3,2,-5)
	robj_loaded_monkey.transform.position = Vec3f(0,2,-5)
	robj_loaded_bannana.transform.position = Vec3f(3,2,-5)

	robj_loaded_cube.transform.position = Vec3f(-3,5,-5)
	robj_loaded_cube.transform.rotation = Vec3f(0,0,0)

	robj_cube_sphere.transform.position = Vec3f(0,0,-3)
	robj_cube_sphere.transform.rotation = Vec3f(0,0,-3)

	robj_plane.transform.rotation = Vec3f(0,-pi*0.5,0)
	robj_plane.transform.position = Vec3f(0,-2,-5)
	robj_plane.transform.scaling = Vec3f(10,10,1)

	robj_iplane.transform.position = Vec3f(-3,0,-3)
	robj_icosphere.transform.position = Vec3f(3,0,-3)

	# cam
	#camera_fps.transform.position = Vec3f(0,0,-2)
	#camera_fps.transform.rotation = Vec3f(0,0,0) #deg2rad(180)

	CameraManager.setPerspective(camera_fps, 60.0f0, Float32(windowsize[1]/windowsize[2]), 0.001f0, 10000.0f0)
	CameraManager.setPerspective(camera_bg, 53.0f0, 1.0f0, 0.001f0, 100.0f0)

	# enable/disable for rendering
	robj_fg.enabled=true
	robj_bg.enabled=true

	#robj_iplane.enabled=false
	#robj_icosphere.enabled=false
	#robj_plane.enabled=false
	#robj_cube.enabled=false

	# transparency order
	RenderManager.renderSortedList(true)
	#robj_model.order=Symbol("-")
	robj_bg.order=Symbol("0")
	robj_cube_sphere.order=Symbol("1")
	robj_plane.order=Symbol("1")
	robj_iplane.order=Symbol("T")
	robj_icosphere.order=Symbol("1")
	robj_cube.order=Symbol("1")
	robj_loaded_cube.order=Symbol("1")
	robj_fg.order=Symbol("Z")

	println("Generate Cubes...")
	max=4
	for z=1:max
		for y=1:max
			for x=1:max
				o=RenderManager.create(Symbol(:CUBE,x+(y-1)*max+(z-1)*max*max))
				o.enabled=false
				#RenderManager.register(o,true,false)
				o.model=model_cube
				o.shaderProgram=shader_texture
				o.textures = textures_cube
				o.camera=camera_fps
				o.transform.position = Vec3f(x*2-max,y*2,z*2)
			end
		end
	end

	global SPOTLIGHT.color = Vec4f(1,0,0,1)
end

function OnReload()
	global WINDOW=this.objref.WINDOW
	global windowsize = Vec2f(WINDOW.size[1],WINDOW.size[2])

	setProgramTick(SCRIPTTICK)
end

function OnInit()
	FILE = Base.@__FILE__
  println("Script: $(basename(FILE)), $(this.args), time: $(mtime(FILE))")
end

max_stars_repo_name: Gilga/JGE

max_issues_repo_name: Gilga/JGE

max_forks_repo_name: Gilga/JGE


text:
"""
    cannibalism(N::UnipartiteNetwork)

This function measure the number of species that practices cannibalism in a given network. A species is cannibal if she interacts with herself.
"""
function cannibalism(N::UnipartiteNetwork)
    sum(N[s,s] for s in species(N))
end


text:
"Generate a noise image for experimentation."
function noisepic_2(A, physwidth, physheight)
    n_rows, n_cols = size(A)
    no = Matrix{Float64}(undef, n_rows, n_cols)
    n_rows = size(no)[1]
    n_cols_in_division = div(size(no)[2], 3)
    length_one_cycle = CUTOFF_23 * DU_26
    dimension_aligned(OU_26 + (-0.4 * physwidth, 0.0m),  OU_26 +  (-0.4 * physwidth, length_one_cycle ))
    n_cycles = physheight / length_one_cycle
    for row in 1:n_rows
        for col in 1:n_cols
            divno = 1 + col ÷ n_cols_in_division
            rate = 1.414 * divno * n_cycles / n_rows
            # Noise is a function generating noisy data,
            # but where the lowest frequency has wavelength 1.41.
            # I.e., one full cycle when x or y varies from 0 to 1.41.
            no[row, col] = if row < div(n_rows, 2)
                noise(row * rate, col * rate)
            else
                noise(row * rate, col * rate, detail = 6, persistence = 0.9)
            end
        end
    end
    normalize_datarange(no)
end
function convolute_image_2(xs, ys, f_xy, n_xy, h, cutoff)
    # Output buffer
    M = Array{Float64}(undef, length(ys), length(xs)) # Image processing convention: a column has a horizontal line of pixels
    # Window, triangle. TODO Experiment with FIR, Hanning
    wf = SVector{NS_2}(range(1.0, 0.1, length = NS_2));
    wb = SVector{NS_2}(range(1.0, 0.1, length = NS_2));
    #x = xs[1]::Quantity{Float64}
    #y = ys[1]
    @assert eltype(xs) <: Quantity{Float64}
    @assert eltype(ys) <: Quantity{Float64}
    @assert n_xy(xs[1], ys[1]) isa Float64
    @assert first(ys) < last(ys)
    # prepare a mutable but fast working buffer forward, and another backward.
    vxf = similar(SVector{NS_2}(fill(xs[1], NS_2)))
    vyf = similar(vxf)
    vxb = similar(vxf)
    vyb = similar(vxf)
    rowsy = length(ys)
    for (i::Int64, x::Quantity{Float64}) in enumerate(xs), (j::Int64, y::Quantity{Float64}) in enumerate(ys)
        vxf[1], vyf[1], vxb[1], vyb[1] = x, y, x, y
        # Find the streamline path
        # Put coordinates to sample in forward and aft buffers (no need to store them in memory really?)
        rk4_steps!(f_xy, vxf, vyf, h)
        rk4_steps!(f_xy, vxb, vyb, -h)
        # Find the intensity for our one pixel.
        pv  = convolute_pixel(wf, wb, vxf, vyf, vxb, vyb, n_xy, h, cutoff)
        # Find the original indexes and update the image matrix
        ii = rowsy + 1 - j
        jj = i
        M[ii, jj] = pv / oneunit(pv)
    end
    normalize_datarange(M)
end

"""
    convolute_pixel(wf, wb, vxf, vyf, vxb, vyb, nxy, h, cutoff))
where
    wf    Window forward: A vector defining a convolution (filter) window, or kernel, with the same length as buffers.
          The second value is used for weighting the second coordinate.
    wb    Window backward: Similar, but for the backwards direction. The second value is used for weighting the second coordinate.
          Hence, it will be identical to wf in many cases.
    vxf   Vector or buffer of x coordinates, forward direction
    vyf   Vector y coordinates, forward.
    vxb   Vector x, backward
    vyb   Vector y, backwards direction
    nxy   Function, noise field
    h     Step length, duration
    cutoff Velocity, substituted to calculate out-of bounds streamlines
"""
function convolute_pixel(wf, wb, vxf, vyf, vxb, vyb, nxy, h, cutoff)
    @assert length(vxf) == length(vyf)
    @assert length(vxb) == length(vyb)
    x0 = vxf[1]
    y0 = vyf[1]
    prevcontrib = nxy(x0, y0) * cutoff * h
    if isnan(x0) || isnan(y0)
        return 0.0 * prevcontrib
    end
    pv = 0.0 * prevcontrib
    for i = 2:length(vxf)
        prevcontrib = _convolve_contrib(prevcontrib, wf, nxy, vxf, vyf, i)
        pv += prevcontrib
    end
    x0 = vxb[1]
    y0 = vyb[1]
    prevcontrib = nxy(x0, y0) * cutoff * h
    if isnan(x0) || isnan(y0)
        return 0.0 * prevcontrib
    end
    for i = 2:length(vxb)
        prevcontrib = _convolve_contrib(prevcontrib, wb, nxy, vxb, vyb, i)
        pv += prevcontrib
    end
    pv
end
"""
    convolute_pixel(w, vx,  vy, nxy, h, cutoff)
where
    w     Window : A vector defining a convolution (filter) window, or kernel, with the same length as buffers.
    vx    Vector or buffer of x coordinates
    vy    Vector y coordinates
    nxy   Function, noise field
    h     Step length, duration
    cutoff Maximum magnitude
"""
function convolute_pixel(w, vx, vy, nxy, h, cutoff)
    @assert length(vx) == length(vy)
    @assert length(w) == length(vy)
    x0 = vx[1]
    y0 = vy[1]
    prevcontrib = nxy(x0, y0) * cutoff * h
    if isnan(x0) || isnan(y0)
        return 0.0 * prevcontrib
    end
    pv = 0.0 * prevcontrib
    for i = 2:length(vx)
        prevcontrib = _convolve_contrib(prevcontrib, w, nxy, vx, vy, i)
        pv += prevcontrib
    end
    pv
end
function _convolve_contrib(prev, w, nxy, vx, vy, i)
    @assert i > 1
    x = vx[i]
    y = vy[i]
    x0 = vx[i - 1]
    y0 = vy[i - 1]
    if isnan(x0) || isnan(y0) || isnan(x) || isnan(y)
        return prev
    end
    k = w[i]
    ds =  hypot(x - x0, y - y0)
    n = nxy(x, y)
    n * k * ds
end


text:
#
using UnPack, BenchmarkTools

mutable struct Field{T}
    u ::Array{T,2} # value
    u1::Array{T,2} # histories
    u2::Array{T,2}
    u3::Array{T,2}
end

#----------------------------------#
mutable struct Gield{T}
    u ::Array{T}
    uh::Vector{Array}
end

#----------------------------------#
function updateHist!(g::Gield)

    for i=length(g.uh):-1:2
        g.uh[i] .= g.uh[i-1]
    end

    g.uh[1] .= g.u
    return
end
#----------------------------------#
function updateHist!(u::AbstractArray)

    for i=length(u):-1:2
        u[i] .= u[i-1]
    end
    return
end
#----------------------------------#
function updateHist!(fld::Field)
    @unpack u,u1,u2,u3 = fld

    u3[:] .= u2[:]
    u2[:] .= u1[:]
    u1[:] .= u[:]

    return
end
#----------------------------------#
n = 10
a = ones(n,n)
b = Array[ (3-i).*a for i in 1:3]
f = Field(3*a,2*a,1*a,0*a)
g = Gield(3*a,b)
t = [3.,2.,1.,0.]

#println("t: $t")
#@btime updateHist!(t)
#println("t: $t")
println("#==================================#")
println("f: $(f.u[1]), $(f.u1[1]), $(f.u2[1]), $(f.u3[1])")
@btime updateHist!(f)
println("f: $(f.u[1]), $(f.u1[1]), $(f.u2[1]), $(f.u3[1])")
println("#==================================#")
println("g: $(g.u[1]), $(g.uh[1][1]), $(g.uh[2][1]), $(g.uh[3][1])")
@btime updateHist!(g)
println("g: $(g.u[1]), $(g.uh[1][1]), $(g.uh[2][1]), $(g.uh[3][1])")
println("#==================================#")
#


text:
using Makie
import Statistics: mean
import GMT: triangulate

"""
    PlotTopography(channelList::Vector{String}, channelValues::Vector{Number}; gridSize = 1000)

Plots a topographical map of the head over the desired points
given by `channelList` and their assigned `channelValues`,
with default grid size being `gridSize = 1000`.

# Example
```jldoctest
channelExamples = String["Fpz", "Fp1", "Fp2", "AF3", "AF4", "AF7", "AF8", "Fz", "F1", "F2", "F3", "F4", "F5", "F6", "F7", "F8", "FCz", "FC1", "FC2", "FC3", "FC4", "FC5", "FC6", "FT7", "FT8", "Cz", "C1", "C2", "C3", "C4", "C5", "C6", "T7", "T8", "CP1", "CP2", "CP3", "CP4", "CP5", "CP6", "TP7", "TP8", "Pz", "P3", "P4", "P5", "P6", "P7", "P8", "POz", "PO3", "PO4", "PO5", "PO6", "PO7", "PO8", "Oz", "O1", "O2"];

PlotTopography(channelExamples, rand(length(channelExamples)))
```
"""
function PlotTopography(channelList::Vector{String}, channelValues::Vector{T}; gridSize = 1000) where T <: Number

    # Basic checking
    length(channelList) == length(channelValues) || error(
        "The lengths of channelList($(length(channelList))) and
         channelValues($(length(channelValues))) should be the SAME.",
    )

    # Retrieve current `channelList` information from the standard channel information
    standardChannelLabel, standardCoordinateX, standardCoordinateY = standardChannels.label, standardChannels.x, standardChannels.y
    channelIndex = indexin(channelList, standardChannelLabel)
    channelXs, channelYs =
        standardCoordinateX[channelIndex], standardCoordinateY[channelIndex]
    maxX, maxY = maximum(channelXs), maximum(channelYs)
    channelPositions = Point2f.(channelXs, channelYs)

    # Extend data beyond the boundary, see:
    # https://www.mathworks.com/matlabcentral/fileexchange/72729-topographic-eeg-meg-plot
    extendData = append_nearest_values(channelXs, channelYs, channelValues)

    # Interpolate data with function `triangualte` from `GMT.jl``,
    # see: `gmthelp(triangulate)` or `?triangulate` for more information.
    interpolateData =
        triangulate(extendData, region = "-1/1/-1/1", inc = string(gridSize, "+n"))
    xs, ys, zs = interpolateData.x, interpolateData.y, interpolateData.z

    # Remove data outside of circle
    zs[abs.(complex.(repeat(xs, outer = gridSize), repeat(ys, inner = gridSize))) .> 1] .= NaN

    # Initianize the figure
    fig = Figure(resolution = (1300, 1200))
    ga = fig[1, 1] = GridLayout()
    gb = fig[1, 2] = GridLayout()

    ax1 = Axis(
        ga[1, 1],
        xzoomlock = true,
        yzoomlock = true,
        xgridvisible = false,
        ygridvisible = false,
        xticksvisible = false,
        yticksvisible = false,
        xticklabelsvisible = false,
        yticklabelsvisible = false,
        bottomspinevisible = false,
        leftspinevisible = false,
        rightspinevisible = false,
        topspinevisible = false,
    )

    # Draw the contour
    contourf!(ax1, xs, ys, zs)

    # Draw the scatter of the electrodes
    scatter!(
        ax1,
        channelPositions,
        color = :black,
        strokecolor = :black,
        strokewidth = 2,
        markersize = 6,
    )

    # Draw the labels of the electrodes
    text!(ax1, channelList, position = channelPositions)

    # Draw the circle around the head
    poly!(
        ax1,
        Circle(Point2f(0, 0), maxX),
        color = :transparent,
        strokecolor = :black,
        strokewidth = 4,
    )

    # Draw the nose
    lines!(ax1, nosePoints(-10, maxX), color = :black, linewidth = 4)
    lines!(ax1, nosePoints(+10, maxX), color = :black, linewidth = 4)

    # Draw the ears
    lines!(ax1, earPoints(-0.75, maxX), color = :black, linewidth = 4)
    lines!(ax1, earPoints(+0.75, maxX), color = :black, linewidth = 4)

    # Add Colorbar
    Colorbar(gb[1, 1], width = 20, ticks = 0:0.1:1)

    # Some decorations
    colsize!(fig.layout, 1, Aspect(1, 1.0))

    save("figure.png", fig, px_per_unit = 1) # output size = 1300 x 1000 pixels

    return fig
end

"""
    append_nearest_values(Xs, Ys, Vals; n = 8, radius = 1.2, k = 4)

Append `n = 8` extra points to the original data at a circle with `radius = 1.2`.
The values of these points are calculated from the nearst `k = 4` points' values
`Vals` at the original coordinates `(Xs, Ys)`.
"""
function append_nearest_values(Xs, Ys, Vals; n = 8, radius = 1.2, k = 4)
    newPoints = Array{Float64,2}(undef, n, 3)
    for i = 1:n
        coordinates = radius * exp(i * 2pi * im / n) # Euler's identity
        distances = coordinates .- complex.(Xs, Ys)
        newVal = mean(Vals[sortperm(abs.(distances)).<= k])
        newPoints[i, :] = [real(coordinates), imag(coordinates), newVal]
    end
    return vcat([Xs Ys Vals], newPoints)
end

"""
    nosePoints(angle, radius)

Generate two points to define the left (`angle > 0`) or right (`angle < 0`) line
to draw the contour of the avator nose, starting from the surface with `radius`.
"""
nosePoints(angle, radius) =
    [radius * Point2f(cosd(90 + angle), sind(90 + angle)), 0.85 * Point2f(0, 1)]

"""
    earPoints(focusX, radius; npoints = 100)

Generate `npoints = 100` points to define the ellipse curve of the left (`FocusX < 0`)
or right (``FocusX > 0`) ear, with the focus point being at `(FocusX, 0)`.
"""
function earPoints(focusX, radius; npoints = 100)
    coordinates = exp.(collect(1:npoints+1) * 2pi * im / npoints) # Euler's identity
    modifyCoordinates =
        complex.(focusX .+ 0.09 .* real(coordinates), 0.18 .* imag(coordinates))
    modifyCoordinates[abs.(modifyCoordinates).< radius] .= NaN
    return Point2f.(real(modifyCoordinates), imag(modifyCoordinates))
end

"""
    standardChannels

Standard 1005 EEG electrode positions. from: https://github.com/sappelhoff/eeg_positions/blob/main/data/Nz-T10-Iz-T9/standard_1005_2D.tsv
"""
const standardChannels = (
    label = String["AF1", "AF10", "AF10h", "AF1h", "AF2", "AF2h", "AF3", "AF3h", "AF4", "AF4h", "AF5", "AF5h", "AF6", "AF6h", "AF7", "AF7h", "AF8", "AF8h", "AF9", "AF9h", "AFF1", "AFF10", "AFF10h", "AFF1h", "AFF2", "AFF2h", "AFF3", "AFF3h", "AFF4", "AFF4h", "AFF5", "AFF5h", "AFF6", "AFF6h", "AFF7", "AFF7h", "AFF8", "AFF8h", "AFF9", "AFF9h", "AFFz", "AFp1", "AFp10", "AFp10h", "AFp1h", "AFp2", "AFp2h", "AFp3", "AFp3h", "AFp4", "AFp4h", "AFp5", "AFp5h", "AFp6", "AFp6h", "AFp7", "AFp7h", "AFp8", "AFp8h", "AFp9", "AFp9h", "AFpz", "AFz", "C1", "C1h", "C2", "C2h", "C3", "C3h", "C4", "C4h", "C5", "C5h", "C6", "C6h", "CCP1", "CCP1h", "CCP2", "CCP2h", "CCP3", "CCP3h", "CCP4", "CCP4h", "CCP5", "CCP5h", "CCP6", "CCP6h", "CCPz", "CP1", "CP1h", "CP2", "CP2h", "CP3", "CP3h", "CP4", "CP4h", "CP5", "CP5h", "CP6", "CP6h", "CPP1", "CPP1h", "CPP2", "CPP2h", "CPP3", "CPP3h", "CPP4", "CPP4h", "CPP5", "CPP5h", "CPP6", "CPP6h", "CPPz", "CPz", "Cz", "F1", "F10", "F10h", "F1h", "F2", "F2h", "F3", "F3h", "F4", "F4h", "F5", "F5h", "F6", "F6h", "F7", "F7h", "F8", "F8h", "F9", "F9h", "FC1", "FC1h", "FC2", "FC2h", "FC3", "FC3h", "FC4", "FC4h", "FC5", "FC5h", "FC6", "FC6h", "FCC1", "FCC1h", "FCC2", "FCC2h", "FCC3", "FCC3h", "FCC4", "FCC4h", "FCC5", "FCC5h", "FCC6", "FCC6h", "FCCz", "FCz", "FFC1", "FFC1h", "FFC2", "FFC2h", "FFC3", "FFC3h", "FFC4", "FFC4h", "FFC5", "FFC5h", "FFC6", "FFC6h", "FFCz", "FFT10", "FFT10h", "FFT7", "FFT7h", "FFT8", "FFT8h", "FFT9", "FFT9h", "FT10", "FT10h", "FT7", "FT7h", "FT8", "FT8h", "FT9", "FT9h", "FTT10", "FTT10h", "FTT7", "FTT7h", "FTT8", "FTT8h", "FTT9", "FTT9h", "Fp1", "Fp1h", "Fp2", "Fp2h", "Fpz", "Fz", "I1", "I1h", "I2", "I2h", "Iz", "LPA", "N1", "N1h", "N2", "N2h", "NAS", "NFp1", "NFp1h", "NFp2", "NFp2h", "NFpz", "Nz", "O1", "O1h", "O2", "O2h", "OI1", "OI1h", "OI2", "OI2h", "OIz", "Oz", "P1", "P10", "P10h", "P1h", "P2", "P2h", "P3", "P3h", "P4", "P4h", "P5", "P5h", "P6", "P6h", "P7", "P7h", "P8", "P8h", "P9", "P9h", "PO1", "PO10", "PO10h", "PO1h", "PO2", "PO2h", "PO3", "PO3h", "PO4", "PO4h", "PO5", "PO5h", "PO6", "PO6h", "PO7", "PO7h", "PO8", "PO8h", "PO9", "PO9h", "POO1", "POO10", "POO10h", "POO1h", "POO2", "POO2h", "POO3", "POO3h", "POO4", "POO4h", "POO5", "POO5h", "POO6", "POO6h", "POO7", "POO7h", "POO8", "POO8h", "POO9", "POO9h", "POOz", "POz", "PPO1", "PPO10", "PPO10h", "PPO1h", "PPO2", "PPO2h", "PPO3", "PPO3h", "PPO4", "PPO4h", "PPO5", "PPO5h", "PPO6", "PPO6h", "PPO7", "PPO7h", "PPO8", "PPO8h", "PPO9", "PPO9h", "PPOz", "Pz", "RPA", "T10", "T10h", "T7", "T7h", "T8", "T8h", "T9", "T9h", "TP10", "TP10h", "TP7", "TP7h", "TP8", "TP8h", "TP9", "TP9h", "TPP10", "TPP10h", "TPP7", "TPP7h", "TPP8", "TPP8h", "TPP9", "TPP9h", "TTP10", "TTP10h", "TTP7", "TTP7h", "TTP8", "TTP8h", "TTP9", "TTP9h"],
    x = [-0.1025, 0.5878, 0.5021, -0.0512, 0.1025, 0.0512, -0.2067, -0.1543, 0.2067, 0.1543, -0.3143, -0.26, 0.3143, 0.26, -0.427, -0.3699, 0.427, 0.3699, -0.5878, -0.5021, -0.1207, 0.7071, 0.6039, -0.0602, 0.1207, 0.0602, -0.2444, -0.182, 0.2444, 0.182, -0.3743, -0.3084, 0.3743, 0.3084, -0.5138, -0.4426, 0.5138, 0.4426, -0.7071, -0.6039, 0.0, -0.0803, 0.454, 0.3878, -0.0401, 0.0803, 0.0401, -0.1615, -0.1207, 0.1615, 0.1207, -0.2443, -0.2027, 0.2443, 0.2027, -0.3299, -0.2867, 0.3299, 0.2867, -0.454, -0.3878, 0.0, 0.0, -0.1584, -0.0787, 0.1584, 0.0787, -0.3249, -0.2401, 0.3249, 0.2401, -0.5095, -0.4142, 0.5095, 0.4142, -0.157, -0.078, 0.157, 0.078, -0.3219, -0.2379, 0.3219, 0.2379, -0.5042, -0.4102, 0.5042, 0.4102, 0.0, -0.1527, -0.0759, 0.1527, 0.0759, -0.3126, -0.2313, 0.3126, 0.2313, -0.4881, -0.3978, 0.4881, 0.3978, -0.1455, -0.0724, 0.1455, 0.0724, -0.2969, -0.22, 0.2969, 0.22, -0.4612, -0.377, 0.4612, 0.377, 0.0, 0.0, 0.0, -0.1349, 0.809, 0.691, -0.0672, 0.1349, 0.0672, -0.2744, -0.2038, 0.2744, 0.2038, -0.4234, -0.3474, 0.4234, 0.3474, -0.5879, -0.5032, 0.5879, 0.5032, -0.809, -0.691, -0.1527, -0.0759, 0.1527, 0.0759, -0.3126, -0.2313, 0.3126, 0.2313, -0.4881, -0.3978, 0.4881, 0.3978, -0.157, -0.078, 0.157, 0.078, -0.3219, -0.2379, 0.3219, 0.2379, -0.5042, -0.4102, 0.5042, 0.4102, 0.0, 0.0, -0.1455, -0.0724, 0.1455, 0.0724, -0.2969, -0.22, 0.2969, 0.22, -0.4612, -0.377, 0.4612, 0.377, 0.0, 0.891, 0.761, -0.6474, -0.5509, 0.6474, 0.5509, -0.891, -0.761, 0.9511, 0.8123, -0.691, -0.5852, 0.691, 0.5852, -0.9511, -0.8123, 0.9877, 0.8436, -0.7176, -0.6059, 0.7176, 0.6059, -0.9877, -0.8436, -0.2245, -0.1137, 0.2245, 0.1137, 0.0, 0.0, -0.309, -0.1564, 0.309, 0.1564, 0.0, -1.0, -0.309, -0.1564, 0.309, 0.1564, 0.0, -0.2639, -0.1336, 0.2639, 0.1336, 0.0, 0.0, -0.2245, -0.1137, 0.2245, 0.1137, -0.2639, -0.1336, 0.2639, 0.1336, 0.0, 0.0, -0.1349, 0.809, 0.691, -0.0672, 0.1349, 0.0672, -0.2744, -0.2038, 0.2744, 0.2038, -0.4234, -0.3474, 0.4234, 0.3474, -0.5879, -0.5032, 0.5879, 0.5032, -0.809, -0.691, -0.1025, 0.5878, 0.5021, -0.0512, 0.1025, 0.0512, -0.2067, -0.1543, 0.2067, 0.1543, -0.3143, -0.26, 0.3143, 0.26, -0.427, -0.3699, 0.427, 0.3699, -0.5878, -0.5021, -0.0803, 0.454, 0.3878, -0.0401, 0.0803, 0.0401, -0.1615, -0.1207, 0.1615, 0.1207, -0.2443, -0.2027, 0.2443, 0.2027, -0.3299, -0.2867, 0.3299, 0.2867, -0.454, -0.3878, 0.0, 0.0, -0.1207, 0.7071, 0.6039, -0.0602, 0.1207, 0.0602, -0.2444, -0.182, 0.2444, 0.182, -0.3743, -0.3084, 0.3743, 0.3084, -0.5138, -0.4426, 0.5138, 0.4426, -0.7071, -0.6039, 0.0, 0.0, 1.0, 1.0, 0.8541, -0.7266, -0.6128, 0.7266, 0.6128, -1.0, -0.8541, 0.9511, 0.8123, -0.691, -0.5852, 0.691, 0.5852, -0.9511, -0.8123, 0.891, 0.761, -0.6474, -0.5509, 0.6474, 0.5509, -0.891, -0.761, 0.9877, 0.8436, -0.7176, -0.6059, 0.7176, 0.6059, -0.9877, -0.8436],
    y = [0.5139, 0.809, 0.691, 0.5105, 0.5139, 0.5105, 0.5274, 0.5194, 0.5274, 0.5194, 0.5513, 0.5379, 0.5513, 0.5379, 0.5879, 0.5678, 0.5879, 0.5678, 0.809, 0.691, 0.4195, 0.7071, 0.6039, 0.4155, 0.4195, 0.4155, 0.4361, 0.4263, 0.4361, 0.4263, 0.4661, 0.4493, 0.4661, 0.4493, 0.5138, 0.4874, 0.5138, 0.4874, 0.7071, 0.6039, 0.4142, 0.6148, 0.891, 0.761, 0.6133, 0.6148, 0.6133, 0.621, 0.6174, 0.621, 0.6174, 0.6317, 0.6258, 0.6317, 0.6258, 0.6474, 0.6388, 0.6474, 0.6388, 0.891, 0.761, 0.6128, 0.5095, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -0.0804, -0.0791, -0.0804, -0.0791, -0.0857, -0.0825, -0.0857, -0.0825, -0.096, -0.0901, -0.096, -0.0901, -0.0787, -0.1616, -0.1591, -0.1616, -0.1591, -0.1718, -0.1657, -0.1718, -0.1657, -0.1912, -0.1802, -0.1912, -0.1802, -0.2445, -0.2412, -0.2445, -0.2412, -0.2587, -0.2503, -0.2587, -0.2503, -0.2852, -0.2702, -0.2852, -0.2702, -0.2401, -0.1584, 0.0, 0.3302, 0.5878, 0.5021, 0.3262, 0.3302, 0.3262, 0.3467, 0.3369, 0.3467, 0.3369, 0.3771, 0.3599, 0.3771, 0.3599, 0.427, 0.3992, 0.427, 0.3992, 0.5878, 0.5021, 0.1616, 0.1591, 0.1616, 0.1591, 0.1718, 0.1657, 0.1718, 0.1657, 0.1912, 0.1802, 0.1912, 0.1802, 0.0804, 0.0791, 0.0804, 0.0791, 0.0857, 0.0825, 0.0857, 0.0825, 0.096, 0.0901, 0.096, 0.0901, 0.0787, 0.1584, 0.2445, 0.2412, 0.2445, 0.2412, 0.2587, 0.2503, 0.2587, 0.2503, 0.2852, 0.2702, 0.2852, 0.2702, 0.2401, 0.454, 0.3878, 0.3299, 0.3048, 0.3299, 0.3048, 0.454, 0.3878, 0.309, 0.2639, 0.2245, 0.2057, 0.2245, 0.2057, 0.309, 0.2639, 0.1564, 0.1336, 0.1137, 0.1036, 0.1137, 0.1036, 0.1564, 0.1336, 0.691, 0.7176, 0.691, 0.7176, 0.7266, 0.3249, -0.9511, -0.9877, -0.9511, -0.9877, -1.0, 0.0, 0.9511, 0.9877, 0.9511, 0.9877, 1.0, 0.8123, 0.8436, 0.8123, 0.8436, 0.8541, 1.0, -0.691, -0.7176, -0.691, -0.7176, -0.8123, -0.8436, -0.8123, -0.8436, -0.8541, -0.7266, -0.3302, -0.5878, -0.5021, -0.3262, -0.3302, -0.3262, -0.3467, -0.3369, -0.3467, -0.3369, -0.3771, -0.3599, -0.3771, -0.3599, -0.427, -0.3992, -0.427, -0.3992, -0.5878, -0.5021, -0.5139, -0.809, -0.691, -0.5105, -0.5139, -0.5105, -0.5274, -0.5194, -0.5274, -0.5194, -0.5513, -0.5379, -0.5513, -0.5379, -0.5879, -0.5678, -0.5879, -0.5678, -0.809, -0.691, -0.6148, -0.891, -0.761, -0.6133, -0.6148, -0.6133, -0.621, -0.6174, -0.621, -0.6174, -0.6317, -0.6258, -0.6317, -0.6258, -0.6474, -0.6388, -0.6474, -0.6388, -0.891, -0.761, -0.6128, -0.5095, -0.4195, -0.7071, -0.6039, -0.4155, -0.4195, -0.4155, -0.4361, -0.4263, -0.4361, -0.4263, -0.4661, -0.4493, -0.4661, -0.4493, -0.5138, -0.4874, -0.5138, -0.4874, -0.7071, -0.6039, -0.4142, -0.3249, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -0.309, -0.2639, -0.2245, -0.2057, -0.2245, -0.2057, -0.309, -0.2639, -0.454, -0.3878, -0.3299, -0.3048, -0.3299, -0.3048, -0.454, -0.3878, -0.1564, -0.1336, -0.1137, -0.1036, -0.1137, -0.1036, -0.1564, -0.1336]
);


text:
"""
    ht_p2_p1

The famous driven-cavity benchmark is solved here with Hood-Taylor combination
of quadratic and linear triangles.

The formulation is the one derived in Reddy, Introduction to the finite element
method, 1993. Page 486 ff.
"""
module ht_p2_p1

using LinearAlgebra
using StaticArrays
using MeshCore: nrelations, nentities
using MeshSteward: T6block, T6toT3
using MeshSteward: Mesh, attach!, baseincrel, boundary
using MeshSteward: vselect, geometry, summary
using MeshSteward: vtkwrite
using Elfel.RefShapes: manifdim, manifdimv
using Elfel.FElements: FEH1_T6, FEH1_T3, refshape, Jacobian
using Elfel.FESpaces: FESpace, ndofs, setebc!, nunknowns, doftype
using Elfel.FESpaces: numberfreedofs!, numberdatadofs!, numberdofs!
using Elfel.FESpaces: scattersysvec!, makeattribute, gathersysvec!, edofcompnt
using Elfel.FESpaces: highestfreedofnum, highestdatadofnum
using Elfel.FEIterators: FEIterator, ndofsperel, elnodes, eldofs
using Elfel.FEIterators: jacjac
using Elfel.QPIterators: QPIterator, bfun, bfungrad, weight
using Elfel.Assemblers: SysmatAssemblerSparse, start!, finish!, assemble!
using Elfel.Assemblers: SysvecAssembler
using Elfel.LocalAssemblers: LocalMatrixAssembler, LocalVectorAssembler, init!
using UnicodePlots

mu = 0.25 # dynamic viscosity
A = 1.0 # length of the side of the square
N = 100;# number of subdivisions along the sides of the square domain

function genmesh()
    # Hood-Taylor pair of meshes is needed
    # This mesh will be for the velocities
    vmesh = Mesh()
    attach!(vmesh, T6block(A, A, N, N), "velocity")
    # This mesh will be used for the pressures. Notice that it needs to be
    # "compatible" with the velocity mesh in the sense that they need to share
    # the nodes at the corners of the triangles.
    pmesh = Mesh()
    attach!(pmesh, T6toT3(baseincrel(vmesh, "velocity")), "pressure")
    return vmesh, pmesh
end

function assembleK(Uxh, Uyh, Ph, tndof, mu)
    function integrateK!(ass, elits, qpits, mu)
        uxnedof, uynedof, pnedof = ndofsperel.(elits)
        kuxux = LocalMatrixAssembler(uxnedof, uxnedof, 0.0)
        kuyuy = LocalMatrixAssembler(uynedof, uynedof, 0.0)
        kuxuy = LocalMatrixAssembler(uxnedof, uynedof, 0.0)
        kuxp = LocalMatrixAssembler(uxnedof, pnedof, 0.0)
        kuyp = LocalMatrixAssembler(uynedof, pnedof, 0.0)
        for el in zip(elits...)
            uxel, uyel, pel = el
            init!(kuxux, eldofs(uxel), eldofs(uxel))
            init!(kuyuy, eldofs(uyel), eldofs(uyel))
            init!(kuxuy, eldofs(uxel), eldofs(uyel))
            init!(kuxp, eldofs(uxel), eldofs(pel))
            init!(kuyp, eldofs(uyel), eldofs(pel))
            for qp in zip(qpits...)
                uxqp, uyqp, pqp = qp
                Jac, J = jacjac(pel, pqp)
                JxW = J * weight(pqp)
                gradNp = bfungrad(pqp, Jac)
                gradNux = bfungrad(uxqp, Jac)
                gradNuy = bfungrad(uyqp, Jac)
                Np = bfun(pqp)
                for j in 1:uxnedof, i in 1:uxnedof
                    kuxux[i, j] += (mu * JxW) * (2 * gradNux[i][1] * gradNux[j][1] + gradNux[i][2] * gradNux[j][2])
                end
                for j in 1:uynedof, i in 1:uynedof
                    kuyuy[i, j] += (mu * JxW) * (gradNuy[i][1] * gradNuy[j][1] + 2 * gradNuy[i][2] * gradNuy[j][2])
                end
                for j in 1:uynedof, i in 1:uxnedof
                    kuxuy[i, j] += (mu * JxW) * (gradNux[i][1] * gradNuy[j][2])
                end
                for j in 1:pnedof, i in 1:uxnedof
                    kuxp[i, j] += (-JxW) * (gradNux[i][1] * Np[j])
                end
                for j in 1:pnedof, i in 1:uynedof
                    kuyp[i, j] += (-JxW) * (gradNuy[i][2] * Np[j])
                end
            end
            assemble!(ass, kuxux)
            assemble!(ass, kuxuy)
            assemble!(ass, transpose(kuxuy))
            assemble!(ass, kuyuy)
            assemble!(ass, kuxp)
            assemble!(ass, transpose(kuxp))
            assemble!(ass, kuyp)
            assemble!(ass, transpose(kuyp))
        end
        return ass
    end

    elits = (FEIterator(Uxh), FEIterator(Uyh), FEIterator(Ph))
    qargs = (kind = :default, npts = 3,)
    qpits = (QPIterator(Uxh, qargs), QPIterator(Uyh, qargs), QPIterator(Ph, qargs))
    ass = SysmatAssemblerSparse(0.0)
    start!(ass, tndof, tndof)
    @time integrateK!(ass, elits, qpits, mu)
    return finish!(ass)
end

function solve!(U, K, F, nu)
    @time KT = K * U
    @time U[1:nu] = K[1:nu, 1:nu] \ (F[1:nu] - KT[1:nu])
end

function run()
    vmesh, pmesh = genmesh()
    # Velocity spaces
    Uxh = FESpace(Float64, vmesh, FEH1_T6(), 1)
    Uyh = FESpace(Float64, vmesh, FEH1_T6(), 1)
    locs = geometry(vmesh)
    inflate = A / N / 100
    # Part of the boundary that is immovable
    boxes = [[0.0 A 0.0 0.0], [0.0 0.0 0.0 A], [A A 0.0 A]]
    for box in boxes
        vl = vselect(locs; box = box, inflate = inflate)
        for i in vl
            setebc!(Uxh, 0, i, 1, 0.0)
            setebc!(Uyh, 0, i, 1, 0.0)
        end
    end
    # The lid
    uxbar = 1.0
    box = [0.0 A A A]
    vl = vselect(locs; box = box, inflate = inflate)
    for i in vl
        setebc!(Uxh, 0, i, 1, uxbar)
        setebc!(Uyh, 0, i, 1, 0.0)
    end
    # Pressure space
    Ph = FESpace(Float64, pmesh, FEH1_T3(), 1)
    setebc!(Ph, 0, 1, 1, 0.0)
    # Number the degrees of freedom
    numberdofs!([Uxh, Uyh, Ph])
    @show tndof = ndofs(Uxh) + ndofs(Uyh) + ndofs(Ph)
    @show tnunk = nunknowns(Uxh) + nunknowns(Uyh) + nunknowns(Ph)
    # Assemble the coefficient matrix
    K = assembleK(Uxh, Uyh, Ph, tndof, mu)
    p = spy(K, canvas = DotCanvas)
    display(p)
    # Solve the system
    U = fill(0.0, tndof)
    gathersysvec!(U, [Uxh, Uyh, Ph])
    F = fill(0.0, tndof)
    solve!(U, K, F, tnunk)
    scattersysvec!([Uxh, Uyh, Ph], U)
    # Postprocessing
    makeattribute(Ph, "p", 1)
    makeattribute(Uxh, "ux", 1)
    makeattribute(Uyh, "uy", 1)
    vtkwrite("ht_p2_p1-p", baseincrel(pmesh), [(name = "p",), ])
    vtkwrite("ht_p2_p1-v", baseincrel(vmesh), [(name = "ux",), (name = "uy",)])
    true
end

end

ht_p2_p1.run()


text:
using ProjectTemplate
using Test

@testset "ProjectTemplate.jl" begin
    # Write your own tests here.
end


text:
using Highway
using JLD

Tf = 14400; T = 12; N = 9
t1 = 10800

Lanes = 2; Lanes_t = 3
A = fld(Tf-t1, T)+1
Sections = [Highway.Sections2_CR[i] for i in 1:2:length(Highway.Sections2_CR)]
S = length(Sections)

J_in = linspace(0.2/Lanes, 0.6/Lanes, 9); J = length(J_in)

peeds_Transition = Array{Float64}(S, Lanes, J)
Densities_Transition = Array{Float64}(S, Lanes, J)
Fluxes_Transition = Array{Float64}(S, Lanes, J)

D_Speeds_Transition = Array{Float64}(S, Lanes, J)
D_Densities_Transition = Array{Float64}(S, Lanes, J)
D_Fluxes_Transition = Array{Float64}(S, Lanes, J)

Fluxes_Transition1 = zeros(S, Lanes, J)
Fluxes_Transition2 = zeros(S, Lanes, J)

Speeds_Transition1 = zeros(S, Lanes, J)
Speeds_Transition2 = zeros(S, Lanes, J)

Densities_Transition1 = zeros(S, Lanes, J)
Densities_Transition2 = zeros(S, Lanes, J)

D_Fluxes_Transition1 = zeros(S, Lanes, J)
D_Fluxes_Transition2 = zeros(S, Lanes, J)

D_Speeds_Transition1 = zeros(S, Lanes, J)
D_Speeds_Transition2 = zeros(S, Lanes, J)

D_Densities_Transition1 = zeros(S, Lanes, J)
D_Densities_Transition2 = zeros(S, Lanes, J)

jldopen("S2-Data-Express.jld", "w") do file
  for  (num_j, j) in enumerate(J_in)

    v0 = zeros(S, Lanes, A)
    d0 = zeros(S, Lanes, A)
    f0 = zeros(S, Lanes, A)

    v1 = zeros(S, Lanes, A)
    d1 = zeros(S, Lanes, A)
    f1 = zeros(S, Lanes, A)

    v2 = zeros(S, Lanes, A)
    d2 = zeros(S, Lanes, A)
    f2 = zeros(S, Lanes, A)

    println("Start simulations J_in = $j")
    for n = 1:N

      Sense_2, S2 = Highway.Sense2()

      flujo_local0 = zeros(S, Lanes, A)
      velocidad_local_promedio0 = zeros(S, Lanes, A)

      flujo_local1 = zeros(S, Lanes, A)
      velocidad_local_promedio1 = zeros(S, Lanes, A)

      flujo_local2 = zeros(S, Lanes, A)
      velocidad_local_promedio2 = zeros(S, Lanes, A)

  	  Arr_Speeds0 = Array(Any, S, Lanes, A)
  	  Arr_Speeds1 = Array(Any, S, Lanes, A)
  	  Arr_Speeds2 = Array(Any, S, Lanes, A)

      for t = 1:A, k = 1:Lanes, s = 1:S
        Arr_Speeds0[s, k, t] = Int8[]
        Arr_Speeds1[s, k, t] = Int8[]
        Arr_Speeds2[s, k, t] = Int8[]
      end

      println("Allocation of arrays finished, start of simulation $n")

      for t = 0:Tf-1

        Merge_Left_Right!(S2[2], S2[1], 2)
        Merge_Right_Left!(S2[1], S2[2], 1)

        Merge_Left_Right!(S2[3], S2[2], 3)
        Merge_Right_Left!(S2[2], S2[3], 2)

        f = j*Lanes/P_ramp_S2[1][1]

        if rand() < j
          Insert_Vehicle!(S2[2].highway, 2t+2, f*P_ramp_S2[1][2]/Lanes)
        end
        if rand() < j
          Insert_Vehicle!(S2[3].highway, 2t+1, f*P_ramp_S2[1][2]/Lanes)
        end

        for (num_ramp, (x0, lramp)) in enumerate(S2_in_ramp_ampli)
            Ramp!(x0, lramp, f*P_ramp_S2[1+num_ramp][1], f*P_ramp_S2[1+num_ramp][2], S2[2], 1)
            num_ramp = x0 = lramp = 0
        end

        for k = 1:Lanes_t
          AccelerationNoiseS2(S2[k].highway)
          DecelerationMove(S2[k])
          if t > t1 && k > 1
            Measure_Fluxes!(S1[k].highway, div(t-t1+1, T)+1, Sections, T,
                        slice(flujo_local0, :, k-1, :), slice(flujo_local1, :, k-1, :),
                        slice(flujo_local2, :, k-1, :))
            Measure_Speeds!(S1[k].highway, div(t-t1+1, T)+1, Sections, T,
                        slice(Arr_Speeds0, :, k-1, :), slice(Arr_Speeds1, :, k-1, :),
                        slice(Arr_Speeds2, :, k-1, :))
          end
        end
      end

      ##### Management of Data
      for t = 1:A, k = 1:Lanes, s = 1:S
        velocidad_local_promedio0[s, k, t] = mean(Arr_Speeds0[s, k, t])
        velocidad_local_promedio1[s, k, t] = mean(Arr_Speeds1[s, k, t])
        velocidad_local_promedio2[s, k, t] = mean(Arr_Speeds2[s, k, t])
      end

      Arr_Speeds0 = Arr_Speeds1 = Arr_Speeds2 = 0

	    Measure_Densities!(f0, v0, d0, A, Lanes, S, N, flujo_local0, velocidad_local_promedio0)
	    Measure_Densities!(f1, v1, d1, A, Lanes, S, N, flujo_local1, velocidad_local_promedio1)
	    Measure_Densities!(f2, v2, d2, A, Lanes, S, N, flujo_local2, velocidad_local_promedio2)

      println("ending of simulation $n")
      flujo_local0 = velocidad_local_promedio0 = flujo_local1 = 0
      velocidad_local_promedio1 = flujo_local2 = velocidad_local_promedio2 = 0
    end

  	writing_Arrays!(Lanes, S, num_j, Fluxes_Transition, Fluxes_Transition1, Fluxes_Transition2,
    			 D_Fluxes_Transition, D_Fluxes_Transition1, D_Fluxes_Transition2, f0, f1, f2)

  	writing_Arrays!(Lanes, S, num_j, Speeds_Transition, Speeds_Transition1, Speeds_Transition2,
    			 D_Speeds_Transition, D_Speeds_Transition1, D_Speeds_Transition2, v0, v1, v2)

  	writing_Arrays!(Lanes, S, num_j, Densities_Transition, Densities_Transition1,
    			 Densities_Transition2, D_Densities_Transition, D_Densities_Transition1,
    			 D_Densities_Transition2, d0, d1, d2)

    f0 = f1 = f2 = v0 = v1 = v2 = d0 = d1 = d2= 0
    println("Ends simulation J_in =  $j")
  end
  println("All simulations over, arrays writing starts")
  write(file, "Fluxes_Transition", Fluxes_Transition)
  write(file, "Fluxes_Transition1", Fluxes_Transition1)
  write(file, "Fluxes_Transition2", Fluxes_Transition2)
  write(file, "D_Fluxes_Transition", D_Fluxes_Transition)
  write(file, "D_Fluxes_Transition1", D_Fluxes_Transition1)
  write(file, "D_Fluxes_Transition2", D_Fluxes_Transition2)

  write(file, "Speeds_Transition", Speeds_Transition)
  write(file, "Speeds_Transition1", Speeds_Transition1)
  write(file, "Speeds_Transition2", Speeds_Transition2)
  write(file, "D_Speeds_Transition", D_Speeds_Transition)
  write(file, "D_Speeds_Transition1", D_Speeds_Transition1)
  write(file, "D_Speeds_Transition2", D_Speeds_Transition2)

  write(file, "Densities_Transition", Densities_Transition)
  write(file, "Densities_Transition1", Densities_Transition1)
  write(file, "Densities_Transition2", Densities_Transition2)
  write(file, "D_Densities_Transition", D_Densities_Transition)
  write(file, "D_Densities_Transition1", D_Densities_Transition1)
  write(file, "D_Densities_Transition2", D_Densities_Transition2)
end
println("The end")

max_stars_repo_name: LeonardoCastro/Cuernavaca

max_issues_repo_name: LeonardoCastro/Cuernavaca

max_forks_repo_name: LeonardoCastro/Cuernavaca


text:
module ExtrusionPolytopesTests

using Test
using Gridap.Helpers
using Gridap.TensorValues
using Gridap.Arrays
using Gridap.Fields
using Gridap.ReferenceFEs

@test num_faces(VERTEX) == 1
@test num_dims(VERTEX) == 0

p = Polytope(HEX_AXIS, HEX_AXIS)
test_polytope(p,optional=true)

@test get_faces(p,0,0) == [[1], [2], [3], [4]]
@test get_faces(p,1,0) == [[1, 2], [3, 4], [1, 3], [2, 4]]
@test get_faces(p,2,0) == [[1, 2, 3, 4]]
@test get_faces(p,0,1) == [[1, 3], [1, 4], [2, 3], [2, 4]]
@test get_faces(p,1,1) == [[1], [2], [3], [4]]
@test get_faces(p,2,1) == [[1, 2, 3, 4]]
@test get_faces(p,0,2) == [[1], [1], [1], [1]]
@test get_faces(p,1,2) == [[1], [1], [1], [1]]
@test get_faces(p,2,2) == [[1]]
@test get_facedims(p) == [0,0,0,0,1,1,1,1,2]
@test get_offsets(p) == [0,4,8]
@test num_facets(p) == 4
@test num_edges(p) == 4

r = Point{2,Float64}[(1, 0), (1, 0), (0, 1), (0, 1)]
@test get_edge_tangent(p) == r

r = [
  [1, 2, 3, 4], [1, 3, 2, 4], [2, 1, 4, 3], [2, 4, 1, 3],
  [3, 1, 4, 2], [3, 4, 1, 2], [4, 2, 3, 1], [4, 3, 2, 1]]
@test get_vertex_permutations(p) == r

r = Point{2,Float64}[(0, -1), (0, 1), (-1, 0), (1, 0)]
@test get_facet_normal(p) == r

@test num_faces(p) == 9
@test num_dims(p) == 2

p = Polytope(TET_AXIS, TET_AXIS, TET_AXIS)
test_polytope(p,optional=true)

@test num_faces(p) == 15
@test num_dims(p) == 3

x = Point{3,Float64}[(0, 0, 0), (1, 0, 0), (0, 1, 0), (0, 0, 1)]
@test get_vertex_coordinates(p) == x

p = SEGMENT
test_polytope(p,optional=true)
@test get_vertex_coordinates(p) == VectorValue{1,Float64}[(0),(1)]
@test get_edge_tangent(p) == VectorValue{1,Float64}[(1)]
@test get_facet_normal(p) == VectorValue{1,Float64}[(-1),(1)]
perm = get_vertex_permutations(p)
@test perm == [[1, 2], [2, 1]]

p = VERTEX
test_polytope(p,optional=true)
@test get_vertex_coordinates(p) == VectorValue{0,Float64}[()]
@test get_edge_tangent(p) == VectorValue{0,Float64}[]
@test get_facet_normal(p) == VectorValue{0,Float64}[]
perm = get_vertex_permutations(p)
@test perm == [[1]]

test_polytope(TRI,optional=true)
test_polytope(QUAD,optional=true)
test_polytope(TET,optional=true)
test_polytope(HEX,optional=true)

perm = get_vertex_permutations(TRI)
@test perm == [[1, 2, 3], [1, 3, 2], [2, 1, 3], [2, 3, 1], [3, 1, 2], [3, 2, 1]]

perm = get_vertex_permutations(QUAD)
@test isa(perm,Vector{Vector{Int}})

@test num_facets(SEGMENT) == 2
@test num_facets(TRI) == 3
@test num_facets(QUAD) == 4
@test num_facets(TET) == 4
@test num_facets(HEX) == 6

@test num_vertices(PYRAMID) == 5
@test num_edges(PYRAMID) == 8
@test num_facets(PYRAMID) == 5

@test is_simplex(TRI) == true
@test is_n_cube(TRI) == false

@test is_simplex(WEDGE) == false
@test is_n_cube(WEDGE) == false

@test is_simplex(QUAD) == false
@test is_n_cube(QUAD) == true

@test is_simplex(VERTEX) == true
@test is_n_cube(VERTEX) == true

d = 1
reffaces = get_reffaces(Polytope{d},WEDGE)
iface_to_ftype = get_face_type(WEDGE,d)
@test length(reffaces) == 1
@test iface_to_ftype == [1, 1, 1, 1, 1, 1, 1, 1, 1]

@test get_face_vertex_permutations(QUAD,1) == [[[1,2],[2,1]],[[1,2],[2,1]],[[1,2],[2,1]],[[1,2],[2,1]]]
@test get_face_vertex_permutations(QUAD,2) == [[[1,2,3,4],[1,3,2,4],[2,1,4,3],[2,4,1,3],[3,1,4,2],[3,4,1,2],[4,2,3,1],[4,3,2,1]]]

r = [
  [[1]],[[1]],[[1]],[[1]],
  [[1,2],[2,1]],[[1,2],[2,1]],[[1,2],[2,1]],[[1,2],[2,1]],
  [[1,2,3,4],[1,3,2,4],[2,1,4,3],[2,4,1,3],[3,1,4,2],[3,4,1,2],[4,2,3,1],[4,3,2,1]]]
@test get_face_vertex_permutations(QUAD) == r

r = [[[1,2],[2,1]],[[1,2],[2,1]],[[1,2],[2,1]]]
@test get_face_vertex_permutations(TRI,1) == r

@test get_face_vertex_permutations(HEX,3) == [[[1, 2, 3, 4, 5, 6, 7, 8]]]
@test get_face_vertex_permutations(TET,3) == [[[1, 2, 3, 4]]]

v,p = simplexify(QUAD)
@test p == TRI

v,p = simplexify(TRI)
@test p == TRI

v,p = simplexify(HEX)
@test p == TET

v,p = simplexify(TET)
@test p == TET

v,p = simplexify(SEGMENT)
@test p == SEGMENT
@test is_simplex(p)

v,p = simplexify(VERTEX)
@test p == VERTEX
@test is_simplex(p)

#using Gridap.Fields
#
#using .ReferenceFEs: _polytopenfaces
#using .ReferenceFEs: NFace
#using .ReferenceFEs: DFace
#using .ReferenceFEs: _polytopemesh
#using .ReferenceFEs: _dimfrom_fs_dimto_fs
#using .ReferenceFEs: _nface_to_nfacedim
#using .ReferenceFEs: _nfaces_vertices
#using .ReferenceFEs: _face_normals
#using .ReferenceFEs: _edge_tangents
#using .ReferenceFEs: _admissible_permutations
#using .ReferenceFEs: _admissible_permutations_n_cube
#
#using Profile
#using ProfileView
#
#
#extrusion = Point{0,Int}()
#p = DFace(extrusion)
#perms = _admissible_permutations(p)
#display(perms)
#
#extrusion = Point(HEX_AXIS)
#p = DFace(extrusion)
#perms = _admissible_permutations(p)
#display(perms)
#
#extrusion = Point(TET_AXIS)
#p = DFace(extrusion)
#perms = _admissible_permutations(p)
#display(perms)
#
#extrusion = Point(TET_AXIS,TET_AXIS)
#p = DFace(extrusion)
#perms = _admissible_permutations(p)
#display(perms)
#
#extrusion = Point(HEX_AXIS,HEX_AXIS)
#p = DFace(extrusion)
#perms = _admissible_permutations(p)
#display(perms)
#
#extrusion = Point(TET_AXIS, TET_AXIS, TET_AXIS)
#p = DFace(extrusion)
#perms = _admissible_permutations(p)
#display(perms)
#
#extrusion = Point(HEX_AXIS, HEX_AXIS, HEX_AXIS)
#p = DFace(extrusion)
#perms = _admissible_permutations(p)
#display(perms)
#
#
#f = DFace{2}(p,21)
#
#f = DFace{3}(p,27)

end # module


text:
# Chain rules to allow differentiable simulations

@non_differentiable check_force_units(args...)
@non_differentiable find_neighbors(args...)
@non_differentiable DistanceVecNeighborFinder(args...)
@non_differentiable run_loggers!(args...)
@non_differentiable visualize(args...)

function ChainRulesCore.rrule(T::Type{<:SVector}, vs::Number...)
    Y = T(vs...)
    function SVector_pullback(Ȳ)
        return NoTangent(), Ȳ...
    end
    return Y, SVector_pullback
end

function ChainRulesCore.rrule(T::Type{<:Atom}, vs...)
    Y = T(vs...)
    function Atom_pullback(Ȳ)
        return NoTangent(), Ȳ.index, Ȳ.charge, Ȳ.mass, Ȳ.σ, Ȳ.ϵ, Ȳ.solute
    end
    return Y, Atom_pullback
end

function ChainRulesCore.rrule(T::Type{<:SpecificInteraction}, vs...)
    Y = T(vs...)
    function SpecificInteraction_pullback(Ȳ)
        return NoTangent(), Ȳ...
    end
    return Y, SpecificInteraction_pullback
end

function ChainRulesCore.rrule(T::Type{<:PairwiseInteraction}, vs...)
    Y = T(vs...)
    function PairwiseInteraction_pullback(Ȳ)
        return NoTangent(), getfield.((Ȳ,), fieldnames(T))...
    end
    return Y, PairwiseInteraction_pullback
end

function ChainRulesCore.rrule(T::Type{<:SpecificForce2Atoms}, vs...)
    Y = T(vs...)
    function SpecificForce2Atoms_pullback(Ȳ)
        return NoTangent(), Ȳ.f1, Ȳ.f2
    end
    return Y, SpecificForce2Atoms_pullback
end

function ChainRulesCore.rrule(T::Type{<:SpecificForce3Atoms}, vs...)
    Y = T(vs...)
    function SpecificForce3Atoms_pullback(Ȳ)
        return NoTangent(), Ȳ.f1, Ȳ.f2, Ȳ.f3
    end
    return Y, SpecificForce3Atoms_pullback
end

function ChainRulesCore.rrule(T::Type{<:SpecificForce4Atoms}, vs...)
    Y = T(vs...)
    function SpecificForce4Atoms_pullback(Ȳ)
        return NoTangent(), Ȳ.f1, Ȳ.f2, Ȳ.f3, Ȳ.f4
    end
    return Y, SpecificForce4Atoms_pullback
end

function ChainRulesCore.rrule(::typeof(sparsevec), is, vs, l)
    Y = sparsevec(is, vs, l)
    @views function sparsevec_pullback(Ȳ)
        return NoTangent(), nothing, Ȳ[is], nothing
    end
    return Y, sparsevec_pullback
end

function ChainRulesCore.rrule(::typeof(accumulateadd), x)
    Y = accumulateadd(x)
    function accumulateadd_pullback(Ȳ)
        return NoTangent(), reverse(accumulate(+, reverse(Ȳ)))
    end
    return Y, accumulateadd_pullback
end

function ChainRulesCore.rrule(::typeof(unsafe_getindex), arr, inds)
    Y = unsafe_getindex(arr, inds)
    function unsafe_getindex_pullback(Ȳ)
        dx = Zygote._zero(arr, eltype(Ȳ))
        dxv = @view dx[inds]
        dxv .= Zygote.accum.(dxv, Zygote._droplike(Ȳ, dxv))
        return NoTangent(), Zygote._project(arr, dx), nothing
    end
    return Y, unsafe_getindex_pullback
end

# Not faster on CPU
function ChainRulesCore.rrule(::typeof(getindices_i), arr::CuArray, neighbors)
    Y = getindices_i(arr, neighbors)
    @views @inbounds function getindices_i_pullback(Ȳ)
        return NoTangent(), accumulate_bounds(Ȳ, neighbors.atom_bounds_i), nothing
    end
    return Y, getindices_i_pullback
end

function ChainRulesCore.rrule(::typeof(getindices_j), arr::CuArray, neighbors)
    Y = getindices_j(arr, neighbors)
    @views @inbounds function getindices_j_pullback(Ȳ)
        return NoTangent(), accumulate_bounds(Ȳ[neighbors.sortperm_j], neighbors.atom_bounds_j), nothing
    end
    return Y, getindices_j_pullback
end

# Required for SVector gradients in RescaleThermostat
function ChainRulesCore.rrule(::typeof(sqrt), x::Real)
    Y = sqrt(x)
    function sqrt_pullback(Ȳ)
        return NoTangent(), sum(Ȳ * inv(2 * Y))
    end
    return Y, sqrt_pullback
end

function ChainRulesCore.rrule(::typeof(reinterpret),
                                ::Type{T},
                                arr::SVector{D, T}) where {D, T}
    Y = reinterpret(T, arr)
    function reinterpret_pullback(Ȳ::Vector{T})
        return NoTangent(), NoTangent(), SVector{D, T}(Ȳ)
    end
    return Y, reinterpret_pullback
end

function ChainRulesCore.rrule(::typeof(reinterpret),
                                ::Type{T},
                                arr::AbstractArray{SVector{D, T}}) where {D, T}
    Y = reinterpret(T, arr)
    function reinterpret_pullback(Ȳ::Vector{T})
        return NoTangent(), NoTangent(), SVector{D, T}.(eachcol(reshape(Ȳ, D, length(Ȳ) ÷ D)))
    end
    return Y, reinterpret_pullback
end

function ChainRulesCore.rrule(::typeof(sum_svec), arr::AbstractArray{SVector{D, T}}) where {D, T}
    Y = sum_svec(arr)
    function sum_svec_pullback(Ȳ::SVector{D, T})
        return NoTangent(), zero(arr) .+ (Ȳ,)
    end
    return Y, sum_svec_pullback
end


text:
module Architectures

export AbstractArchitecture, AbstractMultiArchitecture
export CPU, GPU, MultiGPU
export device, device_event, architecture, array_type, arch_array, unified_array, device_copy_to!

using CUDA
using KernelAbstractions
using CUDAKernels
using Adapt
using OffsetArrays

# Adapt CUDAKernels to multiple devices by splitting stream pool
import CUDAKernels: next_stream

if CUDA.has_cuda_gpu()
using CUDAKernels: STREAM_GC_LOCK

    DEVICE_FREE_STREAMS = Tuple(CUDA.CuStream[] for dev in 1:length(CUDA.devices()))
    DEVICE_STREAMS      = Tuple(CUDA.CuStream[] for dev in 1:length(CUDA.devices()))
    const DEVICE_STREAM_GC_THRESHOLD = Ref{Int}(16)

    function next_stream()
        lock(STREAM_GC_LOCK) do
            handle = CUDA.device().handle + 1
            if !isempty(DEVICE_FREE_STREAMS[handle])
                return pop!(DEVICE_FREE_STREAMS[handle])
            end

            if length(DEVICE_STREAMS[handle]) > DEVICE_STREAM_GC_THRESHOLD[]
                for stream in DEVICE_STREAMS[handle]
                    if CUDA.query(stream)
                        push!(DEVICE_FREE_STREAMS[handle], stream)
                    end
                end
            end

            if !isempty(DEVICE_FREE_STREAMS[handle])
                return pop!(DEVICE_FREE_STREAMS[handle])
            end
            stream = CUDA.CuStream(flags = CUDA.STREAM_NON_BLOCKING)
            push!(DEVICE_STREAMS[handle], stream)
            return stream
        end
    end
end

"""
    AbstractArchitecture

Abstract supertype for architectures supported by Oceananigans.
"""
abstract type AbstractArchitecture end

"""
    AbstractMultiArchitecture

Abstract supertype for Distributed architectures supported by Oceananigans.
"""
abstract type AbstractMultiArchitecture <: AbstractArchitecture end

"""
    CPU <: AbstractArchitecture

Run Oceananigans on one CPU node. Uses multiple threads if the environment
variable `JULIA_NUM_THREADS` is set.
"""
struct CPU <: AbstractArchitecture end

"""
    GPU <: AbstractArchitecture

Run Oceananigans on a single NVIDIA CUDA GPU.
"""
struct GPU <: AbstractArchitecture end

#####
##### These methods are extended in Distributed.jl
#####

device(::CPU) = KernelAbstractions.CPU()
device(::GPU) = CUDAKernels.CUDADevice()

architecture() = nothing
architecture(::Number) = nothing
architecture(::Array) = CPU()
architecture(::CuArray) = GPU()
architecture(a::SubArray) = architecture(parent(a))
architecture(a::OffsetArray) = architecture(parent(a))

"""
    child_architecture(arch)

Return `arch`itecture of child processes.
On single-process, non-distributed systems, return `arch`.
"""
child_architecture(arch) = arch

array_type(::CPU) = Array
array_type(::GPU) = CuArray

arch_array(::CPU, a::Array)   = a
arch_array(::CPU, a::CuArray) = Array(a)
arch_array(::GPU, a::Array)   = CuArray(a)
arch_array(::GPU, a::CuArray) = a

arch_array(arch, a::AbstractRange) = a
arch_array(arch, a::OffsetArray) = OffsetArray(arch_array(arch, a.parent), a.offsets...)
arch_array(arch, ::Nothing)   = nothing
arch_array(arch, a::Number)   = a
arch_array(arch, a::Function) = a

unified_array(::CPU, a) = a
unified_array(::GPU, a) = a

function unified_array(::GPU, arr::AbstractArray)
    buf = Mem.alloc(Mem.Unified, sizeof(arr))
    vec = unsafe_wrap(CuArray{eltype(arr),length(size(arr))}, convert(CuPtr{eltype(arr)}, buf), size(arr))
    finalizer(vec) do _
        Mem.free(buf)
    end
    copyto!(vec, arr)
    return vec
end

## Only for contiguous data!! (i.e. only if the offset for pointer(dst::CuArrat, offset::Int) is 1)
@inline function device_copy_to!(dst::CuArray, src::CuArray; async::Bool = false)
    n = length(src)
    context!(context(src)) do
        GC.@preserve src dst begin
            unsafe_copyto!(pointer(dst, 1), pointer(src, 1), n; async)
        end
    end
    return dst
end

@inline device_copy_to!(dst::Array, src::Array; kw...) = Base.copyto!(dst, src)

device_event(arch) = Event(device(arch))

@inline unsafe_free!(a::CuArray) = CUDA.unsafe_free!(a)
@inline unsafe_free!(a)          = nothing

end # module