text:
function genfunc(
    funcname::Symbol, args::Vector{Symbol}, lhsind::Vector{Any}, rhs::Expr, opt::Ref{Expr}
)
    ex = if isempty(lhsind)
        :($funcname[] := $rhs)
    else
        :($funcname[$(lhsind...)] := $rhs)
    end

    if isassigned(opt)
        return :(@inline $funcname($(args...)) = @tensoropt $(opt[]) $ex)
    else
        return :(@inline $funcname($(args...)) = @tensor $ex)
    end
end

function genfrule(
    funcname::Symbol, args::Vector{Symbol}, lhsind::Vector{Any}, rhs::Expr, opt::Ref{Expr}
)
    ȧrgs, ṙhss = Symbol[], Expr[]
    for arg in args
        @gensym ȧrg
        ṙhs = make_scalar_first(make_only_product(rhs, arg))

        ṙhs = MacroTools.prewalk(ṙhs) do x  # assume to match only once
            @capture(x, $arg) ? ȧrg : x
        end

        push!(ȧrgs, ȧrg)
        push!(ṙhss, ṙhs)
    end
    ṙhs = :(+$(ṙhss...))

    @gensym val v̇al
    lhs = isempty(lhsind) ? :($v̇al[]) : :($v̇al[$(lhsind...)])

    pushbody = if isassigned(opt)
        :(@tensoropt $(opt[]) $lhs := $ṙhs)
    else
        :(@tensor $lhs := $ṙhs)
    end

    return quote
        function ChainRulesCore.frule((_, $(ȧrgs...)), ::typeof($funcname), $(args...))
            ($(ȧrgs...),) = ChainRulesCore.unthunk.(($(ȧrgs...),))
            $val = $(funcname)($(args...))
            $pushbody
            return ($val, $v̇al::typeof($val))
        end
    end
end

function genrrule(
    funcname::Symbol,
    args::Vector{Symbol},
    lhsind::Vector{Any},
    rhs::Expr,
    opt::Ref{Expr},
    isconjs::Vector{Bool},
    indsall::Vector{Union{Symbol,Expr}},
)
    indsall = [lhsind; indsall]

    @gensym Δlhs
    Δlhsex = isempty(lhsind) ? :(first($Δlhs)) : :($Δlhs[$(lhsind...)])

    projargs, projargexs = Symbol[], Expr[]
    ∂args, ∂argexs, ∂argdefs = Symbol[], Expr[], Expr[]
    for (arg, isconj) in zip(args, isconjs)
        @gensym ∂arg projarg
        push!(∂args, ∂arg)
        push!(projargs, projarg)
        push!(projargexs, :(ProjectTo($arg)))

        rhsarg = make_scalar_first(make_only_product(rhs, arg))

        ind = Ref{Vector{Any}}()
        ∂exrhs = MacroTools.prewalk(rhsarg) do x  # assume to match only once
            if @capture(x, $arg[_ind__])
                ind[] = _ind
                isconj ? Δlhsex : :(conj($Δlhsex))
            elseif @capture(x, $arg)
                isconj ? Δlhsex : :(conj($Δlhsex))
            else
                x
            end
        end
        @assert (isassigned(ind) && !isempty(ind[])) || !isassigned(ind)
        istensor = isassigned(ind)

        shouldtr = istensor ? ind[] ≠ unique(ind[]) : false
        indtr = Union{Symbol,Expr}[]

        if shouldtr
            δs = Expr[]
            for (k, i) in enumerate(ind[])
                if i ∉ indtr
                    push!(indtr, i)
                else
                    j = Expr(Symbol("'"), i)
                    while j ∈ [indsall; indtr]
                        j = Expr(Symbol("'"), j)
                    end
                    push!(indtr, j)
                    push!(
                        δs,
                        :(Array{eltype($arg)}(I, size($arg, $k), size($arg, $k))[$i, $j]),
                    )
                end
            end
            ∂exrhs = :(*($∂exrhs, $(δs...)))
        elseif istensor
            append!(indtr, ind[])
        end

        ∂exrhs = isconj ? ∂exrhs : :(conj($∂exrhs))

        @gensym ∂fn
        ∂exdef = if istensor
            if isassigned(opt)
                :(function $∂fn($Δlhs, $(args...))
                    @tensoropt $(opt[]) _[$(indtr...)] := $∂exrhs
                end)
            else
                :($∂fn($Δlhs, $(args...)) = @tensor _[$(indtr...)] := $∂exrhs)
            end
        else
            if isassigned(opt)
                :(function $∂fn($Δlhs, $(args...))
                    return first(@tensoropt $(opt[]) _[] := $∂exrhs)
                end)
            else
                :($∂fn($Δlhs, $(args...)) = first(@tensor _[] := $∂exrhs))
            end
        end
        push!(∂argdefs, ∂exdef)
        if istensor
            @gensym ∂fnadd!! inplaced
            ∂exadd! = if isassigned(opt)
                :(function $∂fnadd!!($inplaced, $Δlhs, $(args...))
                    @tensoropt $(opt[]) $inplaced[$(indtr...)] += $∂exrhs
                end)
            else
                :(function $∂fnadd!!($inplaced, $Δlhs, $(args...))
                    @tensor $inplaced[$(indtr...)] += $∂exrhs
                end)
            end
            push!(∂argdefs, ∂exadd!)
        end

        ∂exarg = if istensor
            :(
                $∂arg = InplaceableThunk(
                    $inplaced -> $∂fnadd!!($inplaced, $Δlhs, $(args...)),
                    Thunk(() -> $projarg($∂fn($Δlhs, $(args...)))),
                )
            )
        else
            :($∂arg = Thunk(() -> $∂fn($Δlhs, $(args...))))
        end

        push!(∂argexs, ∂exarg)
    end

    @gensym lhs pullback
    return quote
        function ChainRulesCore.rrule(::typeof($funcname), $(args...))
            $lhs = $(funcname)($(args...))
            ($(projargs...),) = ($(projargexs...),)
            function $(pullback)($Δlhs)
                return $(Expr(
                    :block,
                    :($Δlhs = Array(ChainRulesCore.unthunk($Δlhs))),
                    ∂argdefs...,
                    ∂argexs...,
                    :(return (ChainRulesCore.NoTangent(), $(∂args...))),
                ))
            end
            return ($lhs, $pullback)
        end
    end
end


text:
using Redis
using Blio

"""
    redis2grh(redis, key)

Return a GUPPI RAW header for HPGUPPI DAQ status buffer stored in redis hash
having key `key`.
"""
function redis2grh(redis, key)
    redishash = hgetall(redis, key)
    # Currently, Blio can only create a GuppiRaw.Header from a file, so we
    # have to create a GUPPI RAW header in a temporary file using key-value
    # pairs from the redis hash.
    mktemp() do path, grhio
        for (k,v) in redishash
            write(grhio, rpad(k,  8)) # 8 + 2 + 70 = 80
            write(grhio, "= ")
            write(grhio, rpad(v, 70))
        end
        write(grhio, rpad("END", 80))
        # Blio requires minimum file size :P
        write(grhio, " " ^ (80*256))
        seekstart(grhio)
        read(grhio, GuppiRaw.Header, skip_padding=false)
    end
end

"""
    hostinst2grh(redis, hostinst; domain="bluse")

Return a GUPPI RAW header for HPGUPPI DAQ status buffer stored in redis hash for
given `domain` and `hostinst`.  Retrieves redis hash for
`\$domain://\$hostinst/status`.
"""
function hostinst2grh(redis, hostinst; domain="bluse")
    redis2grh(redis, "$domain://$hostinst/status")
end

"""
    host2grh(redis, host, instance=0; domain="bluse")

Return a GUPPI RAW header for HPGUPPI DAQ status buffer stored in redis hash for
given `domain`, `host`, and `instance`.  Retrieves redis hash for
`\$domain://\$host/\$instance/status`.
"""
function host2grh(redis, hostname, instance=0; domain="bluse")
    hostinst2grh(redis, "$hostname/$instance", domain=domain)
end

"""
    subarray2grh(redis, subarray; domain="bluse")

Return a GUPPI RAW header for HPGUPPI DAQ status buffer stored in redis hash for
a host/instance that is present in `coordinator:allocated_hosts:\$subarray`.
First host/instance in the returned list for which a corresponding status hash
exists will be used.
"""
function subarray2grh(redis, subarray; domain="bluse")
    hostinsts = lrange(redis, "coordinator:allocated_hosts:$subarray", 0, -1)
    for hi in hostinsts
        grh = hostinst2grh(redis, hi, domain=domain)
        isempty(grh) || return grh
    end
    # Nothing found
    return GuppiRaw.Header()
end

ext: jl


text:
using GnuplotScripting
using Test

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


text:
capa = 10
c = Channel(capa)

function f(c)
    for i in 1:10
        put!(c, i)
    end
end

function g(c)
    for i in c
        sleep(0.5)
        println("From g: ", i)
    end
    println("Finish")
end

@async f(c)
@async g(c)


text:
#___  ____       ___ _____   _   _       _ _
#|  \/  | |     |_  /  __ \ | | | |     (_| |
#| .  . | |       | | /  \/ | | | |_ __  _| |_ ___
#| |\/| | |       | | |     | | | | '_ \| | __/ _ \
#| |  | | |___/\__/ | \__/\ | |_| | | | | | || (_) |
#_______\_____\____/ \____/  _____|_|___|_____\_____ _____ _____
#| ___ \        (_)         | | \ \ / / / __  |  _  / __  |  _  |
#| |_/ _ __ ___  _  ___  ___| |_ \ V /  `' / /| |/' `' / /| |/' |
#|  __| '__/ _ \| |/ _ \/ __| __|/   \    / / |  /| | / / |  /| |
#| |  | | | (_) | |  __| (__| |_/ /^\ \ ./ /__\ |_/ ./ /__\ |_/ /
#\_|  |_|  \___/| |\___|\___|\__\/   \/ \_____/\___/\_____/\___/
#              _/ |
#             |__/
#
# This code is part of the proposal of the team "MLJC UniTo" - University of Turin
# for "ProjectX 2020" Climate Change for AI.
# The code is licensed under MIT 3.0
# Please read readme or comments for credits and further information.

# Compiler: Julia 1.5

# Short description of this file: Numerical Resolution of Euler System

using ModelingToolkit
using Plots, PyPlot
using DifferentialEquations, DiffEqOperators, LinearAlgebra, OrdinaryDiffEq, RecursiveArrayTools

@parameters x y η t θ
@variables u1(..), u2(..), u3(..), u4(..), u5(..), u6(..), u7(..)
@derivatives Dt'~t
@derivatives Dx'~x
@derivatives Dy'~y
@derivatives Dη'~η
@derivatives Dxx''~x
@derivatives Dyy''~y

#functions
u  = u1(x,y,η,t,θ)
v  = u2(x,y,η,t,θ)
w  = u3(x,y,η,t,θ)
ω  = u4(x,y,η,t,θ)
θm = u5(x,y,η,t,θ)
qm = u6(x,y,η,t,θ)
z  = u7(x,y,η,t,θ)

#physical constants
g  = 9.81 #constant
γ  = 1.4 #air cp/cv=1.4
ρd = 1 #dry air density
p0 = 100000 #reference pressure 10^5 Pa
Rd = 8.31 #gas constant for dry air
ηc = 0.2 #page 8 of Advanced Research WRF V4((http://dx.doi.org/10.5065/1dfh-6p97)

#simulation inputs
sum_of_qs = 0.1 #read from file
#dx =
#dy =
#dz =
Δx = 1 #dummy
Δy = 1 #dummy

#scale factors
dist = x^2 + y^2 #dummy for distance on earth
Dxd  = 2*x #Dx(dist)
Dyd  = 2*y #Dy(dist)
mx   = Δx/dist
my   = Δy/dist

#definitions
αd = 1/ρd
p  = p0*(Rd*θm/(p0*αd))^γ
c1 = 2*ηc*ηc/((1 - ηc)*(1 - ηc)*(1 - ηc))
c2 = -ηc*(4 + ηc + ηc*ηc)/((1 - ηc)*(1 - ηc)*(1 - ηc))
c3 = 2*(1 + ηc + ηc*ηc)/((1 - ηc)*(1 - ηc)*(1 - ηc))
c4 = -(1 + ηc)/((1 - ηc)*(1 - ηc)*(1 - ηc))
B  = c1 + c2*η + c3*η*η + c4*η*η*η
ps = 10 #dummy
pt = 1 #dummy
pd = B*(ps - pt) + (η - B)*(p0 - pt) + pt
μd = (ps - pt)*(c2 + 2*c3*η + 3*c4*η*η) + (p0 - pt)*(1 - c2 - 2*c3*η - 3*c4*η*η) #[μd = Dη(pd)]
D2pd = (ps - pt)*(2*c3 + 6*c4*η) + (p0 - pt)*(-2*c3 - 6*c4*η) # (maybe find η as a function of pd)

#ω  = Dt(η)
U  = μd*u/my
V  = μd*v/mx
W  = μd*w/my
Ω  = μd*ω/my
Θm = μd*θm
Qm = μd*qm
ϕ  = g*z
α  = αd/(1 + sum_of_qs)

#geographic parametrization
Ωe = 1 #dummy #angular rotation rate of earth
ψ  = 1 #dummy #latitude on earth
re = 1000 #dummy #radius of earth
αr = 1 #dummy #local rotation angle between y-axis and the meridians
f  = 2*Ωe*sin(ψ)
e  = 2*Ωe*cos(ψ)

#external forcings
Dymx = -Δx/(dist*dist)*Dyd
Dxmy = -Δy/(dist*dist)*Dxd

# Initial and boundary conditions
bcs =       [u1(x,y,η,0,θ) ~ (x^2 + y^2 + η^2),
             u2(x,y,η,0,θ) ~ (x^2 + y^2 + η^2),
             u3(x,y,η,0,θ) ~ (x^2 + y^2 + η^2),
             u4(x,y,η,0,θ) ~ (x^2 + y^2 + η^2),
             u5(x,y,η,0,θ) ~ (x^2 + y^2 + η^2),
             u6(x,y,η,0,θ) ~ (x^2 + y^2 + η^2),
             u7(x,y,η,0,θ) ~ (x^2 + y^2 + η^2)]

PhyBcs =     []

# Space and time domains
domains = [x ∈ IntervalDomain(0.0,1.0),
           y ∈ IntervalDomain(0.0,1.0),
           η ∈ IntervalDomain(0.0,1.0),
           t ∈ IntervalDomain(0.0,1.0)]# Discretization

# Discretization
dx = 0.1; dy= 0.1; dη = 0.1; dt = 0.1

const N_ = 10

const X = reshape([i for i in 1:10 for j in 1:10 for k in 1:10], N_,N_,N_)
const Y = reshape([j for i in 1:10 for j in 1:10 for k in 1:10], N_,N_,N_)

const H = reshape([k for i in 1:10 for j in 1:10 for k in 1:10], N_,N_,N_)

const Mx = Tridiagonal([1.0 for i in 1:N_-1],[-2.0 for i in 1:N_],[1.0 for i in 1:N_-1])
const My = copy(Mx)
const Mη = copy(Mx)

# 1st order differentiation tensors
MMx = zeros(N_,N_,N_)

for i in 1:N_
    MMx[i, 1:end,1:end] = Mx
end

MMy = zeros(N_,N_,N_)

for i in 1:N_
    MMy[1:end, i, 1:end] = My
end

MMη = zeros(N_,N_,N_)

for i in 1:N_
    MMη[1:end, 1:end, i] = Mη
end

# Boundary conditions on derivatives, enforcing Neumann BC
Mx[2,1] = 2.0
Mx[end-1,end] = 2.0
My[1,2] = 2.0
My[end,end-1] = 2.0
Mη[1,2] = 2.0
Mη[end,end-1] = 2.0

const MMxu_1 = zeros(N_,N_,N_);
const MMyu_2 = zeros(N_,N_,N_);
const MMyu_1 = zeros(N_,N_,N_);
const MMηu_1 = zeros(N_,N_,N_);
const MMηu_4 = zeros(N_,N_,N_);
const MMxu_5 = zeros(N_,N_,N_);
const MMηu_5 = zeros(N_,N_,N_);
const MMxu_7 = zeros(N_,N_,N_);
const MMxu_2 = zeros(N_,N_,N_);
const MMηu_2 = zeros(N_,N_,N_);
const MMyu_5 = zeros(N_,N_,N_);
const MMyu_7 = zeros(N_,N_,N_);
const MMxu_3 = zeros(N_,N_,N_);
const MMyu_3 = zeros(N_,N_,N_);
const MMηu_3 = zeros(N_,N_,N_);

MMxu_1

# initial values
u0 = zeros(N_,N_,N_,7)

function ff(du,u,p,t)
  u_1 = @view  u[:,:,:,1]
  u_2 = @view  u[:,:,:,2]
  u_3 = @view  u[:,:,:,3]
  u_4 = @view  u[:,:,:,4]
  u_5 = @view  u[:,:,:,5]
  u_6 = @view  u[:,:,:,6]
  u_7 = @view  u[:,:,:,7]

  du_1 = @view du[:,:, :, 1]
  du_2 = @view du[:,:, :, 2]
  du_3 = @view du[:,:, :, 3]
  du_4 = @view du[:,:, :, 4]
  du_5 = @view du[:,:, :, 5]
  du_6 = @view du[:,:, :, 6]
  du_7 = @view du[:,:, :, 7]

  mul!(MMxu_1,MMx,u_1)
  mul!(MMyu_2,MMy,u_2)
  mul!(MMyu_1,MMy,u_1)
  mul!(MMηu_1,MMη,u_1)
  mul!(MMηu_4,MMη,u_4)
  mul!(MMxu_5,MMx,u_5)
  mul!(MMηu_5,MMη,u_5)
  mul!(MMxu_7,MMx,u_7)
  mul!(MMxu_2,MMx,u_2)
  mul!(MMηu_2,MMη,u_2)
  mul!(MMyu_5,MMy,u_5)
  mul!(MMyu_7,MMy,u_7)
  mul!(MMxu_3,MMx,u_3)
  mul!(MMyu_3,MMy,u_3)
  mul!(MMηu_3,MMη,u_3)
  mul!(MMηu_7,MMη,u_7)
  mul!(MMxu_6,MMx,u_6)
  mul!(MMyu_6,MMy,u_6)
  mul!(MMηu_6,MMη,u_6)

  #= VARIABLES MAPPING
  u  = u1(x,y,η,t,θ)
  v  = u2(x,y,η,t,θ)
  w  = u3(x,y,η,t,θ)
  ω  = u4(x,y,η,t,θ)
  θm = u5(x,y,η,t,θ)
  qm = u6(x,y,η,t,θ)
  z  = u7(x,y,η,t,θ)
  =#

  U  = μd*u_1/my
  V  = μd*u_2/mx
  W  = μd*u_3/my
  Ω  = μd*u_4/my
  Θm = μd*u_5
  Qm = μd*u_6
  ϕ  = g*u_7

  divVu_1 = μd*u_1/Δy*(u_1*Dxd + 2*dist*MMxu_1) + μd/Δx*(u_1*u_2*Dyd + dist*u_1*MMyu_2 + dist*u_2*MMyu_1) + dist/Δy*(D2pd*u_1*u_4 + μd*MMηu_1*u_4 + μd*u_1*MMηu_4)
  divVu_2 = μd/Δy*(u_1*u_2*Dxd + dist*u_1*MMxu_2 + dist*u_2*MMxu_1) + μd/Δx*(u_2*u_2*Dyd + dist*2*u_2*MMyu_2) + 1/my*(D2pd*u_3*u_2 + μd*MMηu_4*u_2 + μd*u_4*MMηu_2)
  divVu_3 = μd/Δy*(Dxd*u_1*u_3 + dist*MMxu_1*u_3 + dist*u_1*MMxu_3) + μd/Δx*(Dyd*u_2*u_3 + dist*MMyu_2*u_3 + dist*u_2*MMyu_3) + 1/my*(D2pd*u_4*u_3 + μd*MMηu_4*u_3 + μd*u_4*MMηu_3)
  divVu_5 = μd/Δy*(Dxd*u_1*u_5 + dist*MMxu_1*u_5 + dist*u_1*MMxu_5) + μd/Δx*(Dyd*u_2*u_5 + dist*MMyu_2*u_5 + dist*u_2*MMyu_5) + 1/my*(D2pd*u_4*u_5 + μd*MMηu_4*u_5 + μd*u_4*MMηu_5)
  divVu_6 = μd/Δy*(Dxd*u_1*u_6 + dist*MMxu_1*u_6 + dist*u_1*MMxu_6) + μd/Δx*(Dyd*u_2*u_6 + dist*MMyu_2*u_6 + dist*u_2*MMyu_6) + 1/my*(D2pd*u_4*u_6 + μd*MMηu_4*u_6 + μd*u_4*MMηu_6)

  innerV∇ϕ = U*g*MMxu_7 + V*g*MMyu_7 + Ω*g*MMηu_7

  FU  =  (mx/my)*(f + u_1*(my/mx)*Dymx - u_2*Dxmy)*V - (u_1/re + e*cos(αr))*W
  FV  = -(my/mx)*((f + u_1*(my/mx)*Dymx - u_2*Dxmy*U + (u_2/re - e*sin(αr))*W))
  FW  = e*(U*cos(αr) - (mx/my)*V*sin(αr) + 1/re*(u_1*U + (mx/my)*u_2*V))
  FΘm = 0 #dummy
  FQm = 0 #dummy

  @. du_1 = - my/μd*(divVu_1 + μd*α*p0*(Rd/(p0*αd))^γ*γ*u_5^(γ-1)*MMxu_5 + (α/αd)*p0*(Rd/(p0*αd))^γ*γ*u_5^(γ-1)*MMηu_5*g*MMxu_7 - FU)
  @. du_2 = - mx/μd*(divVu_2 + μd*α*p0*(Rd/(p0*αd))^γ*γ*u_5^(γ-1)*MMyu_5 + (α/αd)*p0*(Rd/(p0*αd))^γ*γ*u_5^(γ-1)*MMηu_5*g*MMyu_7 - FV)
  @. du_3 = - my/μd*(divVu_3 - g*((α/αd)*p0*(Rd/(p0*αd))^γ*γ*u_5^(γ-1)*MMηu_5 - μd) - FW)

  @. du_5 = - 1/μd*(divVu_5 - FΘm)    #coupling with SFIRE
  @. divV = 0
  @. du_7 = - 1/g*(1.0/μd)*(innerV∇ϕ - g*W)
  @. du_6 = - 1/μd*(divVu_6 - FQm)
end

prob = ODEProblem(ff,u0,(0.0,10.0))
sol = solve(prob,ROCK2(),progress=true,save_everystep=false,save_start=false)
#=
pde_system = PDESystem(eqs,bcs,domains,[x,y,η,t],[u1,u2,u3,u4,u5,u6,u7])

OrdinaryDiffEq.solve(pde_system)
#prob = OrdinaryDiffEq.PDEProblem(pde_system)

xs,ys,ηs,ts = [domain.domain.lower:dx:domain.domain.upper for (dx,domain) in zip([dx,dy,dη,dt],domains)]
#u_predict = [reshape([first(phi([x,y,η,t],res.minimizer)) for x in xs  for y in ys for η in ηs], (length(xs),length(ys),length(ηs)))  for t in ts ]
u_predict  = [reshape([phi([x,y,η,t],res.minimizer)[i] for x in xs for y in ys for η in ηs for t in ts],
              (length(xs),length(ys),length(ηs),length(ts))) for i in 1:7]

u_predict
=#


text:
using ACE.Random: rand_radial, rand_sphere
using ACE: ACEConfig, State

"""
`ZμRnYlm_1pbasis` : utility function to quickly generate a
`Zμ * Rn * Ylmn` 1-particle basis.

todo - finish the docs.
"""
function ZμRnYlm_1pbasis(; init = true, species = nothing, maxdeg = nothing,
                           maxL = maxdeg,
                           Bsel = ACE.SimpleSparseBasis(1, maxdeg),
                           kwargs...)
   RnYlm = ACE.Utils.RnYlm_1pbasis(; maxdeg=maxdeg, maxL=maxL, Bsel = Bsel, kwargs...)
   Zμ = Species1PBasis(species)
   B1p = Zμ * RnYlm
   if init
      ACE.init1pspec!(B1p, Bsel)
   end
   return B1p
end

_rand_atstate(mu0, Zμ, Rn) =
      AtomState{Float64}(mu = rand(Zμ),
                         mu0 = mu0,
                         rr = rand_radial(Rn) * rand_sphere() )

function rand_ACEConfig(B1p, Nat::Integer)
   @assert Set((:μ, :n, :l, :m)) == Set(ACE.symbols(B1p))
   Rn = B1p.bases[2]
   @assert Rn isa ACE.Rn1pBasis
   Zμ = B1p.bases[1]
   @assert Zμ isa Species1PBasis

   mu0 = rand(Zμ)
   Xs = [ _rand_atstate(mu0, Zμ, Rn) for _ = 1:Nat ]
   return ACEConfig(Xs)
end

"""
`PopZμRnYlm_1pbasis` : utility function to quickly generate a
`Pop * Zμ * Rn * Ylmn` 1-particle basis.
"""
function PopZμRnYlm_1pbasis(; init = true, species = nothing, maxdeg = nothing,
                           maxL = maxdeg,
                           Bsel = ACE.SimpleSparseBasis(1, maxdeg),
                           kwargs...)
   RnYlm = ACE.Utils.RnYlm_1pbasis(; maxdeg=maxdeg, maxL=maxL, Bsel = Bsel, kwargs...)
   Zμ = Species1PBasis(species)
   Pop = Pop1PBasis()
   B1p = Pop * Zμ * RnYlm
   if init
      ACE.init1pspec!(B1p, Bsel)
   end
   return B1p
end

_rand_atstate(mu0, Zμ, Rn, Pop) =
     PopAtomState{Float64}(mu = rand(Zμ), mu0 = mu0,
                  rr = rand_radial(Rn) * rand_sphere(),
                  pop = rand() )

function rand_ACEConfig_pop(B1p, Nat::Integer)
   @assert Set((:μ, :n, :l, :m, :P)) == Set(ACE.symbols(B1p))
   Rn = B1p.bases[3]
   @assert Rn isa ACE.Rn1pBasis
   Zμ = B1p.bases[2]
   @assert Zμ isa Species1PBasis
   Pop = B1p.bases[1]
   @assert Pop isa Pop1PBasis

   mu0 = rand(Zμ)
   Xs = [ _rand_atstate(mu0, Zμ, Rn, Pop) for _ = 1:Nat ]
   return ACEConfig(Xs)
end


text:
# same as `sum.jl`, but using the OrcJIT to compile the code and executing it via ccall.

using Test

using LLVM

if length(ARGS) == 2
    x, y = parse.([Int32], ARGS[1:2])
else
    x = Int32(1)
    y = Int32(2)
end

function codegen!(mod::LLVM.Module, name, tm)
    ctx = context(mod)
    param_types = [LLVM.Int32Type(ctx), LLVM.Int32Type(ctx)]
    ret_type = LLVM.Int32Type(ctx)

    triple!(mod, triple(tm))

    ft = LLVM.FunctionType(ret_type, param_types)
    sum = LLVM.Function(mod, name, ft)

    # generate IR
    Builder(ctx) do builder
        entry = BasicBlock(sum, "entry"; ctx)
        position!(builder, entry)

        tmp = add!(builder, parameters(sum)[1], parameters(sum)[2], "tmp")
        ret!(builder, tmp)
    end

    verify(mod)

    ModulePassManager() do pm
        add_library_info!(pm, triple(mod))
        add_transform_info!(pm, tm)
        run!(pm, mod)
    end

    verify(mod)
end

JIT = nothing

if LLVM.has_orc_v2()
    lljit = LLJIT(;tm=JITTargetMachine())

    ts_mod = ThreadSafeModule("jit")
    name = "sum_orc.jl"
    tm = JITTargetMachine()
    ts_mod(m->codegen!(m, name, tm))
    dispose(tm)

    jd = JITDylib(lljit)
    add!(lljit, jd, ts_mod)
    addr = lookup(lljit, name)

    @eval call_sum(x, y) = ccall($(pointer(addr)), Int32, (Int32, Int32), x, y)

    finalizer(LLVM.dispose, lljit)
    JIT = lljit
else
    Context() do ctx
        # Setup jit
        tm = JITTargetMachine()

        orc = OrcJIT(tm)
        register!(orc, GDBRegistrationListener())

        mod = LLVM.Module("jit"; ctx)
        name = mangle(orc, "sum_orc.jl")
        codegen!(mod, name, tm)

        # For debugging:
        #   asm = String(convert(Vector{UInt8}, emit(tm, mod, LLVM.API.LLVMAssemblyFile)))
        #   write(stdout, asm)

        jitted_mod = compile!(orc, mod)

        addr = address(orc, name)
        addr2 = addressin(orc, jitted_mod, name)
        @test addr == addr2
        @test addr.ptr != 0

        unregister!(orc, GDBRegistrationListener())

        # For debugging:
        #   ccall(:jl_breakpoint, Cvoid, (Any,), pointer(addr))
        # Then in GDB
        #   b *(*(uint64_t*)v)
        @eval call_sum(x, y) = ccall($(pointer(addr)), Int32, (Int32, Int32), x, y)

        JIT = orc
    end
end

GC.@preserve JIT begin
    @test call_sum(x, y) == x + y
end


text:
cd("/Users/adrienpajon/Documents/all.nosync/git_repo/pulse_analysis/2020-11-12 xp pulse all seuil/Julia_analysis")

using Glob
using JLD

include("functions/searchAndSave.jl")

#######
encoderFolderPath="../measures/pneu/encoder_file/"
encoderFilename="*.txt"
encoderSaveFilename="imports/import_filename_encoder_pneu.jld"
encoderSaveKey="filenameEncoder"
searchAndSave(encoderFolderPath,encoderFilename,encoderSaveFilename,encoderSaveKey)

#######
IBSFolderPath="../measures/pneu/IBS/"
IBSFilename="*.CSV"
IBSSaveFilename="imports/import_filename_IBS_pneu.jld"
IBSSaveKey="filenameIBS"
searchAndSave(IBSFolderPath,IBSFilename,IBSSaveFilename,IBSSaveKey)

#######
XsensorFolderPath="../measures/pneu/Xsensor/"
XsensorFilename="*.csv"
XsensorSaveFilename="imports/import_filename_Xsensor_pneu.jld"
XsensorSaveKey="filenameXsensor"
searchAndSave(XsensorFolderPath,XsensorFilename,XsensorSaveFilename,XsensorSaveKey)

#######
MTiXsensFolderPath="../measures/pneu/MTiXsens/"
NoMTiXsens="00B41DF6.txt"
#dossier = SHOU R - torse = PROP 1 - siege = uARM R - cuisse = uLEG R - tete = lLEG R
filenameMTiXsensID = (dossier="00B41DF6",torse="00B41E68",siege="00B41E6D",cuisse="00B41E77",tete="00B41E86")
MTiXsensSaveFilename="imports/import_filename_MTiXsens_pneu.jld"
MTiXsensSaveKey=["filenameMTiXsens" "filenameMTiXsensID" "sensorXsensNames"]
searchAndSaveMTiXsens(MTiXsensFolderPath,
                        NoMTiXsens,
                        filenameMTiXsensID,
                        MTiXsensSaveFilename,
                        MTiXsensSaveKey)


text:
################################################
# Approximate info-reward functions
################################################

function negentropy(p::AbstractArray{T}) where T<:Real
    s = zero(T)
    z = zero(T)
    for i = 1:length(p)
        @inbounds pi = p[i]
        if pi > z
            s += pi * log2(pi)
        end
    end
    return log2(length(p)) + s
end

using Distances, Plots; gr()
x = 0:0.01:1
y1 = [negentropy([x[i], 1 - x[i]]) for i = 1:length(x)]
y2 = [cityblock([x[i], 1 - x[i]], [0.5, 0.5]) for i = 1:length(x)]
y3 = [euclidean([x[i], 1 - x[i]], [0.5, 0.5]) for i = 1:length(x)]
y4 = [norm([x[i], 1 - x[i]] - [0.5, 0.5], Inf) for i = 1:length(x)]
plot(x, y1, label = "Negative Entropy", legend = :top,
    xlab = "Belief, P(State = 1)",
    ylab = "Belief Reward",
    linestyle = :solid,
    linewidth = 2)
plot!(x, y2, label = "1-Norm", linestyle = :dash,
    linecolor = :red, linealpha = 1, linewidth = 2)
plot!(x, y3, label = "2-Norm", linestyle = :dot,
    linecolor = :red, linealpha = 1, linewidth = 2)
plot!(x, y4, label = "Infinity-Norm", linestyle = :dashdot,
    linecolor = :red, linealpha = 0.6, linewidth = 2)
fn = string("belief_reward_function.pdf")
path = joinpath(homedir(),".julia\\v0.6\\RobustInfoPOMDP\\data\\figures\\",fn)
savefig(path)

x = 1:0.1:1
y = 1:0.1:1
f(x,y) = negentropy([x - y])
X = repmat(x',length(y),1)
Y = repmat(y,1,length(x))
Z = map(f,X,Y)
p1 = contour(x,y,f,fill=true)
p2 = contour(x,y,Z)
plot(p1,p2)

using ForwardDiff
dne(x) = ForwardDiff.gradient(negentropy, x)

function ev(i::Int, len::Int)
    v = zeros(len)
    v[i] = 1.0
    return v
end

e1 = ev(1,4)
negentropy(e1)
negentropy(c)
negentropy([0.2, 0.25, 0.25, 0.3])
dne([0.2, 0.25, 0.25, 0.3])

ints = 0.1:0.1:0.9
n = length(ints)
alphas = Vector{Vector{Float64}}(n)
for (i, int) in enumerate(ints)
    x = [int, 1 - int]
    m = dne(x)
    y = negentropy(x)
    b = y - dot(m,x)
    a1 = dot(m, [0,1]) + b
    a2 = dot(m, [1,0]) + b
    alphas[i] = [a1, a2]
end

plot(x, y1, linecolor = :black, legend = :bottomright)
plot([0,1], alphas, linealpha = 0.2, linecolor = :black)

y

dot(m,[0,1]) + b
dot(m,[1,0]) + b

ax = [-2.32, 0.848]
plot!([0,1],ax)

using RPOMDPModels
ip = CyberIPOMDP()

rmax = -Inf
b = fill(1/27, 27)
for α in ip.inforeward
    rmax = max(rmax, dot(α, b))
end
rmax

ip.inforeward

# complicated info reward function
ra21 = [1.0, 0.0]
ra22 = [0.9, 0.1]
ra23 = [0.8, 0.2]
ra24 = [0.7, 0.3]
ralphas2 = [ra21, ra22, ra23, ra24, 1-ra24, 1-ra23, 1-ra22, 1-ra21]
plot([0,1], ralphas2)

using Plots; gr()
using Distances

function tinfo(rinfo::Vector{Vector{Float64}}, b::Vector{Float64})
    rmax = -Inf
    for α in rinfo
        rmax = max(rmax, dot(α, b))
    end
    rmax
end

# simple info reward function
irs = [[1.0, -1/3, -1/3, -1/3],
         [-1/3, 1.0, -1/3, -1/3],
         [-1/3, -1/3, 1.0, -1/3],
         [-1/3, -1/3, -1/3, 1.0]]

c = fill(1/4, 4)
e1 = ev(1,4)
tinfo(irs, c)

# complicated info reward function
vhi = -10.0
vlo = 0.1
irc = [[1.0, vhi, vhi, vhi],
        [vhi, 1.0, vhi, vhi],
        [vhi, vhi, 1.0, vhi],
        [vhi, vhi, vhi, 1.0],
        [vlo, -vlo/3, -vlo/3, -vlo/3],
        [-vlo/3, vlo, -vlo/3, -vlo/3],
        [-vlo/3, -vlo/3, vlo, -vlo/3],
        [-vlo/3, -vlo/3, -vlo/3, vlo]]

tinfo(irc, c)
tinfo(irc, e1)
tinfo(irs, [0.5, 0.5, 0.0, 0.0])
tinfo(irc, [0.5, 0.5, 0.0, 0.0])
negentropy([0.5, 0.5, 0.0, 0.0])
norm([0.5, 0.5, 0.0, 0.0] - c, 1)
x = 0:0.001:1.0
y0 = [norm([xi, 1-xi, 0.0, 0.0] - c, 1) for xi in x]
y1 = [negentropy([xi, 1-xi, 0.0, 0.0]) for xi in x]
y2 = [tinfo(irs, [xi, 1-xi, 0.0, 0.0]) for xi in x]
y3 = [tinfo(irc, [xi, 1-xi, 0.0, 0.0]) for xi in x]
plot(x, y2, label = "Simple", legend = :bottomleft,
    xlab = "Belief, P(State = Bad | Position = 1)",
    ylab = "Expected Total Discounted Reward")
plot!(x, y3, label = "Complex")
plot!(x, y2)
plot!(x, y3)

max_stars_repo_name: ajkeith/RobustInfoPOMDP

max_issues_repo_name: ajkeith/RobustInfoPOMDP

max_forks_repo_name: ajkeith/RobustInfoPOMDP


text:
# # Selecting on custom events

# Julio.jl is built on top of Reagents.jl, a framework for writing programs with
# complex nonblocking algorithms and synchronizations.  In fact, `Julia.select`
# is a thin wrapper on top of the choice combinator `|` defined in Reagents.jl.
# Therefore, it is possible to define custom synchronization events using
# Reagents.jl API.

using Reagents
using Reagents: CAS, Computed, Return, Read, PostCommit
using Julio
using Julio: Events
using Test

# Let us define a "leaky" "broadcasting" channel; i.e., "broadcasting" in the
# sense that a `put!` can be `take!`n by multiple tasks and "leaky" in the sense
# that the items will be lost if there are no receiver tasks executing `take!`
# while the sender is executing `put!`.

struct BroadcastChannel{T,Receivers}
    eltype::Val{T}
    lck::Julio.Lock
    receivers::Receivers
end

function BroadcastChannel{T}() where {T}
    receivers = typeof(Julio.Promise{T}())[]
    lck = Julio.Lock()
    return BroadcastChannel(Val(T), lck, receivers)
end

# (We use a lock-based implementation to keep the example simple. It should also
# be possible to use some nonblocking algorithms.)
#
# The receiver requests an item simply by posting a `Julio.Promise`:
#
# (TODO: make it work without touching `Julio.Internal`)

Base.take!(bc::BroadcastChannel{T}) where {T} = Julio.Internal.apply(take!, bc)::T
function Julio.Internal.event(::typeof(take!), bc::BroadcastChannel{T}) where {T}
    (; lck, receivers) = bc
    p = Julio.Promise{T}()
    lock(lck) do
        filter!(isopen, receivers)
        push!(receivers, p)
    end
    return Reagents.WithNack() do nack
        # Cleanup in case this event is not selected:
        Reagents.dissolve(nack ⨟ PostCommit(_ -> close(p)); once = true)
        Return(nothing)
    end ⨟ Events.fetch(p)
end

# An item is sent to all the receivers registered at the time `put!` is called:

function Base.put!(bc::BroadcastChannel{T}, x) where {T}
    (; lck, receivers) = bc
    x = convert(T, x)
    lock(lck) do
        for p in receivers
            Julio.tryput!(p, x)  # tryput! instead of put! to ignore closed promises
        end
        empty!(receivers)
    end
end

#src # TODO: Running `unlock` potentially in other task is bad...?
#src function Julio.Internal.event(::typeof(put!), bc::BroadcastChannel{T}, x) where {T}
#src     (; lck, receivers) = bc
#src     x = convert(T, x)
#src     return Events.lock(lck) ⨟ PostCommit() do _
#src         for p in receivers
#src             Julio.tryput!(p, x)  # ignore closed promises
#src         end
#src         empty!(receivers)
#src         unlock(lck)
#src     end
#src end

# Demo:

function test_broadcastchannel(; ntasks = 4)
    bc = BroadcastChannel{Int}()
    done = Julio.Promise{Nothing}()
    Julio.withtaskgroup() do tg
        tasks = map(1:ntasks) do _
            Julio.spawn!(tg) do
                local items = Int[]
                while true
                    Julio.select(
                        (fetch, done) => Returns(true),
                        (take!, bc) => x -> begin
                            push!(items, x)
                            false
                        end,
                    ) && break
                end
                return items
            end
        end
        try
            for x in 1:2^10
                put!(bc, x)
            end
        finally
            done[] = nothing
        end
        for t in tasks
            items = fetch(t)
            @test all(>(0), diff(items))
        end
        return map(length ∘ fetch, tasks)  #src
    end
end


text:
# A minimalist parser combinator framework. Grammars may be specified using
# strings, regular expressions and the operators & and | (the formner with a
# couple of variants). For example, a simple subset of arithmetic might be
# described as follows:
#
#   number    = r"\d+"
#   product() = number  & ( "*" & product | ϵ )
#   sum()     = product & ( "+" & sum     | ϵ )
#
# Here ϵ represents the empty string. The sequential composition (&) skips
# over any whitespace; there is a version (&̲, i.e., &\underbar) that does not
# skip over any whitespace. This is useful for more specific whitespace
# matching, either for tight juxtaposition, for requiring whitespace (using
# the sp parser), or for matching horizontal/non-breaking space (nbsp if space
# is required, or nbspᵋ if it's optional).
#
# Note that product and sum are specified as functions rather than variables,
# because they access parsers that are not yet available. Such functions are
# automatically called to retrieve the actual parsers to use.
#
# Normally, matched regular expressions become part of the parse result, while
# strings are discarded; the parse result may be altered with the semantic
# action combinator (^), which wraps the result of a subparser, using a
# function or constructor. For example, one might specify
#
#   number ^ s->parse(Int, s)
#
# as part of a rule, to get an actual integer in the parse result, rather than
# its string representation. This functionality can be used to build up custom
# parse trees.

import Base: cmp, push!

const Wrapper = Union{Function, Type}
abstract type OpPolicy end

abstract type Cmb end
const Parser = Union{Regex, AbstractString, Function, Cmb}

# Parsing context
mutable struct Ctx
    grm     :: Parser
    ops     :: Union{OpPolicy, Nil}
    str     :: AbstractString
    idx     :: Int
    log     :: Vector{AbstractString}
end
Ctx(grm, ops, str) = Ctx(grm, ops, str, 1, [])

rewind!(c::Ctx, i::Int) = (c.idx = i; c)
addmsg!(c::Ctx, m) = push!(c.log, m)

function missed!(c::Ctx, p)
    s = c.str
    n = lastindex(s)
    i = min(c.idx, n)
    j = min(nextind(s, min(c.idx + 20, n)), n)
    txt = s[i:j]
    j < n && (txt = string(txt, "..."))
    msg = "$p at index $(c.idx): \"$txt\""
    addmsg!(c, msg)
end

clrmsg!(c::Ctx)    = empty!(c.log)
errmsg(c::Ctx)     = string("syntax error\n",
                     "The following attempted matches failed:\n",
                     join(unique(c.log), "\n"))

# Delegating any special handling of binary operators (used, e.g., in custom
# precedence-based parsing):
isop(x, c::Ctx) = isop(x, c.ops)
node(op, lhs, rhs, c::Ctx) = node(op, lhs, rhs, c.ops)

cmp(a, b, c::Ctx) = cmp(a, b, c.ops)
push!(c::Ctx, lt::Pair) = push!(c.ops, lt)

# Sequential parser combinator
struct Seq <: Cmb
    sub     :: Vector{Parser}
end

# Disjunctive parser combinator
struct Alt <: Cmb
    sub     :: Vector{Parser}
end

# Semantic action combinator: Pass result through callable wrapper
struct Wrap <: Cmb
    p       :: Parser
    w       :: Wrapper
end

function tryparse(parser::Parser, oppolicy::OpPolicy, input::AbstractString;
                  force=false)
    c = Ctx(parser, oppolicy, input)
    res = parse(parser, c)
    if res ≢ nil
        @assert length(res) == 1
        res = first(res)
    end
    (res ≡ nil && force) ? error(join(Parsers.errmsg(c))) : res
end
parse(parser::Parser, oppolicy::OpPolicy, input::AbstractString) =
    tryparse(parser, oppolicy, input, force=true)

parse(f::Function, c::Ctx) = parse(f(), c)

# The parse functions return lists of subexpressions, and this list is spliced
# into the current sequential results, resulting in somewhat flattened parse
# trees. Single subexpressions should thus be returned as singleton lists, and
# dropped results as empty lists.
#
# Wrappers used with the Wrap combinator can return nothing to trigger Wrap to
# return an empty list (rather than, say, a singleton list containing an empty
# list).

# Match a string a the current index, and drop the result
function parse(p::AbstractString, c::Ctx)
    if !startswith(SubString(c.str, c.idx), p)
        missed!(c, p)
        return nil
    end
    c.idx = c.idx + sizeof(p)
    [] # String matches are dropped
end

# Match a regex at the current index, and return the result
function parse(p::Regex, c::Ctx)
    m = match(p, c.str, c.idx)
    if m ≡ nil || m.offset ≠ c.idx
        missed!(c, p)
        return nil
    end
    c.idx = c.idx + sizeof(m.match)
    [m.match]
end

# Parse with subparsers in sequence and return concatenate results
function parse(p::Seq, c::Ctx)
    res = []
    for s in p.sub
        r = parse(s, c)
        r ≡ nil && return nil
        append!(res, r)
    end
    return res
end

# Try parsing with each subparser in turn; return the first successful result
function parse(p::Alt, c::Ctx)
    i = c.idx
    for s in p.sub
        r = parse(s, c)
        r ≢ nil && (clrmsg!(c); return r)
        rewind!(c, i)
    end
    return nil
end

# Parse with subparser and wrap the result using the wrapper
function parse(p::Wrap, c::Ctx)
    r = parse(p.p, c)
    r ≡ nil && return nil
    w = p.w(r, c)
    w ≡ nil ? [] : w isa AbstractVector ? w : [w]
end

## BNF-like syntactic sugar for grammar specifications

# Utility functions for flattening/n-ary parsing of the grammar specification
# itself: Wrap parsers of different kinds in a list, but splice in parsers of
# the same kind (e.g., sequential parsers being sequentially composed).
inseq(x::Seq) = x.sub
inseq(x)      = Parser[x]

inalt(x::Alt) = x.sub
inalt(x)      = Parser[x]

# Wrap a parser; parse result is passed through a function or constructor
(^)(p::Parser, w::Wrapper) = Wrap(p, w)

# Add a descriptor, for use in error messages
(..)(s::AbstractString, p::Parser) = Desc(p, s)

# Use as wrapper, to drop result
drop(x, c) = nil

# Some useful constants (dropped from parse result)
const ϵ        =  ""                         # Empty string (\epsilon)
const sp       =  r"\s+"           ^ drop    # Whitespace
const nbsp     =  r"[^\S\r\n]+"    ^ drop    # Non-breaking space
const spᵋ      =  r"\s*"           ^ drop    # Optional space
const nbspᵋ    =  r"[^\S\r\n]*"    ^ drop    # Optional horizontal space
const bof      =  r"^"             ^ drop    # Beginning of file
const eof      =  r"$"             ^ drop    # End of file
const eol      =  r"([\r\n]|\r\n)" ^ drop    # End of line

# Sequential parsing, whitespace permitted
(&)(l::Parser, r::Parser) = Seq([inseq(l)..., spᵋ, inseq(r)...])

# Sequential parsing, no whitespace permitted (&\underbar)
(&̲)(l::Parser, r::Parser) = Seq([inseq(l)...,      inseq(r)...])

# Alternative parsing
(|)(l::Parser, r::Parser) = Alt([inalt(l)...,      inalt(r)...])


text:
module BinnedPool

# binned memory pool allocator
#
# the core design is a pretty simple:
# - bin allocations into multiple pools according to their size (see `poolidx`)
# - when requested memory, check the pool for unused memory, or allocate dynamically
# - conversely, when released memory, put it in the appropriate pool for future use

using ..CUDA
import ..CUDA: actual_alloc, actual_free, PerDevice, initialize!

using ..PoolUtils

## tunables

const MAX_POOL = 2^27 # 128 MiB

const USAGE_WINDOW = 5

# min and max time between successive background task iterations.
# when the pool usages don't change, scan less regularly.
#
# together with USAGE_WINDOW, this determines how long it takes for objects to get reclaimed
const MIN_DELAY = 1.0
const MAX_DELAY = 5.0

## infrastructure

# TODO: npools
const pool_lock = ReentrantLock()
const pools_used = PerDevice{Vector{Set{Block}}}((dev)->Vector{Set{Block}}())
const pools_avail = PerDevice{Vector{Vector{Block}}}((dev)->Vector{Vector{Block}}())

poolidx(n) = Base.ceil(Int, Base.log2(n))+1
poolsize(idx) = 2^(idx-1)

@assert poolsize(poolidx(MAX_POOL)) <= MAX_POOL "MAX_POOL cutoff should close a pool"

function create_pools(dev, idx)
  if length(pools_avail[dev]) >= idx
    # fast-path without taking a lock
    return
  end

  @lock pool_lock begin
    while length(pools_avail[dev]) < idx
      push!(pools_used[dev], Set{Block}())
      push!(pools_avail[dev], Vector{Block}())
    end
  end
end

## pooling

const freed_lock = NonReentrantLock()
const freed = PerDevice{Vector{Block}}((dev)->Vector{Block}())

# reclaim unused buffers
function reclaim(dev, target_bytes::Int=typemax(Int))
  pool_repopulate(dev)

  @lock pool_lock begin
    @pool_timeit "reclaim" begin
      freed_bytes = 0

      # process pools in reverse, to discard largest buffers first
      for pid in reverse(1:length(pools_avail[dev]))
        bytes = poolsize(pid)
        avail = pools_avail[dev][pid]

        bufcount = length(avail)
        for i in 1:bufcount
          block = pop!(avail)

          actual_free(dev, block)

          freed_bytes += bytes
          if freed_bytes >= target_bytes
            return freed_bytes
          end
        end
      end

      return freed_bytes
    end
  end
end

# repopulate the "available" pools from the list of freed blocks
function pool_repopulate(dev)
  blocks = @lock freed_lock begin
    isempty(freed[dev]) && return
    blocks = Set(freed[dev])
    empty!(freed[dev])
    blocks
  end

  @lock pool_lock begin
    for block in blocks
      pid = poolidx(sizeof(block))

      @inbounds used = pools_used[dev][pid]
      @inbounds avail = pools_avail[dev][pid]

      # mark the buffer as available
      delete!(used, block)
      push!(avail, block)
    end
  end

  return
end

function alloc(dev, bytes)
  if bytes <= MAX_POOL
    pid = poolidx(bytes)
    create_pools(dev, pid)
    bytes = poolsize(pid)
  else
    pid = -1
  end

  block = nothing

  # NOTE: checking the pool is really fast, and not included in the timings
  @lock pool_lock begin
    if pid != -1 && !isempty(pools_avail[dev][pid])
      block = pop!(pools_avail[dev][pid])
    end
  end

  if block === nothing
    @pool_timeit "0. repopulate" begin
      pool_repopulate(dev)
    end

    @lock pool_lock begin
      if pid != -1 && !isempty(pools_avail[dev][pid])
        block = pop!(pools_avail[dev][pid])
      end
    end
  end

  if block === nothing
    @pool_timeit "1. try alloc" begin
      block = actual_alloc(dev, bytes)
    end
  end

  if block === nothing
    @pool_timeit "2a. gc (incremental)" begin
      GC.gc(false)
    end

    @pool_timeit "2b. repopulate" begin
      pool_repopulate(dev)
    end

    @lock pool_lock begin
      if pid != -1 && !isempty(pools_avail[dev][pid])
        block = pop!(pools_avail[dev][pid])
      end
    end
  end

  # TODO: we could return a larger allocation here, but that increases memory pressure and
  #       would require proper block splitting + compaction to be any efficient.

  if block === nothing
    @pool_timeit "3. reclaim" begin
      reclaim(dev, bytes)
    end

    @pool_timeit "4. try alloc" begin
      block = actual_alloc(dev, bytes)
    end
  end

  if block === nothing
    @pool_timeit "5a. gc (full)" begin
      GC.gc(true)
    end

    @pool_timeit "5b. repopulate" begin
      pool_repopulate(dev)
    end

    @lock pool_lock begin
      if pid != -1 && !isempty(pools_avail[dev][pid])
        block = pop!(pools_avail[dev][pid])
      end
    end
  end

  if block === nothing
    @pool_timeit "6. reclaim" begin
      reclaim(dev, bytes)
    end

    @pool_timeit "7. try alloc" begin
      block = actual_alloc(dev, bytes)
    end
  end

  if block === nothing
    @pool_timeit "8. reclaim everything" begin
      reclaim(dev, typemax(Int))
    end

    @pool_timeit "9. try alloc" begin
      block = actual_alloc(dev, bytes, true)
    end
  end

  if block !== nothing && pid != -1
    @lock pool_lock begin
      @inbounds used = pools_used[dev][pid]
      @inbounds avail = pools_avail[dev][pid]

      # mark the buffer as used
      push!(used, block)
    end
  end

  return block
end

function free(dev, block)
  # was this a pooled buffer?
  bytes = sizeof(block)
  if bytes > MAX_POOL
    actual_free(dev, block)
    return
  end

  # we don't do any work here to reduce pressure on the GC (spending time in finalizers)
  # and to simplify locking (preventing concurrent access during GC interventions)
  @spinlock freed_lock begin
    push!(freed[dev], block)
  end
end

function init()
  initialize!(freed, ndevices())

  initialize!(pools_used, ndevices())
  initialize!(pools_avail, ndevices())
end

function cached_memory(dev=device())
  sz = @lock freed_lock mapreduce(sizeof, +, freed[dev]; init=0)
  @lock pool_lock for (pid, pl) in enumerate(pools_avail[dev])
    bytes = poolsize(pid)
    sz += bytes * length(pl)
  end
  return sz
end

end


text:
#**************************************************************************************
# Output_Interface.jl
# =============== part of the GeoEfficiency.jl package.
#
# all the output either to the console or to the csv files is handled here.
#
#**************************************************************************************

const resultsfolder = "results";
const resultdir 	= joinpath(datadir, resultsfolder);		isdir(resultdir) 		|| mkdir(resultdir)
const resultdir_pnt = joinpath(resultdir, "Point");			isdir(resultdir_pnt) 	|| mkdir(resultdir_pnt)
const resultdir_nonPnt = joinpath(resultdir, "non-Point");	isdir(resultdir_nonPnt) || mkdir(resultdir_nonPnt)
countDetectors = 1;

"""
	calc(Detector::RadiationDetector = detectorFactory())

calculate the Geometrical Efficiency of the detector `Detector` and display it on the `console`.
If no detector is supplied it ask for a detector from the `console`.
Also prompt the user to input a source via the `console`.
"""
function calc(Detector::RadiationDetector = detectorFactory())
	global countDetectors
	aPnt, srcRadius, srcLength = source()
	print_with_color(:yellow,"\n\<$(countDetectors)\> $(id(Detector))")
	println("\n - Source(", id(aPnt), ", srcRadius=",srcRadius, ", srcLength=", srcLength, ")")
	try
		println("\n - The Detector Geometrical Efficiency = ", geoEff(Detector, aPnt, srcRadius, srcLength))

	catch err
		println(err)
		input("\n\t To Procced Press any Button")
		return nothing

	end #try
	countDetectors += 1
	print_with_color(:red, repeat(" =",36),"\n\n")
	return nothing

end #function

"""

	calcN()

calculate and display the Geometrical Efficiency.
Prompt the user to input a `detector` and a `source` from the `console`.
Prompt the user `repeatedly` until it exit (give a choice to use the same detector or a new detector).
"""
function calcN()
	Detector = detectorFactory()
	while (true)

		calc(Detector)
		res = input("""\n
    	I- To continue make a choice:
			> using the same detector Press 'd'|'D'
			> using a new detector Press 'n'|'N'\n
    	II- To quit just press return\n
			\n\tyour Choice: """, :blue)|> lowercase;
		if res == "n"
            Detector = detectorFactory()

		elseif res == "d"
            continue

		else
			break

		end #if

	end #while
	return nothing
end #function

function writecsv_head(fname::AbstractString, a, head=[])
	open(fname, "w") do io
		writedlm(io, head, ',')
		writedlm(io, a, ',')
	end
end

"""

	batch()

provide batch calculation of the Geometrical Efficiency based on the data provided from the csv files located in `$(datafolder)`.

Results are saved on a `csv file` named after the detector located in `$(resultdir)`, also a log of the results are displayed on the `console`.

Throw an error if the source location is inappropriate.
\n*****
"""
function batch()
	batch(read_batch_info()...)
end #function

"""

	batch( Detector_info_array::Array{Float64,2},
		srcHeights_array::Array{Float64,1},
			srcRhos_array::Array{Float64,1}=[0.0],
			srcRadii_array::Array{Float64,1}=[0.0],
			srcLengths_array::Array{Float64,1}=[0.0],
			ispoint::Bool=true)

provide batch calculation of the Geometricel efficiecny for each detector in the `Detector_info_array` after applying detectorFactory() to each raw.

A set of sources is constructed of every valid combination of parameter in the `srcRhos_array`, `srcRadii_array`, `srcLengths_array` with conjunction with `ispoint`.

If `ispoint` is true the source type is a point source and the parameters in srcRadii_array , srcLengths_array is completely ignored.
If `ispoint` is false the parameters in srcRhos_array is completely ignored.

Results are saved to a csv file named after the detector located in `$(resultdir)`, also a log of the results are displayed on the `console`.
\n*****
"""
function batch(	Detector_info_array::Array{Float64,2},
				srcHeights_array::Array{Float64,1},
				srcRhos_array::Array{Float64,1}=[0.0],
				srcRadii_array::Array{Float64,1}=[0.0],
				srcLengths_array::Array{Float64,1}=[0.0],
				ispoint::Bool=true)

	return  batch(	getDetectors(Detector_info_array),
					srcHeights_array,
					srcRhos_array,
					srcRadii_array,
					srcLengths_array,
					ispoint)
end #function

"""

	 batch( Detector::RadiationDetector,
			srcHeights_array::Array{Float64,1},
			srcRhos_array::Array{Float64,1}=[0.0],
			srcRadii_array::Array{Float64,1}=[0.0],
			srcLengths_array::Array{Float64,1}=[0.0],
			ispoint::Bool=true)

provide batch calculation of the Geometricel efficiecny for the detector `Detector` .

A set of sources is constructed of every valid combination of parameter in the `srcRhos_array`, `srcRadii_array`, `srcLengths_array` with conjunction with `ispoint`.

If `ispoint` is true the source type is a point source and the parameters in srcRadii_array , srcLengths_array is completely ignored.

If `ispoint` is false the parameters in srcRhos_array is completely ignored.

Results are saved to a csv file named after the detector located in `$(resultdir)`, also a log of the results are displayed on the `console`.

Return a tuple of the `detector array` and the `results array`. the `results array` has coloulmns `Height`, `Rho`, `GeoEfficiency` in the  case of a point while coloulmns `AnchorHeight`, `AnchorRho`, `srcRadius`, `srcLength`, `GeoEfficiency` for non-point sources.
\n*****
"""
function batch(	Detector::RadiationDetector,
				srcHeights_array::Array{Float64,1},
				srcRhos_array::Array{Float64,1}=[0.0],
				srcRadii_array::Array{Float64,1}=[0.0],
				srcLengths_array::Array{Float64,1}=[0.0],
				ispoint::Bool=true)

	return _batch(Val{ispoint},
				Detector::RadiationDetector,
				srcHeights_array,
				srcRhos_array,
				srcRadii_array,
				srcLengths_array
				)
end #function

"""

	 batch( Detectors_array::Array{RadiationDetector,1},
			srcHeights_array::Array{Float64,1},
			srcRhos_array::Array{Float64,1}=[0.0],
			srcRadii_array::Array{Float64,1}=[0.0],
			srcLengths_array::Array{Float64,1}=[0.0],
			ispoint::Bool=true)

provide batch calculation of the Geometricel efficiecny for each detector in the `Detectors_array`.

A set of sources is constructed of every valid combination of parameter in the `srcRhos_array`, `srcRadii_array`, `srcLengths_array` with conjunction with `ispoint`.

If `ispoint` is true the source type is a point source and the parameters in srcRadii_array , srcLengths_array is completely ignored.

If `ispoint` is false the parameters in srcRhos_array is completely ignored.

Results are saved to a csv file named after the detector located in `$(resultdir)`, also a log of the results are displayed on the `console`.
\n*****
"""
function batch(	Detectors_array::Array{RadiationDetector,1},
				srcHeights_array::Array{Float64,1},
				srcRhos_array::Array{Float64,1}=[0.0],
				srcRadii_array::Array{Float64,1}=[0.0],
				srcLengths_array::Array{Float64,1}=[0.0],
				ispoint::Bool=true)

	for Detector = Detectors_array
		batch(Detector,
			srcHeights_array,
			srcRhos_array,
			srcRadii_array,
			srcLengths_array,
			ispoint)

	end # Detectors_array

	input("\n\t the program had termiate, To Exit Press any Button")
	nothing

end #function

"""
batch calclulation for point sources.
"""
function _batch(::Type{Val{true}},
				Detector::RadiationDetector,
				srcHeights_array::Array{Float64,1},
				srcRhos_array::Array{Float64,1}=[0.0],
				srcRadii_array::Array{Float64,1}=[0.0],
				srcLengths_array::Array{Float64,1}=[0.0]
				)

	global countDetectors
	aPnt::Point = Point(0.0, 0.0)
	results::Array{Float64,2} = Array(Float64,(0,0)); out_results::Array{Float64,1} = Float64[];
	cellLabel = "\n\<$(countDetectors)\>$(id(Detector))"
	for srcHeight = srcHeights_array

		aPnt.Height = srcHeight		#setHeight!(aPnt, srcHeight)
		for srcRho = srcRhos_array
			aPnt.Rho = srcRho		#setRho!(aPnt, srcRho)

			try
				x = geoEff(Detector, aPnt)
				push!(out_results, aPnt.Height, aPnt.Rho, x)

			catch err
				isa(err, AssertionError) && break
				rethrow()

			end #try

			print_with_color(:yellow,cellLabel)
			println("\n - Source: ", id(aPnt))
			println("\n - The Detector Geometrical Efficiency = ", out_results[end])
			print_with_color(:red, repeat(" =",36),"\n\n")

		end #for_Rho

	end #for_Height

	results = reshape(out_results, 3, Int(length(out_results)/3)) |> transpose
	println("INFO: Saving <$countDetectors> to '$(id(Detector)).csv'......")
	try
		writecsv_head(joinpath(resultdir_pnt,  "$(id(Detector)).csv"), results, ["Height", "Rho", "GeoEfficiency"]')

	catch err
		warn("'.$(id(Detector)).csv': can't be created, results saved in an alternative file")
		writecsv_head(joinpath(resultdir_pnt, "_$(id(Detector)).csv"), results, ["Height", "Rho", "GeoEfficiency"]')

	end #try
	countDetectors += 1
	return (Detector, results)

end #function

"""
batch calclulation for non-point sources.
"""
function _batch(::Type{Val{false}},
				Detector::RadiationDetector,
				srcHeights_array::Array{Float64,1},
				srcRhos_array::Array{Float64,1}=[0.0],
				srcRadii_array::Array{Float64,1}=[0.0],
				srcLengths_array::Array{Float64,1}=[0.0];
				)

	global countDetectors
	aPnt::Point = Point(0.0, 0.0)
	results::Array{Float64,2} = Array(Float64,(0,0)); out_results::Array{Float64,1} = Float64[];
	cellLabel = "\n\<$(countDetectors)\>$(id(Detector))"
	for srcHeight = srcHeights_array

		aPnt.Height = srcHeight		#setHeight!(aPnt, srcHeight)
		for srcRho = srcRhos_array

			aPnt.Rho = srcRho		#setRho!(aPnt, srcRho)
			for  srcLength = srcLengths_array;

				for srcRadius = srcRadii_array

					try
						x = geoEff(Detector, aPnt, srcRadius , srcLength)
						push!(out_results, aPnt.Height, aPnt.Rho, srcRadius , srcLength, x)

					catch err
						isa(err, AssertionError) && @goto(Next_Height)
						rethrow()

					end #try

					print_with_color(:yellow,cellLabel)
					println("\n - Source[Anchor_", id(aPnt), ", srcRadius=",srcRadius, ", srcLength=", srcLength, "]")
					println("\n - The Detector Geometrical Efficiency = ", out_results[end])
					print_with_color(:red, repeat(" =",36),"\n\n")

				end #for_srcRadius

			@label(Next_srcLength)
			end #for_srcLength

		end #for_Rho

	@label(Next_Height)
	end #for_Height

	results = reshape(out_results, 5, Int(length(out_results)/5)) |> transpose
	println("INFO: Saving <$countDetectors> to '$(id(Detector)).csv'......")
	try err
		writecsv_head(joinpath(resultdir_nonPnt, "$(id(Detector)).csv"), results, ["AnchorHeight", "AnchorRho", "srcRadius", "srcLength", "GeoEfficiency"]')

	catch
		warn("'$(id(Detector)).csv': can't be created, results saved in an alternative file")
		writecsv_head(joinpath(resultdir_nonPnt, "_$(id(Detector)).csv"), results, ["AnchorHeight", "AnchorRho", "srcRadius", "srcLength", "GeoEfficiency"]')

	end #try
	countDetectors += 1;
	return (Detector, results)

end #function


text:
struct Account <: Value

    """
  1: dr drAccount cr crAccount
  increase drAccount , increase crAccount
  infer: there must be an incoming cash in-flow

  2: dr crAccount cr drAccount
  decrease drAccount , decrease crAccount
  infer: cash outflow

  """
    @enum Dr 1 2 # 1 for debit , 2 for Credit
    @enum cr 1 2
    #@enum accountOp(@enum Dr 1 2, @enum Cr 1 2)
    #@enum accountOp(@enum Dr 1, @enum Cr 2)
    ops = (+, -, *, /)
    """
    any Account can genrally add  & remove an mount (currency-agnostic)

    """
    id::Int64
    name::String
    bookValue::Double # these traits are not shared amongst all items in Accounting directory
    #marketValue:Double # viable, for publically traded items

    function add(amount)

        bookValue = bookValue + amount
        return bookValue
    end

    function subtract(amount)
        bookValue > amount ? bookValue = bookValue - amount : false

        bookValue = bookValue - amount, true = bookValue > amount || false
        return bookValue
    end

    drAcc

    function lookup(op)
        # vector syntax is discontinued # op={+, -, *,/}
        flag = false
        for i in enumerate(ops)
            op == ops[i] ? flag = true & break : break

            if (Condition = (op == ops[i]))
                flag = true
                ops[i]
                break
            end
        end
        return flag, ops[i]

    end

op =  lookup(op)

end

max_stars_repo_name: adamwillisXanax/accountingJulia

max_issues_repo_name: adamwillisXanax/accountingJulia

max_forks_repo_name: adamwillisXanax/accountingJulia


text:
@testset "Traits" begin
  @testset "Domain" begin
    # dummy type implementing the Domain trait
    struct DummyDomain{Dim,T} <: Domain{Dim,T}
      origin::Point{Dim,T}
    end
    function TopOptDash.element(domain::DummyDomain{Dim,T}, ind::Int) where {Dim,T}
      c = domain.origin + Vec(ntuple(i->T(ind), Dim))
      r = one(T)
      Ball(c, r)
    end
    TopOptDash.nelements(d::DummyDomain) = 3

    # basic properties
    dom = DummyDomain(P2(0,0))
    @test embeddim(dom) == 2
    @test coordtype(dom) == T

    # indexable/iterable interface
    dom = DummyDomain(P2(0,0))
    @test dom[begin] == Ball(P2(1,1), T(1))
    @test dom[end]   == Ball(P2(3,3), T(1))
    @test eltype(dom) <: Ball{2,T}
    @test length(dom) == 3
    @test collect(dom) == [Ball(P2(i,i), T(1)) for i in 1:3]

    # coordinates of centroids
    dom = DummyDomain(P2(1,1))
    pts = centroid.(Ref(dom), 1:3)
    @test pts == P2[(2,2), (3,3), (4,4)]

    dom = DummyDomain(P2(0,0))
    @test sprint(show, dom) == "3 DummyDomain{2,$T}"
    if T == Float32
      @test sprint(show, MIME"text/plain"(), dom) == "3 DummyDomain{2,Float32}\n  └─Ball{2,Float32}(Point(1.0f0, 1.0f0), 1.0))\n  └─Ball{2,Float32}(Point(2.0f0, 2.0f0), 1.0))\n  └─Ball{2,Float32}(Point(3.0f0, 3.0f0), 1.0))"
    elseif T == Float64
      @test sprint(show, MIME"text/plain"(), dom) == "3 DummyDomain{2,Float64}\n  └─Ball{2,Float64}(Point(1.0, 1.0), 1.0))\n  └─Ball{2,Float64}(Point(2.0, 2.0), 1.0))\n  └─Ball{2,Float64}(Point(3.0, 3.0), 1.0))"
    end
  end

  @testset "Data" begin
    dummydata(domain, table) = DummyData(domain, Dict(paramdim(domain) => table))
    dummymeta(domain, table) = meshdata(domain, Dict(paramdim(domain) => table))

    for (dummy, DummyType) in [(dummydata, DummyData), (dummymeta, MeshData)]
      # fallback constructor with spatial table
      dom = CartesianGrid{T}(2,2)
      tab = dummydata(dom, (a=[1,2,3,4], b=[5,6,7,8]))
      dat = DummyType(tab)
      @test domain(dat) == domain(tab)
      @test values(dat) == values(tab)

      # equality of data sets
      data₁ = dummy(CartesianGrid{T}(2,2), (a=[1,2,3,4], b=[5,6,7,8]))
      data₂ = dummy(CartesianGrid{T}(2,2), (a=[1,2,3,4], b=[5,6,7,8]))
      data₃ = dummy(PointSet(rand(P2,4)), (a=[1,2,3,4], b=[5,6,7,8]))
      @test data₁ == data₂
      @test data₁ != data₃
      @test data₂ != data₃

      # Tables interface
      dom = CartesianGrid{T}(2,2)
      dat = dummy(dom, (a=[1,2,3,4], b=[5,6,7,8]))
      @test Tables.istable(dat)
      @test Tables.rowaccess(dat)
      rows = Tables.rows(dat)
      schema = Tables.schema(rows)
      @test schema.names == (:a,:b,:geometry)
      @test schema.types == (Int, Int, Quadrangle{2,T,Vector{P2}})
      @test collect(rows) == [
        (a=1, b=5, geometry=dom[1]),
        (a=2, b=6, geometry=dom[2]),
        (a=3, b=7, geometry=dom[3]),
        (a=4, b=8, geometry=dom[4])
      ]
      @test collect(Tables.columns(dat)) == [
        [1,2,3,4],
        [5,6,7,8],
        [dom[1],dom[2],dom[3],dom[4]]
      ]
      @test Tables.materializer(dat) <: DummyType

      # Query interface
      dom = CartesianGrid{T}(2,2)
      dat = dummy(dom, (a=[1,2,3,4], b=[5,6,7,8]))
      new = dat |> @mutate(geometry=centroid(_.geometry)) |> DummyType
      @test domain(new) isa PointSet
      @test values(new) == values(dat)

      # column interface
      data = dummy(CartesianGrid{T}(2,2), (a=[1,2,3,4], b=[5,missing,7,8]))
      @test data[:a] == data["a"] == data.a == [1,2,3,4]
      @test isequal(data[:b], [5,missing,7,8])
      @test data[:geometry] == collect(CartesianGrid{T}(2,2))
      @test data[:geometry] == data["geometry"] == data.geometry
      @test_throws ErrorException data[:c]

      # variables interface
      data = dummy(CartesianGrid{T}(2,2), (a=[1,2,3,4], b=[5,missing,7,8]))
      @test variables(data) == (Variable(:a,Int), Variable(:b,Int))
      @test name.(variables(data)) == (:a,:b)
      @test mactype.(variables(data)) == (Int,Int)
      data = dummy(PointSet(rand(P2,4)), (a=[1,2,3,4], b=[5,6,7,8]))
      @test asarray(data, :a) == asarray(data, "a") == [1,2,3,4]
      @test asarray(data, :b) == asarray(data, "b") == [5,6,7,8]
      data = dummy(CartesianGrid{T}(2,2), (a=[1,2,3,4], b=[5,6,7,8]))
      @test asarray(data, :a) == asarray(data, "a") == [1 3; 2 4]
      @test asarray(data, :b) == asarray(data, "b") == [5 7; 6 8]

      data = dummy(CartesianGrid{T}(2,2), (a=[1,2,3,4], b=[5,6,7,8]))
      @test sprint(show, data) == "4 $DummyType{2,$T}"
      if T == Float32
        @test sprint(show, MIME"text/plain"(), data) == "4 $DummyType{2,Float32}\n  variables (rank 2)\n    └─a (Int64)\n    └─b (Int64)\n  domain: 2×2 CartesianGrid{2,Float32}"
      elseif T == Float64
        @test sprint(show, MIME"text/plain"(), data) == "4 $DummyType{2,Float64}\n  variables (rank 2)\n    └─a (Int64)\n    └─b (Int64)\n  domain: 2×2 CartesianGrid{2,Float64}"
      end

      if visualtests
        data = dummy(CartesianGrid{T}(2,2), (a=[1,2,3,4], b=[5,6,7,8]))
        @test_reference "data/$DummyType-$T.png" plot(data)
        data = dummydata(CartesianGrid{T}(2,2), (c=categorical([1,2,3,4]),))
        @test_reference "data/$DummyType-categorical-$T.png" plot(data)
      end
    end
  end
end


text:
using CartesianDDM
using Test

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


text:
module test_stdlib_printf

using Printf
using Test

@test Printf.@sprintf("%.7f", 1/11025) == "0.0000907"

end # module test_stdlib_printf


text:
import WinPath

@test WinPath.in()
@test WinPath.out()


text:
module RedirectRequest

using URIs
using ..Messages
using ..Pairs: setkv
import ..Header
import ..@debug, ..DEBUG_LEVEL

export redirectlayer, nredirects

"""
    redirectlayer(req) -> HTTP.Response

Redirects the request in the case of 3xx response status.
"""
function redirectlayer(handler)
    return function(req; redirect::Bool=true, redirect_limit::Int=3, forwardheaders::Bool=true, response_stream=nothing, kw...)
        if !redirect || redirect_limit == 0
            # no redirecting
            return handler(req; redirect_limit=redirect_limit, kw...)
        end

        count = 0
        while true
            # Verify the url before making the request. Verification is done in
            # the redirect loop to also catch bad redirect URLs.
            verify_url(req.url)
            res = handler(req; redirect_limit=redirect_limit, kw...)

            if (count == redirect_limit ||  !isredirect(res)
                ||  (location = header(res, "Location")) == "")
                return res
            end

            # follow redirect
            oldurl = req.url
            url = resolvereference(req.url, location)
            req = Request(req.method, resource(url), copy(req.headers), req.body;
                url=url, version=req.version, responsebody=response_stream, parent=res, context=req.context)
            if forwardheaders
                req.headers = filter(req.headers) do (header, _)
                    # false return values are filtered out
                    if header == "Host"
                        return false
                    elseif (header in SENSITIVE_HEADERS && !isdomainorsubdomain(url.host, oldurl.host))
                        return false
                    else
                        return true
                    end
                end
            else
                req.headers = Header[]
            end
            @debug 1 "➡  Redirect: $url"
            count += 1
        end
        @assert false "Unreachable!"
    end
end

function nredirects(req)
    return req.parent === nothing ? 0 : (1 + nredirects(req.parent.request))
end

const SENSITIVE_HEADERS = Set([
    "Authorization",
    "Www-Authenticate",
    "Cookie",
    "Cookie2"
])

function isdomainorsubdomain(sub, parent)
    sub == parent && return true
    endswith(sub, parent) || return false
    return sub[length(sub)-length(parent)] == '.'
end

function verify_url(url::URI)
    if !(url.scheme in ("http", "https", "ws", "wss"))
        throw(ArgumentError("missing or unsupported scheme in URL (expected http(s) or ws(s)): $(url)"))
    end
    if isempty(url.host)
        throw(ArgumentError("missing host in URL: $(url)"))
    end
end

end # module RedirectRequest


text:
using DistributedFactorGraphs
using IncrementalInference
using Manifolds
using StaticArrays
using Test

##

@testset "Test Sphere(2) prior and relative (broken)" begin
##

Base.convert(::Type{<:Tuple}, M::Sphere{2, ℝ}) = (:Euclid, :Euclid)
Base.convert(::Type{<:Tuple}, ::IIF.InstanceType{Sphere{2, ℝ}})  = (:Euclid, :Euclid)

@defVariable Sphere2 Sphere(2) [1.0, 0.0, 0.0]
M = getManifold(Sphere2)
@test M == Sphere(2)
pT = getPointType(Sphere2)
@test pT == Vector{Float64}
pϵ = getPointIdentity(Sphere2)
@test pϵ == [1.0, 0.0, 0.0]

@test is_point(getManifold(Sphere2), getPointIdentity(Sphere2))

fg = initfg()

v0 = addVariable!(fg, :x0, Sphere2)

mp = ManifoldPrior(Sphere(2), SA[1., 0, 0], MvNormal([0.01, 0.01]), DefaultOrthonormalBasis(), ExponentialRetraction())
p = addFactor!(fg, [:x0], mp)

doautoinit!(fg, :x0)

vnd = getVariableSolverData(fg, :x0)
@test_broken all(isapprox.(mean(vnd.val), [1,0,0], atol=0.1))
@test_broken all(is_point.(Ref(M), vnd.val))

v1 = addVariable!(fg, :x1, Sphere2)
mf = ManifoldFactor(Sphere(2), MvNormal([0.1, 0.2], [0.05,0.05]))
f = addFactor!(fg, [:x0, :x1], mf)

##
# Debugging Sphere error
smtasks = Task[]
solveTree!(fg; smtasks) #, verbose=true, recordcliqs=ls(fg))
# hists = fetchCliqHistoryAll!(smtasks);

##
end


text:
# This file is part of Intersections.jl

module Intersections
@inline _get_key(arr::T, i) where T = @inbounds arr[i]

@inline function _swap_items!(arr::T, i, j) where T
    @inbounds tmp = arr[i]
    @inbounds arr[i] = arr[j]
    @inbounds arr[j] = tmp
end

include("search.jl")
include("by.jl")
include("merge.jl")
include("imerge.jl")
include("svs.jl")
include("bk.jl")
end


text:
function iterate(m::Integer)
    while m != 1 && m != 89
        s = 0
        while m > 0 # compute sum of squares of digits
            m, d = divrem(m, 10)
            s += d ^ 2
        end
        m = s
    end
    return m
end
itercount(k::Integer) = count(x -> iterate(x) == 89, 1:k)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
#=
This file is auto-generated. Do not edit.
=#
"""
    mutable struct Transfer <: Service
        name::String
        contributingdevices::Vector{Device}
        timeframe::Float64
        requirement::TimeSeries.TimeArray
        _forecasts::InfrastructureSystems.Forecasts
        internal::InfrastructureSystemsInternal
    end

# Arguments
- `name::String`
- `contributingdevices::Vector{Device}`
- `timeframe::Float64`: the relative saturation timeframe
- `requirement::TimeSeries.TimeArray`
- `_forecasts::InfrastructureSystems.Forecasts`
- `internal::InfrastructureSystemsInternal`: power system internal reference, do not modify
"""
mutable struct Transfer <: Service
    name::String
    contributingdevices::Vector{Device}
    "the relative saturation timeframe"
    timeframe::Float64
    requirement::TimeSeries.TimeArray
    _forecasts::InfrastructureSystems.Forecasts
    "power system internal reference, do not modify"
    internal::InfrastructureSystemsInternal
end

function Transfer(name, contributingdevices, timeframe, requirement, _forecasts=InfrastructureSystems.Forecasts(), )
    Transfer(name, contributingdevices, timeframe, requirement, _forecasts, InfrastructureSystemsInternal())
end

function Transfer(; name, contributingdevices, timeframe, requirement, _forecasts=InfrastructureSystems.Forecasts(), )
    Transfer(name, contributingdevices, timeframe, requirement, _forecasts, )
end

# Constructor for demo purposes; non-functional.

function Transfer(::Nothing)
    Transfer(;
        name="init",
        contributingdevices=[ThermalStandard(nothing)],
        timeframe=0.0,
        requirement=[],
        _forecasts=InfrastructureSystems.Forecasts(),
    )
end

"""Get Transfer name."""
get_name(value::Transfer) = value.name
"""Get Transfer contributingdevices."""
get_contributingdevices(value::Transfer) = value.contributingdevices
"""Get Transfer timeframe."""
get_timeframe(value::Transfer) = value.timeframe
"""Get Transfer requirement."""
get_requirement(value::Transfer) = value.requirement
"""Get Transfer _forecasts."""
get__forecasts(value::Transfer) = value._forecasts
"""Get Transfer internal."""
get_internal(value::Transfer) = value.internal


text:
using LightGraphs
using SuperFastGraphs

g = loadgraph("test/testdata/kite.lg", "kite")
println(topKcc(g, 10))

max_stars_repo_name: mameli/SuperFastGraphs

max_issues_repo_name: mameli/SuperFastGraphs

max_forks_repo_name: mameli/SuperFastGraphs


text:
function run_TGV(parts,n::Int,order::Int,dt::Real,tf::Real)

  # Parameters
  L = 2π
  Re = 800.0
  ν = 1/Re

  # TGV initial solution
  u₀(x) = VectorValue(cos(x[1])*sin(x[2])*sin(x[3]),-sin(x[1])*cos(x[2])*sin(x[3]),0.0)
  p₀(x) = 1/16*(cos(2*x[1])+cos(2*x[2]))*(cos(2*x[3])+2)

  # Discretization
  domain = (0,L,0,L,0,L)
  cells = (n,n,n)
  𝒯 = CartesianDiscreteModel(parts,domain,cells;isperiodic=(true,true,true))

  # Triangulation and Integration measure
  Ω = Interior(𝒯)
  dΩ = Measure(Ω,2*order)

  # FE spaces
  refFEᵤ = ReferenceFE(lagrangian,VectorValue{3,Float64},order)
  refFEₚ = ReferenceFE(lagrangian,Float64,order-1)
  V = TestFESpace(𝒯,refFEᵤ)
  Q = TestFESpace(𝒯,refFEₚ)
  U = TransientTrialFESpace(V)
  P = TrialFESpace(Q)
  Y = MultiFieldFESpace([V,Q])
  X = TransientMultiFieldFESpace([U,P])

  # Explicit FE functions
  global ηₙₕ = interpolate(VectorValue(0.0,0.0,0.0),U(0.0))
  global uₙₕ = interpolate(u₀,U(0.0))
  global fv_u = zero_free_values(U(0.0))

  # Stabilization Parameters
  c₁ = 12.0
  c₂ = 2.0
  cc = 4.0
  h = L/(n*order)
  τₘ = 1/(c₁*ν/h^2 + c₂*(meas∘uₙₕ)/h)
  τc = cc *(h^2/(c₁*τₘ))

  # Weak form
  c(a,u,v) = 0.5*((∇(u)'⋅a)⋅v - u⋅(∇(v)'⋅a))
  res(t,(u,p),(v,q)) = ∫( ∂t(u)⋅v  + c(u,u,v) + ν*(∇(u)⊙∇(v)) - p*(∇⋅v) + (∇⋅u)*q +
                          τₘ*((∇(u)'⋅u - ηₙₕ)⋅(∇(v)'⋅u)) + τc*((∇⋅u)*(∇⋅v)) )dΩ
  jac(t,(u,p),(du,dp),(v,q)) = ∫( c(du,u,v) + c(u,du,v) + ν*(∇(du)⊙∇(v)) - dp*(∇⋅v) + (∇⋅du)*q +
                                  τₘ*((∇(u)'⋅u - ηₙₕ)⋅(∇(v)'⋅du) + (∇(du)'⋅u + ∇(u)'⋅du)⋅(∇(v)'⋅u)) +
                                  τc*((∇⋅du)*(∇⋅v)) )dΩ
  jac_t(t,(u,p),(dut,dpt),(v,q)) = ∫( dut⋅v )dΩ
  op = TransientFEOperator(res,jac,jac_t,X,Y)

  # Orthogonal projection
  a(η,κ) = ∫( τₘ*(η⋅κ) )dΩ
  b(κ) = ∫( τₘ*((∇(uₙₕ)'⋅uₙₕ)⋅κ) )dΩ
  op_proj = AffineFEOperator(a,b,U,V)
  ls = PETScLinearSolver()

  # Transient solver
  xₕ₀ = interpolate([u₀,p₀],X(0.0))
  nls = NLSolver(ls,show_trace=true,method=:newton,iterations=10)
  #nls = PETScNonlinearSolver()
  ode_solver = ThetaMethod(nls,dt,0.5)

  # Solution (lazy)
  xₕₜ = solve(ode_solver,op,xₕ₀,0,tf)

  K = Float64[]
  E = Float64[]
  T = Float64[]

  # Iterate over steps
  createpvd(parts,"results/TGV") do pvd
    for (xₕ,t) in xₕₜ
      println("--------------------")
      println("Start of time: $t")
      uₕ,pₕ=xₕ
      ωₕ = ∇×uₕ
      push!(K,0.5*(∑(∫(uₕ⋅uₕ)dΩ))/(2π^3))
      push!(E,ν*(∑(∫(ωₕ⋅ωₕ)dΩ))/(2π^3))
      push!(T,t)
      println("updating global variables")
      uₙₕ = interpolate!(uₕ,fv_u,U(t))
      ηₙₕ = solve(ls,op_proj)
      pvd[t] = createvtk(Ω,"results/TGV_$t";cellfields=["u"=>uₕ,"p"=>pₕ,"eta"=>ηₙₕ,"w"=>ωₕ],nsubcells=10)
      println("End of time: $t")
    end
 end

 return (T,K,E)
end


text:
"""
    Page

Abstract supertype for web pages and mock web pages.
All Page subtypes are assumed to define a `url::String`.
"""
abstract type Page end

"""
    WebPage(url::AbstractString; selector::Function=clean) <: Page

WebPage to be scanned with `url`.
The `selector` is a function which can be used to scan only a specific region of the page for changes.
By default, only the `body` is taken into account.
To take the full page, pass `selector=identity`.
"""
struct WebPage <: Page
    url::String
    selector::Function

    function WebPage(url::AbstractString; selector=clean)
        return new(string(url)::String, selector)
    end
end

"""
    MockPage(url::AbstractString, html::String; selector::Function=identity) <: Page

Mock web page used for testing.
"""
struct MockPage <: Page
    url::String
    html::String
    selector::Function

    function MockPage(url::AbstractString, html::String; selector::Function=clean)
        return new(string(url)::String, html, selector)
    end
end

"""
    PageScan(page::Page, content::String)

Store scanned `content` of `page`.
"""
struct PageScan
    page::Page
    content::String

    function PageScan(page::Page, content)
        return new(page, string(content)::String)
    end
end

urls(scans::AbstractVector{PageScan}) = [scan.page for scan in scans]

function scan(page::MockPage)::PageScan
    return PageScan(page, page.html)
end

function scan(page::WebPage)::PageScan
    url = string(page.url)::String
    response = get(url)
    content = String(response.body)
    selection = page.selector(content)
    return PageScan(page, selection)
end


text:
@testset "1178.number-of-valid-words-for-each-puzzle.jl" begin
    @test find_num_of_valid_words(
        ["aaaa", "asas", "able", "ability", "actt", "actor", "access"],
        ["aboveyz", "abrodyz", "abslute", "absoryz", "actresz", "gaswxyz"],
    ) == [1, 1, 3, 2, 4, 0]
    @test find_num_of_valid_words(
        ["apple", "pleas", "please"],
        ["aelwxyz", "aelpxyz", "aelpsxy", "saelpxy", "xaelpsy"],
    ) == [0, 1, 3, 2, 0]
end


text:
using Test

@testset "CspExamples.jl" begin
    include("CspExamples.jl")
end


text:
using DashBootstrapComponents

modal = html_div([
    html_div(
        [
            dbc_label("Backdrop:"),
            dbc_radioitems(
                id = "backdrop-selector",
                options = [
                    Dict("label" => "True (default)", "value" => true),
                    Dict("label" => "False", "value" => false),
                    Dict("label" => "'static'", "value" => "static"),
                ],
                inline = true,
                value = true,
            ),
        ],
        className = "mb-2",
    ),
    dbc_button("Open modal", id = "open-backdrop", n_clicks = 0),
    dbc_modal(
        [
            dbc_modalheader(dbc_modaltitle("Header"), close_button = true),
            dbc_modalbody("Change the backdrop of this modal with the radio buttons"),
            dbc_modalfooter(
                dbc_button(
                    "Close",
                    id = "close-backdrop",
                    className = "ms-auto",
                    n_clicks = 0,
                ),
            ),
        ],
        id = "modal-backdrop",
        is_open = false,
    ),
])

callback!(
    app,
    Output("modal-backdrop", "backdrop"),
    Input("backdrop-selector", "value"),
) do backdrop
    return backdrop
end;

callback!(
    app,
    Output("modal-backdrop", "is_open"),
    Input("open-backdrop", "n_clicks"),
    Input("close-backdrop", "n_clicks"),
    State("modal-backdrop", "is_open"),
) do n1, n2, is_open
    return n1 > 0 || n2 > 0 ? is_open == 0 : is_open
end;

max_issues_repo_name: sthagen/dash-bootstrap-components

max_forks_repo_name: sthagen/dash-bootstrap-components


text:
function countInc(a::Array{Int, 1}, left=1, right=length(a))
	if right == left
		return 0
	elseif right-left == 1
		if a[left] > a[right]
			return 1
		else
			return 0
		end
	end
	mid = (right+left) ÷ 2
	return countInc(a, left, mid) + countInc(a, mid, right)
end


text:
nosubp_options = (; solve_subproblems = false) # needed for SOS2 in MOIPajarito

insts = OrderedDict()

insts["test"] = [
    ((2, 2, SOS2(), 21), nosubp_options),
    ((2, 2, CCBounded(), 21), nosubp_options),
    ((2, 2, LogIBBounded(), 21), nosubp_options),
    ((3, 3, SOS2(), 11), nosubp_options),
    ((3, 3, CCBounded(), 11), nosubp_options),
    ((3, 3, LogIBBounded(), 11), nosubp_options),
    # ((1, 2, CCNonneg(), 11), nosubp_options),
    # ((1, 2, LogIBNonneg(), 11), nosubp_options),
]

function ballpacking_insts(pwl::PWLSOS2, max_p::Int)
    ts = [((4, 3, pwl, 7), nosubp_options)] # compile
    for p in 3:max_p
        push!(ts, ((4, 3, pwl, 2^p - 1), nosubp_options))
    end
    return ts
end

insts["sos2"] = ballpacking_insts(SOS2(), 10)
insts["logib"] = ballpacking_insts(LogIBBounded(), 13)
insts["cc"] = ballpacking_insts(CCBounded(), 7)

return (BallPacking, insts)


text:
module TestSets
using Test

using Catlab.Theories, Catlab.CategoricalAlgebra.Sets

# Sets from Julia types
#######################

# Callables.
f = SetFunction(x -> 2x, TypeSet(Int), TypeSet(Int))
g = SetFunction(x -> 3x, TypeSet(Int), TypeSet(Int))
@test dom(f) == TypeSet(Int)
@test codom(f) == TypeSet(Int)
@test f(1) == 2

# Constants.
c = ConstantFunction("foo", TypeSet(Int))
@test dom(c) == TypeSet(Int)
@test codom(c) == TypeSet(String)
@test c(1) == "foo"

# Identities.
X = TypeSet(Int)
@test dom(id(X)) == X
@test codom(id(X)) == X
@test id(X)(1) == 1

# Composition.
h = compose(f,g)
@test h(1) == 6
@test compose(id(dom(f)), f) == f
@test compose(f, id(codom(f))) == f

c = ConstantFunction(5, TypeSet(Int))
@test compose(c, f) == ConstantFunction(f(5), TypeSet(Int))
@test compose(f, c) == c
@test compose(c, c) == c

# Predicated sets
#################

odds = PredicatedSet(Int, isodd)
evens = PredicatedSet(Int, iseven)
@test sprint(show, odds) == "PredicatedSet($(Int), isodd)"

plus_one_to_odd = SetFunction(x -> x+1, odds, evens)
plus_one_to_even = SetFunction(x -> x+1, evens, odds)

@test plus_one_to_odd(1) == 2
@test_throws ErrorException plus_one_to_odd(2) == 3
@test plus_one_to_even(2) == 3
@test_throws ErrorException plus_one_to_even(3) == 4

plus_two_to_odd = compose(plus_one_to_odd, plus_one_to_even)
@test plus_two_to_odd(3) == 5
@test dom(plus_two_to_odd) == odds
@test codom(plus_two_to_odd) == odds
@test_throws ErrorException plus_two_to_odd(2) == 4

# Limits
########

# Terminal object.
@test ob(terminal(TypeSet)) == TypeSet(Nothing)
@test delete(terminal(TypeSet), TypeSet(Int))(3) == nothing

# Binary product.
lim = product(TypeSet(Int), TypeSet(String))
@test ob(lim) == TypeSet(Tuple{Int,String})
π1, π2 = lim
@test (π1((1,"foo")), π2((1,"foo"))) == (1,"foo")

f = SetFunction(length, TypeSet(Vector{String}), TypeSet(Int))
g = SetFunction(v -> sprint(join, v), TypeSet(Vector{String}), TypeSet(String))
@test pair(lim, f, g)(["foo", "bar", "baz"]) == (3, "foobarbaz")

# N-ary product.
lim = product(fill(TypeSet(Int), 3))
@test ob(lim) == TypeSet(Tuple{Int,Int,Int})
π1, π2, π3 = lim
@test (π1((1,2,3)), π2((1,2,3)), π3((1,2,3))) == (1,2,3)

fs = [ SetFunction(x -> x+i, TypeSet(Int), TypeSet(Int)) for i in 1:3 ]
@test pair(lim, fs)(3) == (4,5,6)

end


text:
using Distributions, ZVSimulator

d = TDist(5.)

results = psim_rand_vrf(d, nsets=5)

# writedlm_vrf(results, "./", ' ')
# writedlm_reformat("_means.txt", "./_means_pprint.txt", fill("%.3e", 3))


text:
using FractionalCalculus
using SymbolicUtils
using SpecialFunctions
using Test

@testset "FractionalCalculus" begin
    include("./Derivative.jl")
    include("./Integral.jl")
    include("./SymbolicDiffTests.jl")
    include("./SymbolicIntTests.jl")
end


text:
using Terrain
using Rayly, FixedSizeArrays, ColorTypes
using Images

function shade(i::Intersection)
    f = abs(dot(normal(i), i.ray.dir))
    RGB(f,f,f)
end

background = RGB(0.,0,0)

dd = -0.3

cam = SimpleCamera(Point{3}(0.,-1,-3),
                   Vec{3}(0.,1,1dd),
                   Vec{3}(0.,-1dd,1),
                   1.0, 1.0, 1.0,
                   400, 400)

hf = Terrain.Heightfield(1.,1.,100,100)
Terrain.faultline!(hf, 500, Terrain.sinstep(0.1,0.1))

compress(a::AbstractArray, mi, ma) = (a - minimum(a))*(ma-mi)/(maximum(a)-minimum(a)) + mi

hf2 = Heightfield(hf)
hf2.y = compress(hf2.y, -0.2, 0.2)
tris = Terrain.triangle_list(hf2)

acc = ListAccelerator()
add_tris!(acc, tris...)

@time img = render(cam) do ray::Ray
    hit = Rayly.intersect(acc, ray)
    hit != nothing ? shade(hit) : background
end

save_clamped("simple.png", img)


text:
# ---
# title: 1377. Frog Position After T Seconds
# id: problem1377
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Depth-first Search
# link: <https://leetcode.com/problems/frog-position-after-t-seconds/description/>
# hidden: true
# ---
#
# Given an undirected tree consisting of `n` vertices numbered from `1` to `n`.
# A frog starts jumping from **vertex 1**. In one second, the frog jumps from
# its current vertex to another **unvisited** vertex if they are directly
# connected. The frog can not jump back to a visited vertex. In case the frog
# can jump to several vertices, it jumps randomly to one of them with the same
# probability. Otherwise, when the frog can not jump to any unvisited vertex, it
# jumps forever on the same vertex.
#
# The edges of the undirected tree are given in the array `edges`, where
# `edges[i] = [ai, bi]` means that exists an edge connecting the vertices `ai`
# and `bi`.
#
# _Return the probability that after`t` seconds the frog is on the vertex
# `target`._
#
#
#
# **Example 1:**
#
# ![](https://assets.leetcode.com/uploads/2020/02/20/frog_2.png)
#
#
#
#     Input: n = 7, edges = [[1,2],[1,3],[1,7],[2,4],[2,6],[3,5]], t = 2, target = 4
#     Output: 0.16666666666666666
#     Explanation: The figure above shows the given graph. The frog starts at vertex 1, jumping with 1/3 probability to the vertex 2 after **second 1** and then jumping with 1/2 probability to vertex 4 after **second 2**. Thus the probability for the frog is on the vertex 4 after 2 seconds is 1/3 * 1/2 = 1/6 = 0.16666666666666666.
#
#
# **Example 2:**
#
# **![](https://assets.leetcode.com/uploads/2020/02/20/frog_3.png)**
#
#
#
#     Input: n = 7, edges = [[1,2],[1,3],[1,7],[2,4],[2,6],[3,5]], t = 1, target = 7
#     Output: 0.3333333333333333
#     Explanation: The figure above shows the given graph. The frog starts at vertex 1, jumping with 1/3 = 0.3333333333333333 probability to the vertex 7 after **second 1**.
#
#
# **Example 3:**
#
#
#
#     Input: n = 7, edges = [[1,2],[1,3],[1,7],[2,4],[2,6],[3,5]], t = 20, target = 6
#     Output: 0.16666666666666666
#
#
#
#
# **Constraints:**
#
#   * `1 <= n <= 100`
#   * `edges.length == n - 1`
#   * `edges[i].length == 2`
#   * `1 <= ai, bi <= n`
#   * `1 <= t <= 50`
#   * `1 <= target <= n`
#   * Answers within `10-5` of the actual value will be accepted as correct.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
using TransformVariables,
      LogDensityProblems,
      DynamicHMC,
      Statistics,
      ForwardDiff
import LogDensityProblems: ADgradient

export dynamicHMC

"""
    dynamicHMC(
        rng::AbstractRNG,
        m::ConditionalModel,
        _data,
        N::Int = 1000;
        method = logdensity,
        ad_backend = Val(:ForwardDiff),
        reporter = DynamicHMC.NoProgressReport(),
        kwargs...)

Draw `N` samples from the posterior distribution of parameters defined in Soss model `m`, conditional on `_data`. Samples are drawn using Hamiltonial Monte Carlo (HMC) from the `DynamicHMC.jl` package.

This function is essentially a wrapper around `DynamicHMC.mcmc_with_warmup()`. Arguments `reporter`, `ad_backend` [DynamicHMC docs here](https://tamaspapp.eu/DynamicHMC.jl/stable/interface/#DynamicHMC.mcmc_with_warmup))

## Arguments
* `rng`: Random number generator.
* `m`: Soss model.
* `_data`: `NamedTuple` of data to condition on.

## Keyword Arguments
*  `N = 1000`: Number of samples to draw.
*  `method = logdensity`: How to compute the log-density. Options are `logdensity` (delegates to `logdensity` of each component) or `codegen` (symbolic simplification and code generation).
*  `ad_backend = Val(:ForwardDiff)`: Automatic differentiation backend.
*  `reporter = DynamicHMC.NoProgressReport()`: Specify logging during sampling. Default: do not log progress.
*  `kwargs`: Additional keyword arguments passed to core sampling function `DynamicHMC.mcmc_with_warmup()`.

Returns an Array of `Namedtuple` of length `N`. Each entry in the array is a sample of parameters indexed by the parameter symbol.

## Example

```jldoctest
using StableRNGs
rng = StableRNG(42);

m = @model x begin
    β ~ Normal()
    yhat = β .* x
    y ~ For(eachindex(x)) do j
        Normal(yhat[j], 2.0)
    end
end

x = randn(rng, 3);
truth = [-0.41, 1.21, 0.11];

post = dynamicHMC(rng, m(x=x), (y=truth,));
E_β = mean(getfield.(post, :β))

println("Posterior mean β: " * string(round(E_β, digits=2)))

# output
Posterior mean β: 0.25
```

"""
function dynamicHMC(
    rng::AbstractRNG,
    m::ConditionalModel,
    N::Int = 1000;
    # method = logpdf,
    ad_backend = Val(:ForwardDiff),
    reporter = DynamicHMC.NoProgressReport(),
    kwargs...,
)

    M = getmoduletypencoding(m)

    ℓ = if haskey(kwargs, :ℓ)
        codegen(m; ℓ = kwargs[:ℓ])
    else
        (a, o, pars) -> _logdensity(M, Model(m), a, o, pars)
    end

    _argvals = argvals(m)
    _obs = observations(m)

    logp(pars) = ℓ(_argvals, _obs, pars)

    t = xform(m)
    P = LogDensityProblems.TransformedLogDensity(t, logp)
    ∇P = LogDensityProblems.ADgradient(ad_backend, P)

    results = DynamicHMC.mcmc_with_warmup(
        rng,
        ∇P,
        N;
        reporter = reporter
    )
    T = typeof(t(zeros(t.dimension)))

    x = TupleArray{T,1}(undef, N)

    for j in 1:N
        @inbounds x[j] = TransformVariables.transform(t, results.chain[j])
    end

    return x
    # samples = TransformVariables.transform.(t, results.chain)
    # return samples
end

function dynamicHMC(
    rng::AbstractRNG,
    m::ConditionalModel,
    ::Val{Inf};
    method = logdensity,
    ad_backend = Val(:ForwardDiff),
    reporter = DynamicHMC.NoProgressReport(),
    kwargs...,
)
    _data = m.obs
    ℓ(pars) = logdensity(m, merge(pars, _data), method)
    t = xform(m, _data)
    P = LogDensityProblems.TransformedLogDensity(t, ℓ)
    ∇P = LogDensityProblems.ADgradient(ad_backend, P)

    results = DynamicHMC.mcmc_keep_warmup(
        rng,
        ∇P,
        0;
        reporter = reporter,
        kwargs...,
    )
    steps = DynamicHMC.mcmc_steps(
        results.sampling_logdensity,
        results.final_warmup_state,
    )
    return results, steps
end

function dynamicHMC(m::ConditionalModel, args...; kwargs...)
    return dynamicHMC(Random.GLOBAL_RNG, m, args...; kwargs...)
end

# using ResumableFunctions

# export stream

# @resumable function stream(
#     rng::AbstractRNG,
#     f::typeof(dynamicHMC),
#     m::ConditionalModel,
#     _data::NamedTuple,
# )
#     t = xform(m, _data)
#     (results, steps) = dynamicHMC(rng, m, _data, Val(Inf))
#     Q = results.final_warmup_state.Q
#     while true
#         Q, tree_stats = DynamicHMC.mcmc_next_step(steps, Q)
#         @yield (merge(t(Q.q), (_ℓ = Q.ℓq,)), tree_stats)
#     end
# end

# function stream(
#     f::typeof(dynamicHMC),
#     m::ConditionalModel,
#     _data::NamedTuple;
#     kwargs...,
# )
#     return stream(Random.GLOBAL_RNG, f, m, _data; kwargs...)
# end


text:
using HTTP, availability, Sockets, JSON, Random, JuliaDB

const AVAILABILITY_ROUTER = HTTP.Router()

supplyDataset = load("temp.dat")

function toSupplyArray(res)
        list = [];

        for i in res
                push!(list, availability.Supply(i.itemId, i.uom, i.nodeId, i.type, i.quantity))
        end

        return list
end

function getSupply(req::HTTP.Request)
        # @show req

        itemId = HTTP.URIs.splitpath(req.target)[2]
        # @show "itemId: $(itemId)"

        result = filter(r -> isequal(r.itemId, itemId), supplyDataset)
        # @show result

        out = Dict(
                "totalRecords" => length(result),
                "totalQuantity" => ((t -> select(t, :quantity)) ∘ (t -> filter(r -> isequal(r.itemId, itemId), t)))(result) |> sum,
                "supplies" => toSupplyArray(result))

        if(length(result) < 1)
                status = 404
        else
                status = 200
        end

        response = HTTP.Response(status, JSON.json(out))
        HTTP.setheader(response, "Content-Type" => "application/json")
        HTTP.setheader(response, "myHeader" => "headerValue1")

        return response

end
function createSupply(req::HTTP.Request)
        @show req
        println("Supply created")
        HTTP.Response(200)
end

HTTP.@register(AVAILABILITY_ROUTER, "GET", "/supplies/*", getSupply)
HTTP.@register(AVAILABILITY_ROUTER, "POST", "/supplies", createSupply)

# @async HTTP.serve(AVAILABILITY_ROUTER, Sockets.localhost, 8081)
HTTP.serve(AVAILABILITY_ROUTER, Sockets.localhost, 8081)

max_stars_repo_name: kawadesumit/availability

max_issues_repo_name: kawadesumit/availability

max_forks_repo_name: kawadesumit/availability


text:
"""
    Octic{T}

The scalar conservation law with flux
``
  \\partial_t u + \\partial_x u^8 = 0
``
in one space dimensions using `T` as scalar type.
"""
struct Octic{T} <: ScalarBalanceLaw{T,1} end

function Octic(T=Float64)
  Octic{T}()
end

function Base.show(io::IO, model::Octic{T}) where {T}
  print(io, "Scalar conservation law {T=", T, "} with flux f(u) = u^8")
end

"""
    flux(u, model::Octic)

Compute the flux of `u` for `model`.
"""
@inline flux(u, model::Octic) = ((u^2)^2)^2

"""
    speed(u::Real, model::Octic)

Compute the speed f'(`u`) for `model`.
"""
@inline speed(u, model::Octic) = 8*(u^2)^2*u^3

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

"""
    (::GodunovFlux)(uₗ::T, uᵣ::T, model::Octic{T}) where {T}

Compute Godunov's flux between `uₗ` and `uᵣ` for `model`.
"""
function (::GodunovFlux)(uₗ::T, uᵣ::T, model::Octic{T}) where {T}
    if uₗ < uᵣ
      if uₗ < 0 && 0 < uᵣ
        zero(T)
      else
        min(flux(uₗ, model), flux(uᵣ, model))
      end
    else
      max(flux(uₗ, model), flux(uᵣ, model))
    end
end

"""
    (::EnergyConservativeFlux)(uₗ::T, uᵣ::T, model::Octic{T}) where {T}

Compute the energy (L₂ entropy) conservative flux between `uₗ` and `uᵣ` for `model`.
"""
function (flux::EnergyConservativeFlux)(uₗ::T, uᵣ::T, model::Octic{T}) where T<:Real
    ( ((uₗ^2)^2)^2 + (uₗ^2)^2*uₗ^3*uᵣ + (uₗ^3)^2*uᵣ^2 + uₗ^2*uₗ^3*uᵣ^3 + (uₗ^2)^2*(uᵣ^2)^2
    + uₗ^3*uᵣ^2*uᵣ^3 + uₗ^2*(uᵣ^3)^2 + uₗ*(uᵣ^2)^2*uᵣ^3 + ((uᵣ^2)^2)^2 ) / 9
end

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

@inline function add_numerical_fluxes_inner_loop2!(du, fluxes, u, balance_law::Octic,
                                                    Nx, basis, boundaries_included::Val{false},
                                                    jacx, parallel)
    Pp1 = length(basis.nodes)
    Rl = basis.interp_left
    Rr = basis.interp_right
    one_9 = 1 / 9
    two_9 = 2 / 9

    # add numerical fluxes
    @inbounds for ix in Base.OneTo(Nx)
        Rul = zero(eltype(u))
        Rur = zero(eltype(u))
        for nx in Base.OneTo(Pp1)
            Rul += Rl[nx]*u[nx,ix]
            Rur += Rr[nx]*u[nx,ix]
        end
        Ru2l = zero(eltype(u))
        Ru2r = zero(eltype(u))
        for nx in Base.OneTo(Pp1)
            tmp = u[nx,ix]^2
            Ru2l += Rl[nx]*tmp
            Ru2r += Rr[nx]*tmp
        end
        Ru3l = zero(eltype(u))
        Ru3r = zero(eltype(u))
        for nx in Base.OneTo(Pp1)
            tmp = u[nx,ix]^3
            Ru3l += Rl[nx]*tmp
            Ru3r += Rr[nx]*tmp
        end
        Ru4l = zero(eltype(u))
        Ru4r = zero(eltype(u))
        for nx in Base.OneTo(Pp1)
            tmp = (u[nx,ix]^2)^2
            Ru4l += Rl[nx]*tmp
            Ru4r += Rr[nx]*tmp
        end
        Ru5l = zero(eltype(u))
        Ru5r = zero(eltype(u))
        for nx in Base.OneTo(Pp1)
            tmp = u[nx,ix]^2*u[nx,ix]^3
            Ru5l += Rl[nx]*tmp
            Ru5r += Rr[nx]*tmp
        end
        Ru6l = zero(eltype(u))
        Ru6r = zero(eltype(u))
        for nx in Base.OneTo(Pp1)
            tmp = (u[nx,ix]^3)^2
            Ru6l += Rl[nx]*tmp
            Ru6r += Rr[nx]*tmp
        end
        Ru7l = zero(eltype(u))
        Ru7r = zero(eltype(u))
        for nx in Base.OneTo(Pp1)
            tmp = (u[nx,ix]^3)^2*u[nx,ix]
            Ru7l += Rl[nx]*tmp
            Ru7r += Rr[nx]*tmp
        end
        Ru8l = zero(eltype(u))
        Ru8r = zero(eltype(u))
        for nx in Base.OneTo(Pp1)
            tmp = ((u[nx,ix]^2)^2)^2
            Ru8l += Rl[nx]*tmp
            Ru8r += Rr[nx]*tmp
        end

        for nx in Base.OneTo(Pp1)
            du[nx,ix] += (  (fluxes[ix  ]   - two_9 * Ru8l
                                            - two_9 * u[nx,ix] * Ru7l
                                            - two_9 * u[nx,ix]^2 * Ru6l
                                            - two_9 * u[nx,ix]^3 * Ru5l
                                            - one_9 * ((Rul^2)^2)^2
                                            + one_9 * u[nx,ix] * Rul^3*(Rul^2)^2
                                            - one_9 * Ru2l * (Rul^3)^2
                                            + one_9 * u[nx,ix]^2 * (Rul^3)^2
                                            - one_9 * Ru3l * Rul^2*Rul^3
                                            + one_9 * u[nx,ix]^3 * Rul^2*Rul^3
                                            - one_9 * Ru4l * (Rul^2)^2
                                            + one_9 * u[nx,ix] * Ru4l * Rul^3
                                            - one_9 * Ru4l * Ru2l * Rul^2
                                            + one_9 * u[nx,ix]^2 * Ru4l * Rul^2
                                            - one_9 * Ru4l * Ru3l * Rul
                                            + one_9 * u[nx,ix]^3 * Ru4l * Rul
                                            - one_9 * Ru4l^2) * Rl[nx]
                          - (fluxes[ix+1]   - two_9 * Ru8r
                                            - two_9 * u[nx,ix] * Ru7r
                                            - two_9 * u[nx,ix]^2 * Ru6r
                                            - two_9 * u[nx,ix]^3 * Ru5r
                                            - one_9 * ((Rur^2)^2)^2
                                            + one_9 * u[nx,ix] * Rur^3*(Rur^2)^2
                                            - one_9 * Ru2r * (Rur^3)^2
                                            + one_9 * u[nx,ix]^2 * (Rur^3)^2
                                            - one_9 * Ru3r * Rur^2*Rur^3
                                            + one_9 * u[nx,ix]^3 * Rur^2*Rur^3
                                            - one_9 * Ru4r * (Rur^2)^2
                                            + one_9 * u[nx,ix] * Ru4r * Rur^3
                                            - one_9 * Ru4r * Ru2r * Rur^2
                                            + one_9 * u[nx,ix]^2 * Ru4r * Rur^2
                                            - one_9 * Ru4r * Ru3r * Rur
                                            + one_9 * u[nx,ix]^3 * Ru4r * Rur
                                            - one_9 * Ru4r^2) * Rr[nx]
                            ) * jacx / basis.weights[nx]
        end
    end

    nothing
end


text:
using Test, Pkg
using Base: UUID
using Pkg.Types: VersionRange
using Pkg.Operations: load_package_data_raw
import Pkg.TOML

@testset "(Registry|Package|Versions|Deps|Compat).toml" begin
    reg = TOML.parsefile("Registry.toml")
    reguuids = Set{UUID}(UUID(x) for x in keys(reg["packages"]))
    stdlibuuids = Set{UUID}(x for x in keys(Pkg.Types.gather_stdlib_uuids()))
    alluuids = reguuids ∪ stdlibuuids

    # Test that each entry in Registry.toml has a corresponding Package.toml
    # at the expected path with the correct uuid and name
    for (uuid, data) in reg["packages"]
        # Package.toml testing
        pkg = TOML.parsefile(abspath(data["path"], "Package.toml"))
        @test UUID(uuid) == UUID(pkg["uuid"])
        @test data["name"] == pkg["name"]
        @test haskey(pkg, "repo")

        # Versions.toml testing
        vers = TOML.parsefile(abspath(data["path"], "Versions.toml"))
        for (v, data) in vers
            @test VersionNumber(v) isa VersionNumber
            @test haskey(data, "git-tree-sha1")
        end

        # Deps.toml testing
        depsfile = abspath(data["path"], "Deps.toml")
        if isfile(depsfile)
            deps = TOML.parsefile(depsfile)
            # Require all deps to exist in the General registry or be a stdlib
            depuuids = Set{UUID}(UUID(x) for (_, d) in deps for (_, x) in d)
            @test depuuids ⊆ alluuids
            # Test that the way Pkg loads this data works
            @test Pkg.Operations.load_package_data_raw(UUID, depsfile) isa
                Dict{Pkg.Types.VersionRange,Dict{String,UUID}}
        end

        # Compat.toml testing
        compatfile = abspath(data["path"], "Compat.toml")
        if isfile(compatfile)
            compat = TOML.parsefile(compatfile)
            # Test that all names with compat is a dependency
            compatnames = Set{String}(x for (_, d) in compat for (x, _) in d)
            if !(isempty(compatnames) || (length(compatnames)== 1 && "julia" in compatnames))
                depnames = Set{String}(x for (_, d) in TOML.parsefile(depsfile) for (x, _) in d)
                push!(depnames, "julia") # All packages has an implicit dependency on julia
                @assert compatnames ⊆ depnames
            end
            # Test that the way Pkg loads this data works
            @test Pkg.Operations.load_package_data_raw(Pkg.Types.VersionSpec, compatfile) isa
                Dict{Pkg.Types.VersionRange,Dict{String,Pkg.Types.VersionSpec}}
        end
    end
end

max_stars_repo_name: lruthotto/General

max_issues_repo_name: lruthotto/General

max_forks_repo_name: lruthotto/General


text:
module GeoMakie

using Reexport
using GeometryTypes, EarCut, Colors, MakieLayout, AbstractPlotting, ImageMagick, GeoJSON, Glob

using MakieLayout: Formatting

@reexport using Colors, AbstractPlotting
@reexport using GeoInterface, Proj4, NASAEarthObservations

import AbstractPlotting: convert_arguments, convert_attribute
using AbstractPlotting.DocStringExtensions

# fix conflicts
const AbstractGeometry = GeometryTypes.AbstractGeometry
const Point = AbstractPlotting.Point
const attributes = AbstractPlotting.attributes
const transform! = Proj4.transform!
const volume = AbstractPlotting.volume

include("conversions.jl")
include("projections.jl")
include("utils.jl")
include("tick_compat.jl")
include("data.jl")

include("recipes/geoaxis.jl")
include("recipes/stock.jl")

# export recipes

export geoaxis, geoaxis!

# export projections

export LonLat, PlateCarree, Mercator, TransverseMercator, UTM, AlbersEqualArea, Robinson, WinkelTripel

# export utilities

export rec_project, rec_point, xygrid

function __init__()
    Base.delete_method(methods(coastlines!).ms[6])
    Base.delete_method(methods(earth!).ms[6])
end

end # module

## TODOS
# TODO note you need to return a tuple rom conv-arg