text:
# This file shows an example use of Pyramids.jl to implement an algorithm that requires Complex Steerable Pyramids.
#
# This is an implementation of "Phase-Based Frame Interpolation for Video," by Meyer et. al., from CVPR 2015. [1]
#
# [1] https://www.disneyresearch.com/publication/phasebased/

using Images
using ColorTypes
using Pyramids
using Interpolations

type PhasePyramid <: PyramidType
end

function shift_correction(phi_delta::ImagePyramid; shift_limit=0.5)
    corrected_phase_delta = ImagePyramid(phi_delta)

    for level = corrected_phase_delta.num_levels:-1:1
        for orientation = 1:corrected_phase_delta.num_orientations
            phi_l = subband(corrected_phase_delta, level, orientation=orientation)
            dims_l = size(phi_l)
            phi_limit = (π / (phi_delta.scale^(phi_delta.num_levels - level))) * shift_limit

            if level < phi_delta.num_levels
                phi_lp1 = subband(corrected_phase_delta, level+1, orientation=orientation)
                dims_lp1 = size(phi_lp1)
                phi_lp1_itp = interpolate(phi_lp1, BSpline(Linear()), OnGrid())

                for r = 1:dims_l[1]
                    for c = 1:dims_l[2]

                        r_lp1 = ((r-1)/(dims_l[1]-1)) * (dims_lp1[1]-1) + 1
                        c_lp1 = ((c-1)/(dims_l[2]-1)) * (dims_lp1[2]-1) + 1

                        # unwrap
                        while phi_l[r,c] > (phi_lp1_itp[r_lp1,c_lp1]/phi_delta.scale + π)
                            phi_l[r,c] -= 2*π
                        end

                        while phi_l[r,c] < (phi_lp1_itp[r_lp1,c_lp1]/phi_delta.scale - π)
                            phi_l[r,c] += 2*π
                        end

                        confidence = atan2(sin(phi_l[r,c] - phi_lp1_itp[r_lp1,c_lp1]/phi_delta.scale), cos(phi_l[r,c] - phi_lp1_itp[r_lp1,c_lp1]/phi_delta.scale))

                        if (abs(confidence) > π/2)
                            phi_l[r,c] = phi_lp1_itp[r_lp1,c_lp1]/phi_delta.scale
                        end

                        if abs(phi_l[r,c]) > phi_limit
                            phi_l[r,c] = phi_lp1_itp[r_lp1,c_lp1]/phi_delta.scale
                        end

                    end
                end
            else
                phi_l[abs.(phi_l) .> phi_limit] = 0
            end

            update_subband!(corrected_phase_delta, level, phi_l, orientation=orientation)
        end
    end

    return corrected_phase_delta
end

function adjust_phase(phase_delta::ImagePyramid, corrected_phase_delta::ImagePyramid)
    adjusted_phase = ImagePyramid(phase_delta)

    update_subband!(adjusted_phase, 0, adjust_phase(subband(phase_delta, 0), subband(corrected_phase_delta, 0)))

    for l = 1:pyramid1.num_levels
        for o = 1:pyramid1.num_orientations
            update_subband!(adjusted_phase, l, adjust_phase(subband(phase_delta, l, orientation=o), subband(corrected_phase_delta, l, orientation=o)), orientation=o)
        end
    end

    return adjusted_phase
end

function adjust_phase(phase_delta::Array, corrected_phase_delta)
    adjusted_phase_delta = copy(phase_delta)

    for n = 1:length(adjusted_phase_delta)
        while adjusted_phase_delta[n] > (corrected_phase_delta[n] + π)
            adjusted_phase_delta[n] -= 2*π
        end

        while adjusted_phase_delta[n] < (corrected_phase_delta[n] - π)
            adjusted_phase_delta[n] += 2*π
        end
    end

    return adjusted_phase_delta
end

function blend_and_interpolate(pyramid1::ImagePyramid, pyramid2::ImagePyramid, phase_delta::ImagePyramid, alpha)
    blended_pyramid = ImagePyramid(pyramid1)

    for l = 1:pyramid1.num_levels
        for o = 1:pyramid1.num_orientations
            new_band = ((1 - alpha) * abs.(subband(pyramid1, l, orientation=o)) + alpha * abs.(subband(pyramid2, l, orientation=o))) .* exp(complex(0,1) * (angle.(subband(pyramid1, l, orientation=o)) + alpha * subband(phase_delta, l, orientation=o)))
            update_subband!(blended_pyramid, l, new_band, orientation=o)
        end
    end

    # blend low frequency residuals
    update_subband!(blended_pyramid, pyramid1.num_levels+1, subband(pyramid1, pyramid1.num_levels+1) * (1 - alpha) + alpha * subband(pyramid2, pyramid2.num_levels+1))

    # choose one high frequency residual
    if (alpha > 0.5)
        update_subband!(blended_pyramid, 0, subband(pyramid2, 0))
    else
        update_subband!(blended_pyramid, 0, subband(pyramid1, 0))
    end

    return blended_pyramid

end

function phase_difference(pyramid1::ImagePyramid, pyramid2::ImagePyramid)
    phase_bands = Dict{Int, Union{Array,Dict{Int, Array}}}()

    phase_bands[0] = angle.(subband(pyramid2, 0)) - angle.(subband(pyramid1, 0))

    for l = 1:pyramid1.num_levels
        this_band = Dict{Int, Array}()

        for o = 1:pyramid1.num_orientations
            this_band[o] = angle.(subband(pyramid2, l, orientation=o)) - angle.(subband(pyramid1, l, orientation=o))
        end

        phase_bands[l] = copy(this_band)
    end

    return ImagePyramid(phase_bands, pyramid1.scale, PhasePyramid(), pyramid1.num_levels, pyramid1.num_orientations)
end

println("Loading images")

# convert to LAB, extract the L channel as a 1D array and process it
# update documentation

im1 = load("frame_0.jpg")
im1 = channelview(Lab.(im1))

im2 = load("frame_1.jpg")
im2 = channelview(Lab.(im2))

L1 = im1[1,:,:]
L2 = im2[1,:,:]

println("Converting images to complex steerable pyramids")

pyramid1 = ImagePyramid(L1, ComplexSteerablePyramid(), scale=0.5^0.25)
pyramid2 = ImagePyramid(L2, ComplexSteerablePyramid(), scale=0.5^0.25)

phase_delta = phase_difference(pyramid1, pyramid2)
corrected_phase_delta = shift_correction(phase_delta, shift_limit=1)
adjusted_phase_delta = adjust_phase(phase_delta, corrected_phase_delta)

for alpha = 0:0.2:1.0
    println("Generating image with alpha $(alpha)")
    newpyr = blend_and_interpolate(pyramid1, pyramid2, adjusted_phase_delta, alpha)
    newim = toimage(newpyr)

    newLabIm = zeros(im1)
    newLabIm[1,:,:] = newim
    newLabIm[2,:,:] = (1 - alpha) * im1[2,:,:] + alpha * im2[2,:,:]
    newLabIm[3,:,:] = (1 - alpha) * im1[3,:,:] + alpha * im2[3,:,:]

    newLabIm = colorview(Lab, newLabIm)
    rgbIm = RGB.(newLabIm)
    save("interpolated_frame_$(alpha).png", rgbIm)
end


text:
@testset "Evaluation Metrics" begin
    @testset "Rouge" begin
        candidate_sentence = "Brazil, Russia, China and India are growing nations"
        candidate_summary =  "Brazil, Russia, China and India are growing nations. They are all an important part of BRIC as well as regular part of G20 summits."

        reference_sentences = ["Brazil, Russia, India and China are growing nations", "Brazil and India are two of the developing nations that are part of the BRIC"]
        reference_summaries = ["Brazil, Russia, India and China are the next big poltical powers in the global economy. Together referred to as BRIC(S) along with South Korea.", "Brazil, Russia, India and China are together known as the  BRIC(S) and have been invited to the G20 summit."]

        @test rouge_n(reference_summaries, candidate_summary, 1, avg=true) >= 0.505
        @test rouge_n(reference_summaries, candidate_summary, 2, avg=true) >= 0.131

        @test rouge_n(reference_sentences, candidate_sentence, 2, avg=true) >= 0.349
        @test rouge_n(reference_sentences, candidate_sentence, 1, avg=true) >= 0.666

        @test rouge_l_summary(reference_summaries, candidate_summary, 8, true) >= 0.4256
    end
end

ext: jl


text:
dir = "../examples"
files = [
         "afn.jl",
        ]
for file in files
    @testset "$file" begin include("$dir/$file") end
end


text:
using Flux
println("neural_adapter_tests")
using DiffEqFlux
println("Starting Soon!")
using ModelingToolkit
using DiffEqBase
using Test, NeuralPDE
println("Starting Soon!")
using GalacticOptim
using Optim
using Quadrature,Cubature, Cuba
using QuasiMonteCarlo
using SciMLBase
import ModelingToolkit: Interval, infimum, supremum
using Statistics

using Random
Random.seed!(100)

cb = function (p,l)
    println("Current loss is: $l")
    return false
end

## Example, 2D Poisson equation with Neural adapter
println("Example, 2D Poisson equation with Neural adapter")
@parameters x y
@variables u(..)
Dxx = Differential(x)^2
Dyy = Differential(y)^2

# 2D PDE
eq  = Dxx(u(x,y)) + Dyy(u(x,y)) ~ -sin(pi*x)*sin(pi*y)

# Initial and boundary conditions
bcs = [u(0,y) ~ 0.0, u(1,y) ~ -sin(pi*1)*sin(pi*y),
       u(x,0) ~ 0.0, u(x,1) ~ -sin(pi*x)*sin(pi*1)]
# Space and time domains
domains = [x ∈ Interval(0.0,1.0),
           y ∈ Interval(0.0,1.0)]
quadrature_strategy = NeuralPDE.QuadratureTraining(reltol=1e-2,abstol=1e-2,
                                                   maxiters =50, batch=100)
inner = 8
af = Flux.tanh
chain1 = Chain(Dense(2,inner,af),
               Dense(inner,inner,af),
               Dense(inner,1))
initθ = Float64.(DiffEqFlux.initial_params(chain1))
discretization = NeuralPDE.PhysicsInformedNN(chain1,
                                             quadrature_strategy;
                                             init_params = initθ)

@named pde_system = PDESystem(eq,bcs,domains,[x,y],[u(x,y)])
prob = NeuralPDE.discretize(pde_system,discretization)
sym_prob = NeuralPDE.symbolic_discretize(pde_system,discretization)
res = GalacticOptim.solve(prob, BFGS();  maxiters=2000)
phi = discretization.phi

inner_ = 10
af = Flux.tanh
chain2 = FastChain(FastDense(2,inner_,af),
                   FastDense(inner_,inner_,af),
                   FastDense(inner_,inner_,af),
                   FastDense(inner_,1))

initθ2 =Float64.(DiffEqFlux.initial_params(chain2))

function loss(cord,θ)
    chain2(cord,θ) .- phi(cord,res.minimizer)
end

grid_strategy = NeuralPDE.GridTraining(0.05)
quadrature_strategy = NeuralPDE.QuadratureTraining(reltol=1e-5,abstol=1e-5,
                                                   maxiters =50, batch=100)
stochastic_strategy = NeuralPDE.StochasticTraining(400)
quasirandom_strategy = NeuralPDE.QuasiRandomTraining(400,resampling =false,minibatch = 200)
quasirandom_strategy_resampling = NeuralPDE.QuasiRandomTraining(250)

strategies1 = [grid_strategy,quadrature_strategy]

reses_1 = map(strategies1) do strategy_
    println("Neural adapter Poisson equation, strategy: $(nameof(typeof(strategy_)))")
    prob_ = NeuralPDE.neural_adapter(loss,initθ2,pde_system, strategy_)
    res_ = GalacticOptim.solve(prob_, ADAM(0.01); maxiters=8000)
    prob_ = remake(prob_,u0=res_.minimizer)
    res_ = GalacticOptim.solve(prob_, BFGS(); maxiters=200)
end
strategies2 = [stochastic_strategy,quasirandom_strategy]# quasirandom_strategy_resampling]
reses_2 = map(strategies2) do strategy_
    println("Neural adapter Poisson equation, strategy: $(nameof(typeof(strategy_)))")
    prob_ = NeuralPDE.neural_adapter(loss,initθ2,pde_system, strategy_)
    res_ = GalacticOptim.solve(prob_, ADAM(0.01); maxiters=8000)
    prob_ = remake(prob_,u0=res_.minimizer)
    res_ = GalacticOptim.solve(prob_, BFGS(); maxiters=200)
end
reses_ = [reses_1;reses_2;]

discretizations = map(res_ -> NeuralPDE.PhysicsInformedNN(chain2,
                                                          grid_strategy;
                                                          init_params = res_.minimizer),reses_)

probs = map(discret -> NeuralPDE.discretize(pde_system,discret),discretizations)
phis = map(discret -> discret.phi ,discretizations)

xs,ys = [infimum(d.domain):0.01:supremum(d.domain) for d in domains]
analytic_sol_func(x,y) = (sin(pi*x)*sin(pi*y))/(2pi^2)

u_predict = reshape([first(phi([x,y],res.minimizer)) for x in xs for y in ys],(length(xs),length(ys)))

u_predicts = map(zip(phis,reses_)) do (phi_,res_)
     reshape([first(phi_([x,y],res_.minimizer)) for x in xs for y in ys],(length(xs),length(ys)))
end

u_real = reshape([analytic_sol_func(x,y) for x in xs for y in ys], (length(xs),length(ys)))

@test u_predict ≈ u_real rtol = 0.1
@test u_predicts[1] ≈ u_real rtol = 0.2
map(u_predicts[2:end]) do upred
    @test upred ≈ u_real rtol = 0.3
end

#using Plots
# i=3
# diff_u = abs.(u_predict .- u_real)
# diff_u_ = abs.(u_predicts[i] .- u_real)
# p2 = plot(xs, ys, u_predict, linetype=:contourf,title = "predict");
# p1 = plot(xs, ys, u_real, linetype=:contourf,title = "analytic");
# p5 = plot(xs, ys, diff_u,linetype=:contourf,title = "error");
# p3 = plot(xs, ys, u_predicts[i],linetype=:contourf,title = "predict_");
# p6 = plot(xs, ys, diff_u_,linetype=:contourf,title = "error_");
# plot(p2,p1,p5,p3,p6)

## Example, 2D Poisson equation, domain decomposition
println("Example, 2D Poisson equation, domain decomposition")
@parameters x y
@variables u(..)
Dxx = Differential(x)^2
Dyy = Differential(y)^2

eq  = Dxx(u(x,y)) + Dyy(u(x,y)) ~ -sin(pi*x)*sin(pi*y)

bcs = [u(0,y) ~ 0.0, u(1,y) ~ -sin(pi*1)*sin(pi*y),
       u(x,0) ~ 0.0, u(x,1) ~ -sin(pi*x)*sin(pi*1)]

# Space
x_0 = 0.0
x_end = 1.0
x_domain = Interval(x_0, x_end)
y_domain = Interval(0.0, 1.0)
domains = [x ∈ x_domain,
           y ∈ y_domain]

count_decomp = 10

# Neural network
af = Flux.tanh
inner = 12
chains = [FastChain(FastDense(2, inner, af), FastDense(inner, inner, af), FastDense(inner, 1)) for _ in 1:count_decomp]
initθs = map(c -> Float64.(c), DiffEqFlux.initial_params.(chains))

xs_ = infimum(x_domain):1/count_decomp:supremum(x_domain)
xs_domain = [(xs_[i], xs_[i+1]) for i in 1:length(xs_)-1]
domains_map = map(xs_domain) do (xs_dom)
    x_domain_ = Interval(xs_dom...)
    domains_ = [x ∈ x_domain_,
                y ∈ y_domain]
end

analytic_sol_func(x,y) = (sin(pi*x)*sin(pi*y))/(2pi^2)
function create_bcs(x_domain_,phi_bound)
    x_0, x_e =  x_domain_.left, x_domain_.right
    if x_0 == 0.0
        bcs = [u(0,y) ~ 0.0,
               u(x_e,y) ~ analytic_sol_func(x_e,y),
               u(x,0) ~ 0.0,
               u(x,1) ~ -sin(pi*x)*sin(pi*1)]
        return bcs
    end
    bcs = [u(x_0,y) ~ phi_bound(x_0,y),
           u(x_e,y) ~ analytic_sol_func(x_e,y),
           u(x,0) ~ 0.0,
           u(x,1) ~ -sin(pi*x)*sin(pi*1)]
    bcs
end

reses = []
phis = []
pde_system_map = []

for i in 1:count_decomp
    println("decomposition $i")
    domains_ = domains_map[i]
    phi_in(cord) = phis[i-1](cord,reses[i-1].minimizer)
    # phi_bound(x,y) = if (x isa Matrix)  phi_in(vcat(x, fill(y,size(x)))) else  phi_in(vcat(fill(x,size(y)),y)) end
    phi_bound(x,y) = phi_in(vcat(x,y))
    @register_symbolic phi_bound(x,y)
    global phi_bound
    Base.Broadcast.broadcasted(::typeof(phi_bound), x,y) = phi_bound(x,y)
    bcs_ = create_bcs(domains_[1].domain, phi_bound)
    @named pde_system_ = PDESystem(eq, bcs_, domains_, [x, y], [u(x,y)])
    push!(pde_system_map,pde_system_)
    strategy = NeuralPDE.GridTraining([0.1/count_decomp, 0.1])

    discretization = NeuralPDE.PhysicsInformedNN(chains[i], strategy; init_params=initθs[i])

    prob = NeuralPDE.discretize(pde_system_,discretization)
    symprob = NeuralPDE.symbolic_discretize(pde_system_,discretization)
    res_ = GalacticOptim.solve(prob, BFGS(), maxiters=1500)
    @show res_.minimum
    phi = discretization.phi
    push!(reses, res_)
    push!(phis, phi)
end

# function plot_(i)
#     xs, ys = [infimum(d.domain):dx:supremum(d.domain) for (dx,d) in zip([0.001,0.01], domains_map[i])]
#     u_predict = reshape([first(phis[i]([x,y],reses[i].minimizer)) for x in xs for y in ys],(length(xs),length(ys)))
#     u_real = reshape([analytic_sol_func(x,y) for x in xs for y in ys], (length(xs),length(ys)))
#     diff_u = abs.(u_predict .- u_real)
#     p1 = plot(xs, ys, u_real, linetype=:contourf,title = "analytic");
#     p2 = plot(xs, ys, u_predict, linetype=:contourf,title = "predict");
#     p3 = plot(xs, ys, diff_u,linetype=:contourf,title = "error");
#     plot(p1,p2,p3)
# end
# ps =[plot_(i) for i in 1:count_decomp]

function compose_result(dx)
    u_predict_array = Float64[]
    diff_u_array = Float64[]
    ys = infimum(domains[2].domain):dx:supremum(domains[2].domain)
    xs_ = infimum(x_domain):dx:supremum(x_domain)
    xs = collect(xs_)
    function index_of_interval(x_)
        for (i,x_domain) in enumerate(xs_domain)
            if x_<= x_domain[2] && x_>= x_domain[1]
                return i
            end
        end
    end
    for x_ in xs
        i = index_of_interval(x_)
        u_predict_sub = [first(phis[i]([x_,y],reses[i].minimizer)) for y in ys]
        u_real_sub = [analytic_sol_func(x_,y)  for y in ys]
        diff_u_sub = u_predict_sub .- u_real_sub
        append!(u_predict_array,u_predict_sub)
        append!(diff_u_array,diff_u_sub)
    end
    xs,ys = [infimum(d.domain):dx:supremum(d.domain) for d in domains]
    u_predict = reshape(u_predict_array,(length(xs),length(ys)))
    diff_u = reshape(diff_u_array, (length(xs),length(ys)))
    u_predict, diff_u
end
dx= 0.01
u_predict, diff_u = compose_result(dx)

inner_ = 18
af = Flux.tanh
chain2 = FastChain(FastDense(2,inner_,af),
                   FastDense(inner_,inner_,af),
                   FastDense(inner_,inner_,af),
                   FastDense(inner_,inner_,af),
                   FastDense(inner_,1))

initθ2 =Float64.(DiffEqFlux.initial_params(chain2))

@named pde_system = PDESystem(eq, bcs, domains, [x, y], [u(x,y)])

losses = map(1:count_decomp) do i
    loss(cord,θ) = chain2(cord,θ) .- phis[i](cord,reses[i].minimizer)
end

prob_ = NeuralPDE.neural_adapter(losses,initθ2, pde_system_map,NeuralPDE.GridTraining([0.1/count_decomp,0.1]))
res_ = GalacticOptim.solve(prob_, BFGS(); maxiters=2000)
@show res_.minimum
prob_ = NeuralPDE.neural_adapter(losses,res_.minimizer, pde_system_map, NeuralPDE.GridTraining(0.01))
res_ = GalacticOptim.solve(prob_, BFGS(); maxiters=1000)
@show res_.minimum

parameterless_type_θ = DiffEqBase.parameterless_type(initθ2)
phi_ = NeuralPDE.get_phi(chain2,parameterless_type_θ)

xs,ys = [infimum(d.domain):dx:supremum(d.domain) for d in domains]
u_predict_ = reshape([first(phi_([x,y],res_.minimizer)) for x in xs for y in ys],(length(xs),length(ys)))
u_real = reshape([analytic_sol_func(x,y) for x in xs for y in ys], (length(xs),length(ys)))
diff_u_ = u_predict_ .- u_real

@test u_predict ≈ u_real rtol = 0.1
@test u_predict_ ≈ u_real rtol = 0.1

# p1 = plot(xs, ys, u_predict, linetype=:contourf,title = "predict 1");
# p2 = plot(xs, ys, u_predict_,linetype=:contourf,title = "predict 2");
# p3 = plot(xs, ys, u_real, linetype=:contourf,title = "analytic");
# p4 = plot(xs, ys, diff_u,linetype=:contourf,title = "error 1");
# p5 = plot(xs, ys, diff_u_,linetype=:contourf,title = "error 2");
# plot(p1,p2,p3,p4,p5)


text:
module Examples

using ..Garlic

include("aob.jl")
include("rcp.jl")

end


text:
function filter_stdlib_list(stdlib_list::AbstractVector{StdlibInfo}, stdlibs_to_include::AbstractString)::Vector{StdlibInfo}
    _stdlibs_to_include_string = strip(stdlibs_to_include)
    if ( isempty(_stdlibs_to_include_string) ) || ( lowercase(_stdlibs_to_include_string) == "all" )
        return stdlib_list
    else
        _stdlibs_to_include_vector = convert(Vector{String}, strip.(split(_stdlibs_to_include_string, ",")))::Vector{String}
        return filter_stdlib_list(stdlib_list, _stdlibs_to_include_vector)
    end
end

function filter_stdlib_list(stdlib_list::AbstractVector{<:StdlibInfo}, stdlibs_to_include::AbstractVector{<:AbstractString})::Vector{StdlibInfo}
    filtered_stdlib_list = Vector{StdlibInfo}(undef, 0)
    for stdlib in stdlib_list
        if stdlib.name in stdlibs_to_include
            push!(filtered_stdlib_list, stdlib)
        end
    end
    return filtered_stdlib_list
end


text:
using AutocorrelationTime
using Base.Test

# write your own tests here
@test 0.4 < integrated_autocorrelation_time(rand(100000)) < 0.6 #Random numbers have integrated autocorrelation time 1/2


text:
using FastGroupBy, DataBench
svec = DataBench.gen_string_vec_var_len(1_000_000, 24);
julia_lsd_radixsort_elapsed = @elapsed FastGroupBy.radixsort_lsd16!(svec)
@assert issorted(svec)

svec = DataBench.gen_string_vec_var_len(1_000_000, 24);
julia_lsd_radixsort_elapsed = @elapsed FastGroupBy.radixsort_lsd24!(svec)
@assert issorted(svec)

svec = DataBench.gen_string_vec_var_len(10_000_000, 24);
julia_lsd_radixsort_elapsed = @elapsed FastGroupBy.radixsort_lsd16!(svec)
@assert issorted(svec)

svec = DataBench.gen_string_vec_var_len(10_000_000, 24);
julia_lsd_radixsort_elapsed = @elapsed FastGroupBy.radixsort_lsd24!(svec)
@assert issorted(svec)

svec = DataBench.gen_string_vec_var_len(100_000_000, 24);
julia_lsd_radixsort_elapsed = @elapsed FastGroupBy.radixsort_lsd16!(svec)
@assert issorted(svec)

svec = DataBench.gen_string_vec_var_len(100_000_000, 24);
julia_lsd_radixsort_elapsed = @elapsed FastGroupBy.radixsort_lsd24!(svec)
@assert issorted(svec)

svec = DataBench.gen_string_vec_var_len(1_000_000, 32);
julia_lsd_radixsort_elapsed = @elapsed FastGroupBy.radixsort_lsd16!(svec)
@assert issorted(svec)

svec = DataBench.gen_string_vec_var_len(1_000_000, 32);
julia_lsd_radixsort_elapsed = @elapsed FastGroupBy.radixsort_lsd24!(svec)
@assert issorted(svec)

svec = DataBench.gen_string_vec_var_len(1_000_000, 32);
julia_lsd_radixsort_elapsed = @elapsed FastGroupBy.radixsort_lsd32!(svec)
@assert issorted(svec)

svec = DataBench.gen_string_vec_var_len(10_000_000, 32);
julia_lsd_radixsort_elapsed = @elapsed FastGroupBy.radixsort_lsd16!(svec)
@assert issorted(svec)

svec = DataBench.gen_string_vec_var_len(10_000_000, 32);
julia_lsd_radixsort_elapsed = @elapsed FastGroupBy.radixsort_lsd24!(svec)
@assert issorted(svec)

svec = DataBench.gen_string_vec_var_len(10_000_000, 32);
julia_lsd_radixsort_elapsed = @elapsed FastGroupBy.radixsort_lsd32!(svec)
@assert issorted(svec)

svec = DataBench.gen_string_vec_var_len(100_000_000, 32);
julia_lsd_radixsort_elapsed = @elapsed FastGroupBy.radixsort_lsd16!(svec)
@assert issorted(svec)

svec = DataBench.gen_string_vec_var_len(100_000_000, 32);
julia_lsd_radixsort_elapsed = @elapsed FastGroupBy.radixsort_lsd24!(svec)
@assert issorted(svec)

svec = DataBench.gen_string_vec_var_len(100_000_000, 32);
julia_lsd_radixsort_elapsed = @elapsed FastGroupBy.radixsort_lsd32!(svec)
@assert issorted(svec)

svec = DataBench.gen_string_vec_var_len(10_000_000, 64);
julia_lsd_radixsort_elapsed = @elapsed FastGroupBy.radixsort_lsd16!(svec)
@assert issorted(svec)

svec = DataBench.gen_string_vec_var_len(10_000_000, 64);
julia_lsd_radixsort_elapsed = @elapsed FastGroupBy.radixsort_lsd24!(svec)
@assert issorted(svec)

svec = DataBench.gen_string_vec_var_len(10_000_000, 64);
julia_lsd_radixsort_elapsed = @elapsed FastGroupBy.radixsort_lsd32!(svec)
@assert issorted(svec)

svec = DataBench.gen_string_vec_var_len(100_000_000, 64);
julia_lsd_radixsort_elapsed = @elapsed FastGroupBy.radixsort_lsd16!(svec)
@assert issorted(svec)

svec = DataBench.gen_string_vec_var_len(100_000_000, 64);
julia_lsd_radixsort_elapsed = @elapsed FastGroupBy.radixsort_lsd24!(svec)
@assert issorted(svec)

svec = DataBench.gen_string_vec_var_len(100_000_000, 64);
julia_lsd_radixsort_elapsed = @elapsed FastGroupBy.radixsort_lsd32!(svec)
@assert issorted(svec)


text:
#####
##### Grabbing properties for thermodynamic variables
#####

function thermodynamic_field(model)
    model.thermodynamic_variable isa Energy  && return model.tracers.ρe
    model.thermodynamic_variable isa Entropy && return model.tracers.ρs
end

function intermediate_thermodynamic_field(model)
    model.thermodynamic_variable isa Energy  && return model.time_stepper.intermediate_fields.tracers.ρe
    model.thermodynamic_variable isa Entropy && return model.time_stepper.intermediate_fields.tracers.ρs
end

#####
##### Kernels for diagnosing thermodynamic variables
#####

@kernel function compute_velocities!(velocities, grid, momenta, total_density)
    i, j, k = @index(Global, NTuple)

    @inbounds velocities.u[i, j, k] = momenta.ρu[i, j, k] / ℑxᶠᵃᵃ(i, j, k, grid, total_density)
    @inbounds velocities.v[i, j, k] = momenta.ρv[i, j, k] / ℑyᵃᶠᵃ(i, j, k, grid, total_density)
    @inbounds velocities.w[i, j, k] = momenta.ρw[i, j, k] / ℑzᵃᵃᶠ(i, j, k, grid, total_density)
end

@kernel function compute_p_over_ρ!(p_over_ρ, grid, thermodynamic_variable, gases, gravity, total_density, momenta, tracers)
    i, j, k = @index(Global, NTuple)

    @inbounds p_over_ρ[i, j, k] = diagnose_p_over_ρ(i, j, k, grid, thermodynamic_variable, gases, gravity, total_density, momenta, tracers)
end

@kernel function compute_temperature!(temperature, grid, thermodynamic_variable, gases, gravity, total_density, momenta, tracers)
    i, j, k = @index(Global, NTuple)

    @inbounds temperature[i, j, k] = diagnose_temperature(i, j, k, grid, thermodynamic_variable, gases, gravity, total_density, momenta, tracers)
end

function compute_temperature!(model)
    temperature = intermediate_thermodynamic_field(model)

    temperature, total_density, momenta, tracers =
        datatuples(temperature, model.total_density, model.momenta, model.tracers)

    compute_temperature_event =
        launch!(model.architecture, model.grid, :xyz, compute_temperature!,
                temperature, model.grid, model.thermodynamic_variable, model.gases,
                model.gravity, total_density, momenta, tracers,
                dependencies=Event(device(model.architecture)))

    wait(device(model.architecture), compute_temperature_event)

    return temperature
end

#####
##### CFL
#####

function cfl(model, Δt)
    # We will store the velocities in the time stepper's intermediate fields.
    velocities = (
        u = model.time_stepper.intermediate_fields.ρu,
        v = model.time_stepper.intermediate_fields.ρv,
        w = model.time_stepper.intermediate_fields.ρw
    )

    velocities, momenta, total_density =
        datatuples(velocities, model.momenta, model.total_density)

    velocity_event =
        launch!(model.architecture, model.grid, :xyz, compute_velocities!,
                velocities, model.grid, momenta, total_density,
                dependencies=Event(device(model.architecture)))

    wait(device(model.architecture), velocity_event)

    u_max = maximum(velocities.u.parent)
    v_max = maximum(velocities.v.parent)
    w_max = maximum(velocities.w.parent)

    Δx, Δy, Δz = model.grid.Δx, model.grid.Δy, model.grid.Δz
    CFL = Δt / min(Δx/u_max, Δy/v_max, Δz/w_max)

    return CFL
end

#####
##### Acoustic CFL
#####

function acoustic_cfl(model, Δt)
    # We will store p/ρ in the time stepper's intermediate field for the thermodynamic variable.
    p_over_ρ = intermediate_thermodynamic_field(model)

    p_over_ρ, total_density, momenta, tracers =
        datatuples(p_over_ρ, model.total_density, model.momenta, model.tracers)

    compute_p_over_ρ_event =
        launch!(model.architecture, model.grid, :xyz, compute_p_over_ρ!,
                p_over_ρ, model.grid, model.thermodynamic_variable, model.gases,
                model.gravity, total_density, momenta, tracers,
                dependencies=Event(device(model.architecture)))

    wait(device(model.architecture), compute_p_over_ρ_event)

    p_over_ρ_max = maximum(p_over_ρ.parent)

    γ_max = maximum(gas.cₚ / gas.cᵥ for gas in model.gases)
    c_max = √(γ_max * p_over_ρ_max)

    Δx, Δy, Δz = model.grid.Δx, model.grid.Δy, model.grid.Δz
    acoustic_CFL = Δt / min(Δx/c_max, Δy/c_max, Δz/c_max)

    return acoustic_CFL
end


text:
module ED

using ..Positions
using ..SimLib
using ..SimLib: Maybe, FArray

using Arpack
using Distributed
using FilteredMatrices
import JLD2
using KrylovKit
using LinearAlgebra
using Printf: @sprintf
using Reexport
using SharedArrays
using XXZNumerics
using SpinSymmetry

export EDDataDescriptor, EDDerivedDataDescriptor, EDDerivedData, Full, Sparse, ShiftInvert, ShiftInvertARPACK, POLFED, ed_size, run_ed
export Serial, Threaded, Parallel

include("task.jl")
include("models/models.jl")
include("diagonalizationtype.jl")
include("ed_datadescriptor.jl")
include("ed_derived_descriptor.jl")
include("runmode.jl")
include("ed_core.jl")

@reexport using .EDModels

end #module


text:
using PrettyTables

"""
match_to_id :
Match each row of a sample to a leaf id (key)
"""
function match_to_id(tree,X::Array{Float64,2})
    n_samples, n_features = size(X)
    id_vec = Array{Int64,1}()

    for i in 1:n_samples
        node = tree
        while node.is_leaf == false
            if X[i,node.feature] < node.threshold
                node = node.l
            else
                node = node.r
            end
        end
        push!(id_vec,node.id)
    end
    return(id_vec)
end

"""
recursive_action :
Perform a recursive action (e.g. print) over the nodes of a tree
"""
function recursive_action(tree,action::Function)
    if tree.is_leaf == false
        action(tree)
        recursive_action(tree.l,action)
        recursive_action(tree.r,action)
    else
        action(tree)
    end
end

"""
remove_region :
Remove region index information from a node
"""
function remove_region!(node)
    node.n = size(node.region)[1]
    node.region = nothing
end

"""
tree_to_array :
convert a tree to an array
"""
function tree_to_array(tree,parent=-1,tree_tab=Array{Union{Int64,Float64,Bool},2}(undef,0,8))
    if tree.is_leaf == false
        tree_tab = vcat(Union{Int64,Float64,Bool}[tree.id parent tree.is_leaf tree.value tree.feature tree.threshold tree.depth tree.n],
                    tree_to_array(tree.l,tree.id,tree_tab),
                    tree_to_array(tree.r,tree.id,tree_tab))
    else
        tree_tab = Union{Int64,Float64,Bool}[tree.id parent tree.is_leaf tree.value NaN NaN tree.depth tree.n]
    end
    return(tree_tab)
end

"""
print a tree

"""
function print(
    tree
    ;digits=2)
    tree = tree.tree
    header = ["id", "parent","leaf","value","feature","threshold","depth","n"]
    tree_tab = tree_to_array(tree)
    tree_tab = tree_tab[sortperm(tree_tab[:,7]),:]
    formatter = Dict(4 => (v,i)-> round(v;digits=digits),
                    6 => (v,i)-> round(v;digits=digits))
    pretty_table(tree_tab,header;formatter=formatter)
end

"""
node_id :
Generate a unique ID for each node
"""
function node_id(ids::Array{Int64,1})
    if (length(ids) == 0)
        return(0)
    else
        return(1 + maximum(ids))
    end
end


text:
import StaticArrays: SizedVector
export PauliString

# TODO: expand to Clifford?
struct PauliString{N,BT<:ConstantGate{1},VT<:SizedVector{N,BT}} <: CompositeBlock{N}
    blocks::VT
    function PauliString(blocks::SizedVector{N,BT}) where {N,BT<:ConstantGate{1}}
        @warn "`PauliString` will be moved to `YaoExtensions.jl` in the next release."
        new{N,BT,typeof(blocks)}(blocks)
    end
end

# NOTE: PauliString has a fixed size `N`, thus by default, it should use
#      SizedVector, or this block could be actually not correct.

"""
    PauliString(xs::PauliGate...)

Create a `PauliString` from some Pauli gates.

# Example

```jldoctest; setup=:(using YaoBlocks)
julia> PauliString(X, Y, Z)
┌ Warning: `PauliString` will be moved to `YaoExtensions.jl` in the next release.
└ @ YaoBlocks ~/.julia/dev/YaoBlocks/src/composite/pauli_strings.jl:8
nqubits: 3
PauliString
├─ X gate
├─ Y gate
└─ Z gate
```
"""
PauliString(xs::PauliGate...) = PauliString(SizedVector{length(xs),PauliGate}(xs))

"""
    PauliString(list::Vector)

Create a `PauliString` from a list of Pauli gates.

# Example

```jldoctest; setup=:(using YaoBlocks)
julia> PauliString([X, Y, Z])
┌ Warning: `PauliString` will be moved to `YaoExtensions.jl` in the next release.
└ @ YaoBlocks ~/.julia/dev/YaoBlocks/src/composite/pauli_strings.jl:8
nqubits: 3
PauliString
├─ X gate
├─ Y gate
└─ Z gate
```
"""
function PauliString(xs::Vector)
    for each in xs
        if !(each isa PauliGate)
            error("expect pauli gates")
        end
    end
    return PauliString(SizedVector{length(xs),PauliGate}(xs))
end

subblocks(ps::PauliString) = ps.blocks
chsubblocks(pb::PauliString, blocks::Vector) = PauliString(blocks)
chsubblocks(pb::PauliString, it) = PauliString(collect(it))

occupied_locs(ps::PauliString) = (findall(x -> !(x isa I2Gate), ps.blocks)...,)

cache_key(ps::PauliString) = map(cache_key, ps.blocks)

YaoBase.ishermitian(::PauliString) = true
YaoBase.isreflexive(::PauliString) = true
YaoBase.isunitary(::PauliString) = true

Base.copy(ps::PauliString) = PauliString(copy(ps.blocks))
Base.getindex(ps::PauliString, x) = getindex(ps.blocks, x)
Base.lastindex(ps::PauliString) = lastindex(ps.blocks)
Base.iterate(ps::PauliString) = iterate(ps.blocks)
Base.iterate(ps::PauliString, st) = iterate(ps.blocks, st)
Base.length(ps::PauliString) = length(ps.blocks)
Base.eltype(ps::PauliString) = eltype(ps.blocks)
Base.eachindex(ps::PauliString) = eachindex(ps.blocks)
Base.getindex(ps::PauliString, index::Union{UnitRange,Vector}) =
    PauliString(getindex(ps.blocks, index))
function Base.setindex!(ps::PauliString, v::PauliGate, index::Union{Int})
    ps.blocks[index] = v
    return ps
end

function Base.:(==)(lhs::PauliString{N}, rhs::PauliString{N}) where {N}
    (length(lhs.blocks) == length(rhs.blocks)) && all(lhs.blocks .== rhs.blocks)
end

xgates(ps::PauliString{N}) where {N} =
    RepeatedBlock{N}(X, (findall(x -> x isa XGate, (ps.blocks...,))...,))
ygates(ps::PauliString{N}) where {N} =
    RepeatedBlock{N}(Y, (findall(x -> x isa YGate, (ps.blocks...,))...,))
zgates(ps::PauliString{N}) where {N} =
    RepeatedBlock{N}(Z, (findall(x -> x isa ZGate, (ps.blocks...,))...,))

function apply!(reg::AbstractRegister, ps::PauliString)
    _check_size(reg, ps)
    for pauligates in [xgates, ygates, zgates]
        blk = pauligates(ps)
        apply!(reg, blk)
    end
    return reg
end

function mat(::Type{T}, ps::PauliString) where {T}
    return mat(T, xgates(ps)) * mat(T, ygates(ps)) * mat(T, zgates(ps))
end


text:
function vtuwrite(filename, cgnodes, cgelcon, cgcells, celltype, scalars, vectors, fields)

filename = filename * ".vtu";

# get dimensions
npoints,nd = size(cgnodes);
ncells,nve = size(cgcells);

# Preparing output for CG
if nd==2
    cgnodes = hcat(cgnodes,zeros(npoints,1));
end
outputCG = zeros(npoints,3);

# Get endianness
if ENDIAN_BOM==0x04030201
    byte_order = "LittleEndian";
elseif ENDIAN_BOM==0x01020304
    byte_order = "BigEndian";
else
    error("Endian is not valid");
end

# binary format
formattype = "appended";

# float and integer byte size
fbytesize = 4; # Float32 format
ibytesize = 4; # Int32   format
obytesize = 8; # Int64   format

nsc = length(scalars);
if nsc>1
    if nsc==2
        sidx = [1];
        sidy = [2];
    else
        sidx = collect(1:2:(nsc-1));
        sidy = collect(2:2:nsc);
    end
else
    sidx = [];
    sidy = [];
end
nvt = length(vectors);
if nvt>1
    if nvt==2
        vidx = [1];
        vidy = [2];
    else
        vidx = collect(1:2:(nvt-1));
        vidy = collect(2:2:nvt);
    end
else
    vidx = [];
    vidy = [];
end

# Open VTK output file
fid = open(filename, "w");
# VTK DataFile Version
mystr = "<?xml version=\"1.0\"?>\n";
mystr = mystr * "<VTKFile type=\"UnstructuredGrid\" version=\"1.0\" byte_order=\"" * byte_order * "\" header_type=\"UInt64\">\n";
mystr = mystr * "  <UnstructuredGrid>\n";
mystr = mystr * "    <Piece NumberOfPoints=\"" * string(npoints) * "\"" *  " NumberOfCells=\"" * string(ncells) * "\">\n";

###########################################################################
#                      1 - WRITE VTU METADATA                             #
#    (Write all dataset description for the VTU binary-appended format    #
###########################################################################
offset = 0;
if (length(sidx)>0) || (length(vidx)>0)
    mystr = mystr * "      <PointData Scalars=\"scalars\">\n";
end
for ii = 1:length(sidx)
    title = scalars[sidx[ii]];
    mystr = mystr * "        <DataArray type=\"Float32\" Name=\"" * title * "\" Format=\"" * formattype * "\" offset=\"" * string(offset) * "\">\n";
    mystr = mystr * "        </DataArray>\n";
    offset = offset + npoints*fbytesize + obytesize;
end

for ii = 1:length(vidx)
    title = vectors[vidx[ii]];
    mystr = mystr * "        <DataArray type=\"Float32\" Name=\"" * title * "\" NumberOfComponents=\"3" * "\" Format=\"" * formattype * "\" offset=\"" * string(offset) * "\">\n";
    mystr = mystr * "        </DataArray>\n";
    offset = offset + 3*npoints*fbytesize + obytesize;
end
if (length(sidx)>0) || (length(vidx)>0)
    mystr = mystr * "      </PointData>\n";
end
mystr = mystr * "      <Points>\n";
mystr = mystr * "        <DataArray type=\"Float32\" Name=\"" * "points" * "\" NumberOfComponents=\"3" * "\" Format=\"" * formattype * "\" offset=\"" * string(offset) * "\">\n";
mystr = mystr * "        </DataArray>\n";
mystr = mystr * "      </Points>\n";
offset = offset + 3*npoints*fbytesize + obytesize;
mystr = mystr * "      <Cells>\n";
mystr = mystr * "        <DataArray type=\"Int32\" Name=\"" * "connectivity" * "\" Format=\"" * formattype * "\" offset=\"" * string(offset) * "\">\n";
mystr = mystr * "        </DataArray>\n";
offset = offset + ncells*nve*ibytesize + obytesize;
mystr = mystr * "        <DataArray type=\"Int32\" Name=\"" * "offsets" * "\" Format=\"" * formattype * "\" offset=\"" * string(offset) * "\">\n";
mystr = mystr * "        </DataArray>\n";
offset = offset + ncells*ibytesize + obytesize;
mystr = mystr * "        <DataArray type=\"UInt8\" Name=\"" * "types" * "\" Format=\"" * formattype * "\" offset=\"" * string(offset) * "\">\n";
mystr = mystr * "        </DataArray>\n";
mystr = mystr * "      </Cells>\n";
mystr = mystr * "    </Piece>\n";
mystr = mystr * "  </UnstructuredGrid>\n";
mystr = mystr * "  <AppendedData encoding=\"raw\">\n"
mystr = mystr * "   _";

write(fid,mystr);

###########################################################################
#                   2 - WRITE UNSTRUCTURED DATA                           #
#        (Vectors and Scalars datasets all projected on CG mesh           #
###########################################################################
# This part writes the dataset attributes for the mesh nodes or elements.
# It is customary to write the dataset before the mesh description.

#  Write all the scalar dataset first
for ii = 1:length(sidx)
    output = fields[:,scalars[sidy[ii]],:];#varargin{ sidx(ii) + 2};
    outputCG[cgelcon[:],1] = output[:];
    write(fid, Int64(npoints*fbytesize));
    write(fid, Float32.(outputCG[:,1]));
end

#  Write all the vector datasets then
for ii = 1:length(vidx)
    tmp = fields[:,vectors[vidy[ii]][1],:];
    output = zeros(length(tmp),3);
    output[:,1] = tmp[:];
    tmp = fields[:,vectors[vidy[ii]][2],:];
    output[:,2] = tmp[:];
    if nd==3
        tmp = fields[:,vectors[vidy[ii]][3],:];
        output[:,3] = tmp[:];
    end
    # Convert to CG Field
    outputCG[cgelcon[:],1] = output[:,1];
    outputCG[cgelcon[:],2] = output[:,2];
    outputCG[cgelcon[:],3] = output[:,3];
    # For a more accurate conversion, use the slower
    write(fid, Int64(3*npoints*fbytesize));
    write(fid, Float32.(outputCG'));
end

###########################################################################
#                3 - WRITE UNSTRUCTURED_GRID MESH                         #
#       (Mesh nodes coordinates connectivities, offset ant type)          #
###########################################################################
# Write 3D coordinates of mesh points
output = cgnodes';
write(fid, Int64(3*npoints*fbytesize));
write(fid, Float32.(output));

# Write (sub)-cells connectivities
output = cgcells';
write(fid, Int64(nve*ncells*ibytesize));
write(fid, Int32.(output));

# Write (sub)-cells offsets
output = collect(nve:nve:nve*ncells);
write(fid, Int64(ncells*ibytesize));
write(fid, Int32.(output));

# Write (sub)-cells type
output = celltype*ones(Int,1,ncells);
write(fid, Int64(ncells));
write(fid, UInt8.(output));

write(fid,"  </AppendedData>\n");
write(fid,"</VTKFile>\n");
close(fid);

end

max_stars_repo_name: glwagner/Exasim

max_issues_repo_name: glwagner/Exasim

max_forks_repo_name: glwagner/Exasim


text:
import LinearAlgebra: norm, axpy!

abstract type AbstractCeedVector end

struct CeedVectorActive <: AbstractCeedVector end
Base.getindex(::CeedVectorActive) = C.CEED_VECTOR_ACTIVE[]

struct CeedVectorNone <: AbstractCeedVector end
Base.getindex(::CeedVectorNone) = C.CEED_VECTOR_NONE[]

mutable struct CeedVector <: AbstractCeedVector
    ref::RefValue{C.CeedVector}
    arr::Union{Nothing,AbstractArray}
    CeedVector(ref::Ref{C.CeedVector}) = new(ref, nothing)
end

"""
    CeedVector(c::Ceed, len::Integer)

Creates a `CeedVector` of given length.
"""
function CeedVector(c::Ceed, len::Integer)
    ref = Ref{C.CeedVector}()
    C.CeedVectorCreate(c[], len, ref)
    obj = CeedVector(ref)
    finalizer(obj) do x
        # ccall(:jl_safe_printf, Cvoid, (Cstring, Cstring), "Finalizing %s.\n", repr(x))
        destroy(x)
    end
    return obj
end
destroy(v::CeedVector) = C.CeedVectorDestroy(v.ref) # COV_EXCL_LINE
Base.getindex(v::CeedVector) = v.ref[]

Base.summary(io::IO, v::CeedVector) = print(io, length(v), "-element CeedVector")
function Base.show(io::IO, ::MIME"text/plain", v::CeedVector)
    summary(io, v)
    println(io, ":")
    witharray_read(v, MEM_HOST) do arr
        Base.print_array(io, arr)
    end
end
Base.show(io::IO, v::CeedVector) = witharray_read(a -> show(io, a), v, MEM_HOST)

function Base.length(::Type{T}, v::CeedVector) where {T}
    len = Ref{C.CeedInt}()
    C.CeedVectorGetLength(v[], len)
    return T(len[])
end

Base.ndims(::CeedVector) = 1
Base.ndims(::Type{CeedVector}) = 1
Base.axes(v::CeedVector) = (Base.OneTo(length(v)),)
Base.size(v::CeedVector) = (length(Int, v),)
Base.length(v::CeedVector) = length(Int, v)

"""
    setvalue!(v::CeedVector, val::Real)

Set the [`CeedVector`](@ref) to a constant value.
"""
setvalue!(v::CeedVector, val::Real) = C.CeedVectorSetValue(v[], val)
"""
    setindex!(v::CeedVector, val::Real)
    v[] = val

Set the [`CeedVector`](@ref) to a constant value, synonymous to [`setvalue!`](@ref).
"""
Base.setindex!(v::CeedVector, val::Real) = setvalue!(v, val)

"""
    norm(v::CeedVector, ntype::NormType)

Return the norm of the given [`CeedVector`](@ref).

The norm type can either be specified as one of `NORM_1`, `NORM_2`, `NORM_MAX`.
"""
function norm(v::CeedVector, ntype::NormType)
    nrm = Ref{CeedScalar}()
    C.CeedVectorNorm(v[], ntype, nrm)
    nrm[]
end

"""
    norm(v::CeedVector, p::Real)

Return the norm of the given [`CeedVector`](@ref), see [`norm(::CeedVector,
::NormType)`](@ref).

`p` can have value 1, 2, or Inf, corresponding to `NORM_1`, `NORM_2`, and `NORM_MAX`,
respectively.
"""
function norm(v::CeedVector, p::Real)
    if p == 1
        ntype = NORM_1
    elseif p == 2
        ntype = NORM_2
    elseif isinf(p)
        ntype = NORM_MAX
    else
        error("norm(v::CeedVector, p): p must be 1, 2, or Inf")
    end
    norm(v, ntype)
end

"""
    reciprocal!(v::CeedVector)

Set `v` to be equal to its elementwise reciprocal.
"""
reciprocal!(v::CeedVector) = C.CeedVectorReciprocal(v[])

"""
    setarray!(v::CeedVector, mtype::MemType, cmode::CopyMode, arr)

Set the array used by a [`CeedVector`](@ref), freeing any previously allocated array if
applicable. The backend may copy values to a different [`MemType`](@ref). See also
[`syncarray!`](@ref) and [`takearray!`](@ref).

!!! warning "Avoid OWN_POINTER CopyMode"
    The [`CopyMode`](@ref) `OWN_POINTER` is not suitable for use with arrays that are
    allocated by Julia, since those cannot be properly freed from libCEED.
"""
function setarray!(v::CeedVector, mtype::MemType, cmode::CopyMode, arr)
    C.CeedVectorSetArray(v[], mtype, cmode, arr)
    if cmode == USE_POINTER
        v.arr = arr
    end
end

"""
    syncarray!(v::CeedVector, mtype::MemType)

Sync the [`CeedVector`](@ref) to a specified [`MemType`](@ref). This function is used to
force synchronization of arrays set with [`setarray!`](@ref). If the requested memtype is
already synchronized, this function results in a no-op.
"""
syncarray!(v::CeedVector, mtype::MemType) = C.CeedVectorSyncArray(v[], mtype)

"""
    takearray!(v::CeedVector, mtype::MemType)

Take ownership of the [`CeedVector`](@ref) array and remove the array from the
[`CeedVector`](@ref). The caller is responsible for managing and freeing the array. The
array is returns as a `Ptr{CeedScalar}`.
"""
function takearray!(v::CeedVector, mtype::MemType)
    ptr = Ref{Ptr{CeedScalar}}()
    C.CeedVectorTakeArray(v[], mtype, ptr)
    v.arr = nothing
    ptr[]
end

# Helper function to parse arguments of @witharray and @witharray_read
function witharray_parse(assignment, args)
    if !Meta.isexpr(assignment, :(=))
        error("@witharray must have first argument of the form v_arr=v") # COV_EXCL_LINE
    end
    arr = assignment.args[1]
    v = assignment.args[2]
    mtype = MEM_HOST
    sz = :((length($(esc(v))),))
    body = args[end]
    for i = 1:length(args)-1
        a = args[i]
        if !Meta.isexpr(a, :(=))
            error("Incorrect call to @witharray or @witharray_read") # COV_EXCL_LINE
        end
        if a.args[1] == :mtype
            mtype = a.args[2]
        elseif a.args[1] == :size
            sz = esc(a.args[2])
        end
    end
    arr, v, sz, mtype, body
end

"""
    @witharray(v_arr=v, [size=(dims...)], [mtype=MEM_HOST], body)

Executes `body`, having extracted the contents of the [`CeedVector`](@ref) `v` as an array
with name `v_arr`. If the [`memory type`](@ref MemType) `mtype` is not provided, `MEM_HOST`
will be used. If the size is not specified, a flat vector will be assumed.

# Examples
Negate the contents of `CeedVector` `v`:
```
@witharray v_arr=v v_arr .*= -1.0
```
"""
macro witharray(assignment, args...)
    arr, v, sz, mtype, body = witharray_parse(assignment, args)
    quote
        arr_ref = Ref{Ptr{C.CeedScalar}}()
        C.CeedVectorGetArray($(esc(v))[], $(esc(mtype)), arr_ref)
        try
            $(esc(arr)) = UnsafeArray(arr_ref[], Int.($sz))
            $(esc(body))
        finally
            C.CeedVectorRestoreArray($(esc(v))[], arr_ref)
        end
    end
end

"""
    @witharray_read(v_arr=v, [size=(dims...)], [mtype=MEM_HOST], body)

Same as [`@witharray`](@ref), but provides read-only access to the data.
"""
macro witharray_read(assignment, args...)
    arr, v, sz, mtype, body = witharray_parse(assignment, args)
    quote
        arr_ref = Ref{Ptr{C.CeedScalar}}()
        C.CeedVectorGetArrayRead($(esc(v))[], $(esc(mtype)), arr_ref)
        try
            $(esc(arr)) = UnsafeArray(arr_ref[], Int.($sz))
            $(esc(body))
        finally
            C.CeedVectorRestoreArrayRead($(esc(v))[], arr_ref)
        end
    end
end

"""
    setindex!(v::CeedVector, v2::AbstractArray)
    v[] = v2

Sets the values of [`CeedVector`](@ref) `v` equal to those of `v2` using broadcasting.
"""
Base.setindex!(v::CeedVector, v2::AbstractArray) = @witharray(a = v, a .= v2)

"""
    CeedVector(c::Ceed, v2::AbstractVector; mtype=MEM_HOST, cmode=COPY_VALUES)

Creates a new [`CeedVector`](@ref) using the contents of the given vector `v2`. By default,
the contents of `v2` will be copied to the new [`CeedVector`](@ref), but this behavior can
be changed by specifying a different `cmode`.
"""
function CeedVector(c::Ceed, v2::AbstractVector; mtype=MEM_HOST, cmode=COPY_VALUES)
    v = CeedVector(c, length(v2))
    setarray!(v, mtype, cmode, v2)
    v
end

"""
    Vector(v::CeedVector)

Create a new `Vector` by copying the contents of `v`.
"""
function Base.Vector(v::CeedVector)
    v2 = Vector{CeedScalar}(undef, length(v))
    @witharray_read(a = v, v2 .= a)
end

"""
    witharray(f, v::CeedVector, mtype=MEM_HOST)

Calls `f` with an array containing the data of the `CeedVector` `v`, using [`memory
type`](@ref MemType) `mtype`.

Because of performance issues involving closures, if `f` is a complex operation, it may be
more efficient to use the macro version `@witharray` (cf. the section on "Performance of
captured variable" in the [Julia
documentation](https://docs.julialang.org/en/v1/manual/performance-tips) and related [GitHub
issue](https://github.com/JuliaLang/julia/issues/15276).

# Examples

Return the sum of a vector:
```
witharray(sum, v)
```
"""
function witharray(f, v::CeedVector, mtype::MemType=MEM_HOST)
    arr_ref = Ref{Ptr{C.CeedScalar}}()
    C.CeedVectorGetArray(v[], mtype, arr_ref)
    arr = UnsafeArray(arr_ref[], (length(v),))
    res = try
        f(arr)
    finally
        C.CeedVectorRestoreArray(v[], arr_ref)
    end
    return res
end

"""
    witharray_read(f, v::CeedVector, mtype::MemType=MEM_HOST)

Same as [`witharray`](@ref), but with read-only access to the data.

# Examples

Display the contents of a vector:
```
witharray_read(display, v)
```
"""
function witharray_read(f, v::CeedVector, mtype::MemType=MEM_HOST)
    arr_ref = Ref{Ptr{C.CeedScalar}}()
    C.CeedVectorGetArrayRead(v[], mtype, arr_ref)
    arr = UnsafeArray(arr_ref[], (length(v),))
    res = try
        f(arr)
    finally
        C.CeedVectorRestoreArrayRead(v[], arr_ref)
    end
    return res
end

"""
    scale!(v::CeedVector, a::Real)

Overwrite `v` with `a*v` for scalar `a`. Returns `v`.
"""
function scale!(v::CeedVector, a::Real)
    C.CeedVectorScale(v[], a)
    return v
end

"""
    axpy!(a::Real, x::CeedVector, y::CeedVector)

Overwrite `y` with `x*a + y`, where `a` is a scalar. Returns `y`.

!!! warning "Different argument order"
    In order to be consistent with `LinearAlgebra.axpy!`, the arguments are passed in order: `a`,
    `x`, `y`. This is different than the order of arguments of the C function `CeedVectorAXPY`.
"""
function axpy!(a::Real, x::CeedVector, y::CeedVector)
    C.CeedVectorAXPY(y[], a, x[])
    return y
end

"""
    pointwisemult!(w::CeedVector, x::CeedVector, y::CeedVector)

Overwrite `w` with `x .* y`. Any subset of x, y, and w may be the same vector. Returns `w`.
"""
function pointwisemult!(w::CeedVector, x::CeedVector, y::CeedVector)
    C.CeedVectorPointwiseMult(w[], x[], y[])
    return w
end

max_stars_repo_name: AdelekeBankole/libCEED

max_issues_repo_name: AdelekeBankole/libCEED

max_forks_repo_name: AdelekeBankole/libCEED


text:
using GeometricIntegrators
using ChargedParticleDynamics.GuidingCenter4d.TokamakFastSurface
using GeometricExamples

const Δt     = 5.
const ntime  = 5000
const nx     = 200
const ny     = 200
const nsave  = 10
const nplot  = 5
const run_id = "poincare_2nd_fast_firk_dt5"

include("guiding_center_4d_settings_firk.jl")

tableau_list = get_tableau_list_firk()

pinv = guiding_center_4d_ode_poincare_invariant_2nd(Δt, nx, ny, ntime, nsave)

include("guiding_center_4d_fast_poincare_invariant_2nd.jl")


text:
list_test_str = split(
"2199943210
3987894921
9856789892
8767896789
9899965678",
'\n'
)

list_str = let input = @__DIR__() * "/day9_input.txt"
    isfile(input) ? readlines(input) : list_test_str
end

# Test function

function testit(name, result, check::Function)
    if check(result)
        @info "$name successful, with result $result"
    else
        @warn "$name failed, with result $result"
    end
end

# Part 1

function part1(list)
    m = fill(9, length(list)+2, length(list[1])+2)
    m[2:end-1, 2:end-1] .= parse.(Int, permutedims(reduce(hcat, collect.(list))))
    res = 0
    h, w = size(m)
    for i in 2:h-1
        for j in 2:w-1
            if all(p -> m[i, j] < m[i+first(p), j+last(p)], ((1, 0), (-1, 0), (0, 1),  (0, -1)))
                res += 1 + m[i, j]
            end
        end
    end
    return res
end

testit("Test Part 1", part1(list_test_str), ==(15))
testit("Challenge Part 1", part1(list_str), ==(591))

# Part 2

function part2(list)
    function poke!(c, i, j, id)
        if c[i, j] == 0
            c[i, j] = id
            for (u,v) in ((1, 0), (-1, 0), (0, 1), (0, -1))
                poke!(c, i+u, j+v, id)
            end
        end
    end
    c = fill(-1, length(list)+2, length(list[1])+2)
    c[2:end-1, 2:end-1] .= -(parse.(Int, permutedims(reduce(hcat, collect.(list)))) .== 9)
    h, w = size(c)
    id = 0
    for i in 2:h-1
        for j in 2:w-1
            if c[i, j] == 0
                id += 1
                poke!(c, i, j, id)
            end
        end
    end
    return prod(sort([count(c .== j) for j in 1:maximum(c)], rev=true)[1:3])
end

testit("Test Part 2", part2(list_test_str), ==(1134))
testit("Challenge Part 2", part2(list_str), ==(1113424))


text:
# This file was generated, do not modify it. # hide
features, targets = MLJ.@load_boston
features = DataFrames.DataFrame(features)
@show size(features)
@show targets[1:3]
first(features, 3) |> MLJ.pretty


text:
abstract type AbstractResult end # Common for all result types, e.g., SimResult and FreqRespResult

## SimResult: the output of lsim etc. ==========================================
abstract type AbstractSimResult <: AbstractResult end # Result of a time-domain simulation

struct SimResult{Ty, Tt, Tx, Tu, Ts} <: AbstractSimResult # Result of lsim
    y::Ty
    t::Tt
    x::Tx
    u::Tu
    sys::Ts
end

# To emulate, e.g., lsim(sys, u)[1] -> y
function Base.getindex(r::SimResult, i::Int)
    return getfield(r, i)
end

function Base.getindex(r::SimResult, v::AbstractVector)
    return getfield.((r,), v)
end

# to allow destructuring, e.g., y,t,x = lsim(sys, u)
# This performs explicit iteration in the type domain to ensure inferability
Base.iterate(r::SimResult)              = (r.y, Val(:t))
Base.iterate(r::SimResult, ::Val{:t})   = (r.t, Val(:x))
Base.iterate(r::SimResult, ::Val{:x})   = (r.x, Val(:u))
Base.iterate(r::SimResult, ::Val{:u})   = (r.u, Val(:done))
Base.iterate(r::SimResult, ::Val{:done}) = nothing

function Base.getproperty(r::SimResult, s::Symbol)
    s ∈ fieldnames(SimResult) && (return getfield(r, s))
    if s === :nx
        return size(r.x, 1)
    elseif s === :nu
        return size(r.u, 1)
    elseif s === :ny
        return size(r.y, 1)
    else
        throw(ArgumentError("Unsupported property $s"))
    end
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

struct ConnectionMonitorsApi <: SwaggerApi
    client::Swagger.Client
end

"""

Create or update a connection monitor.
Param: resourceGroupName::String (required)
Param: networkWatcherName::String (required)
Param: connectionMonitorName::String (required)
Param: parameters::ConnectionMonitor (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: ConnectionMonitorResult
"""
function _swaggerinternal_connectionMonitorsCreateOrUpdate(_api::ConnectionMonitorsApi, resourceGroupName::String, networkWatcherName::String, connectionMonitorName::String, parameters, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "PUT", ConnectionMonitorResult, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Network/networkWatchers/{networkWatcherName}/connectionMonitors/{connectionMonitorName}", ["azure_auth"], parameters)
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "networkWatcherName", networkWatcherName)  # type String
    Swagger.set_param(_ctx.path, "connectionMonitorName", connectionMonitorName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function connectionMonitorsCreateOrUpdate(_api::ConnectionMonitorsApi, resourceGroupName::String, networkWatcherName::String, connectionMonitorName::String, parameters, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_connectionMonitorsCreateOrUpdate(_api, resourceGroupName, networkWatcherName, connectionMonitorName, parameters, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function connectionMonitorsCreateOrUpdate(_api::ConnectionMonitorsApi, response_stream::Channel, resourceGroupName::String, networkWatcherName::String, connectionMonitorName::String, parameters, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_connectionMonitorsCreateOrUpdate(_api, resourceGroupName, networkWatcherName, connectionMonitorName, parameters, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""

Deletes the specified connection monitor.
Param: resourceGroupName::String (required)
Param: networkWatcherName::String (required)
Param: connectionMonitorName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: Nothing
"""
function _swaggerinternal_connectionMonitorsDelete(_api::ConnectionMonitorsApi, resourceGroupName::String, networkWatcherName::String, connectionMonitorName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "DELETE", Nothing, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Network/networkWatchers/{networkWatcherName}/connectionMonitors/{connectionMonitorName}", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "networkWatcherName", networkWatcherName)  # type String
    Swagger.set_param(_ctx.path, "connectionMonitorName", connectionMonitorName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function connectionMonitorsDelete(_api::ConnectionMonitorsApi, resourceGroupName::String, networkWatcherName::String, connectionMonitorName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_connectionMonitorsDelete(_api, resourceGroupName, networkWatcherName, connectionMonitorName, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function connectionMonitorsDelete(_api::ConnectionMonitorsApi, response_stream::Channel, resourceGroupName::String, networkWatcherName::String, connectionMonitorName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_connectionMonitorsDelete(_api, resourceGroupName, networkWatcherName, connectionMonitorName, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""

Gets a connection monitor by name.
Param: resourceGroupName::String (required)
Param: networkWatcherName::String (required)
Param: connectionMonitorName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: ConnectionMonitorResult
"""
function _swaggerinternal_connectionMonitorsGet(_api::ConnectionMonitorsApi, resourceGroupName::String, networkWatcherName::String, connectionMonitorName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", ConnectionMonitorResult, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Network/networkWatchers/{networkWatcherName}/connectionMonitors/{connectionMonitorName}", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "networkWatcherName", networkWatcherName)  # type String
    Swagger.set_param(_ctx.path, "connectionMonitorName", connectionMonitorName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function connectionMonitorsGet(_api::ConnectionMonitorsApi, resourceGroupName::String, networkWatcherName::String, connectionMonitorName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_connectionMonitorsGet(_api, resourceGroupName, networkWatcherName, connectionMonitorName, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function connectionMonitorsGet(_api::ConnectionMonitorsApi, response_stream::Channel, resourceGroupName::String, networkWatcherName::String, connectionMonitorName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_connectionMonitorsGet(_api, resourceGroupName, networkWatcherName, connectionMonitorName, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""

Lists all connection monitors for the specified Network Watcher.
Param: resourceGroupName::String (required)
Param: networkWatcherName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: ConnectionMonitorListResult
"""
function _swaggerinternal_connectionMonitorsList(_api::ConnectionMonitorsApi, resourceGroupName::String, networkWatcherName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", ConnectionMonitorListResult, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Network/networkWatchers/{networkWatcherName}/connectionMonitors", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "networkWatcherName", networkWatcherName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function connectionMonitorsList(_api::ConnectionMonitorsApi, resourceGroupName::String, networkWatcherName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_connectionMonitorsList(_api, resourceGroupName, networkWatcherName, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function connectionMonitorsList(_api::ConnectionMonitorsApi, response_stream::Channel, resourceGroupName::String, networkWatcherName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_connectionMonitorsList(_api, resourceGroupName, networkWatcherName, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""

Query a snapshot of the most recent connection states.
Param: resourceGroupName::String (required)
Param: networkWatcherName::String (required)
Param: connectionMonitorName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: ConnectionMonitorQueryResult
"""
function _swaggerinternal_connectionMonitorsQuery(_api::ConnectionMonitorsApi, resourceGroupName::String, networkWatcherName::String, connectionMonitorName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "POST", ConnectionMonitorQueryResult, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Network/networkWatchers/{networkWatcherName}/connectionMonitors/{connectionMonitorName}/query", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "networkWatcherName", networkWatcherName)  # type String
    Swagger.set_param(_ctx.path, "connectionMonitorName", connectionMonitorName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function connectionMonitorsQuery(_api::ConnectionMonitorsApi, resourceGroupName::String, networkWatcherName::String, connectionMonitorName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_connectionMonitorsQuery(_api, resourceGroupName, networkWatcherName, connectionMonitorName, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function connectionMonitorsQuery(_api::ConnectionMonitorsApi, response_stream::Channel, resourceGroupName::String, networkWatcherName::String, connectionMonitorName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_connectionMonitorsQuery(_api, resourceGroupName, networkWatcherName, connectionMonitorName, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""

Starts the specified connection monitor.
Param: resourceGroupName::String (required)
Param: networkWatcherName::String (required)
Param: connectionMonitorName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: Nothing
"""
function _swaggerinternal_connectionMonitorsStart(_api::ConnectionMonitorsApi, resourceGroupName::String, networkWatcherName::String, connectionMonitorName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "POST", Nothing, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Network/networkWatchers/{networkWatcherName}/connectionMonitors/{connectionMonitorName}/start", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "networkWatcherName", networkWatcherName)  # type String
    Swagger.set_param(_ctx.path, "connectionMonitorName", connectionMonitorName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function connectionMonitorsStart(_api::ConnectionMonitorsApi, resourceGroupName::String, networkWatcherName::String, connectionMonitorName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_connectionMonitorsStart(_api, resourceGroupName, networkWatcherName, connectionMonitorName, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function connectionMonitorsStart(_api::ConnectionMonitorsApi, response_stream::Channel, resourceGroupName::String, networkWatcherName::String, connectionMonitorName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_connectionMonitorsStart(_api, resourceGroupName, networkWatcherName, connectionMonitorName, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""

Stops the specified connection monitor.
Param: resourceGroupName::String (required)
Param: networkWatcherName::String (required)
Param: connectionMonitorName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: Nothing
"""
function _swaggerinternal_connectionMonitorsStop(_api::ConnectionMonitorsApi, resourceGroupName::String, networkWatcherName::String, connectionMonitorName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "POST", Nothing, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Network/networkWatchers/{networkWatcherName}/connectionMonitors/{connectionMonitorName}/stop", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "networkWatcherName", networkWatcherName)  # type String
    Swagger.set_param(_ctx.path, "connectionMonitorName", connectionMonitorName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function connectionMonitorsStop(_api::ConnectionMonitorsApi, resourceGroupName::String, networkWatcherName::String, connectionMonitorName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_connectionMonitorsStop(_api, resourceGroupName, networkWatcherName, connectionMonitorName, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function connectionMonitorsStop(_api::ConnectionMonitorsApi, response_stream::Channel, resourceGroupName::String, networkWatcherName::String, connectionMonitorName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_connectionMonitorsStop(_api, resourceGroupName, networkWatcherName, connectionMonitorName, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""

Update tags of the specified connection monitor.
Param: resourceGroupName::String (required)
Param: networkWatcherName::String (required)
Param: connectionMonitorName::String (required)
Param: parameters::TagsObject (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: ConnectionMonitorResult
"""
function _swaggerinternal_connectionMonitorsUpdateTags(_api::ConnectionMonitorsApi, resourceGroupName::String, networkWatcherName::String, connectionMonitorName::String, parameters, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "PATCH", ConnectionMonitorResult, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Network/networkWatchers/{networkWatcherName}/connectionMonitors/{connectionMonitorName}", ["azure_auth"], parameters)
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "networkWatcherName", networkWatcherName)  # type String
    Swagger.set_param(_ctx.path, "connectionMonitorName", connectionMonitorName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function connectionMonitorsUpdateTags(_api::ConnectionMonitorsApi, resourceGroupName::String, networkWatcherName::String, connectionMonitorName::String, parameters, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_connectionMonitorsUpdateTags(_api, resourceGroupName, networkWatcherName, connectionMonitorName, parameters, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function connectionMonitorsUpdateTags(_api::ConnectionMonitorsApi, response_stream::Channel, resourceGroupName::String, networkWatcherName::String, connectionMonitorName::String, parameters, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_connectionMonitorsUpdateTags(_api, resourceGroupName, networkWatcherName, connectionMonitorName, parameters, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

export connectionMonitorsCreateOrUpdate, connectionMonitorsDelete, connectionMonitorsGet, connectionMonitorsList, connectionMonitorsQuery, connectionMonitorsStart, connectionMonitorsStop, connectionMonitorsUpdateTags


text:
# Project on adaptive evolvability as first described in cce/Adaptive_evolvability9_6_21.docx.

# For each of number_sources starting phenotypes choose number_target random phenotypes within
#     distance_limit of the source phenotype.
# For each starting phenotype and each target phenotype do number_runs epochal evolutions with
#     a maximum of step_limit steps and fitness defined by hamming distance from target.
# Record the number of successes and the averge number of steps for successful runs
# max_steps is the maximum number of steps on the evolution to find the starting chromosome for a run
function compare_evolvabilities( p::Parameters, source_list::GoalList, number_targets::Int64,
    distance_limit::Float64, step_limit::Int64, number_runs::Int64;
    csvfile::String="", evo_df_file::String="", evo_df_column::Symbol=:d_evolvability,
    iteration_limit::Int64=200, max_steps::Int64=200000, distance_function::Function=hamming_distance )
  evo_df = DataFrame()
  if length(evo_df_file) > 0
    try
      evo_df = read_dataframe(evo_df_file)
    catch
      error("reading dataframe from ",evo_df_file," failed")
    end
    @assert size(evo_df)[1] == length(source_list)  # Make sure goal lists are compatible
  end
  df = DataFrame()
  df.source = Goal[]
  #df.numinputs = Int64[]
  #df.numoutputs = Int64[]
  #df.numgates= Int64[]
  #df.numlevelsback = Int64[]
  df.numtargets = Int64[]
  df.dist_limit = Float64[]
  df.step_limit = Int64[]
  df.number_runs = Int64[]
  df.avg_num_successes = Float64[]
  df.avg_steps = Float64[]
  println("size(df): ",size(df))
  # list of parameter values that go into the dataframe
  #df_params_list = Any[p.numinputs, p.numoutputs, p.numinteriors, p.numlevelsback, number_targets, distance_limit, step_limit, number_runs]
  #println("df_params_list: ",df_params_list)
  #for src in source_list
  df_rows = pmap( src->run_compare_evolvabilities( p, src, number_targets, distance_limit, step_limit, number_runs,
        iteration_limit=iteration_limit, max_steps=max_steps, distance_function=distance_function ), source_list )
  for df_row in df_rows
    push!(df,df_row)
  end
  df.evolvability = evo_df[:,evo_df_column]
  if length(csvfile) > 0
    open( csvfile, "w" ) do f
      hostname = chomp(open("/etc/hostname") do f read(f,String) end)
      println(f,"# date and time: ",Dates.now())
      println(f,"# host: ",hostname," with ",nprocs()-1,"  processes: " )
      println(f,"# funcs: ", Main.CGP.default_funcs(p.numinputs))
      print_parameters(f,p,comment=true)
      CSV.write( f, df, append=true, writeheader=true )
    end
  end
  df
end

# Helper function for the above function for a single src goal.
function run_compare_evolvabilities( p::Parameters, src::Goal, number_targets::Int64,
  distance_limit::Float64, step_limit::Int64, number_runs::Int64;
  iteration_limit::Int64=200, max_steps::Int64=200000, distance_function::Function=hamming_distance )
  println("src: ",src)
  df_params_list =
    Any[
      #p.numinputs,
      #p.numoutputs,
      #p.numinteriors,
      #p.numlevelsback,
      number_targets,
      distance_limit,
      step_limit,
      number_runs
    ]
  iteration = 0
  target_list = Goal[]
  while iteration < iteration_limit && length(target_list) < number_targets
    target = randgoal( p.numinputs, p.numoutputs )
    #hd = hamming_distance( src, target, p.numinputs )
    hd = distance_function( src, target, p.numinputs )
    #println("target: ",target,"  hd: ",hd)
    if hd <= distance_limit && target != src
      push!(target_list,target)
      println("target: ",target,"  hd: ",hd,"  iteration: ",iteration)
    end
    iteration += 1
  end
  if iteration == iteration_limit
    error("target list iteration exceeded iteration limit of ",iteration_limit," for source src: ",src)
  end
  println("target_list: ",target_list)
  successes = 0
  sum_steps = 0
  for target in target_list
    ch = random_chromosome( p )
    (src_ch,steps) = neutral_evolution(ch, src, max_steps )
    while iteration < iteration_limit && steps == max_steps
      (src_ch,steps) = neutral_evolution(ch, src, max_steps )
      iteration += 1
    end
    if iteration == iteration_limit
      error("iteration to src exceeded iteration limit of ",iteration_limit," for source src: ",src," and target: ",target)
    end
    @assert  steps < max_steps
    for run = 1:number_runs
      (target_ch,steps) = neutral_evolution(src_ch, target, step_limit )
      successes += steps < step_limit ? 1 : 0
      sum_steps += steps < step_limit ? steps : 0
    end
  end
  avg_successes = successes/length(target_list)/number_runs
  avg_steps = sum_steps/successes
  df_row = vcat([src],df_params_list,[avg_successes,avg_steps])
  return df_row
end

# Doesn't make these assignments global
function test_inputs()
  p = Parameters(3,1,7,4)
   #source_list=randgoallist( 10, p.numinputs, p.numoutputs )
   source_list=map(x->[MyInt(x)],collect(0:255))
   number_targets=3
   distance_limit=0.3
   step_limit=200
   number_runs=10
   evo_df_file="../data/9_3_21/robust_evo_by_walks9_3_21B.csv"
end


text:
@testset "JJ terms" begin
    @testset "jj coupling of equivalent electrons" begin
        # Table 4.5, Cowan 1981
        foreach([
            (ro"1s",ro"2p-") => [(0,2) => 0, 1 => 1//2],
            (ro"2p",ro"3d-") => [(0,4) => 0, (1,3) => 3//2, 2 => (0,2)],
            (ro"3d",ro"4f-") => [(0,6) => 0, (1,5) => 5//2, (2,4) => (0,2,4),
                                3 => (3//2,5//2,9//2)],
            (ro"4f",ro"5g-") => [(0,8) => 0, (1,7) => 7//2, (2,6) => (0,2,4,6),
                                (3,5) => (3//2,5//2,7//2,9//2,11/2,15//2),
                                4 => (0,2,2,4,4,5,6,8)],
            (ro"5g",ro"6h-") => [(0,10) => 0, (1,9) => 9//2, (2,8) => (0,2,4,6,8),
                                (3,7) => (3//2,5//2,7//2,9//2,9//2,11//2,13//2,15//2,17//2,21//2),
                                (4,6) => (0,0,2,2,3,4,4,4,5,6,6,6,7,8,8,9,10,12),
                                5 => (1//2,3//2,5//2,5//2,7//2,7//2,9//2,9//2,9//2,11//2,11//2,13//2,13//2,15//2,15//2,17//2,17//2,19//2,21//2,25//2)]
        ]) do (orbs,wsj)
            foreach(orbs) do orb
                foreach(wsj) do (ws,j)
                    js = map(HalfInteger, isa(j, Number) ? [j] : j)
                    foreach(ws) do w
                        # Should actually test without unique
                        @test unique(intermediate_terms(orb,w)) == unique(js)
                    end
                end
            end
        end
    end
end


text:
module SparseMatrixAssemblers

using Test
using Gridap.Arrays
using Gridap.TensorValues
using Gridap.ReferenceFEs
using Gridap.Geometry
using Gridap.Integration
using Gridap.Fields
using Gridap.Algebra
using SparseArrays
using Gridap.FESpaces
using Gridap.CellData

domain =(0,1,0,1)
partition = (2,2)
model = CartesianDiscreteModel(domain,partition)

reffe = ReferenceFE(lagrangian,Float64,1)
V = FESpace(model,reffe,dirichlet_tags=[1,2,3,4,6,5])
U = V

b(x) = x[2]

v = get_cell_shapefuns(V)
u = get_cell_shapefuns_trial(U)

degree = 2
trian = get_triangulation(model)
quad = CellQuadrature(trian,degree)

btrian = BoundaryTriangulation(model)
bquad = CellQuadrature(btrian,degree)

b0trian = Triangulation(btrian,Int[])
b0quad = CellQuadrature(b0trian,degree)

cellmat = integrate(∇(v)⊙∇(u),quad)
cellvec = integrate(v⊙b,quad)
cellmatvec = pair_arrays(cellmat,cellvec)
cellids = collect(1:num_cells(trian))

bcellmat = integrate(v*u,bquad)
bcellvec = integrate(v*3,bquad)
bcellmatvec = pair_arrays(bcellmat,bcellvec)
bcellids = btrian.glue.face_to_cell

b0cellmat = integrate(v*u,b0quad)
b0cellvec = integrate(v*3,b0quad)
b0cellmatvec = pair_arrays(b0cellmat,b0cellvec)
b0cellids = get_cell_to_bgcell(b0trian)
@test length(b0cellids) == 0

term_to_cellmat = [cellmat, bcellmat, b0cellmat]
term_to_cellvec = [cellvec, bcellvec, b0cellvec]
term_to_cellids = [cellids, bcellids, b0cellids]
term_to_cellmatvec = [ cellmatvec, bcellmatvec, b0cellmatvec ]

mtypes = [
  SparseMatrixCSC,
  SparseMatrixCSR,
  SparseMatrixCSR{0},
  SparseMatrixCSR{1},
  SparseMatrixCSR{0,Float64, Int},
  SparseMatrixCSR{1,Float64, Int},
  SymSparseMatrixCSR,
  SymSparseMatrixCSR{0},
  SymSparseMatrixCSR{1},
  SymSparseMatrixCSR{0, Float64, Int},
  SymSparseMatrixCSR{1, Float64, Int}]

for T in mtypes

  matvecdata = ( term_to_cellmatvec , term_to_cellids, term_to_cellids)
  matdata = (term_to_cellmat,term_to_cellids,term_to_cellids)
  vecdata = (term_to_cellvec,term_to_cellids)
  data = (matvecdata,matdata,vecdata)

  assem = SparseMatrixAssembler(T,Vector{Float64},U,V)
  test_sparse_matrix_assembler(assem,matdata,vecdata,data)

  matdata = ([cellmat],[cellids],[cellids])
  vecdata = ([cellvec],[cellids])

  mat = assemble_matrix(assem,matdata)
  vec = assemble_vector(assem,vecdata)

  x = mat \ vec

  assemble_matrix!(mat,assem,matdata)
  assemble_vector!(vec,assem,vecdata)

  x2 = mat \ vec
  @test x ≈ x2

  @test vec ≈ [0.0625, 0.125, 0.0625]
  @test mat[1, 1]  ≈  1.333333333333333
  @test mat[2, 1]  ≈ -0.33333333333333
  @test mat[1, 2]  ≈ -0.33333333333333
  @test mat[2, 2]  ≈ 2.666666666666666
  @test mat[3, 2]  ≈ -0.33333333333333
  @test mat[2, 3]  ≈ -0.33333333333333
  @test mat[3, 3]  ≈ 1.333333333333333

  data = (([cellmatvec],[cellids],[cellids]),([],[],[]),([],[]))
  mat, vec = allocate_matrix_and_vector(assem,data)
  assemble_matrix_and_vector!(mat,vec,assem,data)
  assemble_matrix_and_vector!(mat,vec,assem,data)

  @test vec ≈ [0.0625, 0.125, 0.0625]
  @test mat[1, 1]  ≈  1.333333333333333
  @test mat[2, 1]  ≈ -0.33333333333333

  x3 = mat \ vec
  @test x ≈ x3

  mat, vec = assemble_matrix_and_vector(assem,data)

  x4 = mat \ vec
  @test x ≈ x4

  @test vec ≈ [0.0625, 0.125, 0.0625]
  @test mat[1, 1]  ≈  1.333333333333333
  @test mat[2, 1]  ≈ -0.33333333333333

end

strian = SkeletonTriangulation(model)
squad = CellQuadrature(strian,degree)

scellmat = integrate(jump(v)*u.⁻,squad)
@test isa(scellmat[1],BlockArrayCoo)
@test is_zero_block(scellmat[1],1,1)
@test is_zero_block(scellmat[1],2,1)
@test is_nonzero_block(scellmat[1],1,2)
@test is_nonzero_block(scellmat[1],2,2)

scellvec = integrate(mean(v*3),squad)
@test isa(scellvec[1],BlockArrayCoo)
scellmatvec = pair_arrays(scellmat,scellvec)
scellids = get_cell_to_bgcell(strian)
zh = zero(V)
scellvals = get_cell_dof_values(zh,scellids)
scellmatvec = attach_dirichlet(scellmatvec,scellvals)

matvecdata = ([scellmatvec],[scellids],[scellids])
matdata = ([scellmat],[scellids],[scellids])
vecdata = ([scellvec],[scellids])
data = (matvecdata,matdata,vecdata)

assem = SparseMatrixAssembler(U,V)
test_sparse_matrix_assembler(assem,matdata,vecdata,data)

A = assemble_matrix(assem,matdata)
@test A == zeros(num_free_dofs(V),num_free_dofs(U))

# Allowing a non-blocked matrix in a context where you would expect a blocked one.
scellmat = map(i->zeros(size(i)),scellmat)
matdata = ([scellmat],[scellids],[scellids])
A = assemble_matrix(assem,matdata)
@test A == zeros(num_free_dofs(V),num_free_dofs(U))

end # module


text:
# Julia wrapper for header: /home/osboxes/.julia/v0.6/Sundials/deps/usr/include/sunmatrix/sunmatrix_band.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

function SUNBandMatrix(N::sunindextype, mu::sunindextype, ml::sunindextype, smu::sunindextype)
    ccall((:SUNBandMatrix, libsundials_sunmatrixband), SUNMatrix, (sunindextype, sunindextype, sunindextype, sunindextype), N, mu, ml, smu)
end

function SUNBandMatrix_Print(A::SUNMatrix, outfile)
    ccall((:SUNBandMatrix_Print, libsundials_sunmatrixband), Void, (SUNMatrix, Ptr{Cint}), A, outfile)
end

function SUNBandMatrix_Rows(A::SUNMatrix)
    ccall((:SUNBandMatrix_Rows, libsundials_sunmatrixband), sunindextype, (SUNMatrix,), A)
end

function SUNBandMatrix_Columns(A::SUNMatrix)
    ccall((:SUNBandMatrix_Columns, libsundials_sunmatrixband), sunindextype, (SUNMatrix,), A)
end

function SUNBandMatrix_LowerBandwidth(A::SUNMatrix)
    ccall((:SUNBandMatrix_LowerBandwidth, libsundials_sunmatrixband), sunindextype, (SUNMatrix,), A)
end

function SUNBandMatrix_UpperBandwidth(A::SUNMatrix)
    ccall((:SUNBandMatrix_UpperBandwidth, libsundials_sunmatrixband), sunindextype, (SUNMatrix,), A)
end

function SUNBandMatrix_StoredUpperBandwidth(A::SUNMatrix)
    ccall((:SUNBandMatrix_StoredUpperBandwidth, libsundials_sunmatrixband), sunindextype, (SUNMatrix,), A)
end

function SUNBandMatrix_LDim(A::SUNMatrix)
    ccall((:SUNBandMatrix_LDim, libsundials_sunmatrixband), sunindextype, (SUNMatrix,), A)
end

function SUNBandMatrix_Data(A::SUNMatrix)
    ccall((:SUNBandMatrix_Data, libsundials_sunmatrixband), Ptr{realtype}, (SUNMatrix,), A)
end

function SUNBandMatrix_Cols(A::SUNMatrix)
    ccall((:SUNBandMatrix_Cols, libsundials_sunmatrixband), Ptr{Ptr{realtype}}, (SUNMatrix,), A)
end

function SUNBandMatrix_Column(A::SUNMatrix, j::sunindextype)
    ccall((:SUNBandMatrix_Column, libsundials_sunmatrixband), Ptr{realtype}, (SUNMatrix, sunindextype), A, j)
end

function SUNMatGetID_Band(A::SUNMatrix)
    ccall((:SUNMatGetID_Band, libsundials_sunmatrixband), SUNMatrix_ID, (SUNMatrix,), A)
end

function SUNMatClone_Band(A::SUNMatrix)
    ccall((:SUNMatClone_Band, libsundials_sunmatrixband), SUNMatrix, (SUNMatrix,), A)
end

function SUNMatDestroy_Band(A::SUNMatrix)
    ccall((:SUNMatDestroy_Band, libsundials_sunmatrixband), Void, (SUNMatrix,), A)
end

function SUNMatZero_Band(A::SUNMatrix)
    ccall((:SUNMatZero_Band, libsundials_sunmatrixband), Cint, (SUNMatrix,), A)
end

function SUNMatCopy_Band(A::SUNMatrix, B::SUNMatrix)
    ccall((:SUNMatCopy_Band, libsundials_sunmatrixband), Cint, (SUNMatrix, SUNMatrix), A, B)
end

function SUNMatScaleAdd_Band(c::realtype, A::SUNMatrix, B::SUNMatrix)
    ccall((:SUNMatScaleAdd_Band, libsundials_sunmatrixband), Cint, (realtype, SUNMatrix, SUNMatrix), c, A, B)
end

function SUNMatScaleAddI_Band(c::realtype, A::SUNMatrix)
    ccall((:SUNMatScaleAddI_Band, libsundials_sunmatrixband), Cint, (realtype, SUNMatrix), c, A)
end

function __SUNMatMatvec_Band(A::SUNMatrix, x::N_Vector, y::N_Vector)
    ccall((:SUNMatMatvec_Band, libsundials_sunmatrixband), Cint, (SUNMatrix, N_Vector, N_Vector), A, x, y)
end

function SUNMatMatvec_Band(A, x, y)
    __x = convert(NVector, x)
    __y = convert(NVector, y)
    __SUNMatMatvec_Band(A, convert(N_Vector, __x), convert(N_Vector, __y))
end

function SUNMatSpace_Band(A::SUNMatrix, lenrw, leniw)
    ccall((:SUNMatSpace_Band, libsundials_sunmatrixband), Cint, (SUNMatrix, Ptr{Clong}, Ptr{Clong}), A, lenrw, leniw)
end
# Julia wrapper for header: /home/osboxes/.julia/v0.6/Sundials/deps/usr/include/sunmatrix/sunmatrix_dense.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

function SUNDenseMatrix(M::sunindextype, N::sunindextype)
    ccall((:SUNDenseMatrix, libsundials_sunmatrixdense), SUNMatrix, (sunindextype, sunindextype), M, N)
end

function SUNDenseMatrix_Print(A::SUNMatrix, outfile)
    ccall((:SUNDenseMatrix_Print, libsundials_sunmatrixdense), Void, (SUNMatrix, Ptr{Cint}), A, outfile)
end

function SUNDenseMatrix_Rows(A::SUNMatrix)
    ccall((:SUNDenseMatrix_Rows, libsundials_sunmatrixdense), sunindextype, (SUNMatrix,), A)
end

function SUNDenseMatrix_Columns(A::SUNMatrix)
    ccall((:SUNDenseMatrix_Columns, libsundials_sunmatrixdense), sunindextype, (SUNMatrix,), A)
end

function SUNDenseMatrix_LData(A::SUNMatrix)
    ccall((:SUNDenseMatrix_LData, libsundials_sunmatrixdense), sunindextype, (SUNMatrix,), A)
end

function SUNDenseMatrix_Data(A::SUNMatrix)
    ccall((:SUNDenseMatrix_Data, libsundials_sunmatrixdense), Ptr{realtype}, (SUNMatrix,), A)
end

function SUNDenseMatrix_Cols(A::SUNMatrix)
    ccall((:SUNDenseMatrix_Cols, libsundials_sunmatrixdense), Ptr{Ptr{realtype}}, (SUNMatrix,), A)
end

function SUNDenseMatrix_Column(A::SUNMatrix, j::sunindextype)
    ccall((:SUNDenseMatrix_Column, libsundials_sunmatrixdense), Ptr{realtype}, (SUNMatrix, sunindextype), A, j)
end

function SUNMatGetID_Dense(A::SUNMatrix)
    ccall((:SUNMatGetID_Dense, libsundials_sunmatrixdense), SUNMatrix_ID, (SUNMatrix,), A)
end

function SUNMatClone_Dense(A::SUNMatrix)
    ccall((:SUNMatClone_Dense, libsundials_sunmatrixdense), SUNMatrix, (SUNMatrix,), A)
end

function SUNMatDestroy_Dense(A::SUNMatrix)
    ccall((:SUNMatDestroy_Dense, libsundials_sunmatrixdense), Void, (SUNMatrix,), A)
end

function SUNMatZero_Dense(A::SUNMatrix)
    ccall((:SUNMatZero_Dense, libsundials_sunmatrixdense), Cint, (SUNMatrix,), A)
end

function SUNMatCopy_Dense(A::SUNMatrix, B::SUNMatrix)
    ccall((:SUNMatCopy_Dense, libsundials_sunmatrixdense), Cint, (SUNMatrix, SUNMatrix), A, B)
end

function SUNMatScaleAdd_Dense(c::realtype, A::SUNMatrix, B::SUNMatrix)
    ccall((:SUNMatScaleAdd_Dense, libsundials_sunmatrixdense), Cint, (realtype, SUNMatrix, SUNMatrix), c, A, B)
end

function SUNMatScaleAddI_Dense(c::realtype, A::SUNMatrix)
    ccall((:SUNMatScaleAddI_Dense, libsundials_sunmatrixdense), Cint, (realtype, SUNMatrix), c, A)
end

function __SUNMatMatvec_Dense(A::SUNMatrix, x::N_Vector, y::N_Vector)
    ccall((:SUNMatMatvec_Dense, libsundials_sunmatrixdense), Cint, (SUNMatrix, N_Vector, N_Vector), A, x, y)
end

function SUNMatMatvec_Dense(A, x, y)
    __x = convert(NVector, x)
    __y = convert(NVector, y)
    __SUNMatMatvec_Dense(A, convert(N_Vector, __x), convert(N_Vector, __y))
end

function SUNMatSpace_Dense(A::SUNMatrix, lenrw, leniw)
    ccall((:SUNMatSpace_Dense, libsundials_sunmatrixdense), Cint, (SUNMatrix, Ptr{Clong}, Ptr{Clong}), A, lenrw, leniw)
end
# Julia wrapper for header: /home/osboxes/.julia/v0.6/Sundials/deps/usr/include/sunmatrix/sunmatrix_sparse.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

function __SUNSparseMatrix(M::sunindextype, N::sunindextype, NNZ::sunindextype, sparsetype::Cint)
    ccall((:SUNSparseMatrix, libsundials_sunmatrixsparse), SUNMatrix, (sunindextype, sunindextype, sunindextype, Cint), M, N, NNZ, sparsetype)
end

function SUNSparseMatrix(M, N, NNZ, sparsetype)
    __SUNSparseMatrix(M, N, NNZ, convert(Cint, sparsetype))
end

function __SUNSparseFromDenseMatrix(A::SUNMatrix, droptol::realtype, sparsetype::Cint)
    ccall((:SUNSparseFromDenseMatrix, libsundials_sunmatrixdense), SUNMatrix, (SUNMatrix, realtype, Cint), A, droptol, sparsetype)
end

function SUNSparseFromDenseMatrix(A, droptol, sparsetype)
    __SUNSparseFromDenseMatrix(A, droptol, convert(Cint, sparsetype))
end

function __SUNSparseFromBandMatrix(A::SUNMatrix, droptol::realtype, sparsetype::Cint)
    ccall((:SUNSparseFromBandMatrix, libsundials_sunmatrixband), SUNMatrix, (SUNMatrix, realtype, Cint), A, droptol, sparsetype)
end

function SUNSparseFromBandMatrix(A, droptol, sparsetype)
    __SUNSparseFromBandMatrix(A, droptol, convert(Cint, sparsetype))
end

function SUNSparseMatrix_Realloc(A::SUNMatrix)
    ccall((:SUNSparseMatrix_Realloc, libsundials_sunmatrixsparse), Cint, (SUNMatrix,), A)
end

function SUNSparseMatrix_Print(A::SUNMatrix, outfile)
    ccall((:SUNSparseMatrix_Print, libsundials_sunmatrixsparse), Void, (SUNMatrix, Ptr{Cint}), A, outfile)
end

function SUNSparseMatrix_Rows(A::SUNMatrix)
    ccall((:SUNSparseMatrix_Rows, libsundials_sunmatrixsparse), sunindextype, (SUNMatrix,), A)
end

function SUNSparseMatrix_Columns(A::SUNMatrix)
    ccall((:SUNSparseMatrix_Columns, libsundials_sunmatrixsparse), sunindextype, (SUNMatrix,), A)
end

function SUNSparseMatrix_NNZ(A::SUNMatrix)
    ccall((:SUNSparseMatrix_NNZ, libsundials_sunmatrixsparse), sunindextype, (SUNMatrix,), A)
end

function SUNSparseMatrix_NP(A::SUNMatrix)
    ccall((:SUNSparseMatrix_NP, libsundials_sunmatrixsparse), sunindextype, (SUNMatrix,), A)
end

function SUNSparseMatrix_SparseType(A::SUNMatrix)
    ccall((:SUNSparseMatrix_SparseType, libsundials_sunmatrixsparse), Cint, (SUNMatrix,), A)
end

function SUNSparseMatrix_Data(A::SUNMatrix)
    ccall((:SUNSparseMatrix_Data, libsundials_sunmatrixsparse), Ptr{realtype}, (SUNMatrix,), A)
end

function SUNSparseMatrix_IndexValues(A::SUNMatrix)
    ccall((:SUNSparseMatrix_IndexValues, libsundials_sunmatrixsparse), Ptr{sunindextype}, (SUNMatrix,), A)
end

function SUNSparseMatrix_IndexPointers(A::SUNMatrix)
    ccall((:SUNSparseMatrix_IndexPointers, libsundials_sunmatrixsparse), Ptr{sunindextype}, (SUNMatrix,), A)
end

function SUNMatGetID_Sparse(A::SUNMatrix)
    ccall((:SUNMatGetID_Sparse, libsundials_sunmatrixsparse), SUNMatrix_ID, (SUNMatrix,), A)
end

function SUNMatClone_Sparse(A::SUNMatrix)
    ccall((:SUNMatClone_Sparse, libsundials_sunmatrixsparse), SUNMatrix, (SUNMatrix,), A)
end

function SUNMatDestroy_Sparse(A::SUNMatrix)
    ccall((:SUNMatDestroy_Sparse, libsundials_sunmatrixsparse), Void, (SUNMatrix,), A)
end

function SUNMatZero_Sparse(A::SUNMatrix)
    ccall((:SUNMatZero_Sparse, libsundials_sunmatrixsparse), Cint, (SUNMatrix,), A)
end

function SUNMatCopy_Sparse(A::SUNMatrix, B::SUNMatrix)
    ccall((:SUNMatCopy_Sparse, libsundials_sunmatrixsparse), Cint, (SUNMatrix, SUNMatrix), A, B)
end

function SUNMatScaleAdd_Sparse(c::realtype, A::SUNMatrix, B::SUNMatrix)
    ccall((:SUNMatScaleAdd_Sparse, libsundials_sunmatrixsparse), Cint, (realtype, SUNMatrix, SUNMatrix), c, A, B)
end

function SUNMatScaleAddI_Sparse(c::realtype, A::SUNMatrix)
    ccall((:SUNMatScaleAddI_Sparse, libsundials_sunmatrixsparse), Cint, (realtype, SUNMatrix), c, A)
end

function __SUNMatMatvec_Sparse(A::SUNMatrix, x::N_Vector, y::N_Vector)
    ccall((:SUNMatMatvec_Sparse, libsundials_sunmatrixsparse), Cint, (SUNMatrix, N_Vector, N_Vector), A, x, y)
end

function SUNMatMatvec_Sparse(A, x, y)
    __x = convert(NVector, x)
    __y = convert(NVector, y)
    __SUNMatMatvec_Sparse(A, convert(N_Vector, __x), convert(N_Vector, __y))
end

function SUNMatSpace_Sparse(A::SUNMatrix, lenrw, leniw)
    ccall((:SUNMatSpace_Sparse, libsundials_sunmatrixsparse), Cint, (SUNMatrix, Ptr{Clong}, Ptr{Clong}), A, lenrw, leniw)
end


text:
# This file is auto-generated by AWSMetadata.jl
using AWS
using AWS.AWSServices: appconfigdata
using AWS.Compat
using AWS.UUIDs

"""
    get_latest_configuration(configuration_token)
    get_latest_configuration(configuration_token, params::Dict{String,<:Any})

Retrieves the latest deployed configuration. This API may return empty Configuration data
if the client already has the latest version. See StartConfigurationSession to obtain an
InitialConfigurationToken to call this API.  Each call to GetLatestConfiguration returns a
new ConfigurationToken (NextPollConfigurationToken in the response). This new token MUST be
provided to the next call to GetLatestConfiguration when polling for configuration updates.
To avoid excess charges, we recommend that you include the ClientConfigurationVersion value
with every call to GetConfiguration. This value must be saved on your client. Subsequent
calls to GetConfiguration must pass this value by using the ClientConfigurationVersion
parameter.

# Arguments
- `configuration_token`: Token describing the current state of the configuration session.
  To obtain a token, first call the StartConfigurationSession API. Note that every call to
  GetLatestConfiguration will return a new ConfigurationToken (NextPollConfigurationToken in
  the response) and MUST be provided to subsequent GetLatestConfiguration API calls.

"""
function get_latest_configuration(
    configuration_token; aws_config::AbstractAWSConfig=global_aws_config()
)
    return appconfigdata(
        "GET",
        "/configuration",
        Dict{String,Any}("configuration_token" => configuration_token);
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function get_latest_configuration(
    configuration_token,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return appconfigdata(
        "GET",
        "/configuration",
        Dict{String,Any}(
            mergewith(
                _merge,
                Dict{String,Any}("configuration_token" => configuration_token),
                params,
            ),
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    start_configuration_session(application_identifier, configuration_profile_identifier, environment_identifier)
    start_configuration_session(application_identifier, configuration_profile_identifier, environment_identifier, params::Dict{String,<:Any})

Starts a configuration session used to retrieve a deployed configuration. See the
GetLatestConfiguration API for more details.

# Arguments
- `application_identifier`: The application ID or the application name.
- `configuration_profile_identifier`: The configuration profile ID or the configuration
  profile name.
- `environment_identifier`: The environment ID or the environment name.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"RequiredMinimumPollIntervalInSeconds"`: The interval at which your client will poll for
  configuration. If provided, the service will throw a BadRequestException if the client
  polls before the specified poll interval. By default, client poll intervals are not
  enforced.
"""
function start_configuration_session(
    ApplicationIdentifier,
    ConfigurationProfileIdentifier,
    EnvironmentIdentifier;
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return appconfigdata(
        "POST",
        "/configurationsessions",
        Dict{String,Any}(
            "ApplicationIdentifier" => ApplicationIdentifier,
            "ConfigurationProfileIdentifier" => ConfigurationProfileIdentifier,
            "EnvironmentIdentifier" => EnvironmentIdentifier,
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function start_configuration_session(
    ApplicationIdentifier,
    ConfigurationProfileIdentifier,
    EnvironmentIdentifier,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return appconfigdata(
        "POST",
        "/configurationsessions",
        Dict{String,Any}(
            mergewith(
                _merge,
                Dict{String,Any}(
                    "ApplicationIdentifier" => ApplicationIdentifier,
                    "ConfigurationProfileIdentifier" => ConfigurationProfileIdentifier,
                    "EnvironmentIdentifier" => EnvironmentIdentifier,
                ),
                params,
            ),
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end