text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "Antic"
version = v"0.2.2"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://github.com/wbhart/antic/archive/0.2.2.tar.gz",
                  "9e8d518fc9b27807d5bfd173cc81d779a6a1e07ae37bc59427c5223283e9e4d4")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
cd antic*/
if [[ ${target} == *musl* ]]; then
   export CFLAGS=-D_GNU_SOURCE=1;
elif [[ ${target} == *mingw* ]]; then
   sed -i -e "s#/lib\>#/$(basename ${libdir})#g" configure
   extraflags=--build=MINGW${nbits};
fi
./configure --prefix=$prefix --disable-static --enable-shared --with-gmp=$prefix --with-mpfr=$prefix --with-flint=$prefix ${extraflags}
make -j${nproc}
make install LIBDIR=$(basename ${libdir})
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()

# The products that we will ensure are always built
products = [
    LibraryProduct("libantic", :libantic)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="FLINT_jll", uuid="e134572f-a0d5-539d-bddf-3cad8db41a82"))
    Dependency("GMP_jll", v"6.1.2"),
    Dependency("MPFR_jll", v"4.0.2"),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: mortenpi/Yggdrasil

max_issues_repo_name: mortenpi/Yggdrasil

max_forks_repo_name: mortenpi/Yggdrasil


text:
using JSON
using Optim

include("equations.jl")

"""
The parameter, which needs to be predicted.
"""
@enum Predict TEMPERATURE PRESSURE

"""
    EquilibriumPoint(eq, predict, phase; max_P=nothing, max_T=nothing,
                     tol=nothing, incr=nothing, incr_min=nothing,
                     incr_den=nothing)

Contains the data that necessary for predicting an equilibrium point using an
equation of state.
"""
type EquilibriumPoint

    eq::Equation
    predict::Predict
    phase::Phase
    max_P::Float64
    max_T::Float64
    tol::Float64
    incr::Float64
    incr_min::Float64
    incr_den::Float64

    function EquilibriumPoint(eq, predict, phase; max_P=nothing, max_T=nothing,
                              tol=nothing, incr=nothing, incr_min=nothing,
                              incr_den=nothing)

        max_P == nothing && (max_P = max(eq.mix.Pcr...) * 1.1)
        max_T == nothing && (max_T = max(eq.mix.Tcr...) * 1.3)
        tol == nothing && (tol = 0.001)
        if predict == PRESSURE
            incr == nothing && (incr = 101325.0)
            incr_min == nothing && (incr_min = 100.0)
            incr_den == nothing && (incr_den = 100.0)
        else
            incr == nothing && (incr = 10.0)
            incr_min == nothing && (incr_min = 0.05)
            incr_den == nothing && (incr_den = 2.0)
        end

        return new(eq, predict, phase, max_P, max_T, tol, incr,
                   incr_min, incr_den)
    end
end

"""
    CriticalPoint(ep; T=nothing, P=nothing, incr=nothing, incr_min=nothing,
                  incr_den=nothing, corr_coef=nothing)

Contains the data that necessary for predicting an critical point using an
equation of state.
"""
type CriticalPoint

    ep::EquilibriumPoint
    T::Float64
    P::Float64
    incr::Float64
    incr_min::Float64
    incr_den::Float64
    corr_coef::Float64

    function CriticalPoint(ep; T=nothing, P=nothing, incr=nothing,
                           incr_min=nothing, incr_den=nothing,
                           corr_coef=nothing)
        T == nothing && (T = sum(ep.eq.mix.Tb .* ep.eq.mix.x))
        P == nothing && (P = 101325.0)
        if ep.predict == PRESSURE
            incr == nothing && (incr = 10.0)
            incr_min == nothing && (incr_min = 0.05)
            incr_den == nothing && (incr_den = 14.14)
            corr_coef == nothing && (corr_coef = 0.9)
        else
            incr == nothing && (incr = 202650.0)
            incr_min == nothing && (incr_min = 3000.0)
            incr_den == nothing && (incr_den = 8.21)
            corr_coef == nothing && (corr_coef = 0.99)
        end

        return new(ep, T, P, incr, incr_min, incr_den, corr_coef)

    end
end

################################################################################
## Equilibrium Point
################################################################################

"""
    solve(ep::EquilibriumPoint, T::Float64, P::Float64)

Solves an equation of state and returns an equilibrium points for the given
properties.
"""
function solve(eps_::Array{EquilibriumPoint, 1}, T::Float64, P::Float64)
	function app(a::Array{EquilibriumData, 1}, b::Array{EquilibriumData, 1})
        append!(a, b)
        return a
    end
    eds = @parallel app for ep ∈ eps_
        [solve(ep, T, P)]
    end
    return eds
end

"""
    solve(ep::EquilibriumPoint, T::Float64, P::Float64)

Solves an equation of state and returns an equilibrium point for the given
properties.
"""
function solve(ep::EquilibriumPoint, T::Float64, P::Float64)
    incr = ep.incr
    ed = predict!(ep, T, P)
    ep.incr = incr
    return ed
end

"""
    predict!(ep::EquilibriumPoint, T::Float64, P::Float64)

Solves an equation of state and returns an equilibrium point for the given
properties. But it chenges the increment.
"""
function predict!(ep::EquilibriumPoint, T::Float64, P::Float64)

    s = 0
    if ep.predict == PRESSURE
        crit = ep.max_P
        slv = (x) -> solve(ep.eq, T, x)
        adj_s = (x) -> x
        val = P
    else
        crit = ep.max_T
        slv = (x) -> solve(ep.eq, x, P)
        adj_s = (x) -> 1.1 - x / 10.0
        val = T
    end
    Rval = val
    n_iter = 0
    while true
        ed = slv(val)
        ep.incr < ep.incr_min && return nothing
        if ed == nothing
            if val > crit
                ep.incr /= ep.incr_den
                val = Rval
            end
            val += ep.incr
        else
            if ep.phase == VAPOR
                s = adj_s(update_x!(ep, ed))
                val /= s
				n_iter += 1
                if n_iter > 1000 ep.tol *= 2 end
            else
                s = adj_s(update_y!(ep, ed))
                val *= s
				n_iter += 1
                if n_iter > 1000 ep.tol *= 2 end
            end
        end
		abs(1 - s) < ep.tol && return ed
    end
end

"""
    update_x!(ep, ed)

Updates liquid fraction in the given mixture and returns the S coefficient for
the further computations.
"""
function update_x!(ep, ed)
    k = ed.ϕv ./ ed.ϕl
    ky = k .* ed.y
    s = sum(ky)
    ep.eq.mix.x = ky / s
    return s
end

"""
    update_y!(ep, ed)

Updates vapor fraction in the given mixture and returns the S coefficient for
the further computations.
"""
function update_y!(ep, ed)
    k = ed.ϕl ./ ed.ϕv
    kx = k .* ed.x
    s = sum(kx)
    ep.eq.mix.y = kx / s
    return s
end

################################################################################
## Critical Point
################################################################################

"""
    solve(cp::CriticalPoint)

Solves an equation of state and returns an array of equilibrium points, where
last of them a critical, for the given properties.
"""
function solve(cp::CriticalPoint)
    return predict!(deepcopy(cp))
end

"""
    predict!(cp::CriticalPoint)

Solves an equation of state and returns an an array of equilibrium points, where
last of them a critical, for the given properties. But is changes the original
properties.
"""
function predict!(cp::CriticalPoint)

    eds::Array{EquilibriumData, 1}, ed = [], 1
    if cp.ep.predict == PRESSURE
        upd1! = (cp::CriticalPoint) -> cp.T -= cp.incr
        upd2! = (cp::CriticalPoint) -> cp.T += cp.incr
        upd3! = (cp::CriticalPoint) -> cp.P *= cp.corr_coef
    else
        upd1! = (cp::CriticalPoint) -> cp.P -= cp.incr
        upd2! = (cp::CriticalPoint) -> cp.P += cp.incr
        upd3! = (cp::CriticalPoint) -> cp.T *= cp.corr_coef
    end

    while ed != nothing
        ed = solve(cp.ep, cp.T, cp.P)
        if ed != nothing
            cp.T, cp.P = ed.T, ed.P
            push!(eds, ed)
            upd3!(cp)
        elseif cp.incr / cp.incr_den > cp.incr_min
            ed = 1
            upd1!(cp)
            cp.incr /= cp.incr_den
        end
        upd2!(cp)
    end

    return eds
end

"""
    solve(cps::Array{CriticalPoint, 1})

Solves multiple critical points in parallel (if multiple threads is defined).
"""
function solve(cps::Array{CriticalPoint, 1})
    function app(a::Array{Array{EquilibriumData, 1}, 1},
                b::Array{Array{EquilibriumData, 1}, 1})
        append!(a, b)
        return a
    end
    edss = @parallel app for cp ∈ cps
        [solve(cp)]
    end
    return edss
end

"""
    fit(cps::Array{CriticalPoint, 1},
        k_bounds::Tuple{Float64,Float64}=(-1.0,1.0); Tcr_exp=nothing,
        Pcr_exp=nothing, Vcr_exp=nothing, w=[1.0,1.0,1.0], opt_set...)

Optimizes binary interaction parameter for given experimental critical data.
You may provide additional arguments to the optimization function (from Optim
using package) using keyword arguments.
"""
function fit(cps::Array{CriticalPoint, 1},
             k_bounds::Tuple{Float64,Float64}=(-0.1,0.1); Tcr_exp=nothing,
             Pcr_exp=nothing, Vcr_exp=nothing, w=[1.0,1.0,1.0], opt_set...)

    prms = [p for p ∈ [Tcr_exp, Pcr_exp, Vcr_exp] if p != nothing]
    if length(prms) == 0
        error("You have to specify at least one of the critical parametes.")
    end
    function to_opt(k::Float64)
        for cp ∈ cps
            cp.ep.eq.mix.k = [0.0 k; k 0.0]
        end
        edss = solve(cps)
        Tcr_calc::Array{Float64, 1} = [x[end].T for x ∈ edss]
        Pcr_calc::Array{Float64, 1} = [x[end].P for x ∈ edss]
        Vcr_calc::Array{Float64, 1} = [x[end].Vl for x ∈ edss]
        s::Float64 = 0.0
        for (i, (e, c)) ∈ enumerate(zip([Tcr_exp, Pcr_exp, Vcr_exp],
                                      [Tcr_calc, Pcr_calc, Vcr_calc]))
            if e != nothing
                s += sum(w[i] * ((e - c) ./ e).^2)
            end
        end
        return s
    end
    return optimize(to_opt, k_bounds[1], k_bounds[2]; show_trace=true,
                    opt_set...)
end

################################################################################
## PostProcessing Functions
################################################################################

"""
    dump_json(eds::Array{EquilibriumData, 1}, filename::AbstractString)

Saves the given array of equilibrium points in json format.
"""
function dump_json(eds::Array{EquilibriumData, 1}, filename::AbstractString)

    typedict(x) = Dict(fn=>getfield(x, fn) for fn ∈ fieldnames(x))
    to_json = [typedict(x) for x ∈ eds]
    open(filename, "w") do file
        txt = JSON.json(to_json)
        write(file, txt)
    end
end

"""
    dump_json(edss::Array{Array{EquilibriumData, 1}, 1},
                   filename::AbstractString)

Saves the given array of arrays of equilibrium points in json format.
The main use is for saving the product of multiple critical point valuations.
"""
function dump_json(edss::Array{Array{EquilibriumData, 1}, 1},
                   filename::AbstractString)

    typedict(x) = Dict(fn=>getfield(x, fn) for fn ∈ fieldnames(x))
    to_json = [[typedict(x) for x ∈ y] for y ∈ edss]
    open(filename, "w") do file
        txt = JSON.json(to_json)
        write(file, txt)
    end
end

ext: jl


text:
function cubemesh(m,n,o,elemtype)
# Generate mesh for unit cube

p = [[x; y; z] for x = 0.0:1.0/Float64(m):1.0, y = 0.0:1.0/Float64(n):1.0, z = 0.0:1.0/Float64(o):1.0][:]
p = hcat(p...);

m = m+1;
n = n+1;
o = o+1;

if elemtype==0
    t0=[[5 6 7 3]; [5 6 2 3]; [5 1 2 3]; [6 8 7 4]; [6 7 4 3]; [6 2 4 3]];

    map=[1 2];
    map=[map  map.+m];
    map=[map map.+m*n];
    t=map[t0];
    t=kron(t,ones(m-1,1))+kron(ones(size(t)),collect(0:m-2));
    t=kron(t,ones(n-1,1))+kron(ones(size(t)),collect(0:n-2).*m);
    t=kron(t,ones(o-1,1))+kron(ones(size(t)),collect(0:o-2).*(m*n));
else
    t = [1 2 m+2 m+1 m*n+1 m*n+2 m*n+m+2 m*n+m+1];
    t=kron(t,ones(o-1,1))+kron(ones(size(t)),collect(0:o-2).*(m*n));
    t=kron(t,ones(n-1,1))+kron(ones(size(t)),collect(0:n-2).*m);
    t=kron(t,ones(m-1,1))+kron(ones(size(t)),collect(0:m-2));
end

t = Int.(t');

return p, t

end

max_stars_repo_name: glwagner/Exasim

max_issues_repo_name: glwagner/Exasim

max_forks_repo_name: glwagner/Exasim


text:
"""

    full_computation()

A test function which sets up an entire scene, with an object and a light source, as well
as an `ImageFilm` camera, renders the image and writes it to disk.

"""
function full_computation()

    # Create material with flat white spectrum
    whitespec = SampledSpectrum(300, 850, [1.0 for i = 1:100])
    redspec = SampledSpectrum(300, 850, collect(range(0.0, 1.0, length=100)))

    CIE_table = make_CIE_table(whitespec)

    mat = Lambert(whitespec)

    # Create sphere shape
    Tsph = translation(0.0, 0.0, 0.0)
    sph = Sphere(1, 1.0, Tsph)

    # Create open cylinder
    Tcyl = translation(2.0, 1.0, 0.0) * scaling(1.0, 1.0, 2.0)
    cyl = Cylinder(Tcyl)

    # Create primitive combining shape and material
    sph_primitive = GeometricPrimitive(sph, mat, nothing, 1)
    cyl_primitive = GeometricPrimitive(cyl, mat, nothing, 2)

    # Create list of primitives
    primitives = GeometricPrimitive[sph_primitive, cyl_primitive]

    # Generate bounding box hierarchy from primitives
    stuff = BVHAccelerator(primitives)

    # Create a distant light source with a single-line 532 nm spectrum
    p_light = DistantLight(redspec, rotation(Z_AXIS, π/8) * rotation(X_AXIS, π/2))

    # Create scene with bounding box hierarchy and light source
    scene = Scene(stuff, LightSource[p_light])

    # Make a camera
    resX = 1024
    resY = 1024
    widthX = 5
    widthY = 5
    F = ImageFilm(resX, resY, 1.0, TriangleFilter(1, 1), CIE_table)
    camera_position = rotation(Y_AXIS, π/2) * translation(0.0, 0.0, -3.0)
    image_cam = OrthographicCamera(camera_position, widthX, widthY, 0.0, 0.0, F)

    # Set up pixel sampler with 9 rays per pixel (3 rays in both X and Y directions)
    sampler = StratifiedSampler(resX, resY, 3)

    # Declare an integrator with maximum ray depth of 1
    whitted = WhittedIntegrator(1)

    # Run the renderer
    @time render(scene, image_cam, whitted, sampler ; debug=true)

    # Write the output
    write_image(F, "test.png")

    return true
end

@testset "Full computation" begin
    @test full_computation()
    @test isfile("test.png")
end

max_stars_repo_name: dronir/GRAYTR

max_issues_repo_name: dronir/GRAYTR

max_forks_repo_name: dronir/GRAYTR


text:
@dim TestDim "Test dimension"

@testset "dims creation macro" begin
    @test name(TestDim) == "Test dimension"
    @test shortname(TestDim) == "TestDim"
    @test val(TestDim(:test)) == :test
    @test metadata(TestDim(1, AllignedGrid(), "metadata")) == "metadata"
    @test units(TestDim) == nothing
    @test label(TestDim) == "Test dimension"
    @test eltype(TestDim(1)) == Int
    @test eltype(TestDim([1,2,3])) == Vector{Int}
    @test length(TestDim(1)) == 1
    @test length(TestDim([1,2,3])) == 3
end

# Basic dim and array initialisation
a = ones(5, 4)
# Must construct with a tuple for dims/refdims
@test_throws MethodError DimensionalArray(a, X((140, 148)))
@test_throws MethodError DimensionalArray(a, (X((140, 148)), Y((2, 11))), Z(1))
da = DimensionalArray(a, (X((140, 148)), Y((2, 11))))

dimz = dims(da)
@test slicedims(dimz, (2:4, 3)) == ((X(LinRange(142,146,3)),), (Y(8.0),))
@test name(dimz) == ("X", "Y")
@test shortname(dimz) == ("X", "Y")
@test units(dimz) == (nothing, nothing)
@test label(dimz) == ("X, Y")

a = [1 2 3 4
     2 3 4 5
     3 4 5 6]
dimz = (X((143, 145), AllignedGrid(), nothing),
        Y((-38, -35), AllignedGrid(), nothing))
da = DimensionalArray(a, dimz)

@testset "dims" begin
    @test dims(dimz) === dimz
    @test dims(dimz, X) === dimz[1]
    @test dims(dimz, Y) === dimz[2]
    @test_throws ArgumentError dims(dimz, Time)
    @test typeof(dims(da)) == Tuple{X{LinRange{Float64},AllignedGrid{Ordered{Forward,Forward},Center,UnknownSampling,Nothing},Nothing},
                                Y{LinRange{Float64},AllignedGrid{Ordered{Forward,Forward},Center,UnknownSampling,Nothing},Nothing}}
end

@testset "arbitrary dim names" begin
    dimz = (Dim{:row}((10, 30)), Dim{:column}((-20, 10)))
    @test name(dimz) == ("Dim row", "Dim column")
    @test shortname(dimz) == ("row", "column")
    @test label(dimz) == ("Dim row, Dim column")
end


text:
ħ = 1.054571800e-34
Planck = 6.626070040e-34
Boltzmann = 1.38064852e-23

_h = 6.62607004
_k = 1.38064852

"""
     temperature_2_β(T)

Convert physical temperature `T` in mK to thermodynamic `β` in the unit of inverse angular frequency, that is to say ``β = ħ/kT``.
"""
function  temperature_2_β(T)
    5*_h/_k/pi/T
end

"""
    temperature_2_freq(T)

Convert temperature from mK to GHz.
"""
function temperature_2_freq(T)
    _k/_h/10*T
end

"""
    β_2_temperature(β)

Convert thermodynamic `β` in the unit of inverse angular frequency to physical temperature `T` in GHz.
"""
function β_2_temperature(β)
    5*_h/_k/pi/β
end

"""
    freq_2_temperature(freq)

Convert frequency in GHz to temperature in mK.
"""
function freq_2_temperature(freq)
    10 * freq * _h / _k
end


text:
"""
    AbstractSusceptibility

- A susceptibility is used to build dispersion models.
- A dispersion model is a collection of susceptibilities of the same type or different types.
- The supported susceptibilities  `LorentzianSusceptibility` `DrudeSusceptibility`
  `PoleResidue`
- Each susceptibility returns a symmetric 3*3 matrix.
- The sigma field is used when the susceptibility is isotropic
- The sigma_diag and sigma_offdiag vectors are used when the susceptibility is anisotropic and can
be used to turn on/off the anisotropic terms.
"""
abstract type AbstractSusceptibility end

@with_kw struct LorentzianSusceptibility <: AbstractSusceptibility
    frequency
    gamma
    sigma = 1.0
    sigma_diag = sigma * [1, 1, 1]
    sigma_offdiag = [0, 0, 0]
end

@with_kw struct DrudeSusceptibility <: AbstractSusceptibility
    frequency
    gamma
    sigma = 1.0
    sigma_diag = sigma * [1, 1, 1]
    sigma_offdiag = [0, 0, 0]
end

@with_kw struct PoleResidueSusceptibility <: AbstractSusceptibility
    a
    c
    sigma = 1.0
    sigma_diag = sigma * [1, 1, 1]
    sigma_offdiag = [0, 0, 0]
end

function eval_susceptibility(s::LorentzianSusceptibility, f)
    @unpack frequency, gamma, sigma, sigma_diag, sigma_offdiag = s
    if iszero(gamma)
        κ = frequency^2 / (frequency^2 - f^2)
    else
        κ = frequency^2 / (frequency^2 - f^2 - im * gamma * f)
    end
    [
        sigma_diag[1] sigma_offdiag[1] sigma_offdiag[2]
        sigma_offdiag[1] sigma_diag[2] sigma_offdiag[3]
        sigma_offdiag[2] sigma_offdiag[3] sigma_diag[3]
    ] * κ

end

function eval_susceptibility(s::DrudeSusceptibility, f)
    @unpack frequency, gamma, sigma, sigma_diag, sigma_offdiag = s
    if iszero(gamma)
        κ = -frequency^2 / (f * (f + im * gamma))
    else
        κ = -frequency^2 / (f * (f))
    end
    [
        sigma_diag[1] sigma_offdiag[1] sigma_offdiag[2]
        sigma_offdiag[1] sigma_diag[2] sigma_offdiag[3]
        sigma_offdiag[2] sigma_offdiag[3] sigma_diag[3]
    ] * κ
end

function eval_susceptibility(s::PoleResidueSusceptibility, f)
    @unpack a, c, sigma, sigma_diag, sigma_offdiag = s
    ω = 2π * f
    κ = c /(im *  ω - a)
    [
        sigma_diag[1] sigma_offdiag[1] sigma_offdiag[2]
        sigma_offdiag[1] sigma_diag[2] sigma_offdiag[3]
        sigma_offdiag[2] sigma_offdiag[3] sigma_diag[3]
    ] * κ
end


text:
############################################################
# joAbstractOperator #######################################
############################################################

# includes
#include("joAbstractOperator/constructors.jl")
include("joAbstractOperator/base_functions.jl")
include("joAbstractOperator/extra_functions.jl")


text:
using Flux, Test
using Flux.CuArrays
using Flux: gpu

@info "Testing GPU Support"

@testset "CuArrays" begin

CuArrays.allowscalar(false)

x = randn(5, 5)
cx = gpu(x)
@test cx isa CuArray

@test Flux.onecold(gpu([1.0, 2.0, 3.0])) == 3

x = Flux.onehotbatch([1, 2, 3], 1:3)
cx = gpu(x)
@test cx isa Flux.OneHotMatrix && cx.data isa CuArray
@test (cx .+ 1) isa CuArray

m = Chain(Dense(10, 5, tanh), Dense(5, 2), softmax)
cm = gpu(m)

@test all(p isa CuArray for p in params(cm))
@test cm(gpu(rand(10, 10))) isa CuArray{Float32,2}

x = [1.,2.,3.]
cx = gpu(x)
@test Flux.crossentropy(x,x) ≈ Flux.crossentropy(cx,cx)
@test Flux.crossentropy(x,x, weight=1.0) ≈ Flux.crossentropy(cx,cx, weight=1.0)
@test Flux.crossentropy(x,x, weight=[1.0;2.0;3.0]) ≈ Flux.crossentropy(cx,cx, weight=cu([1.0;2.0;3.0]))

x = [-1.1491, 0.8619, 0.3127]
y = [1, 1, 0.]
@test Flux.binarycrossentropy.(σ.(x),y) ≈ Array(Flux.binarycrossentropy.(cu(σ.(x)),cu(y)))
@test Flux.logitbinarycrossentropy.(x,y) ≈ Array(Flux.logitbinarycrossentropy.(cu(x),cu(y)))

xs = rand(5, 5)
ys = Flux.onehotbatch(1:5,1:5)
@test collect(cu(xs) .+ cu(ys)) ≈ collect(xs .+ ys)

c = gpu(Conv((2,2),3=>4))
x = gpu(rand(10, 10, 3, 2))
l = c(gpu(rand(10,10,3,2)))
@test gradient(x -> sum(c(x)), x)[1] isa CuArray

c = gpu(CrossCor((2,2),3=>4))
x = gpu(rand(10, 10, 3, 2))
l = c(gpu(rand(10,10,3,2)))
@test gradient(x -> sum(c(x)), x)[1] isa CuArray

end

@testset "onecold gpu" begin
  y = Flux.onehotbatch(ones(3), 1:10) |> gpu;
  @test Flux.onecold(y) isa CuArray
  @test y[3,:] isa CuArray
end

@testset "restructure gpu" begin
  dudt = Dense(1,1) |> gpu
  p,re = Flux.destructure(dudt)
  foo(x) = sum(re(p)(x))
  @test gradient(foo, cu(rand(1)))[1] isa CuArray
end

if CuArrays.has_cudnn()
  @info "Testing Flux/CUDNN"
  include("cudnn.jl")
  include("curnn.jl")
  include("layers.jl")
else
  @warn "CUDNN unavailable, not testing GPU DNN support"
end


text:
####################
# value extraction #
####################

@inline extract_value!(::Type{T}, out::DiffResult, ydual) where {T} =
    DiffResults.value!(d -> value(T,d), out, ydual)
@inline extract_value!(::Type{T}, out, ydual) where {T} = out # ???

@inline function extract_value!(::Type{T}, out, y, ydual) where {T}
    map!(d -> value(T,d), y, ydual)
    copy_value!(out, y)
end

@inline copy_value!(out::DiffResult, y) = DiffResults.value!(out, y)
@inline copy_value!(out, y) = out

###################################
# vector mode function evaluation #
###################################

@generated function dualize(::Type{T}, x::StaticArray) where T
    N = length(x)
    dx = Expr(:tuple, [:(Dual{T}(x[$i], chunk, Val{$i}())) for i in 1:N]...)
    V = StaticArrays.similar_type(x, Dual{T,eltype(x),N})
    return quote
        chunk = Chunk{$N}()
        $(Expr(:meta, :inline))
        return $V($(dx))
    end
end

@inline static_dual_eval(::Type{T}, f, x::StaticArray) where T = f(dualize(T, x))

function vector_mode_dual_eval(f::F, x, cfg::Union{JacobianConfig,GradientConfig}) where {F}
    xdual = cfg.duals
    seed!(xdual, x, cfg.seeds)
    return f(xdual)
end

function vector_mode_dual_eval(f!::F, y, x, cfg::JacobianConfig) where {F}
    ydual, xdual = cfg.duals
    seed!(xdual, x, cfg.seeds)
    seed!(ydual, y)
    f!(ydual, xdual)
    return ydual
end

##################################
# seed construction/manipulation #
##################################

@generated function construct_seeds(::Type{Partials{N,V}}) where {N,V}
    return Expr(:tuple, [:(single_seed(Partials{N,V}, Val{$i}())) for i in 1:N]...)
end

function seed!(duals::AbstractArray{Dual{T,V,N}}, x,
               seed::Partials{N,V} = zero(Partials{N,V})) where {T,V,N}
    duals .= Dual{T,V,N}.(x, Ref(seed))
    return duals
end

function seed!(duals::AbstractArray{Dual{T,V,N}}, x,
               seeds::NTuple{N,Partials{N,V}}) where {T,V,N}
    dual_inds = 1:N
    duals[dual_inds] .= Dual{T,V,N}.(view(x,dual_inds), seeds)
    return duals
end

function seed!(duals::AbstractArray{Dual{T,V,N}}, x, index,
               seed::Partials{N,V} = zero(Partials{N,V})) where {T,V,N}
    offset = index - 1
    dual_inds = (1:N) .+ offset
    duals[dual_inds] .= Dual{T,V,N}.(view(x, dual_inds), Ref(seed))
    return duals
end

function seed!(duals::AbstractArray{Dual{T,V,N}}, x, index,
               seeds::NTuple{N,Partials{N,V}}, chunksize = N) where {T,V,N}
    offset = index - 1
    seed_inds = 1:chunksize
    dual_inds = seed_inds .+ offset
    duals[dual_inds] .= Dual{T,V,N}.(view(x, dual_inds), getindex.(Ref(seeds), seed_inds))
    return duals
end


text:
__precompile__()

module Bio

import BioSequences
import GenomicFeatures
import BioAlignments
import BioStructures
import GeneticVariation
import BioServices
import BioTools

# TODO
#import Phylogenies

const Seq = BioSequences
const Intervals = GenomicFeatures
const Align = BioAlignments
const Structure = BioStructures
const Var = GeneticVariation
const Services = BioServices
const Tools = BioTools

# TODO
#const Phylo = Phylogenies

end  # module Bio


text:
## LSMR works on any matrix with the following methods:
## eltype, size, mul!
## LSMR works on any vector with the following methods:
## eltype, length, rmul!, norm

#############################################################################
##
## Define preconditioned matrix
##
#############################################################################

struct PreconditionedMatrix{TA, Tp, Tx}
    A::TA
    P::Tp
    tmp::Tx
    tmp2::Tx
end

eltype(A::PreconditionedMatrix) = eltype(A.A)
size(A::PreconditionedMatrix, i::Integer) = size(A.A, i)

# MyAdjoint that is a type which is NOT an AbstractMatrix to avoid method ambiguity
struct MyAdjoint{S}
    parent::S
end
adjoint(x::MyAdjoint) = x.parent

Base.adjoint(M::PreconditionedMatrix) = MyAdjoint(M)
function mul!(b, pm::PreconditionedMatrix{TA, Tp, Tx}, a, α::Number, β::Number) where {TA, Tp, Tx}
    ldiv!(pm.tmp, pm.P, a)
    mul!(b, pm.A, pm.tmp, α, β)
    return b
end

function mul!(b, Cpm::MyAdjoint{PreconditionedMatrix{TA, Tp, Tx}}, a, α::Number, β::Number) where {TA, Tp, Tx}
    T = eltype(b)
    pm = adjoint(Cpm)
    β = convert(T, β)
    mul!(pm.tmp, pm.A',  a, one(T), zero(T))
    ldiv!(pm.tmp2, pm.P, pm.tmp)
    if β != one(T)
        if β == zero(T)
            fill!(b, β)
        else
            rmul!(b, β)
        end
    end
    axpy!(α, pm.tmp2, b)
    return b
end

#############################################################################
##
## Define dampened matrix
##
#############################################################################

struct DampenedVector{T1, T2}
    y::T1 # dimension of f(x)
    x::T2 # dimension of x
end
eltype(a::DampenedVector) = promote_type(eltype(a.y), eltype(a.x))
length(a::DampenedVector) = length(a.y) + length(a.x)
function rmul!(a::DampenedVector, α::Number)
    rmul!(a.y, α)
    rmul!(a.x, α)
    return a
end
norm(a::DampenedVector) = sqrt(norm(a.y)^2 + norm(a.x)^2)

struct DampenedMatrix{TA, Tx}
    A::TA
    diagonal::Tx
end
eltype(A::DampenedMatrix) = promote_type(eltype(A.A), eltype(A.diagonal))
function size(A::DampenedMatrix, dim::Integer)
    m, n = size(A.A)
    l = length(A.diagonal)
    dim == 1 ? (m + l) :
    dim == 2 ? n : 1
end
Base.adjoint(M::DampenedMatrix) = MyAdjoint(M)

function mul!(b, mw::DampenedMatrix, a, α::Number, β::Number) where {T}
    if β != 1
        rmul!(b, β)
    end
    mul!(b.y, mw.A, a, α, 1.0)
    map!((z, x, y)-> z + α * x * y, b.x, b.x, a, mw.diagonal)
    return b
end
function mul!(b, Cmw::MyAdjoint{DampenedMatrix{TA, Tx}}, a, α::Number, β::Number) where {TA, Tx}
    T = eltype(b)
    mw = adjoint(Cmw)
    β = convert(T, β)
    if β != one(T)
        if β == zero(T)
            fill!(b, β)
        else
            rmul!(b, β)
        end
    end
    mul!(b, mw.A', a.y, α, one(T))
    map!((z, x, y)-> z + α * x * y, b, b, a.x, mw.diagonal)
    return b
end

#############################################################################
##
## Define inverse preconditioner
##
##############################################################################

struct InverseDiagonal{Tx}
    _::Tx
end
function ldiv!(y, ID::InverseDiagonal, x)
    map!(*, y, x, ID._)
end

#############################################################################
##
## Preconditioner
##
##############################################################################
function getpreconditioner(nls::LeastSquaresProblem, optimizer::Union{Dogleg{LSMR{Nothing, Nothing}}, LevenbergMarquardt{LSMR{Nothing, Nothing}}})
    preconditioner! = function(out, x, J, damp)
        colsumabs2!(out._, J)
        Tout = eltype(out._)
        if damp != 0
            axpy!(one(Tout), damp, out._)
        end
        map!(x -> x > zero(Tout) ? 1 / sqrt(x) : zero(Tout), out._, out._)
        return out
    end
    P = InverseDiagonal(_zeros(nls.x))
    return preconditioner!, P
end

function getpreconditioner(nls::LeastSquaresProblem, optimizer::Union{Dogleg{LSMR}, LevenbergMarquardt{LSMR}})
    return optimizer.solver.preconditioner!, optimizer.solver.P
end

#############################################################################
##
## solve J'J \ J'y
##
## we use LSMR on Ax = y
## with A = J / sqrt(diag(J'J)) (diagonal preconditioner)
##
## LSMR works on any matrix with the following methods:
## eltype, size, mul!
##
##############################################################################

struct LSMRAllocatedSolver{Tx0, Tx1, Tx2, Tx22, Tx3, Tx4, Tx5, Ty} <: AbstractAllocatedSolver
    preconditioner!::Tx0
    P::Tx1
    tmp::Tx2
    tmp2::Tx22
    v::Tx3
    h::Tx4
    hbar::Tx5
    u::Ty
end

function AbstractAllocatedSolver(nls::LeastSquaresProblem, optimizer::Dogleg{LSMR{T1, T2}}) where {T1, T2}
    preconditioner!, P = getpreconditioner(nls, optimizer)
    LSMRAllocatedSolver(preconditioner!, P, _zeros(nls.x), _zeros(nls.x),
        _zeros(nls.x), _zeros(nls.x), _zeros(nls.x), _zeros(nls.y))
end

function ldiv!(x, J, y, A::LSMRAllocatedSolver)
    preconditioner!, P, tmp, tmp2, v, h, hbar, u = A.preconditioner!, A.P, A.tmp, A.tmp2, A.v, A.h, A.hbar, A.u

    # prepare x
    fill!(x, 0)

    # prepare b
    copyto!(u, y)

    # prepare A
    fill!(tmp, 0)
    preconditioner!(P, x, J, 0)
    A = PreconditionedMatrix(J, P, tmp, tmp2)

    # solve
    x, ch = lsmr!(x, A, u, v, h, hbar)
    ldiv!(tmp, P, x)
    copyto!(x, tmp)
    return x, ch.mvps
end

##############################################################################
##
## solve (J'J + diagm(damp)) \ J'y (used in LevenbergMarquardt)
## No need to solve exactly :
## "An Inexact Levenberg-Marquardt Method for Large Sparse Nonlinear Least Squares"
## Weight Holt (1985)
##
## We use LSMR on A x = b with
## A = |J          |  + diagonal preconditioner
##     |diag(damp) |
## b = vcat(y, zeros(damp))
##

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

function getpreconditioner(nls::LeastSquaresProblem, optimizer::LevenbergMarquardt{LSMR{T1, T2}}) where {T1, T2}
    return optimizer.solver.preconditioner!, optimizer.solver.P
end

struct LSMRDampenedAllocatedSolver{Tx0, Tx1, Tx2, Tx22, Tx3, Tx4, Tx5, Tx6, Ty} <: AbstractAllocatedSolver
    preconditioner!::Tx0
    P::Tx1
    tmp::Tx2
    tmp2::Tx22
    v::Tx3
    h::Tx4
    hbar::Tx5
    zerosvector::Tx6
    u::Ty
end

function AbstractAllocatedSolver(nls::LeastSquaresProblem,  optimizer::LevenbergMarquardt{LSMR{T1, T2}}) where {T1, T2}
    preconditioner!, P = getpreconditioner(nls, optimizer)
    LSMRDampenedAllocatedSolver(preconditioner!, P, _zeros(nls.x), _zeros(nls.x), _zeros(nls.x), _zeros(nls.x), _zeros(nls.x),  _zeros(nls.x), _zeros(nls.y))
end

function ldiv!(x, J, y, damp, A::LSMRDampenedAllocatedSolver)
    preconditioner!, P, tmp, tmp2, v, h, hbar, zerosvector, u =
            A.preconditioner!, A.P, A.tmp, A.tmp2, A.v, A.h, A.hbar, A.zerosvector, A.u
    # prepare x
    fill!(x, 0)

    # prepare b
    copyto!(u, y)
    fill!(zerosvector, 0)
    b = DampenedVector(u, zerosvector)

    # prepare A
    fill!(tmp, 0)
    preconditioner!(P, x, J, damp)
    map!(sqrt, damp, damp)
    A = PreconditionedMatrix(DampenedMatrix(J, damp), P, tmp, tmp2)
    # solve
    x, ch = lsmr!(x, A, b, v, h, hbar, btol = 0.5)
    ldiv!(tmp, P, x)
    copyto!(x, tmp)
    return x, ch.mvps
end


text:
@testset "categorical sampler" begin

	cc = dummyChoiceContext()

	@testset "getparams, numparams, paramranges, sample" begin

		@testset "using default construction" begin

			s = DataGenerators.CategoricalSampler(4)

			@testset "numparams and paramranges" begin
			    @test DataGenerators.numparams(s) == 4
			    prs = DataGenerators.paramranges(s)
			    @test typeof(prs) <: Vector{Tuple{Float64,Float64}}
			    @test prs == fill((0.0,1.0), DataGenerators.numparams(s))
			end

			@testset "default params" begin
			    @test DataGenerators.getparams(s) == [0.25, 0.25, 0.25, 0.25]
			end

			@mtestset "default sampling" reps=Main.REPS alpha=Main.ALPHA begin
			    x, trace = DataGenerators.sample(s, (0,1), cc)
			    @test typeof(x) <: Int
			    @mtest_values_are [1,2,3,4] x
		        @mtest_distributed_as Categorical([0.25,0.25,0.25,0.25]) x
			end

		end

		@testset "non-default construction" begin

			s = DataGenerators.CategoricalSampler(5, [0.3,0.2,0.1,0.2,0.2])

			@testset "constructor with params" begin

			    @test getparams(s) == [0.3,0.2,0.1,0.2,0.2]

			    @mtestset "consistent with categorical" reps=Main.REPS alpha=Main.ALPHA begin
	            	x, trace = DataGenerators.sample(s, (0,1), cc)
			        @mtest_distributed_as Categorical([0.3,0.2,0.1,0.2,0.2]) x
			    end

			end

		end

	end

	@testset "setparams" begin

		s = DataGenerators.CategoricalSampler(4)
		prs = DataGenerators.paramranges(s)
		midparams = map(pr->robustmidpoint(pr[1],pr[2]), prs)

		@testset "valid parameters $params" for params in [[0.1,0.2,0.3,0.4], [0.4,0.1,0.1,0.4], [0.3,0.2,0.3,0.2],]

		    DataGenerators.setparams!(s, params)

		    @mtestset "consistent with categorical" reps=Main.REPS alpha=Main.ALPHA begin
	        	x, trace = DataGenerators.sample(s, (0,1), cc)
		        @mtest_distributed_as Categorical(params) x
			end

		end

		@testset "boundary values for parameter index $pidx" for pidx in 1:length(prs)

		    pr = prs[pidx]
		    params = copy(midparams)

			@testset "bound index $bidx" for bidx in 1:2

		        params[pidx] = pr[bidx]
		        DataGenerators.setparams!(s, params)

			    @mtestset "consistent with categorical" reps=Main.REPS alpha=Main.ALPHA begin
		        	x, trace = DataGenerators.sample(s, (0,1), cc)
			        @mtest_distributed_as  Categorical(params./sum(params)) x
				end

				@testset "range check exception" begin
		            params[pidx] = bidx == 1 ? prevfloat(pr[bidx]) : nextfloat(pr[bidx])
		            @test_throws ErrorException DataGenerators.setparams!(s, params)
				end

			end

		end

		@testset "setparams normalises weights" begin

		    DataGenerators.setparams!(s, [0.4, 0.6, 0.7, 0.3])
		    @test getparams(s) == [0.2, 0.3, 0.35, 0.15]

		    @mtestset "consistent with categorical" reps=Main.REPS alpha=Main.ALPHA begin
	        	x, trace = DataGenerators.sample(s, (0,1), cc)
		        @mtest_distributed_as Categorical([0.2,0.3,0.35,0.15]) x
			end

		end

		@testset "setparams adjusts when all weights are zero" begin

		    DataGenerators.setparams!(s, [0.0, 0.0, 0.0, 0.0])
		    @test getparams(s) == [0.25, 0.25, 0.25, 0.25]

		    @mtestset "consistent with categorical" reps=Main.REPS alpha=Main.ALPHA begin
	        	x, trace = DataGenerators.sample(s, (0,1), cc)
		        @mtest_distributed_as Categorical([0.25,0.25,0.25,0.25]) x
			end

		end

		@testset "setparams with wrong number of parameters" begin
		    @test_throws ErrorException DataGenerators.setparams!(s, midparams[1:end-1])
		    @test_throws ErrorException DataGenerators.setparams!(s, [midparams; 0.5])
		end

	end

	@testset "estimateparams" begin

		otherparams = [0.35, 0.15, 0.1, 0.25, 0.15]

		@testset "equal bounds" begin

		    params = [0.2, 0.2, 0.2, 0.2, 0.2]
		    s1 = DataGenerators.CategoricalSampler(5, params)
		    s2 = DataGenerators.CategoricalSampler(5, otherparams)
		    traces = map(1:100) do i
		        x, trace = DataGenerators.sample(s1, (0,1), cc)
		        trace
		    end
		    estimateparams!(s2, traces)

		    @mtestset "consistent with categorical" reps=Main.REPS alpha=Main.ALPHA begin
	        	x, trace = DataGenerators.sample(s2, (0,1), cc)
		        @mtest_distributed_as Categorical(params) x
			end

		end

		@testset "last category never sampled" begin

		    params = [0.3, 0.4, 0.2, 0.1, 0.0]
		    s1 = DataGenerators.CategoricalSampler(5, params)
		    s2 = DataGenerators.CategoricalSampler(5, otherparams)
		    traces = map(1:100) do i
		        x, trace = DataGenerators.sample(s1, (0,1), cc)
		        trace
		    end
		    estimateparams!(s2, traces)

		    @mtestset "consistent with categorical" reps=Main.REPS alpha=Main.ALPHA begin
	        	x, trace = DataGenerators.sample(s2, (0,1), cc)
		        @mtest_distributed_as Categorical(params) x
			end

		end

		@testset "too few traces" begin

		    params = [0.1, 0.2, 0.4, 0.2, 0.1]
		    s1 = DataGenerators.CategoricalSampler(5, params)
		    s2 = DataGenerators.CategoricalSampler(5, otherparams)
		    traces = map(1:0) do i
		        x, trace = DataGenerators.sample(s1, (0,1), cc)
		        trace
		    end
		    estimateparams!(s2, traces)

		    @mtestset "consistent with categorical" reps=Main.REPS alpha=Main.ALPHA begin
	        	x, trace = DataGenerators.sample(s2, (0,1), cc)
		        @mtest_distributed_as Categorical(otherparams) x
			end

		end

	end

end


text:
using ProgressMeter, War

"""
    war_result_vector(r::Vector{Int})::Vector{Int}

Return information about a game of war as a vector with the following components:
* [1] length of the game (number of battles)
* [2] number of single wars
* [3] number of double wars
* [4] number of triple wars
* [5] number of quadruple (or longer) wars
"""
function war_result_vector(r::Vector{Int})::Vector{Int}
    result = zeros(Int, 5)
    dr = abs.(diff(r))
    result[1] = length(dr)

    result[2] = count(x -> 1 < x <= 5, dr)
    result[3] = count(x -> 5 < x <= 9, dr)
    result[4] = count(x -> 9 < x <= 13, dr)
    result[5] = count(x -> 14 < x, dr)

    return result
end

function war_report(r::Vector{Int})
    v = war_result_vector(r)
    println("Number of moves:          $(v[1])")
    println("Number of single wars:    $(v[2])")
    println("Number of double wars:    $(v[3])")
    println("Number of triple wars:    $(v[4])")
    println("Number of super wars:     $(v[5])")
end

"""
    mean_report(reps::Int = 100_000)

Play `reps` games of war and report statistics on the average
number of battles, wars, double wars, and so on.
"""
function mean_report(reps::Int = 100_000)
    v = zeros(Int, 5)
    PM = Progress(reps)
    for k = 1:reps
        v += war_result_vector(play_war())
        next!(PM)
    end
    v /= reps
    println("Number of games played:           $reps")
    println("Average number of turns:          $(v[1])")
    println("Average number of single wars:    $(v[2])")
    println("Average number of double wars:    $(v[3])")
    println("Average number of triple wars:    $(v[4])")
    println("Average number of super wars:     $(v[5])")
end


text:
const SSD1306 = "libssd1306_i2c"
const ssd1306_err_t = Ptr{Cvoid}
const ssd1306_framebuffer_t = Ptr{Cvoid}
const ssd1306_i2c_t = Ptr{Cvoid}

ssd1306_err_create() = @ccall SSD1306.ssd1306_err_create(C_NULL::Ptr{Cvoid})::ssd1306_err_t
ssd1306_framebuffer_create(width, height, err)::ssd1306_framebuffer_t = @ccall SSD1306.ssd1306_framebuffer_create(width::UInt8, height::UInt8, err::Ptr{Cvoid})::Ptr{Cvoid}
ssd1306_framebuffer_draw_line(fbp, x0, y0, x1, y1, color)::Int = @ccall SSD1306.ssd1306_framebuffer_draw_line(fbp::Ptr{Cvoid}, x0::UInt8, y0::UInt8, x1::UInt8, y1::UInt8, color::Bool)::Cint
ssd1306_err_destroy(err) = @ccall SSD1306.ssd1306_err_destroy(err::ssd1306_err_t)::Cvoid
ssd1306_i2c_display_update(oled::ssd1306_i2c_t, fbp::ssd1306_framebuffer_t) = @ccall SSD1306.ssd1306_i2c_display_update(oled::Ptr{Cvoid}, fbp::Ptr{Cvoid})::Cint
ssd1306_i2c_open(filename, daddr, width, height, logerr)::ssd1306_i2c_t = @ccall SSD1306.ssd1306_i2c_open(filename::Cstring, daddr::UInt8, width::UInt8, height::UInt8, logerr::Ptr{Cvoid})::Ptr{Cvoid}
ssd1306_i2c_display_initialize(oled::ssd1306_i2c_t) = @ccall SSD1306.ssd1306_i2c_display_initialize(oled::Ptr{Cvoid})::Cint
ssd1306_i2c_display_clear(oled::ssd1306_i2c_t) = @ccall SSD1306.ssd1306_i2c_display_clear(oled::Ptr{Cvoid})::Cint
ssd1306_framebuffer_destroy(fbp::ssd1306_framebuffer_t) = @ccall SSD1306.ssd1306_framebuffer_destroy(fbp::Ptr{Cvoid})::Cvoid
ssd1306_i2c_close(oled::ssd1306_i2c_t) = @ccall SSD1306.ssd1306_i2c_close(oled::Ptr{Cvoid})::Cvoid

function main()
    oled = ssd1306_i2c_open("/dev/i2c-1", 0x3c, 128, 64, C_NULL)
    err = ssd1306_i2c_display_initialize(oled)
    ssd1306_i2c_display_clear(oled)
    errp = ssd1306_err_create()
    fbp = ssd1306_framebuffer_create(128, 64, errp)
    err = ssd1306_framebuffer_draw_line(fbp, 0, 0, 127, 63, true)
    ssd1306_i2c_display_update(oled, fbp)
    sleep(5)
    ssd1306_framebuffer_destroy(fbp)
    ssd1306_i2c_close(oled)
    ssd1306_err_destroy(errp)
end

main()


text:
using Oceananigans.Operators: Vᶜᶜᶜ
using Oceananigans.Fields: ZeroField

@inline _advective_tracer_flux_x(args...) = advective_tracer_flux_x(args...)
@inline _advective_tracer_flux_y(args...) = advective_tracer_flux_y(args...)
@inline _advective_tracer_flux_z(args...) = advective_tracer_flux_z(args...)

const ZeroU = NamedTuple{(:u, :v, :w), Tuple{ZeroField, ZeroField, ZeroField}}

@inline div_Uc(i, j, k, grid, advection, ::ZeroU, c) = zero(eltype(grid))
@inline div_Uc(i, j, k, grid, advection, U, ::ZeroField) = zero(eltype(grid))

@inline div_Uc(i, j, k, grid, ::Nothing, U, c) = zero(eltype(grid))
@inline div_Uc(i, j, k, grid, ::Nothing, ::ZeroU, c) = zero(eltype(grid))
@inline div_Uc(i, j, k, grid, ::Nothing, U, ::ZeroField) = zero(eltype(grid))

#####
##### Tracer advection operator
#####

"""
    div_uc(i, j, k, grid, advection, U, c)

Calculates the divergence of the flux of a tracer quantity ``c`` being advected by
a velocity field, ``𝛁⋅(𝐯 c)``,

    1/V * [δxᶜᵃᵃ(Ax * u * ℑxᶠᵃᵃ(c)) + δyᵃᶜᵃ(Ay * v * ℑyᵃᶠᵃ(c)) + δzᵃᵃᶜ(Az * w * ℑzᵃᵃᶠ(c))]

which ends up at the location `ccc`.
"""
@inline function div_Uc(i, j, k, grid, advection, U, c)
    1/Vᶜᶜᶜ(i, j, k, grid) * (δxᶜᵃᵃ(i, j, k, grid, _advective_tracer_flux_x, advection, U.u, c) +
                             δyᵃᶜᵃ(i, j, k, grid, _advective_tracer_flux_y, advection, U.v, c) +
                             δzᵃᵃᶜ(i, j, k, grid, _advective_tracer_flux_z, advection, U.w, c))
end


text:
using Test
using KernelAbstractions
using CUDAKernels

import CloudMicrophysics
import CLIMAParameters
import Thermodynamics

const CM = CloudMicrophysics
const CP = CLIMAParameters
const TD = Thermodynamics

struct EarthParameterSet <: CP.AbstractEarthParameterSet end
const prs = EarthParameterSet()

const AM = CloudMicrophysics.AerosolModel
const AA = CloudMicrophysics.AerosolActivation
const CM0 = CloudMicrophysics.Microphysics_0M
const CM1 = CloudMicrophysics.Microphysics_1M
const CP0 = CLIMAParameters.Atmos.Microphysics_0M

const liquid = CM.Microphysics_1M.LiquidType()
const ice = CM.Microphysics_1M.IceType()
const rain = CM.Microphysics_1M.RainType()
const snow = CM.Microphysics_1M.SnowType()

if get(ARGS, 1, "Array") == "CuArray"
    using CUDA
    ArrayType = CUDA.CuArray
    CUDA.allowscalar(false)
    device(::Type{T}) where {T <: CuArray} = CUDADevice()
else
    ArrayType = Array
    device(::Type{T}) where {T <: Array} = CPU()
end

@show ArrayType

@kernel function test_aerosol_activation_kernel!(
    prs,
    output::AbstractArray{FT},
    r,
    stdev,
    N,
    ϵ,
    ϕ,
    M,
    ν,
    ρ,
    κ,
) where {FT}

    i = @index(Group, Linear)

    # atmospheric conditions (taken from aerosol activation tests)
    T::FT = 294.0       # air temperature K
    p::FT = 100000.0    # air pressure Pa
    w::FT = 0.5         # vertical velocity m/s
    p_vs::FT = TD.saturation_vapor_pressure(prs, T, TD.Liquid())
    q_vs::FT = 1 / (1 - CP.Planet.molmass_ratio(prs) * (p_vs - p) / p_vs)
    # water vapor specific humidity (saturated)
    q = TD.PhasePartition(q_vs, FT(0.0), FT(0.0))

    args = (T, p, w, q)

    @inbounds begin
        mode_B = AM.Mode_B(
            FT(r[i]),
            FT(stdev[i]),
            FT(N[i]),
            (FT(1.0),),
            (FT(ϵ[i]),),
            (FT(ϕ[i]),),
            (FT(M[i]),),
            (FT(ν[i]),),
            (FT(ρ[i]),),
            1,
        )
        mode_κ = AM.Mode_κ(
            FT(r[i]),
            FT(stdev[i]),
            FT(N[i]),
            (FT(1.0),),
            (FT(1.0),),
            (FT(M[i]),),
            (FT(κ[i]),),
            1,
        )

        arsl_dst_B = AM.AerosolDistribution((mode_B,))
        arsl_dst_κ = AM.AerosolDistribution((mode_κ,))

        output[1, i] = AA.mean_hygroscopicity_parameter(prs, arsl_dst_B)[1]
        output[2, i] = AA.mean_hygroscopicity_parameter(prs, arsl_dst_κ)[1]

        output[3, i] = AA.total_N_activated(prs, arsl_dst_B, args...)
        output[4, i] = AA.total_N_activated(prs, arsl_dst_κ, args...)

        output[5, i] = AA.total_M_activated(prs, arsl_dst_B, args...)
        output[6, i] = AA.total_M_activated(prs, arsl_dst_κ, args...)
    end
end

@kernel function test_0_moment_micro_kernel!(
    prs,
    output::AbstractArray{FT},
    liquid_frac,
    qc,
    qt,
) where {FT}

    i = @index(Group, Linear)

    @inbounds begin
        ql::FT = qc[i] * liquid_frac[i]
        qi::FT = (1 - liquid_frac[i]) * qc[i]
        q = TD.PhasePartition(FT(qt[i]), ql, qi)

        output[1, i] = CM0.remove_precipitation(prs, q)

        _τ_precip = CP0.τ_precip(prs)
        _qc_0 = CP0.qc_0(prs)

        output[2, i] = -max(0, ql + qi - _qc_0) / _τ_precip
    end
end

@kernel function test_1_moment_micro_accretion_kernel!(
    prs,
    output::AbstractArray{FT},
    ρ,
    qt,
    qi,
    qs,
    ql,
    qr,
) where {FT}

    i = @index(Group, Linear)

    @inbounds begin
        output[1, i] = CM1.accretion(prs, liquid, rain, ql[i], qr[i], ρ[i])
        output[2, i] = CM1.accretion(prs, ice, snow, qi[i], qs[i], ρ[i])
        output[3, i] = CM1.accretion(prs, liquid, snow, ql[i], qs[i], ρ[i])
        output[4, i] = CM1.accretion(prs, ice, rain, qi[i], qr[i], ρ[i])
        output[5, i] = CM1.accretion_rain_sink(prs, qi[i], qr[i], ρ[i])
        output[6, i] =
            CM1.accretion_snow_rain(prs, snow, rain, qs[i], qr[i], ρ[i])
        output[7, i] =
            CM1.accretion_snow_rain(prs, rain, snow, qr[i], qs[i], ρ[i])
    end
end

@kernel function test_1_moment_micro_snow_melt_kernel!(
    prs,
    output::AbstractArray{FT},
    ρ,
    T,
    qs,
) where {FT}

    i = @index(Group, Linear)

    @inbounds begin
        output[i] = CM1.snow_melt(prs, qs[i], ρ[i], T[i])
    end
end

@testset "Aerosol activation kernels" begin
    FT = Float32
    data_length = 2
    output = ArrayType(Array{FT}(undef, 6, data_length))
    fill!(output, -44.0)

    dev = device(ArrayType)
    work_groups = (1,)
    ndrange = (data_length,)

    r = ArrayType([0.243 * 1e-6, 1.5 * 1e-6])
    stdev = ArrayType([1.4, 2.1])
    N = ArrayType([100.0 * 1e6, 1.0 * 1e6])
    ϵ = ArrayType([1.0, 1.0])
    ϕ = ArrayType([1.0, 0.9])
    M = ArrayType([0.132, 0.058443])
    ν = ArrayType([3.0, 2.0])
    ρ = ArrayType([1770.0, 2170.0])
    κ = ArrayType([0.53, 1.12])

    kernel! = test_aerosol_activation_kernel!(dev, work_groups)
    event =
        kernel!(prs, output, r, stdev, N, ϵ, ϕ, M, ν, ρ, κ, ndrange = ndrange)
    wait(dev, event)

    # test if all aerosol activation output is positive
    @test all(Array(output)[:, :] .>= FT(0))
    # test if higroscopicity parameter is the same for κ and B modes
    @test all(isapprox(Array(output)[1, :], Array(output)[2, :], rtol = 0.3))
    # test if the number and mass activated are the same for κ and B modes
    @test all(isapprox(Array(output)[3, :], Array(output)[4, :], rtol = 1e-5))
    @test all(isapprox(Array(output)[5, :], Array(output)[6, :], rtol = 1e-5))
end

@testset "0-moment microphysics kernels" begin
    FT = Float32
    data_length = 3
    output = ArrayType(Array{FT}(undef, 2, data_length))
    fill!(output, -44.0)

    dev = device(ArrayType)
    work_groups = (1,)
    ndrange = (data_length,)

    liquid_frac = ArrayType([0.0, 0.5, 1.0])
    qt = ArrayType([13e-3, 13e-3, 13e-3])
    qc = ArrayType([3e-3, 4e-3, 5e-3])

    kernel! = test_0_moment_micro_kernel!(dev, work_groups)
    event = kernel!(prs, output, liquid_frac, qc, qt, ndrange = ndrange)
    wait(dev, event)

    # test 0-moment rain removal is callable and returns a reasonable value
    @test all(isequal(Array(output)[1, :], Array(output)[2, :]))
end

@testset "1-moment microphysics kernels" begin
    FT = Float32
    data_length = 2
    output = ArrayType(Array{FT}(undef, 7, data_length))
    fill!(output, -44.0)

    dev = device(ArrayType)
    work_groups = (1,)
    ndrange = (data_length,)

    ρ = ArrayType([1.2, 1.2])
    qt = ArrayType([0.0, 20e-3])
    qi = ArrayType([0.0, 5e-4])
    qs = ArrayType([0.0, 5e-4])
    ql = ArrayType([0.0, 5e-4])
    qr = ArrayType([0.0, 5e-4])

    kernel! = test_1_moment_micro_accretion_kernel!(dev, work_groups)
    event = kernel!(prs, output, ρ, qt, qi, qs, ql, qr, ndrange = ndrange)
    wait(dev, event)

    # test 1-moment accretion is callable and returns a reasonable value
    @test all(Array(output)[:, 1] .== FT(0))
    @test Array(output)[1, 2] ≈ 1.4150106417043544e-6
    @test Array(output)[2, 2] ≈ 2.453070979562392e-7
    @test Array(output)[3, 2] ≈ 2.453070979562392e-7
    @test Array(output)[4, 2] ≈ 1.768763302130443e-6
    @test Array(output)[5, 2] ≈ 3.085229094251214e-5
    @test Array(output)[6, 2] ≈ 2.1705865794293408e-4
    @test Array(output)[7, 2] ≈ 6.0118801860768854e-5

    data_length = 3
    output = ArrayType(Array{FT}(undef, 1, data_length))
    fill!(output, -44.0)

    dev = device(ArrayType)
    work_groups = (1,)
    ndrange = (data_length,)

    ρ = ArrayType([1.2, 1.2, 1.2])
    T = ArrayType([273.15 + 2, 273.15 + 2, 273.15 - 2])
    qs = ArrayType([1e-4, 0.0, 1e-4])

    kernel! = test_1_moment_micro_snow_melt_kernel!(dev, work_groups)
    event = kernel!(prs, output, ρ, T, qs, ndrange = ndrange)
    wait(dev, event)

    # test if 1-moment snow melt is callable and returns reasonable values
    @test Array(output)[1] ≈ 9.518235437405256e-6
    @test Array(output)[2] ≈ 0.0
    @test Array(output)[3] ≈ 0.0
end


text:
"""
    Map the pairwise force to a 2d kernel so at best we get one gpu thread
    calculating the force (or the abscence of it) between one pair of particles,
    with all pairs being done simultaneously!
"""
function ForceKernel(Y,F,R,k,N,t)
    index_x = (blockIdx().x-1)*blockDim().x + threadIdx().x # starts at 1 for blockidx !
    stride_x = gridDim().x*blockDim().x

    index_y = (blockIdx().y-1)*blockDim().y + threadIdx().y
    stride_y = gridDim().y*blockDim().y

    for i in index_x:stride_x:N
        for j in index_y:stride_y:N
            if (i == j)
                continue
            end
            r = 2*R # min dist by radii
            @inbounds rx = Y[t-1,j,1]-Y[t-1,i,1] # sep vector x
            @inbounds ry = Y[t-1,j,2]-Y[t-1,i,2] # sep vector y
            d2 = rx*rx+ry*ry
            if (d2 < r*r) # distance test without sqrt
                # apply a force to handle overlap
                d = CUDA.sqrt(d2) # need it now
                mag = -k*(r-d)

                F[i,1] += mag*rx/d
                F[i,2] += mag*ry/d
            end
        end
    end
end

"""
    Do an Euler-Maruyama step with the calculated forces and constants, 1d kernel
    so at best one thread for each particle.
"""
function StepKernel(Y,F,V,v0,Dr,Dt,dt,RAND,N,t)
    index_x = (blockIdx().x-1)*blockDim().x + threadIdx().x # starts at 1 for blockidx !
    stride_x = gridDim().x*blockDim().x

    for i in index_x:stride_x:N
        # angular diffusion
        @inbounds Y[t,i,3] = Y[t-1,i,3] + RAND[t-1,i,3]*CUDA.sqrt(2.0*Dr*dt)
        DT = CUDA.sqrt(2.0*Dt*dt)
        # save increments for boundary kernel
        @inbounds V[t,i,1] = dt * (v0*CUDA.cos(Y[t-1,i,3]) + F[i,1] + RAND[t-1,i,1]*DT)
        @inbounds V[t,i,2] = dt * (v0*CUDA.sin(Y[t-1,i,3]) + F[i,2] + RAND[t-1,i,2]*DT)
        # update (subject to boundary kernel)
        @inbounds Y[t,i,1] = Y[t-1,i,1] + V[t,i,1]
        @inbounds Y[t,i,2] = Y[t-1,i,2] + V[t,i,2]

        # reset force while we are here
        F[i,1] = 0.0
        F[i,2] = 0.0
    end
end

"""
    Keep particles in the box [0,L]x[0,L] by ``elastic collision''
"""
function BoundsKernel(Y,F,R,V,L,N,t)
    index_x = (blockIdx().x-1)*blockDim().x + threadIdx().x # starts at 1 for blockidx !
    stride_x = gridDim().x*blockDim().x

    for i in index_x:stride_x:N
        ux = 0. |> Float32 # be careful this is needed  for atan2...
        uy = 0. |> Float32
        flag = false
        ang = 0. |> Float32
        @inbounds vx = V[t,i,1]
        @inbounds vy = V[t,i,2]

        @inbounds OUTX = Y[t,i,1] < R || Y[t,i,1] > L-R
        @inbounds OUTY = Y[t,i,2] < R || Y[t,i,2] > L-R

        if (OUTX)
            ux = -1.0*vx |> Float32
            ang = CUDA.atan2(vy,ux)
            flag = true
        end

        if (OUTY)
            uy = -1.0*vy |> Float32
            if (flag)
                ang = CUDA.atan2(uy,ux)
            else
                ang = CUDA.atan2(uy,vx)
            end
            flag = true
        end

        if (flag)
            @inbounds Y[t,i,1] += ux
            @inbounds Y[t,i,2] += uy
            @inbounds Y[t,i,3] = ang
        end
    end
end

"""
    Encapsulate GPU kernels in a function, computes data in batches to
        balance host/device memory shunting with memory on device.
"""
function GPUSteps(X,steps::Int=1;R=1.0,v0=1.0,dt=0.00166,Dt=0.0,Dr=0.005,L=10,k=300.,threads=256)
    Y = zeros(1+steps,size(X,1),size(X,2))
    Y[1,:,:] = X
    Y = Y|>cu
    V = zeros(size(Y))|>cu
    # forces
    F = zeros(size(Y,2),2)|>cu

    # curand not quite here yet?? https://discourse.julialang.org/t/how-to-generate-a-random-number-in-cuda-kernel-function/50364
    # will have to precompute random variables
    RAND = zeros(steps,size(X,1),3)
    if Dt > 0
        RAND[:,:,1:2] = randn(steps,size(X,1),2)
    end

    if Dr > 0
        RAND[:,:,3] = randn(steps,size(X,1))
    end
    RAND=RAND|>cu

    for i in 2:steps+1

        CUDA.@sync begin
            @cuda threads=threads blocks=ceil(Int,size(X,1)/threads) ForceKernel(Y,F,R,k,size(X,1),i)
        end

        CUDA.@sync begin
            @cuda threads=threads blocks=ceil(Int,size(X,1)/threads) StepKernel(Y,F,V,v0,Dr,Dt,dt,RAND,size(X,1),i)
        end

        CUDA.@sync begin
            @cuda threads=threads blocks=ceil(Int,size(X,1)/threads) BoundsKernel(Y,F,R,V,L,size(X,1),i)
        end

    end

    return Y[2:end,:,:] |> Array
end

max_stars_repo_name: Jerboa-app/ActiveBrownianParticles

max_issues_repo_name: Jerboa-app/ActiveBrownianParticles

max_forks_repo_name: Jerboa-app/ActiveBrownianParticles


text:
using OrdinaryDiffEq, ParameterizedFunctions, ForwardDiff

f = @ode_def begin
  dx = a*x - b*x*y
  dy = -c*y + x*y
end a b c

p = [1.5,1.0,3.0]
prob = ODEProblem(f,[1.0;1.0],(0.0,10.0),p)
t = 0.0:0.5:10.0

function G(p)
  tmp_prob = remake(prob,u0=convert.(eltype(p),prob.u0),p=p)
  sol = solve(tmp_prob,Vern9(),abstol=1e-14,reltol=1e-14,saveat=t)
  A = convert(Array,sol)
  sum(((1 .- A).^2)./2)
end
# G([1.5,1.0,3.0])
# res2 = ForwardDiff.gradient(G,[1.5,1.0,3.0])

@info("Running forwarddiff")
# G([1.5,1.0,3.0])
es2 = ForwardDiff.gradient(G,[1.5,1.0,3.0])
# using Zygote
# Zygote.gradient(G,[1.5,1.0,3.0])


text:
import Base: adjoint, *, +, -, getindex, hcat, vcat, reshape
export jacobian

###### ' ######
adjoint(L::T) where {T <: AbstractOperator} = AdjointOperator(L)

######+,-######
(+)(L::T) where {T <: AbstractOperator} = L
(-)(L::T) where {T <: AbstractOperator} = Scale(-1.0, L)
(+)(L1::AbstractOperator, L2::AbstractOperator) = Sum(L1,  L2 )
(-)(L1::AbstractOperator, L2::AbstractOperator) = Sum(L1, -L2 )

###### * ######
function (*)(L::AbstractOperator, b::AbstractArray)
  C = codomainType(L)
  if typeof(C) <: Tuple
    y = ArrayPartition(zeros.(codomainType(L), size(L, 1))...)
  else
	  y = zeros(codomainType(L), size(L, 1))
  end
	mul!(y, L, b)
	return y
end

#(*)(L::AbstractOperator, b::Tuple) = (*)(L, ArrayPartition(b...))

*(coeff::T, L::AbstractOperator) where {T<:Number} = Scale(coeff,L)
*(L1::AbstractOperator, L2::AbstractOperator) = Compose(L1,L2)

# getindex
function getindex(A::AbstractOperator,idx...)
	if ndoms(A,2) == 1
		Gout = GetIndex(codomainType(A),size(A,1),idx)
		return Gout*A
	elseif length(idx) == 1  && ndoms(A,2) == length(idx[1])
		return permute(A,idx[1])
	else
		error("cannot split operator of type $(typeof(A))")
	end
end

#get index of HCAT returns HCAT (or Operator)
function getindex(H::HCAT, idx::Union{AbstractArray,Int})

	unfolded = vcat([[i... ] for i in H.idxs]...)
	if length(idx) == length(unfolded)
		return permute(H,idx)
	else
		new_H = ()
		for i in idx
			for ii in eachindex(H.idxs)
				if i in H.idxs[ii]
					if typeof(H.idxs[ii]) <: Int
						new_H = (new_H...,H.A[ii])
					else
					error("cannot split operator")
					end
				end
			end
		end
		return HCAT(new_H,H.buf)
	end
end

#get index of HCAT returns HCAT (or Operator)
function getindex(H::VCAT, idx::Union{AbstractArray,Int})

	unfolded = vcat([[i... ] for i in H.idxs]...)
	if length(idx) == length(unfolded)
		return permute(H,idx)
	else
		new_H = ()
		for i in idx
			for ii in eachindex(H.idxs)
				if i in H.idxs[ii]
					if typeof(H.idxs[ii]) <: Int
						new_H = (new_H...,H.A[ii])
					else
					error("cannot split operator")
					end
				end
			end
		end
		return VCAT(new_H,H.buf)
	end
end

getindex(H::A, idx::Union{AbstractArray,Int}) where {L <: HCAT, D, S, A<: AffineAdd{L,D,S}} =
AffineAdd(getindex(H.A, idx), H.d, S)

# get index of scale
getindex(A::S,idx...) where {T, L, S <:Scale{T,L}} = Scale(A.coeff,A.coeff_conj,getindex(A.A,idx...))

hcat(L::Vararg{AbstractOperator}) = HCAT(L...)
vcat(L::Vararg{AbstractOperator}) = VCAT(L...)

###### reshape ######
reshape(L::A, idx::NTuple{N,Int}) where {N,A<:AbstractOperator} = Reshape(L,idx)
reshape(L::A, idx::Vararg{Int}) where {A<:AbstractOperator} = Reshape(L,idx)

###### jacobian ######
jacobian = Jacobian


text:
function CreateQueue(env::AWSEnv, msg::CreateQueueType)
    sqsresp::SQSResponse = call_sqs(env, "CreateQueue" , msg)
    if  (sqsresp.pd != nothing) && (sqsresp.obj == nothing)
        sqsresp.obj = CreateQueueResponseType(sqsresp.pd)
    end
    sqsresp
end
function CreateQueue(env::AWSEnv; kwargs...)
    msg=CreateQueueType()
    for p in kwargs
        setfield!(msg, p[1], p[2])
    end
    CreateQueue(env, msg)
end
export CreateQueue

function DeleteQueue(env::AWSEnv, msg::DeleteQueueType)
    sqsresp::SQSResponse = call_sqs(env, "DeleteQueue" , msg)
    if  (sqsresp.pd != nothing) && (sqsresp.obj == nothing)
        sqsresp.obj = DeleteQueueResponseType(sqsresp.pd)
    end
    sqsresp
end
function DeleteQueue(env::AWSEnv; kwargs...)
    msg=DeleteQueueType()
    for p in kwargs
        setfield!(msg, p[1], p[2])
    end
    DeleteQueue(env, msg)
end
export DeleteQueue

function GetQueueAttributes(env::AWSEnv, msg::GetQueueAttributesType)
    sqsresp::SQSResponse = call_sqs(env, "GetQueueAttributes" , msg)
    if  (sqsresp.pd != nothing) && (sqsresp.obj == nothing)
        sqsresp.obj = GetQueueAttributesResponseType(sqsresp.pd)
    end
    sqsresp
end
function GetQueueAttributes(env::AWSEnv; kwargs...)
    msg=GetQueueAttributesType()
    for p in kwargs
        setfield!(msg, p[1], p[2])
    end
    GetQueueAttributes(env, msg)
end
export GetQueueAttributes

function GetQueueUrl(env::AWSEnv, msg::GetQueueUrlType)
    sqsresp::SQSResponse = call_sqs(env, "GetQueueUrl" , msg)
    if  (sqsresp.pd != nothing) && (sqsresp.obj == nothing)
        sqsresp.obj = GetQueueUrlResponseType(sqsresp.pd)
    end
    sqsresp
end
function GetQueueUrl(env::AWSEnv; kwargs...)
    msg=GetQueueUrlType()
    for p in kwargs
        setfield!(msg, p[1], p[2])
    end
    GetQueueUrl(env, msg)
end
export GetQueueUrl

function ListQueues(env::AWSEnv, msg::ListQueuesType)
    sqsresp::SQSResponse = call_sqs(env, "ListQueues" , msg)
    if  (sqsresp.pd != nothing) && (sqsresp.obj == nothing)
        sqsresp.obj = ListQueuesResponseType(sqsresp.pd)
    end
    sqsresp
end
function ListQueues(env::AWSEnv; kwargs...)
    msg=ListQueuesType()
    for p in kwargs
        setfield!(msg, p[1], p[2])
    end
    ListQueues(env, msg)
end
export ListQueues

function ReceiveMessage(env::AWSEnv, msg::ReceiveMessageType)
    sqsresp::SQSResponse = call_sqs(env, "ReceiveMessage" , msg, false)
    if  (sqsresp.pd != nothing) && (sqsresp.obj == nothing)
        sqsresp.obj = ReceiveMessageResponseType(sqsresp.pd)
    end
    sqsresp
end
function ReceiveMessage(env::AWSEnv; kwargs...)
    msg=ReceiveMessageType()
    for p in kwargs
        setfield!(msg, p[1], p[2])
    end
    ReceiveMessage(env, msg)
end
export ReceiveMessage

function SendMessage(env::AWSEnv, msg::SendMessageType)
    sqsresp::SQSResponse = call_sqs(env, "SendMessage" , msg, false)
    if  (sqsresp.pd != nothing) && (sqsresp.obj == nothing)
        sqsresp.obj = SendMessageResponseType(sqsresp.pd)
    end
    sqsresp
end
function SendMessage(env::AWSEnv; kwargs...)
    msg=SendMessageType()
    for p in kwargs
        setfield!(msg, p[1], p[2])
    end
    SendMessage(env, msg)
end
export SendMessage

function SetQueueAttributes(env::AWSEnv, msg::SetQueueAttributesType)
    sqsresp::SQSResponse = call_sqs(env, "SetQueueAttributes" , msg)
    if  (sqsresp.pd != nothing) && (sqsresp.obj == nothing)
        sqsresp.obj = SetQueueAttributesResponseType(sqsresp.pd)
    end
    sqsresp
end
function SetQueueAttributes(env::AWSEnv; kwargs...)
    msg=SetQueueAttributesType()
    for p in kwargs
        setfield!(msg, p[1], p[2])
    end
    SetQueueAttributes(env, msg)
end
export SetQueueAttributes

function DeleteMessage(env::AWSEnv, msg::DeleteMessageType)
    sqsresp::SQSResponse = call_sqs(env, "DeleteMessage" , msg)
    if  (sqsresp.pd != nothing) && (sqsresp.obj == nothing)
        sqsresp.obj = DeleteMessageResponseType(sqsresp.pd)
    end
    sqsresp
end
function DeleteMessage(env::AWSEnv; kwargs...)
    msg=DeleteMessageType()
    for p in kwargs
        setfield!(msg, p[1], p[2])
    end
    DeleteMessage(env, msg)
end
export DeleteMessage

#=
ChangeMessageVisibility
ChangeMessageVisibilityBatch
ListDeadLetterSourceQueues
PurgeQueue
AddPermission
RemovePermission
=#

function SendMessageBatch(env::AWSEnv, msg::SendMessageBatchType)
    sqsresp::SQSResponse = call_sqs(env, "SendMessageBatch" , msg)
    if  (sqsresp.pd != nothing) && (sqsresp.obj == nothing)
        sqsresp.obj = SendMessageBatchResponseType(sqsresp.pd)
    end
    sqsresp
end
function SendMessageBatch(env::AWSEnv; kwargs...)
    msg=SendMessageBatchType()
    for p in kwargs
        setfield!(msg, p[1], p[2])
    end
    SendMessageBatch(env, msg)
end
export SendMessageBatch

function DeleteMessageBatch(env::AWSEnv, msg::DeleteMessageBatchType)
    sqsresp::SQSResponse = call_sqs(env, "DeleteMessageBatch" , msg)
    if  (sqsresp.pd != nothing) && (sqsresp.obj == nothing)
        sqsresp.obj = DeleteMessageBatchResponseType(sqsresp.pd)
    end
    sqsresp
end
function DeleteMessageBatch(env::AWSEnv; kwargs...)
    msg=DeleteMessageBatchType()
    for p in kwargs
        setfield!(msg, p[1], p[2])
    end
    DeleteMessageBatch(env, msg)
end
export DeleteMessageBatch

function ChangeMessageVisibility(env::AWSEnv, msg::ChangeMessageVisibilityType)
    sqsresp::SQSResponse = call_sqs(env, "ChangeMessageVisibility" , msg)
    if  (sqsresp.pd != nothing) && (sqsresp.obj == nothing)
        sqsresp.obj = ChangeMessageVisibilityResponseType(sqsresp.pd)
    end
    sqsresp
end
function ChangeMessageVisibility(env::AWSEnv; kwargs...)
    msg=ChangeMessageVisibilityType()
    for p in kwargs
        setfield!(msg, p[1], p[2])
    end
    ChangeMessageVisibility(env, msg)
end
export ChangeMessageVisibility

function ChangeMessageVisibilityBatch(env::AWSEnv, msg::ChangeMessageVisibilityBatchType)
    sqsresp::SQSResponse = call_sqs(env, "ChangeMessageVisibilityBatch" , msg)
    if  (sqsresp.pd != nothing) && (sqsresp.obj == nothing)
        sqsresp.obj = ChangeMessageVisibilityResponseBatchType(sqsresp.pd)
    end
    sqsresp
end
function ChangeMessageVisibilityBatch(env::AWSEnv; kwargs...)
    msg=ChangeMessageVisibilityBatchType()
    for p in kwargs
        setfield!(msg, p[1], p[2])
    end
    ChangeMessageVisibilityBatch(env, msg)
end
export ChangeMessageVisibilityBatch

function ListDeadLetterSourceQueues(env::AWSEnv, msg::ListDeadLetterSourceQueuesType)
    sqsresp::SQSResponse = call_sqs(env, "ListDeadLetterSourceQueues" , msg)
    if  (sqsresp.pd != nothing) && (sqsresp.obj == nothing)
        sqsresp.obj = ListDeadLetterSourceQueuesResponseType(sqsresp.pd)
    end
    sqsresp
end
function ListDeadLetterSourceQueues(env::AWSEnv; kwargs...)
    msg=ListDeadLetterSourceQueuesType()
    for p in kwargs
        setfield!(msg, p[1], p[2])
    end
    ListDeadLetterSourceQueues(env, msg)
end
export ListDeadLetterSourceQueues

function PurgeQueue(env::AWSEnv, msg::PurgeQueueType)
    sqsresp::SQSResponse = call_sqs(env, "PurgeQueue" , msg)
    if  (sqsresp.pd != nothing) && (sqsresp.obj == nothing)
        sqsresp.obj = PurgeQueueResponseType(sqsresp.pd)
    end
    sqsresp
end
function PurgeQueue(env::AWSEnv; kwargs...)
    msg=PurgeQueueType()
    for p in kwargs
        setfield!(msg, p[1], p[2])
    end
    PurgeQueue(env, msg)
end
export PurgeQueue

function AddPermission(env::AWSEnv, msg::AddPermissionType)
    sqsresp::SQSResponse = call_sqs(env, "AddPermission" , msg)
    if  (sqsresp.pd != nothing) && (sqsresp.obj == nothing)
        sqsresp.obj = AddPermissionResponseType(sqsresp.pd)
    end
    sqsresp
end
function AddPermission(env::AWSEnv; kwargs...)
    msg=AddPermissionType()
    for p in kwargs
        setfield!(msg, p[1], p[2])
    end
    AddPermission(env, msg)
end
export AddPermission

function RemovePermission(env::AWSEnv, msg::RemovePermissionType)
    sqsresp::SQSResponse = call_sqs(env, "RemovePermission" , msg)
    if  (sqsresp.pd != nothing) && (sqsresp.obj == nothing)
        sqsresp.obj = RemovePermissionResponseType(sqsresp.pd)
    end
    sqsresp
end
function RemovePermission(env::AWSEnv; kwargs...)
    msg=RemovePermissionType()
    for p in kwargs
        setfield!(msg, p[1], p[2])
    end
    RemovePermission(env, msg)
end
export RemovePermission

ValidRqstMsgs = Dict(
    "CreateQueueType"=>true,
    "GetQueueUrlType"=>true,
    "ListQueuesType"=>true,
    "DeleteQueueType"=>true,
    "DeleteMessageType"=>true,
    "GetQueueAttributesType"=>true,
    "ReceiveMessageType"=>true,
    "SendMessageType"=>true,
    "SetQueueAttributesType"=>true,
    "SendMessageBatchType"=>true,
    "DeleteMessageBatchType"=>true,
    "ChangeMessageVisibilityType"=>true,
    "ChangeMessageVisibilityBatchType"=>true,
    "ListDeadLetterSourceQueuesType"=>true,
    "PurgeQueueType"=>true,
    "AddPermissionType"=>true,
    "RemovePermissionType"=>true
)


text:
using Statistics, StatsBase
using PyCall
xr = pyimport("xarray")
using Elefridge

path = "/Users/milan/cams"
filelist = filter(x->endswith(x,"no2.grib"),readdir(path))
gribfile = xr.open_dataset(joinpath(path,filelist[1]),engine="cfgrib")
no2 = gribfile.no2.values
level = 1:size(no2)[1]

# Float32 information entropy
bitmin = Int(reinterpret(UInt32,minimum(no2)))
bitmax = Int(reinterpret(UInt32,maximum(no2)))
diff = bitmax-bitmin
B = Array{UInt16,1}(undef,diff+1)

for x in no2
    B[reinterpret(UInt32,x)-bitmin+1] += 0x1
end

# throw away all zero entries (don't contribute to entropy)
B = B[B .!= 0x0000]
p = B / sum(B)

H_F32 = entropy(p,2)

# Lin24 information entropy
L = LinQuant24Array(no2)
B = zeros(UInt16,2^24)

for x in L
    B[x+1] += 0x1
end

B = B[B .!= 0x0000]
p = B / sum(B)

H_Lin24 = entropy(p,2)

# Log16 information entropy
L = LogQuant16Array(no2)
B_log = zeros(UInt16,2^16)

for x in L
    B_log[x+1] += 0x1
end

p = B_log / sum(B_log)
H_Log16 = entropy(p,2)

## bitpattern histogram

fig,(ax1,ax2) = subplots(2,1)

ax1.plot(B)
ax2.plot(B_log)

ax1.set_ylabel("N")
ax2.set_ylabel("N")
ax2.set_xlabel("bit pattern")

ax1.set_xticks([1,2^23,2^24])
ax1.set_xticklabels(["0x000000","0x800000","0xffffff"])

ax2.set_xticks([1,2^15,2^16])
ax2.set_xticklabels(["0x0000","0x8000","0xffff"])

ax1.set_title("Bitpattern histogram: 24-bit lin quantisation",loc="left")
ax2.set_title("Bitpattern histogram: 16-bit log quantisation",loc="left")
ax1.set_title("a",loc="right",fontweight="bold")
ax2.set_title("b",loc="right",fontweight="bold")

ax1.text(0.75,0.9,"Entropy = 16.1bit",transform=ax1.transAxes)
ax2.text(0.75,0.9,"Entropy = 15.0bit",transform=ax2.transAxes)

tight_layout()


text:
"""
    BCType

Abstract supertype for boundary condition types.
"""
abstract type BCType end

"""
    Periodic

A type specifying a periodic boundary condition.

A condition may not be specified with a `Periodic` boundary condition.
"""
struct Periodic <: BCType end

"""
    Flux

A type specifying a boundary condition on the flux of a field.
"""
struct Flux <: BCType end

"""
    Gradient

A type specifying a boundary condition on the derivative or gradient of a field. Also
called a Neumann boundary condition.
"""
struct Gradient <: BCType end

"""
    Value

A type specifying a boundary condition on the value of a field. Also called a Dirchlet
boundary condition.
"""
struct Value <: BCType end

"""
    NormalFlow

A type specifying the component of a velocity field normal to a boundary.

Thus `NormalFlow` can only be applied to `u` along x, `v` along y, or `w` along z.
For all other cases --- fields located at (Cell, Cell, Cell), or `u`, `v`,
and `w` in (y, z), (x, z), and (x, y), respectively, either `Value`,
`Gradient`, or `Flux` conditions must be used.

Note that `NormalFlow` differs from a zero `Value` boundary condition:
`Value` imposes values at cell centers, while `NormalFlow` imposes values
_on_ a boundary, at cell faces. Only wall-normal components of the velocity field are defined
on cell faces with respect to the wall-normal direction, and therefore only wall-normal
components of the velocity field are defined on boundaries.
Both tracers and wall-tangential components of velocity fields
are defined at cell centers with respect to the wall-normal direction.
"""
struct NormalFlow <: BCType end


text:
module FeatureCreators

export AbstractFeatureCreator, create_features, feature_size

"""
    AbstractFeatureCreator
An abstract feature creator, for feature transformation from the states.
"""
abstract type AbstractFeatureCreator end

"""
    create_features(fc::AbstractFeatureCreator, s)
Actually create the features
"""
function create_features(fc::AbstractFeatureCreator, s; kwargs...)
    throw("Implement create features for $(typof(fc))")
end

"""
    feature_size(fc::AbstractFeatureCreator)
Get size of feature vector the features assume exists.
"""
function feature_size(fc::AbstractFeatureCreator)
    throw("Implement feature size for $(typof(fc))")
end

export TileCoder
include("FeatureCreators/TileCoder.jl")

export KenervaCoder
include("FeatureCreators/KenervaCoder.jl")

end

max_stars_repo_name: mkschleg/JuliaRL

max_issues_repo_name: mkschleg/JuliaRL

max_forks_repo_name: mkschleg/JuliaRL


text:
"""
        wpm(decisionMat, weights, fns)

Apply WPM (Weighted Product Method) for a given matrix and weights.

# Arguments:
 - `decisionMat::DataFrame`: n × m matrix of objective values for n alterntives and m criteria
 - `weights::Array{Float64, 1}`: m-vector of weights that sum up to 1.0. If the sum of weights is not 1.0, it is automatically normalized.
 - `fns::Array{Function, 1}`: m-vector of functions to be applied on the columns.

# Description
wpm() applies the WPM method to rank n alterntives subject to m criteria which are supposed to be
either maximized or minimized.

# Output
- `::WPMResult`: WPMResult object that holds multiple outputs including scores, rankings, and best index.

# Examples
```julia-repl
julia> decmat = [3        12.5        2        120        14        3;
       5        15        3        110        38        4;
       3        13        2        120        19        3;
       4        14        2        100        31        4;
       3        15        1.5        125        40        4]
5×6 Array{Float64,2}:
 3.0  12.5  2.0  120.0  14.0  3.0
 5.0  15.0  3.0  110.0  38.0  4.0
 3.0  13.0  2.0  120.0  19.0  3.0
 4.0  14.0  2.0  100.0  31.0  4.0
 3.0  15.0  1.5  125.0  40.0  4.0

julia> df = makeDecisionMatrix(decmat)
5×6 DataFrame
 Row │ Crt1     Crt2     Crt3     Crt4     Crt5     Crt6
     │ Float64  Float64  Float64  Float64  Float64  Float64
─────┼──────────────────────────────────────────────────────
   1 │     3.0     12.5      2.0    120.0     14.0      3.0
   2 │     5.0     15.0      3.0    110.0     38.0      4.0
   3 │     3.0     13.0      2.0    120.0     19.0      3.0
   4 │     4.0     14.0      2.0    100.0     31.0      4.0
   5 │     3.0     15.0      1.5    125.0     40.0      4.0

julia> weights = [0.221, 0.159, 0.175, 0.127, 0.117, 0.201];

julia> fns = [maximum, minimum, minimum, maximum, minimum, maximum];

julia> result = wpm(df, weights, fns);

julia> result.scores
5-element Array{Float64,1}:
 0.7975224331331252
 0.7532541470584717
 0.7647463553356331
 0.7873956894790834
 0.7674278741781709

julia> result.bestIndex
1
```

# References
Zavadskas, E. K., Turskis, Z., Antucheviciene, J., & Zakarevicius, A. (2012). Optimization of Weighted Aggregated Sum Product Assessment. Elektronika Ir Elektrotechnika, 122(6), 3-6. https://doi.org/10.5755/j01.eee.122.6.1810
"""
function wpm(decisionMat::DataFrame, weights::Array{Float64,1}, fns::Array{Function,1})::WPMResult

    row, col = size(decisionMat)
    normalizedDecisionMat = similar(decisionMat)
    w = unitize(weights)
    colminmax = zeros(Float64, col)
    @inbounds for i in 1:col
        colminmax[i] = decisionMat[:, i] |> fns[i]
        if fns[i] == maximum
            normalizedDecisionMat[:, i] = decisionMat[:, i] ./ colminmax[i]
        elseif fns[i] == minimum
            normalizedDecisionMat[:, i] = colminmax[i] ./ decisionMat[:, i]
        end
    end
    scoreMat = similar(normalizedDecisionMat)
    for i in 1:col
        scoreMat[:, i] = normalizedDecisionMat[:, i].^w[i]
    end

    scores = zeros(Float64, row)
    for i in 1:row
        scores[i] = prod(scoreMat[i, :])
    end
    rankings = sortperm(scores)

    bestIndex = rankings |> last

    result = WPMResult(
        decisionMat,
        normalizedDecisionMat,
        w,
        scores,
        rankings,
        bestIndex
    )

    return result
end

"""
        wpm(setting)

Apply WPM (Weighted Product Method) for a given matrix and weights.

# Arguments:
 - `setting::MCDMSetting`: MCDMSetting object.

# Description
wpm() applies the WPM method to rank n alterntives subject to m criteria which are supposed to be
either maximized or minimized.

# Output
- `::WPMResult`: WPMResult object that holds multiple outputs including scores, rankings, and best index.
"""
function wpm(setting::MCDMSetting)::WPMResult
    wpm(
        setting.df,
        setting.weights,
        setting.fns
    )
end

max_stars_repo_name: drvinceknight/JMcDM

max_issues_repo_name: drvinceknight/JMcDM

max_forks_repo_name: drvinceknight/JMcDM


text:
# ---
# title: 1352. Product of the Last K Numbers
# id: problem1352
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Array, Design
# link: <https://leetcode.com/problems/product-of-the-last-k-numbers/description/>
# hidden: true
# ---
#
# Implement the class `ProductOfNumbers` that supports two methods:
#
# 1.` add(int num)`
#
#   * Adds the number `num` to the back of the current list of numbers.
#
# 2.` getProduct(int k)`
#
#   * Returns the product of the last `k` numbers in the current list.
#   * You can assume that always the current list has **at least** `k` numbers.
#
# At any time, the product of any contiguous sequence of numbers will fit into a
# single 32-bit integer without overflowing.
#
#
#
# **Example:**
#
#
#
#     **Input**
#     ["ProductOfNumbers","add","add","add","add","add","getProduct","getProduct","getProduct","add","getProduct"]
#     [[],[3],[0],[2],[5],[4],[2],[3],[4],[8],[2]]
#
#     **Output**
#     [null,null,null,null,null,null,20,40,0,null,32]
#
#     **Explanation**
#     ProductOfNumbers productOfNumbers = new ProductOfNumbers();
#     productOfNumbers.add(3);        // [3]
#     productOfNumbers.add(0);        // [3,0]
#     productOfNumbers.add(2);        // [3,0,2]
#     productOfNumbers.add(5);        // [3,0,2,5]
#     productOfNumbers.add(4);        // [3,0,2,5,4]
#     productOfNumbers.getProduct(2); // return 20. The product of the last 2 numbers is 5 * 4 = 20
#     productOfNumbers.getProduct(3); // return 40. The product of the last 3 numbers is 2 * 5 * 4 = 40
#     productOfNumbers.getProduct(4); // return 0. The product of the last 4 numbers is 0 * 2 * 5 * 4 = 0
#     productOfNumbers.add(8);        // [3,0,2,5,4,8]
#     productOfNumbers.getProduct(2); // return 32. The product of the last 2 numbers is 4 * 8 = 32
#
#
#
#
# **Constraints:**
#
#   * There will be at most `40000` operations considering both `add` and `getProduct`.
#   * `0 <= num <= 100`
#   * `1 <= k <= 40000`
#
#
## @lc code=start
using LeetCode

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


text:
module Dash
import HTTP, JSON2, CodecZlib, MD5
using Sockets
const ROOT_PATH = realpath(joinpath(@__DIR__, ".."))
include("exceptions.jl")
include("Components.jl")
include("Front.jl")
include("HttpHelpers/HttpHelpers.jl")
import .Front
using .Components
using .HttpHelpers

export dash, Component, Front, callback!,
enable_dev_tools!, ClientsideFunction,
run_server, PreventUpdate, no_update, @var_str,
Input, Output, State, make_handler, callback_context,
ALL, MATCH, ALLSMALLER

include("Contexts/Contexts.jl")
include("env.jl")
include("utils.jl")
include("app.jl")
include("resources/registry.jl")
include("resources/application.jl")
include("handlers.jl")
include("server.jl")

@doc """
    module Dash

Julia backend for [Plotly Dash](https://github.com/plotly/dash)

# Examples
```julia

using Dash
using DashHtmlComponents
using DashCoreComponents
app = dash(external_stylesheets=["https://codepen.io/chriddyp/pen/bWLwgP.css"]) do
    html_div() do
        dcc_input(id="graphTitle", value="Let's Dance!", type = "text"),
        html_div(id="outputID"),
        dcc_graph(id="graph",
            figure = (
                data = [(x = [1,2,3], y = [3,2,8], type="bar")],
                layout = Dict(:title => "Graph")
            )
        )

    end
end
callback!(app, Output("outputID", "children"), Input("graphTitle","value"), State("graphTitle","type")) do value, type
    "You've entered: '\$(value)' into a '\$(type)' input control"
end
callback!(app, Output("graph", "figure"), Input("graphTitle", "value")) do value
    (
        data = [
            (x = [1,2,3], y = abs.(randn(3)), type="bar"),
            (x = [1,2,3], y = abs.(randn(3)), type="scatter", mode = "lines+markers", line = (width = 4,))
        ],
        layout = (title = value,)
    )
end
run_server(app, HTTP.Sockets.localhost, 8050)
```

""" Dash

end # module


text:
using Combinatorics, Base.Iterators, Knet

gods = (:t,:f,:r)
responses = (:a,:b)
universes = [ (g,a) for g in permutations(gods), a in permutations(responses) ] |> vec # gods and yesno
questions = powerset(permutations(gods) |> collect) |> collect  # is the universe one of these?
questions = product(1:3, questions) |> collect |> vec           # who is asking the question?

function answers(universe, question)
    gods,yesno = universe
    yes,no = yesno
    askto,univset = question
    correct = (gods in univset)
    ans = (gods[askto] == :t ? (correct ? (yes,) : (no,)) :
           gods[askto] == :f ? (correct ? (no,) : (yes,)) :
           responses)
end

strategies = product(questions, questions, questions) |> collect |> vec

function search()
    for s in progress(strategies)
        good(s) && return s
    end
end

function good(s; verbose=false)
    for u in universes
        for q in s
            a = answers(u, q)
        end
    end
end

#     z = map(universes) do u
#         map(s) do q
#             a = answers(u, q)
#         a1 = answers(u,g1,q1)
#         a2 = answers(u,g2,q2)
#         a3 = answers(u,g3,q3)
#         a4 = answers(u,g4,q4)
#         a5 = answers(u,g5,q5)
#         possible_observations = product(a1,a2,a3,a4,a5) |> collect |> vec |> sort
#         (u, possible_observations)
#     end
#     uniq = true
#     for (p1,p2) in combinations(a,2)
#         u1,o1 = p1
#         u2,o2 = p2
#         gods1,yesno1 = u1
#         gods2,yesno2 = u2
#         if gods1 == gods2
#             continue        # no need to distinguish
#         end
#         for x in o1
#             if x in o2
#                 uniq=false  # same observation in two universes
#                 if verbose
#                     println(u1)
#                     println(u2)
#                 end
#                 break
#             end
#         end
#         uniq || break
#     end
#     return uniq
# end

#     for g1 in 1:1,
#         q1 in questions,

#         g2 in 1:2,
#         q2 in questions,
#         g3 in (g1 == g2 ? (1,2) : (1,2,3)),
#         q3 in questions,
#         g4 in (g1 == g2 ? (1,2) : (1,2,3)),
#         q4 in questions

#         for u in universes
#             a1 = answers(u,g1,q1)
#             a2 = answers(u,g2,q2)
#             a12 = sort(vec(collect(product(a1,a2))))
#         end
#     end
# end

#         a1 = [ answers(u,g1,q1) for u in universes ]
#         for g2 in 1:2, q2 in questions
#             a2 = [ answers(u,g2,q2) for u in universes ]
#             a12 = [ sort(vec(collect(product(a1[i],a2[i])))) for i in 1:U ]

#             for g3 in (g1 == g2 ? (1,2) : (1,2,3)), q3 in questions # question if two answers equal
#                 for g4 in (g1 == g2 ? (1,2) : (1,2,3)), q4 in questions # question if two answers are different
#                     a3 = map(a12) do h12
#                         h123 = []
#                         for h in h12
#                             if h[1][end] == h[2][end]
#                                 a3 = [ answers(u,g3,q3)
#                             end
#                         end
#                     end
#                 end
#             end
#         end
#     end
# end

#             a3 = [ answers(u,g3,q3) for u in universes ]

#                         aa = [ sort(vec(collect(product(a1[i],a2[i],a3[i])))) for i in 1:U ]

#                         uniq = true
#                         for i in 1:U-1, j in i+1:U, x in aa[i]
#                             if x in aa[j]
#                                 uniq = false
#                                 break
#                             end
#                         end

#                         if uniq
#                             println(Any[ (g1,q1), (g2,q2), (g3,q3) ])
#                             println(aa)
#                             return
#                         end
#                     end
#                 end
#             end
#         end
#     end
# end

#=
strategies = permutations(questions, 3) |> collect

function responses(s)
    [ (answer(u,1,s[1]), answer(u,2,s[2]), answer(u,3,s[3])) for u in universes ]
end

function goodstrategy(s)
    a = responses(s)
    length(unique(a)) == 12
end

#questions = [ (a,b,c) for a in 1:3, b in 1:3, c in gods ] |> vec # ask `a` if `b` is `c`
#push!(questions, (1,1,:x), (2,2,:x), (3,3,:x)) # ask `a` question whose answer true

# wlog we ask the first question to god1: 9 choices for question
# wlog we ask the second question to god1 or god2 regardless of the answer to the first question: 2x9 choices
# choice of targets: 111, 112, 121, 122, 123
# !!! third question and its target may depend on whether first two answers were same or different !!!
# represent this as a separate fourth question
# using g3/q3 if first two answers were equal, g4/q4 if different

# strategies = product(1:3, 1:3, 1:3, 1:3, 1:3,
#                      questions, questions, questions, questions, questions) |> collect |> vec

# strategies = filter(strategies) do s
#     s[1] <= 1 && s[2] <= 1+s[1] && s[3] <= 1+s[2] && s[4] <= 1+maximum(s[1:3]) && s[5] <= 1+maximum(s[1:4])
# end
# enumerating strategies:

=#

max_stars_repo_name: denizyuret/research

max_issues_repo_name: denizyuret/research

max_forks_repo_name: denizyuret/research


text:
"""
    traverse!(node::Node, f::Function[, args...], <keyword arguments>)
    traverse(node::Node, f::Function[, args...], <keyword arguments>)

Traverse the nodes of a (sub-)tree, given any starting node in the tree, and apply a function
which is either mutating (use `traverse!`) or not (use `traverse`).

# Arguments

- `node::Node`: An MTG node (*e.g.* the whole mtg returned by `read_mtg()`).
- `f::Function`: a function to apply over each node
- `args::Any`: any argument to pass to the function
- <keyword arguments>:

    - `scale = nothing`: The scale to filter-in (i.e. to keep). Usually a Tuple-alike of integers.
    - `symbol = nothing`: The symbol to filter-in. Usually a Tuple-alike of Strings.
    - `link = nothing`: The link with the previous node to filter-in. Usually a Tuple-alike of Char.
    - `filter_fun = nothing`: Any filtering function taking a node as input, e.g. [`isleaf`](@ref).

# Returns

Nothing for `traverse!` because it mutates the (sub-)tree in-place, or an Array of whatever
the function returns for `traverse`.

# Examples

```julia
file = joinpath(dirname(dirname(pathof(MultiScaleTreeGraph))),"test","files","simple_plant.mtg")
mtg = read_mtg(file)
traverse!(mtg, x -> isleaf(x) ? println(x.name," is a leaf") : nothing)
node_5 is a leaf
node_7 is a leaf

# We can also use the `do...end` block notation when we have a complex set of instructions:
traverse!(mtg) do x
    if isleaf(x)
         println(x.name," is a leaf")
    end
end
```
"""
function traverse!(
        node::Node,
        f::Function,
        args...;
        scale = nothing,
        symbol = nothing,
        link = nothing,
        filter_fun = nothing
    )

    if is_filtered(node, scale, symbol, link, filter_fun)
        if !isempty(args)
            f(node, args...)
        else
            f(node)
        end
    end

    if !isleaf(node)
        for chnode in ordered_children(node)
            traverse!(
                chnode,
                f,
                args...;
                scale = scale, symbol = symbol, link = link, filter_fun = filter_fun
            )
        end
    end
end

function traverse!(node::Node, f::Function;
        scale = nothing,
        symbol = nothing,
        link = nothing,
        filter_fun = nothing)

    if is_filtered(node, scale, symbol, link, filter_fun)
        f(node)
    end

    if !isleaf(node)
        for chnode in ordered_children(node)
            traverse!(
                chnode,
                f;
                scale = scale, symbol = symbol, link = link, filter_fun = filter_fun
            )
        end
    end
end

# Non-mutating version:
# Set-up array of value and call the workhorse (traverse_)
function traverse(
        node::Node,
        f::Function,
        args...;
        scale = nothing,
        symbol = nothing,
        link = nothing,
        filter_fun = nothing
    )

    val = []
    # NB: f has to return someting here, if its a mutating function, use traverse!
    traverse_(
        node,
        f,
        val,
        args...;
        scale = scale, symbol = symbol, link = link, filter_fun = filter_fun)
    return val
end

# Actual workhorse:
function traverse_(
        node::Node,
        f::Function,
        val,
        args...;
        scale, symbol, link, filter_fun
    )

    if is_filtered(node, scale, symbol, link, filter_fun)
        push!(val, f(node, args...))
    end

    if !isleaf(node)
        for chnode in ordered_children(node)
            traverse_(
                chnode,
                f,
                val,
                args...;
                scale = scale, symbol = symbol, link = link, filter_fun = filter_fun
            )
        end
    end
end

# Same but with no arguments passed to the function:
function traverse(
        node::Node,
        f::Function;
        scale = nothing,
        symbol = nothing,
        link = nothing,
        filter_fun = nothing
    )
    val = []
    # NB: f has to return someting here, if its a mutating function, use traverse!
    traverse_(node, f, val, scale = scale, symbol = symbol, link = link, filter_fun = filter_fun)
    return val
end

function traverse_(
        node::Node,
        f::Function,
        val;
        scale = nothing,
        symbol = nothing,
        link = nothing,
        filter_fun = nothing
    )

    if is_filtered(node, scale, symbol, link, filter_fun)
        push!(val, f(node))
    end

    if !isleaf(node)
        for chnode in ordered_children(node)
            traverse_(chnode, f, val, scale = scale, symbol = symbol, link = link, filter_fun = filter_fun)
        end
    end
end

# Used for the do...end block notation
function traverse!(
        f::Function,
        node::Node,
        args...;
        scale = nothing,
        symbol = nothing,
        link = nothing,
        filter_fun = nothing
    )
    traverse!(node, f, args...; scale = scale, symbol = symbol, link = link, filter_fun = filter_fun)
end
# Same here but without arguments
function traverse!(
        f::Function,
        node::Node;
        scale = nothing,
        symbol = nothing,
        link = nothing,
        filter_fun = nothing
    )
    traverse!(node, f; scale = scale, symbol = symbol, link = link, filter_fun = filter_fun)
end

# And with the non-mutating version:
function traverse(
        f::Function,
        node::Node,
        args...;
        scale = nothing,
        symbol = nothing,
        link = nothing,
        filter_fun = nothing
    )
    traverse(node, f, args...; scale = scale, symbol = symbol, link = link, filter_fun = filter_fun)
end
# Same here but without arguments
function traverse(
        f::Function,
        node::Node;
        scale = nothing,
        symbol = nothing,
        link = nothing,
        filter_fun = nothing
    )
    traverse(node, f; scale = scale, symbol = symbol, link = link, filter_fun = filter_fun)
end


text:
using Pkg
Pkg.activate(pwd())

# # Julia package
# ## Introduction

using DifferentialEquations

f(u,p,t) = 0.98*u

u0 = 1.0
tspan = (0.0, 1.0)

#+

prob = ODEProblem(f, u0, tspan)
sol = solve(prob)

#+

using Plots

plot(sol; label="")

#+

sol(0.8)

# ### Exercise:
#
# When calling the `solve` function, we can specify the interpolation way. Solve the ODE
# with linear interpolation (`dense=false`) and the Runge-Kutta method of the fourth order
# (`RK4()`). Plot the results and compare them with the default and original solutions.
#
# ---
# ### Solution:

sol2 = solve(prob, dense=false)
sol3 = solve(prob, RK4())

#+

ts = range(tspan...; length = 100)

plot(ts, t->exp(0.98*t), label="True solution", legend=:topleft);
plot!(ts, t->sol(t), label="Default");
plot!(ts, t->sol2(t), label="Linear");
plot!(ts, t->sol3(t), label="Runge-Kutta")

# ---
#
# ## Lorenz system

function lorenz(u, p, t)
    σ, ρ, β = p
    x_t = σ*(u[2]-u[1])
    y_t = u[1]*(ρ-u[3]) - u[2]
    z_t = u[1]*u[2] - β*u[3]
    return [x_t; y_t; z_t]
end

#+

u0 = [1.0; 0.0; 0.0]
p = [10; 28; 8/3]
tspan = (0.0, 100.0)

#+

prob = ODEProblem(lorenz, u0, tspan, p)
sol = solve(prob)

#+

plot(sol)

#+

plt1 = plot(sol, vars=(1,2,3), label="")

#+

plot(sol, vars=(1,2,3), denseplot=false; label="")

#+

traj = hcat(sol.u...)
plot(traj[1,:], traj[2,:], traj[3,:]; label="")

# ### Exercise:
#
# Use the `nextfloat` function to perturb the first parameter of `p` by the smallest
# possible value. Then solve the Lorenz system again and compare the results by plotting
# the two trajectories next to each other.
#
# ---
# ### Solution:

p0 = (nextfloat(p[1]), p[2:end]...)

#+

prob0 = ODEProblem(lorenz, u0, tspan, p0)
sol0 = solve(prob0)

#+

plt0 = plot(sol0, vars=(1,2,3), label="")
plot(plt1, plt0; layout=(1,2))

# ---

hcat(sol(tspan[2]), sol0(tspan[2]))


text:
export TRMA57_Sham_vs_Nwt_λ

function TRMA57_Sham_vs_Nwt_λ(nlp 		:: AbstractNLPModel,
              		 		  nlpstop 	:: NLPStopping;
							  corr_ho :: Symbol = :Shamanskii_MA57,
							  nwt_res_fact = 0.25,
							  λfact = 1.0,
					 		  kwargs...
               		 		  )

	return TRARC(nlp, nlpstop; TR = TrustRegion(10.0),
				  c = Combi(hessian_sparse, PDataMA57{eltype(nlp.meta.x0)}, (x, y, z) -> solve_modelTRDiag_HO_vs_Nwt_λ(x, y, z; ho_correction = corr_ho, nwt_res_fact = nwt_res_fact, λfact = λfact), preprocessMA57, decreaseFact, Tparam{eltype(nlp.meta.x0)}()),
				  kwargs...
				  )
end


text:
#FileIO Interface

fileio_save(f, doc::AbstractDict) = bson(f.filename, doc)
fileio_save(f, args::Pair...) = bson(f.filename, Dict(args))

#This syntax is already in use to work with |> in FileIO
#fileio_save(f; kws...) = bson(f.filename, Dict(kws))

#function fileio_save(f, args...)
#    @assert length(args) % 2 ==0 "Mismatch between labels and data fields"
#    d = Dict(args[1:2:end] .=> args[2:2:end])
#    bson(f.filename, d)
#end

fileio_load(f) = load(f.filename)

function fileio_load(f, args...)
    data = load(f.filename)
    length(args) > 1 && return map( arg -> data[arg], args)
    return data[args[1]]
end


text:
@testset "$TEST $G" begin

f,fio = mktemp()

p1 = readgraph(joinpath(testdir,"testdata","tutte.gml"), G)

@test_throws ErrorException writegraph("file1.ciao", G())
@test_throws ErrorException readgraph("file.ciao", G)
if G <: ANetwork
    @test_throws ErrorException writenetwork("file2.ciao", G())
    @test_throws ErrorException readnetwork("file.ciao", G)
end
@test_throws ErrorException Erdos.NI() # not implemented fallback

# test :gml
p = readgraph(joinpath(testdir,"testdata","tutte.gml"), G)
g = graph(:tutte, G)
@test p == g

gml1 = readgraph(joinpath(testdir,"testdata", "twographs-10-28.gml"), G)
@test nv(gml1) == 10
@test ne(gml1) == 28
@test writegraph(f, gml1, :gml) == 1
gml1a = readgraph(f, :gml, G)
@test gml1a == gml1

gml1 = readgraph(joinpath(testdir,"testdata", "twounnamedgraphs.gml"), G)
@test typeof(gml1) == G
@test nv(gml1) == 4
@test ne(gml1) == 6
@test writegraph(f, gml1, :gml) == 1
gml1a = readgraph(f, :gml, G)
@test gml1a == gml1

g = G(10,0)
@test writegraph(f, g, :gml) == 1
ga = readgraph(f, :gml, G)
@test g == ga

g = G(10,20)
@test writegraph(f, g, :gml) == 1
ga = readgraph(f, :gml, G)
@test g == ga

g = DG(10,20)
@test writegraph(f, g, :gml) == 1
ga = readgraph(f, :gml, G)
@test g == ga

# test :dot
# g = readgraph(joinpath(testdir, "testdata", "twographs.dot"), G)
# @test g == CompleteGraph(6, G)
@test writegraph(f, g, :dot) == 1
# ga = readgraph(f, :dot, G)
# @test g == ga

g = G(10,0)
@test writegraph(f, g, :dot) == 1
# ga = readgraph(f, :dot, G)
# @test g == ga

g = G(10,20)
@test writegraph(f, g, :dot) == 1
# ga = readgraph(f, :dot, G)
# @test g == ga

g = DG(10,20)
@test writegraph(f, g, :dot) == 1
# ga = readgraph(f, :dot, G)
# @test g == ga

#test :net
g10 = CompleteGraph(10, G)
@test typeof(g10) == G

@test writegraph(f, g10, :net) == 1
@test readgraph(f,:net, G) == g10
@test_throws ErrorException readgraph(f, G)

g10 = PathDiGraph(10, DG)
@test typeof(g10) == DG
@test writegraph(f, g10, :net) == 1
@test readgraph(f,:net, G) == g10

g10 = PathDiGraph(10, DG)
@test writegraph(f, g10, :net) == 1
@test readgraph(f,:net, G) == g10

@test writegraph(f, g10, :net, compress=true) == 1
@test readgraph(f,:net, G, compressed=true) == g10

fname = joinpath(f * ".net.gz")
@test writegraph(fname, g10) == 1
@test readgraph(fname, G) == g10
rm(fname)

g10 = readgraph(joinpath(testdir, "testdata", "kinship.net"), G)
@test nv(g10) == 6
@test ne(g10) == 8

g = G(10,0)
@test writegraph(f, g, :net) == 1
ga = readgraph(f, :net, G)
@test g == ga

g = G(10,20)
@test writegraph(f, g, :net) == 1
ga = readgraph(f, :net, G)
@test g == ga

g = DG(10,20)
@test writegraph(f, g, :net) == 1
ga = readgraph(f, :net, G)
@test g == ga

isfile(f) && rm(f)
 end #testset


text:
using Test

include("lasagna.jl")

@testset "preparation time" begin
    @test preptime(2) == 4
    @test preptime(3) == 6
    @test preptime(8) == 16
end

@testset "remaining time" begin
    @test remaining_time(30) == 30
    @test remaining_time(50) == 10
    @test remaining_time(60) == 0
end

@testset "total working time" begin
    @test total_working_time(3, 20) == 26
end


text:
function solve_poisson_pcg( Lmat::SparseMatrixCSC{Float64,Int64},
                            prec::SparseMatrixCSC{Float64,Int64},
                            rho::Array{Float64,1}, NiterMax::Int64;
                            verbose=false, TOL=5.e-10 )
  #
  Npoints = size(rho)[1]
  phi = zeros( Float64, Npoints ) # XXX or use some starting guess
  #
  r = zeros( Float64, Npoints )
  p = zeros( Float64, Npoints )
  z = zeros( Float64, Npoints )
  #
  nabla2_phi = Lmat*phi
  r = rho - nabla2_phi
  z = apply_prec_ilu0( prec, r)
  p = copy(z)

  rsold = dot( r, z )

  for iter = 1 : NiterMax
    #
    nabla2_phi = Lmat*p
    #
    alpha = rsold/dot( p, nabla2_phi )
    #
    phi = phi + alpha * p
    r = r - alpha * nabla2_phi
    z = apply_prec_ilu0( prec, r)
    #
    rsnew = dot( z, r )
    deltars = rsold - rsnew
    #println(typeof(rsnew))
    #println(rsnew)
    if verbose
      @printf("%8d %18.10e\n", iter, sqrt(abs(rsnew)))
    end
    #
    if sqrt(abs(rsnew)) < TOL
      if verbose
        @printf("Convergence achieved in solve_poison_pcg: N, iter: %d %d\n", Npoints, iter)
      end
      break
    end
    #
    p = z + (rsnew/rsold) * p
    #
    rsold = rsnew
  end
  #
  return phi
  #
end # of function