text:
abstract type AbstractRenderer end

struct MakieRenderer <: AbstractRenderer
    scene::AbstractPlotting.Scene
    image::Observables.Observable{Array{ColorTypes.RGBA{Float32}, 2}}
    cudaimage::CuArray{ColorTypes.RGBA{Float32}, 2}
    numiters::Int64 # Should this be here? Where would be better?
    aspectratio::Float64
    height::Int64
    width::Int64
    #camera::AbstractPlotting::Camera2D
end

scene(renderer::MakieRenderer) = renderer.scene

function init_renderer()
    #scene = Makie.Scene()
    height = 360
    aspectratio = 16/9
    width = Int(floor(height * aspectratio))
    numiters = 250
    xrange, yrange, mandelimg = rendermandelbrotimagecuda(0., 0., 0., numiters, height, aspectratio)
    mandelimgcolor = [ColorTypes.RGBA(Float32(color), Float32(color), Float32(color)) for color in mandelimg]
    cudaimage = CUDA.CuArray(mandelimgcolor)

    #cam =

    #Node = Observables.Node
    #camera = AbstractPlotting.Camera2D(AbstractPlotting.Node(FRect(0,0,width,height)), AbstractPlotting.Node(0.0f0),
        #AbstractPlotting.Node(nothing), AbstractPlotting.Node(nothing), AbstractPlotting.Node(0.0f0),
        #AbstractPlotting.Node(Vec{2, Int}(width, height)), AbstractPlotting.Node(false))
    theme = AbstractPlotting.Attributes(show_axis=false, raw=false, scale_plot=true,
                        padding=Point3(0.0f0, 0.0f0, 0.0f0),
                        align= (:middle, :middle),
                        #resolution=(640, 360), #limits=FRect(0, 0, 640, 360),
                        panbutton=nothing, update_limits=true)
    AbstractPlotting.set_theme!(theme)

    scene = Makie.image(mandelimgcolor)
    #cam = AbstractPlotting.cam2d!(scene; padding=0.000, area=AbstractPlotting.Node(FRect(0,0,width,height)))
    #cam = AbstractPlotting.cam2d!(scene; padding=0.000, area=AbstractPlotting.Node(FRect(0,0,width,height)))
    #AbstractPlotting.update_cam!(scene, cam)
    #AbstractPlotting.update!(scene)
    display(scene)
    #AbstractPlotting.update_cam!(scene)
    #AbstractPlotting.update!(scene)
    #display(scene)
    image = scene.plots[1][:image]
    Observables.notify!(image)
    #AbstractPlotting.update_cam!(scene, camera)
    cam = AbstractPlotting.cam2d!(scene; padding=0.000, area=AbstractPlotting.Node(FRect(0,0,width,height)),
        panbutton=nothing, update_limits=true)
    AbstractPlotting.update_cam!(scene, cam)
    MakieRenderer(scene, image, cudaimage, numiters, aspectratio, height, width)
end

function render_game(renderer::MakieRenderer, gamestate::AbstractGameState)
    numtype = Float64
    # First gather the parameters
    # Renderer parameters
    numiters = renderer.numiters
    outimage = renderer.image
    cudaimage = renderer.cudaimage
    aspectratio = renderer.aspectratio
    height = renderer.height
    width = renderer.width

    numthreads = (16, 16)
    numblocks = ceil(Int, width / numthreads[1]), ceil(Int, height / numthreads[2])

    # GameState parameters
    center = preimage_position(player(gamestate))
    centerx, centery = real(center), imag(center) #.+ 0.4 .* randn.()
    #rotationfactor, scalefactor = rotation(player(gamestate)) + 0.5 * randn(), scale(player(gamestate)) + 0.5 * randn()
    rotationfactor, scalefactor = rotation(player(gamestate)), scale(player(gamestate))

    expscale = numtype(2 ^ -scalefactor)
    yrangeextent = numtype(expscale)
    xrangeextent = numtype(expscale * aspectratio)
    xstart = numtype(-xrangeextent / 2)
    ystart = numtype(-yrangeextent / 2)
    cosrot = numtype(cos(rotationfactor))
    sinrot = numtype(sin(rotationfactor))

    #mobius = current_transform(world(gamestate))
    conformal = current_transform(world(gamestate))
    #a, b, c, d = mobius.a, mobius.b, mobius.c, mobius.d

    # then call
    CUDA.@sync begin
        #@cuda threads=numthreads blocks=numblocks mandelbrotandregiongpu!(cudaimage,
            #centerx, xstart, xrangeextent, centery, ystart, yrangeextent,
            #numiters, height, width, cosrot, sinrot, a, b, c, d)
        @cuda threads=numthreads blocks=numblocks mandelbrotandregiongpurational!(cudaimage,
            centerx, xstart, xrangeextent, centery, ystart, yrangeextent,
            numiters, height, width, cosrot, sinrot, conformal)
    end

    CUDA.copyto!(outimage[], cudaimage)
    Observables.notify!(outimage)
    #AbstractPlotting.update_cam!(renderer.scene)
    sleep(0.0001)
    outimage[]
end

max_stars_repo_name: zoemcc/fractrixzeta

max_issues_repo_name: zoemcc/fractrixzeta

max_forks_repo_name: zoemcc/fractrixzeta


text:
"""
    Transform

Abstract supertype for all feature Transforms.
"""
abstract type Transform end

# Make Transforms callable types
(t::Transform)(x; kwargs...) = apply(x, t; kwargs...)

"""
    transform(::T, data)

Defines the feature engineering pipeline for some type `T`, which comprises a collection of
[`Transform`](@ref)s and other steps to be peformed on the `data`.

The idea around a "transform interface” is to make feature transformations composable, i.e.
the output of any one `Transform` should be valid input to another.

Feature engineering pipelines should obey the same principle and it should be trivial to
add/remove `Transform` steps that compose the pipeline without it breaking.

`transform` should be overloaded for custom types `T` that require feature engineering.
The only requirement is that the return of `transform `is itself "transformable", i.e.
calling [`is_transformable`](@ref) on the output returns true.
"""
function transform end

"""
    transform!(::T, data)

Mutating version of [`transform`](@ref).
"""
function transform! end

ext: jl


text:
using CSV
using DataFrames
using Plots
using JLD2
using DiscretisedFluidQueues
import Distributions
include("/Users/anguslewis/Documents/SFFMProject/DFQExamples.jl/plot_utils.jl")
include("/Users/anguslewis/Documents/SFFMProject/DFQExamples.jl/error_metrics.jl")
fun_number = 6
@halfwidth_plot_defaults()

plotlyjs()
function p()

ks_data = CSV.read((@__DIR__)*"/../data/meshs_ks_func_count_$(fun_number).csv",DataFrame)

ks_data.DG_lin_lim = CSV.read((@__DIR__)*"/../data2/meshs_ks_func_count_$(fun_number).csv",DataFrame).DG_limiter

plot()
linestyles_vec = [:solid,:dash,:dashdot, :dot,:dash]
markerstyles_vec = [:cross,:diamond,:circle,:dot,:diamond]
for (c,col) in enumerate(names(ks_data))
    os = collect(1:2:21)
    (c==5)&&(os.+=1)
    plot!(log10.(os),ks_data[:,col],
        label=col,linestyle=linestyles_vec[c],
        marker=markerstyles_vec[c],
        linewidth=2)
end
plot!(xlabel="Dimension")
plot!(ylabel="Error"); error_ticks!(plot!())
plot!(title="KS error - CDF")
# plot!(legend=(0.275,0.625))
plot!()
@add_lines!(ks_data,("Unif","QBDRAP"),@__DIR__)
savefig((@__DIR__)*"/meshs_ks_error_formatted.pdf")
end
p()

function p()
l2_pdf_data = CSV.read((@__DIR__)*"/../data/meshs_l2_pdf_func_count_$(fun_number).csv",DataFrame)
l2_pdf_data.DG_lin_lim = CSV.read((@__DIR__)*"/../data2/meshs_l2_pdf_func_count_$(fun_number).csv",DataFrame).DG_limiter

plot()
linestyles_vec = [:solid,:dash,:dashdot, :dot,:dash]
markerstyles_vec = [:cross,:diamond,:circle,:cross,:diamond]
for (c,col) in enumerate(names(l2_pdf_data))
    os = collect(1:2:21)
    (c==5)&&(os.+=1)
    plot!(log10.(os),l2_pdf_data[:,col],
        label=col,linestyle=linestyles_vec[c],
        marker=markerstyles_vec[c],
        linewidth=2)
end
plot!(xlabel="Dimension")
plot!(ylabel="Error"); error_ticks!(plot!())
plot!(title="L² error - PDF")
# plot!(legend=(0.25,0.625))
# plot!(ylims=ylims(plot!()).+(-0.3,0))
@add_lines!(l2_pdf_data,("Unif","QBDRAP"),@__DIR__)
savefig((@__DIR__)*"/meshs_l2_pdf_error_formatted.pdf")
end
p()

function p()
l1_cdf_data = CSV.read((@__DIR__)*"/../data/meshs_l1_cdf_func_count_$(fun_number).csv",DataFrame)
l1_cdf_data.DG_lin_lim = CSV.read((@__DIR__)*"/../data2/meshs_l1_cdf_func_count_$(fun_number).csv",DataFrame).DG_limiter

@halfwidth_plot_defaults()
plot()
linestyles_vec = [:solid,:dash,:dashdot, :dot,:dash]
markerstyles_vec = [:cross,:diamond,:circle,:cross,:diamond]
for (c,col) in enumerate(names(l1_cdf_data))
    os = collect(1:2:21)
    (c==5)&&(os.+=1)
    plot!(log10.(os),l1_cdf_data[:,col],
        label=col,linestyle=linestyles_vec[c],
        marker=markerstyles_vec[c],
        linewidth=2)
end
plot!(xlabel="Dimension")
plot!(ylabel="Error"); error_ticks!(plot!())
plot!(title="L¹ error - CDF")
# plot!(legend=(0.225,0.625))
plot!()
@add_lines!(l1_cdf_data,("Unif","QBDRAP"),@__DIR__)
savefig((@__DIR__)*"/meshs_l1_cdf_error_formatted.pdf")
end
p()

function p()
@fullwidth_plot_defaults()
dir = @__DIR__()
os = 1:2:7
plot(layout=(5,length(os)))
linestyles_vec = [:solid,:solid,:solid,:solid]
markerstyles_vec = [:cross,:diamond,:circle]
xs = range(eps(),10-eps(),length=200)
for (c,o) in enumerate(os)
    model = BoundedFluidQueue(fill(0.0,2,2),[1.0; 0.0],10.0)

    mesh = DGMesh(0.0:10.0,o)
    dq = DiscretisedFluidQueue(model,mesh)
    jldopen(dir*"/../../../coeffs_matrix.jld2") do f
        coeffs = f["coeff_matrix"][fun_number,c,1]
        dt = SFMDistribution(coeffs,dq)
        plot!(xs,pdf(dt).(xs,1),subplot=c,label=(c==1 ? "DG" : false),
            linestyle=linestyles_vec[1],ylims = (-0.1,1.1),  yticks=false,
            # marker=markerstyles_vec[1],
            linewidth=2,xticks=false)
    end
    plot!(subplot=c, yticks=(c==1 ? (0.1:0.4:1.2) : false),
        ylabel=(c==1 ? "DG" : false), grid=false,
        title = string("Dim. ",o))

    mesh = DGMesh(0.0:(1.0./o):10.0,1)
    dq = DiscretisedFluidQueue(model,mesh)
    jldopen(dir*"/../../../coeffs_matrix.jld2") do f
        coeffs = f["coeff_matrix"][fun_number,c,2]
        dt = SFMDistribution(coeffs,dq)
        plot!(xs,pdf(dt).(xs,1),subplot=c+4,label=(c==1 ? "Unif." : false),
            linestyle=linestyles_vec[2],ylims = (-0.1,1.1),  yticks=false,
            # marker=markerstyles_vec[2],
            linewidth=2,xticks=false)
    end
    (c==1)&&plot!(subplot=c+4, yticks=(c==1 ? (0.1:0.4:1.2) : false),
        ylabel=(c==1 ? "Unif." : false), grid=false)

    mesh = DGMesh(0.0:10.0,o)
    dq = DiscretisedFluidQueue(model,mesh)
    jldopen(dir*"/../../../coeffs_matrix.jld2") do f
        coeffs = f["coeff_matrix"][fun_number,c,3]
        dt = SFMDistribution(coeffs,dq)
        plot!(xs,pdf(dt).(xs,1),subplot=c+8,label=(c==1 ? "DG (limiter)" : false),
            linestyle=linestyles_vec[3],ylims = (-0.1,1.1),  yticks=false,
            # marker=markerstyles_vec[3],
            linewidth=2,xticks=false)
    end
    (c==1)&&plot!(subplot=c+8, yticks=(c==1 ? (0.1:0.4:1.2) : false),
        ylims = (-0.1,1.1),
        ylabel=(c==1 ? "DG-limit" : false), grid=false)

    mesh = FRAPMesh(0.0:10.0,o)
    dq = DiscretisedFluidQueue(model,mesh)
    jldopen(dir*"/../../../coeffs_matrix.jld2") do f
        coeffs = f["coeff_matrix"][fun_number,c,4]
        dt = SFMDistribution(coeffs,dq)
        plot!(xs,pdf(dt).(xs,1),subplot=c+12,label=(c==1 ? "Q-R" : false),
            linestyle=linestyles_vec[3],ylims = (-0.1,1.1), yticks=false,
            # marker=markerstyles_vec[3],
            linewidth=2,xticks=0:2:10, grid=false)
    end
    (c==1)&&plot!(subplot=c+12, yticks=(c==1 ? (0.1:0.4:1.2) : false),
        ylims = (-0.1,1.1),
        ylabel=(c==1 ? "Q-R" : false), grid=false)
    # plot!(xs,x->2*x*(x<0.5) + Float64(x>=0.5),subplot=c,label=(c==1 ? "truth" : false),
    #     linestyle=:dot,
    #     # marker=markerstyles_vec[3],
    #     color=:grey,
    #     linewidth=2)

    # plot!(subplot=c, yticks=(c==1 ? (0.1:0.3:1.2) : false),
    #     ylims = (-0.1,1.1), xticks = 0.0:3:10,
    #     ylabel=(c==1 ? "probability" : false), grid=false,
    #     title = string("Dim. ",o),
    #     legend=:outertopright)
    plot!(legend=false)

    model = BoundedFluidQueue([0 0; 1.0 -1.0],[1.0; 0.0],10.0)

    mesh = DGMesh(0.0:(2/(o+1)):10.0,2)
    dq = DiscretisedFluidQueue(model,mesh)
    thepdf0(z,i) = (norm_pdf(z)/(1-norm_cdf(0)))*(i==1)*(z>=0.0)
    d0 = SFMDistribution(thepdf0,dq)
    dt = integrate_time(d0,build_full_generator(dq),4.0, StableRK4(0.005); limiter=GeneralisedMUSCL)

    plot!(xs,pdf(dt).(xs,1),subplot=c+16,label=(c==1 ? "DG lin lim" : false),
        linestyle=linestyles_vec[3],ylims = (-0.1,1.1), yticks=false,
        # marker=markerstyles_vec[3],
        linewidth=2,xticks=0:2:10, grid=false)
    (c==1)&&plot!(subplot=c+16, yticks=(c==1 ? (0.1:0.4:1.2) : false),
        ylims = (-0.1,1.1),
        ylabel=(c==1 ? "DG lin lim" : false), grid=false)
    plot!(legend=false)
end
plot!()
savefig((@__DIR__)*"/pdf_formatted.pdf")
end
p()

function p()
@halfwidth_plot_defaults()
L1_cell_probs_errors = DataFrame(DG = Float64[], Unif = Float64[], DG_limiter = Float64[], QBDRAP = Float64[], DG_lin_lim = Float64[])
# truth = zeros(10)
# truth[5] = 1.0
norm_cdf(x) = Distributions.cdf(Distributions.Normal(2.5,0.5),x)
i=1
F4(z) = (norm_cdf(z-4)-norm_cdf(0))/(norm_cdf(10)-norm_cdf(0))*(i==1)*(z>=4)*(z<10) + (z>=10)
truth = F4.(1:10)-F4.(0:9)
truth_pm = zeros(3)
for (c_o,o) in enumerate(1:2:21)
    row = zeros(5)
    for (c_m,m) in enumerate(names(L1_cell_probs_errors))
        (c_m<5)&&jldopen((@__DIR__)*"/../../../coeffs_matrix.jld2") do f
            model = BoundedFluidQueue(fill(0.0,2,2),[1.0; 0.0],10.0)

            mesh = (c_m<4) ? DGMesh(0.0:10.0,o) : FRAPMesh(0.0:10.0,o)
            dq = DiscretisedFluidQueue(model,mesh)

            coeffs = f["coeff_matrix"][fun_number,c_o,c_m]
            dt = SFMDistribution(coeffs,dq)

            x_vals = range(0.5,9.5;length=10)

            approx = sum(cell_probs(dt).(x_vals,(1:2)'),dims=2)
            approx_pm = [coeffs[1:N₋(model)];coeffs[end-N₊(model)+1:end]]

            row[c_m] = log10(sum(abs.(approx-truth)) + sum(abs.(approx_pm-truth_pm)))
        end
        (c_m==5)&&begin
            model = BoundedFluidQueue([0 0; 1.0 -1.0],[1.0; 0.0],10.0)

            mesh = DGMesh(0.0:(2/(o+1)):10.0,2)
            dq = DiscretisedFluidQueue(model,mesh)
            thepdf0(z,i) = (norm_pdf(z)/(1-norm_cdf(0)))*(i==1)*(z>=0.0)
            d0 = SFMDistribution(thepdf0,dq)
            dt = integrate_time(d0,build_full_generator(dq),4.0, StableRK4(0.005); limiter=GeneralisedMUSCL)
            coeffs = dt.coeffs

            x_vals = (mesh.nodes[1:end-1]+mesh.nodes[2:end])./2

            approx = sum(reshape(cell_probs(dt).(x_vals,1),Int((o+1)/2),10),dims=1)[:]
            approx_pm = [coeffs[1:N₋(model)];coeffs[end-N₊(model)+1:end]]

            row[c_m] = log10(sum(abs.(approx-truth)) + sum(abs.(approx_pm-truth_pm)))
        end
    end
    push!(L1_cell_probs_errors,row)
end

file = (@__DIR__)*"/../data/meshs_l1_cell_"*"func_count_$(fun_number)"
CSV.write(file*".csv",L1_cell_probs_errors)
q = plot()
linestyles_vec = [:solid,:dash,:dashdot, :dot,:dash]
markerstyles_vec = [:cross,:diamond,:circle,:dot,:diamond]
for (c,reconstruction) in enumerate(names(L1_cell_probs_errors))
    os = collect(1:2:21)
    (c==5)&&(os.+=1)
    plot!(log10.(os),L1_cell_probs_errors[:,reconstruction],label=reconstruction,
        linestyle=linestyles_vec[c],
        marker=markerstyles_vec[c],
        linewidth=2)
end
plot!(xlabel="Dimension", ylabel="Error",
    title="Error between cell probabilities",)
    # legend=(0.225,0.625))
error_ticks!(plot!())
@add_lines!(L1_cell_probs_errors,("Unif","QBDRAP"),@__DIR__)
savefig((@__DIR__)*"/L1_cell_probs.pdf")
end
p()

function p()
    @fullwidth_plot_defaults
    model = BoundedFluidQueue([0 0; 1.0 -1.0],[1.0; 0.0],10.0)

    mesh = DGMesh(0.0:(2/(21+1)):10.0,2)
    dq = DiscretisedFluidQueue(model,mesh)
    thepdf0(z,i) = (norm_pdf(z)/(1-norm_cdf(0)))*(i==1)*(z>=0.0)
    d0 = SFMDistribution(thepdf0,dq)
    dt = integrate_time(d0,build_full_generator(dq),4.0, StableRK4(0.005); limiter=GeneralisedMUSCL)

    x_vals = range(5,8,length=201)

    dg_lin_lim = pdf(dt).(x_vals,1)

    frapmesh = FRAPMesh(0.0:1.0:10.0,21)
    frapdq = DiscretisedFluidQueue(model,frapmesh)
    frapd0 = SFMDistribution(thepdf0,frapdq)
    frapdt = integrate_time(frapd0,build_full_generator(frapdq),4.0, StableRK4(0.005))

    frap = pdf(frapdt).(x_vals,1)

    # linestyles_vec = [:solid,:dash,:dashdot, :dot,:dash]
    # markerstyles_vec = [:cross,:diamond,:circle,:dot,:diamond]

    plot(x_vals,x->thepdf0(x-4.0,1),label="true PDF",linestyle=:solid,color=:grey)
    plot!(x_vals,frap,label="QBDRAP",linestyle=:solid,color=4,linewidth=2)
    plot!(x_vals,dg_lin_lim,label="DG-lin-lim",linestyle=:solid,color=5,linewidth=2)
    plot!(xlabel="x",ylabel="Density",title="PDF")
    savefig((@__DIR__)*"/order_21_PDFs_formatted.pdf")
end
p()


text:
# L0 pseudo-norm (times a constant)

export NormL0

"""
**``L_0`` pseudo-norm**

    NormL0(λ=1)

Returns the function
```math
f(x) = λ\\cdot\\mathrm{nnz}(x)
```
for a nonnegative parameter `λ`.
"""
struct NormL0{R <: Real} <: ProximableFunction
    lambda::R
    function NormL0{R}(lambda::R) where {R <: Real}
        if lambda < 0
            error("parameter λ must be nonnegative")
        else
            new(lambda)
        end
    end
end

NormL0(lambda::R=1) where {R <: Real} = NormL0{R}(lambda)

function (f::NormL0)(x::AbstractArray{T}) where {R, T <: RealOrComplex{R}}
    return f.lambda * R(count(v -> v != 0, x))
end

function prox!(y::AbstractArray{T}, f::NormL0, x::AbstractArray{T}, gamma::Real=1) where {R, T <: RealOrComplex{R}}
    countnzy = R(0)
    gl = gamma * f.lambda
    for i in eachindex(x)
        over = abs(x[i]) > sqrt(2 * gl)
        y[i] = over * x[i]
        countnzy += over
    end
    return f.lambda * countnzy
end

fun_name(f::NormL0) = "weighted L0 pseudo-norm"
fun_dom(f::NormL0) = "AbstractArray{Real}, AbstractArray{Complex}"
fun_expr(f::NormL0) = "x ↦ λ countnz(x)"
fun_params(f::NormL0) = "λ = $(f.lambda)"

function prox_naive(f::NormL0, x::AbstractArray{T}, gamma::Real=1) where {R, T <: RealOrComplex{R}}
    over = abs.(x) .> sqrt(2 * gamma * f.lambda)
    y = x.*over
    return y, f.lambda * R(count(v -> v != 0, y))
end


text:
module MolSimPy

export MatSciPy

include("MatSciPy.jl")

using .MatSciPy
using ASE: ASEAtoms

end # module


text:
@testset "map" begin
    function _mapview(f, d::AbstractDictionary)
        I = keytype(d)
        T = Core.Compiler.return_type(f, Tuple{eltype(d)})

        return MappedDictionary{I, T, typeof(f), Tuple{typeof(d)}}(f, (d,))
    end

    i = HashIndices([1,2,3,4,5])

    @test issetequal(pairs(map(iseven, i)::HashDictionary), [1=>false, 2=>true, 3=>false, 4=>true, 5=>false])
    @test issetequal(pairs(map(isodd, i)::HashDictionary), [1=>true, 2=>false, 3=>true, 4=>false, 5=>true])

    @test issetequal(pairs(_mapview(iseven, i)::Dictionaries.MappedDictionary), [1=>false, 2=>true, 3=>false, 4=>true, 5=>false])
    @test issetequal(pairs(_mapview(isodd, i)::Dictionaries.MappedDictionary), [1=>true, 2=>false, 3=>true, 4=>false, 5=>true])

    d = HashDictionary([1,2,3,4,5], [1,3,2,4,5])

    @test issetequal(pairs(map(iseven, d)::HashDictionary), [1=>false, 2=>false, 3=>true, 4=>true, 5=>false])
    @test issetequal(pairs(map(isodd, d)::HashDictionary), [1=>true, 2=>true, 3=>false, 4=>false, 5=>true])

    @test issetequal(pairs(_mapview(iseven, d)::Dictionaries.MappedDictionary), [1=>false, 2=>false, 3=>true, 4=>true, 5=>false])
    @test issetequal(pairs(_mapview(isodd, d)::Dictionaries.MappedDictionary), [1=>true, 2=>true, 3=>false, 4=>false, 5=>true])
end


text:
##INVESTIGATE DISTRIBUTION OF RE-ABC LIKELIHOOD ESTIMATES
using RareEventABC
using ProgressMeter
using DataFrames
using Distributions
using Plots; pgfplots()
using StatPlots

##GENERATE OBSERVATIONS
srand(1);
d = 25;
yobs = randn(d) .* 3.0;
a = RareEventABC.IIDNormalDetails(yobs, d);

##GENERATE LIKELIHOOD ESTIMATES
schedule = doSMC(a, [3.0], 5.0, nparticles=200, quiet=true).schedule
llexperiment = Dict();
@showprogress for nparticles in 10:10:100
    llest = Float64[doSMC(a, [3.0], schedule, nparticles=nparticles, d=25, quiet=true).llest for i in 1:100]
    llexperiment[nparticles] = llest
end

##QQ PLOTS
theoretical_quantiles = quantile(Normal(), (1:100)/101)
qqdf = DataFrame(samp=Float64[], theoretical=Float64[], nparticles=Int[])
for nparticles in 20:10:50
    ll = llexperiment[nparticles]
    temp = DataFrame(samp=sort(ll), theoretical=theoretical_quantiles, particles=nparticles)
    qqdf = vcat(qqdf, temp)
end

qqdf2 = qqdf[qqdf[:samp] .> -Inf, :] ##Plots.jl doesn't seem to work with infinite values
default(guidefont=font(25), tickfont=font(25), legendfont=font(25), size=(500,500));
p1 = scatter(qqdf2, :theoretical, :samp, group=:particles,
             markershape=:auto,
             xlims=(-2.5,3.0), ##Hack so that legend isn't in the way
             xlabel="Sample quantile", ylabel="Theoretical quantile");
savefig(p1, "iidnormal_qq.pdf")

##COMPARE ADAPTIVE AND NON-ADAPTIVE ALGORITHMS
b = RareEventABC.IIDNormalDetails([0.0], 1);
eps = 0.2;
##trueL = (cdf(n, eps) - cdf(n,-eps)) / (2*eps)

nparticles = 10;

srand(1);
nreps = 10000;
likelihood_adaptive = Float64[doSMC(b, [1.0], eps, nparticles=nparticles, d=1, quiet=true).llest for i in 1:nreps] |> exp;

srand(2);
schedule = doSMC(b, [1.0], eps, nparticles=nparticles, quiet=true).schedule
likelihood_nonadaptive = Float64[doSMC(b, [1.0], schedule, nparticles=nparticles, d=1, quiet=true).llest for i in 1:nreps] |> exp;

mean(likelihood_adaptive), std(likelihood_adaptive)/sqrt(nreps)
mean(likelihood_nonadaptive), std(likelihood_nonadaptive)/sqrt(nreps)


text:
function judipmap(func, iter::UnitRange; on_error = nothing)
    # Ignore map and workers if single source
    length(iter) == 1 && (return [func(iter[1])])
    # Switch to asyncmap if serial julia
    length(default_worker_pool()) < 2 && (return asyncmap(func, iter))
    # Standard pmap if parallel and multiple sources
    return pmap(func, iter; on_error=on_error)
end

function judipmap(func, model::Array{Model, 1}, args::Vararg{Array, N}) where N
    # overload for multiple models
    argout = Array{Any, 1}(undef, length(model))
    @sync for (i, ai) in enumerate(zip(model, args...))
        @async argout[i] = func(ai...)
    end
    return argout
end


text:
function compute_reservation(solution::Solution, mapf::MAPF; agents=1:nb_agents(mapf))
    reservation = Reservation()
    for a in agents
        path = solution[a]
        update_reservation!(reservation, path, mapf::MAPF)
    end
    return reservation
end

function update_reservation!(reservation::Reservation, path::Path, mapf::MAPF)
    for (t, v) in path
        for g in mapf.group_memberships[v]
            for w in mapf.conflict_groups[g]
                push!(reservation, (t, w))
            end
        end
    end
    return nothing
end


text:
include("init.jl")

open("testconfig.json", "w") do f

    write(f, """{
	"name": "EV_Controller",
	"loglevel": 5,
	"coreType": "zmq",
	"timeDelta": 1.0,

	"endpoints":[
	{
		"name": "EV_Controller/EV6",
		"type": "genmessage",
		"global": true

	},
	{
		"name": "EV_Controller/EV5",
		"type": "genmessage",
		"global": true
	},
	{
		"name": "EV_Controller/EV4",
		"type": "genmessage",
		"global": true
	},
	{
		"name": "EV_Controller/EV3",
		"type": "genmessage",
		"global": true
	},
	{
		"name": "EV_Controller/EV2",
		"type": "genmessage",
		"global": true
	},
	{
		"name": "EV_Controller/EV1",
		"type": "genmessage",
		"global": true
	}
	],
	"filters":
	[
	{
	"name":"filterEV6",  //filters can have names (optional)
	"target":"EV_Controller/EV6",  //this is a key field specifying the target endpoint
	//"inputType":"genmessage",  //can trigger some warnings if there is mismatches for custom filters only used if operation is "custom"
	//"outputType":"genmessage",  //this could be useful if the filter actually translates the data and can be used to automatically order filters
	"mode":"source",  //source is the default but can be "source","destination","clone","cloning"
	"operation":"delay", //currently valid operations are "delay","clone","cloning","timedelay","randomdelay","randomdrop","reroute","redirect","custom"
	"properties":  //additional properties for filters are specified in a property array or object if there is just a single one
	{
		"name":"delay",  //A delay filter just has a single property
		"value":600	//times default to seconds though units can also be specified "200 ms" or similar
	}
	},
	{
	"name":"filterEV5",  //filters can have names (optional)
	"target":"EV_Controller/EV5",  //this is a key field specifying the target endpoint
	//"inputType":"genmessage",  //can trigger some warnings if there is mismatches for custom filters only used if operation is "custom"
	//"outputType":"genmessage",  //this could be useful if the filter actually translates the data and can be used to automatically order filters
	"mode":"source",  //source is the default but can be "source","destination","clone","cloning"
	"operation":"delay", //currently valid operations are "delay","clone","cloning","timedelay","randomdelay","randomdrop","reroute","redirect","custom"
	"properties":  //additional properties for filters are specified in a property array or object if there is just a single one
	{
		"name":"delay",  //A delay filter just has a single property
		"value":600	//times default to seconds though units can also be specified "200 ms" or similar
	}
	},
	{
	"name":"filterEV4",  //filters can have names (optional)
	"target":"EV_Controller/EV4",  //this is a key field specifying the target endpoint
	//"inputType":"genmessage",  //can trigger some warnings if there is mismatches for custom filters only used if operation is "custom"
	//"outputType":"genmessage",  //this could be useful if the filter actually translates the data and can be used to automatically order filters
	"mode":"source",  //source is the default but can be "source","destination","clone","cloning"
	"operation":"delay", //currently valid operations are "delay","clone","cloning","timedelay","randomdelay","randomdrop","reroute","redirect","custom"
	"properties":  //additional properties for filters are specified in a property array or object if there is just a single one
	{
		"name":"delay",  //A delay filter just has a single property
		"value":600	//times default to seconds though units can also be specified "200 ms" or similar
	}
	},
	{
	"name":"filterEV3",  //filters can have names (optional)
	"target":"EV_Controller/EV3",  //this is a key field specifying the target endpoint
	//"inputType":"genmessage",  //can trigger some warnings if there is mismatches for custom filters only used if operation is "custom"
	//"outputType":"genmessage",  //this could be useful if the filter actually translates the data and can be used to automatically order filters
	"mode":"source",  //source is the default but can be "source","destination","clone","cloning"
	"operation":"delay", //currently valid operations are "delay","clone","cloning","timedelay","randomdelay","randomdrop","reroute","redirect","custom"
	"properties":  //additional properties for filters are specified in a property array or object if there is just a single one
	{
		"name":"delay",  //A delay filter just has a single property
		"value":600	//times default to seconds though units can also be specified "200 ms" or similar
	}
	},
	{
	"name":"filterEV2",  //filters can have names (optional)
	"target":"EV_Controller/EV2",  //this is a key field specifying the target endpoint
	//"inputType":"genmessage",  //can trigger some warnings if there is mismatches for custom filters only used if operation is "custom"
	//"outputType":"genmessage",  //this could be useful if the filter actually translates the data and can be used to automatically order filters
	"mode":"source",  //source is the default but can be "source","destination","clone","cloning"
	"operation":"delay", //currently valid operations are "delay","clone","cloning","timedelay","randomdelay","randomdrop","reroute","redirect","custom"
	"properties":  //additional properties for filters are specified in a property array or object if there is just a single one
	{
		"name":"delay",  //A delay filter just has a single property
		"value":600	//times default to seconds though units can also be specified "200 ms" or similar
	}
	},
	{
	"name":"filterEV1",  //filters can have names (optional)
	"target":"EV_Controller/EV1",  //this is a key field specifying the target endpoint
	//"inputType":"genmessage",  //can trigger some warnings if there is mismatches for custom filters only used if operation is "custom"
	//"outputType":"genmessage",  //this could be useful if the filter actually translates the data and can be used to automatically order filters
	"mode":"source",  //source is the default but can be "source","destination","clone","cloning"
	"operation":"delay", //currently valid operations are "delay","clone","cloning","timedelay","randomdelay","randomdrop","reroute","redirect","custom"
	"properties":  //additional properties for filters are specified in a property array or object if there is just a single one
	{
		"name":"delay",  //A delay filter just has a single property
		"value":600	//times default to seconds though units can also be specified "200 ms" or similar
	}
	}
	],
	"subscriptions":[
	{
		"key": "IEEE_123_feeder_0/totalLoad",
		"type": "complex",
		"required": true
	},  {
		"key": "IEEE_123_feeder_0/charge_EV6",
		"type": "complex",
		"required": true
	},	{
		"key": "IEEE_123_feeder_0/charge_EV5",
		"type": "complex",
		"required": true
	},	{
		"key": "IEEE_123_feeder_0/charge_EV4",
		"type": "complex",
		"required": true
	},  {
		"key": "IEEE_123_feeder_0/charge_EV3",
		"type": "complex",
		"required": true
	},  {
		"key": "IEEE_123_feeder_0/charge_EV2",
		"type": "complex",
		"required": true
	},  {
		"key": "IEEE_123_feeder_0/charge_EV1",
		"type": "complex",
		"required": true
	}]
}
""")

end

@testset "Combination Federate" begin

    broker = createBroker();

    cfed = h.helicsCreateCombinationFederateFromConfig("testconfig.json");

    @test h.helicsFederateIsValid(cfed)

    @test h.helicsFederateGetEndpointCount(cfed) == 6
    @test h.helicsFederateGetFilterCount(cfed) == 6
    @test h.helicsFederateGetInputCount(cfed) == 7

    ept = h.helicsFederateGetEndpointByIndex(cfed, 0);

    @test h.helicsEndpointGetName(ept) == "EV_Controller/EV6"

    filt = h.helicsFederateGetFilterByIndex(cfed, 3)

    @test h.helicsFilterGetName(filt) == "EV_Controller/filterEV3"

    f = h.helicsFederateGetFilter(cfed, "EV_Controller/filterEV3")
    @test h.helicsFilterGetName(f) == "EV_Controller/filterEV3"

    ipt = h.helicsFederateGetInputByIndex(cfed, 4)
    @test h.helicsSubscriptionGetKey(ipt) == "IEEE_123_feeder_0/charge_EV3"

    h.helicsFederateDestroy(cfed)
    h.helicsBrokerDestroy(broker)

    rm("testconfig.json")

end


text:
function simple_map(pa :: AnArea_Regular)

  slat = minimum(pa.lat )
  nlat = maximum(pa.lat)

  llon = minimum(pa.lon )
  rlon = maximum(pa.lon )
  if (sum( (pa.lon .> 355) & (pa.lon .< 5)) > 0)
    if (sum( (pa.lon .< 185) & (pa.lon .> 175) ) == 0)
      (llon, rlon) = (rlon, llon)
    end
  end
  m = basemap.Basemap(projection="cyl",
  llcrnrlat=slat ,urcrnrlat=nlat ,llcrnrlon=llon,urcrnrlon=rlon,resolution="c")
end
export simple_map

function DrawMap(m; NoLands :: Bool = false)
  if ( NoLands)
    m[:drawcoastlines](color="#555555")
  else
    m[:fillcontinents](color="#555555")
  end
  m[:drawmeridians](0:30:360.0,labels=[0,0,0,1],fontsize=10);
  m[:drawparallels](-90:10.0:90,labels=[1,0,0,0],fontsize=10);
end
export DrawMap

function lonrotate(plon :: Array{Float64,2}, rr)
  if ((plon[1,1] > 180) & (plon[end,1] < 180))
    lon = mean(plon ,2)
    return [rr[lon .<= 180,:] ; rr[lon .> 180,:]]
  else
    return rr
  end
end

function canvas( pa :: AnArea_Regular,xarr :: Array{Float64,2})
  return map(x->lonrotate(pa.plon,x),(pa.plon, pa.plat, xarr))
end

function canvas( pa :: AnArea_Regular,xarr :: Array{Bool,2})
  (plon,plat,xx) = map( x-> lonrotate(plon,x), (pa.plon, pa.plat, xarr))
  return map(plon[xx], plat[xx])
end

function canvas( pa :: AnArea_Regular,xarr :: BitArray{2})
  (plon,plat,xx) = map( x-> lonrotate(plon,x), (pa.plon, pa.plat, xarr))
  return map(plon[xx], plat[xx])
end
function canvas(x :: Array{Float64,1}, y :: Array{Float64,1},  xarr :: Array{Float64,2})
  (px, py) = begin
    p1 = x .+ zeros(y')
    p2 = y .+ zeros(x')
    (p1, p2')
  end
  return (px, py,xarr)
end

function canvas(x :: Array{Float64,1}, y :: Array{Float64,1},  xarr :: Array{Float64})
  xarr1 = squeeze(xarr, tuple(findin(size(xarr),1)...))
  @assert ( ndims(xarr1) == 2)
  return canvas(x, y, xarr1)
end
export canvas

include("levels.jl")


text:
using Gadfly

set_default_plot_size(6inch, 3inch)

x, y, size = rand(10), rand(10), rand(["foo","bar","pooh"], 10)
levels = sort(unique(size))
p1 = plot(x=x, y=y, size=size, Geom.point, Scale.size_discrete(levels=levels))
p2 = plot(x=x, y=y, size=size, Geom.point, Scale.size_discrete2(levels=levels),
    Theme(key_swatch_shape=Shape.circle))

hstack(p1, p2)


text:
using Test
using Distributed, .Threads

if nprocs() == 1 && nthreads() == 1
    # Lazy stuff should also test distributed
    println("Adding a proc")
    addprocs(1)
end

@testset "basics" begin include("basics.jl"); end
@testset "taxi" begin include("taxi.jl"); end


text:
####
#### length calculations
####

"""
$(TYPEDEF)

Counts for contiguous indexes. Used for combinatoric calculations.
"""
struct Counts
    indexes::UnitRange{Int}
    counts::Vector{Int}
end

Base.pairs(C::Counts) = zip(C.indexes, C.counts)

Base.firstindex(C::Counts) = first(C.indexes)

Base.lastindex(C::Counts) = last(C.indexes)

Base.sum(C::Counts) = sum(C.counts)

"""
$(SIGNATURES)

Calculate the counts
```math
C[k] = ∑ C₁[i] ⋅ C₂[j]  \text{where} i + j = k ≤ cap
```
"""
function count_combinations(cap, C1::Counts, C2::Counts)
    i_first = firstindex(C1) + firstindex(C2)
    i_last = min(cap, lastindex(C1) + lastindex(C2))
    c = zeros(Int, i_last - i_first + 1)
    for (i1, c1) in pairs(C1)
        for (i2, c2) in pairs(C2)
            i = i1 + i2
            if i ≤ i_last
                c[i - i_first + 1] += c1 * c2
            else
                break           # no point in continuing inner loop
            end
        end
    end
    Counts(i_first:i_last, c)
end

####
#### capped cartesian indices iterator
####

"""
$(TYPEDEF)

An iterator over the cartesian indexes `1 ≤ i[1] ≤ I[1]`, …, in colum major order, with
`sum(i) ≤ cap`. Iterator returns N-tuples of integers. Useful for traversing Smolyak
combinations.
"""
struct CappedCartesianIndices{N}
    cap::Int
    I::NTuple{N,Int}
    function CappedCartesianIndices(cap::Int, I::NTuple{N,Int}) where {N}
        @argcheck all(I .≥ 1)
        @argcheck cap ≥ N
        new{N}(cap, I)
    end
end

Base.IteratorSize(::Type{<:CappedCartesianIndices}) = Base.HasLength()

function Base.length(iter::CappedCartesianIndices)
    @unpack cap, I = iter
    C = mapfoldl((C1, C2) -> count_combinations(cap, C1, C2), I) do i
        i = min(cap, i)
        Counts(1:i, ones(Int, i))
    end
    sum(C)
end

Base.IteratorEltype(::Type{<:CappedCartesianIndices}) = Base.HasEltype()

Base.eltype(::Type{CappedCartesianIndices{N}}) where {N} = NTuple{N,Int}

@inline function Base.iterate(iter::CappedCartesianIndices{N}) where {N}
    ι = ntuple(_ -> 1, Val{N}())
    ι, ι
end

@inline function Base.iterate(iter::CappedCartesianIndices, ι)
    @unpack cap, I = iter
    valid, ι = __inc(cap, sum(ι), I, ι)
    valid ? (ι, ι) : nothing
end

@inline __inc(::Int, ::Int, ::Tuple{}, ::Tuple{}) = false, ()

@inline function __inc(cap::Int, ∑ι::Int, I::Tuple{Int}, ι::Tuple{Int})
    ι1 = first(ι)
    (ι1 < first(I) && ∑ι < cap), (ι1 + 1, )
end

@inline function __inc(cap, ∑ι, I, ι)
    ι1 = first(ι)
    if (ι1 < first(I) && ∑ι < cap)
        true, (ι1 + 1, Base.tail(ι)...)
    else
        Δ1 = 1 - ι1
        valid, ιtail = __inc(cap - 1, ∑ι + Δ1 - 1, Base.tail(I), Base.tail(ι))
        valid, (1, ιtail...)
    end
end

"""
$(SIGNATURES)

Largest indices for each axis that occur in the iteration.
"""
function largest_indices(iter::CappedCartesianIndices{N}) where N
    @unpack cap, I = iter
    min.(cap - N + 1, I)
end


text:
# ---
# title: 897. Increasing Order Search Tree
# id: problem897
# author: Indigo
# date: 2021-02-01
# difficulty: Easy
# categories: Tree, Depth-first Search, Recursion
# link: <https://leetcode.com/problems/increasing-order-search-tree/description/>
# hidden: true
# ---
#
# Given the `root` of a binary search tree, rearrange the tree in **in-order**
# so that the leftmost node in the tree is now the root of the tree, and every
# node has no left child and only one right child.
#
#
#
# **Example 1:**
#
# ![](https://assets.leetcode.com/uploads/2020/11/17/ex1.jpg)
#
#
#
#     Input: root = [5,3,6,2,4,null,8,1,null,null,null,7,9]
#     Output: [1,null,2,null,3,null,4,null,5,null,6,null,7,null,8,null,9]
#
#
# **Example 2:**
#
# ![](https://assets.leetcode.com/uploads/2020/11/17/ex2.jpg)
#
#
#
#     Input: root = [5,1,7]
#     Output: [1,null,5,null,7]
#
#
#
#
# **Constraints:**
#
#   * The number of nodes in the given tree will be in the range `[1, 100]`.
#   * `0 <= Node.val <= 1000`
#
#
## @lc code=start
using LeetCode

function increasing_bst(root::TreeNode{Int})
    p = TreeNode(0)
    q = p
    function mid_tranv(root::TreeNode{Int})
        (root.left !== nothing) && mid_tranv(root.left)
        q.right = TreeNode(root.val)
        q = q.right
        (root.right !== nothing) && mid_tranv(root.right)
        nothing
    end
    mid_tranv(root)
    p.right
end
## @lc code=end


text:
@testset "892.surface-area-of-3d-shapes.jl" begin
    @test surface_area([[1,2],[3,4]]) == 34
    @test surface_area([[1,0],[0,2]]) == 16
    @test surface_area([[1,1,1],[1,0,1],[1,1,1]]) == 32
    @test surface_area([[2,2,2],[2,1,2],[2,2,2]]) == 46
    @test surface_area([[2]]) == 10
    @test surface_area(hcat([2,2,2],[2,1,2],[2,2,2])) == 46
end


text:
"""
Computes a global solution for a model via backward time iteration.
The time iteration is  applied directly to the decision rule of the model.

If the initial guess for the decision rule is not explicitly provided, the initial guess is provided by `ConstantDecisionRule`.
If the stochastic process for the model is not explicitly provided, the process is taken from the default provided by the model object, `model.exogenous`.

# Arguments
* `model::NumericModel`: Model object that describes the current model environment.
* `process`: The stochastic process associated with the exogenous variables in the model.
* `init_dr`: Initial guess for the decision rule.
# Returns
* `dr`: Solved decision rule.
"""
function time_iteration_direct(model, dprocess::AbstractDiscretizedProcess,
                               grid, init_dr::AbstractDecisionRule;
                               verbose::Bool=true, maxit::Int=100, trace::Bool=false,
                               tol_η::Float64=1e-7)

    # Grid
    endo_nodes = nodes(ListOfPoints,grid)
    N = length(endo_nodes)

    # Discretized exogenous process
    number_of_smooth_drs(dprocess) = max(n_nodes(dprocess), 1)
    n_m = nsd = number_of_smooth_drs(dprocess)

    p = SVector(model.calibration[:parameters]...)

    # initial guess for controls
    x0 = [init_dr(i, endo_nodes) for i=1:nsd]

    ti_trace = trace ? IterationTrace([x0]) : nothing

    n_x = length(model.calibration[:controls])
    complementarities = true
    if complementarities == true
        x_lb = [controls_lb(model, node(Point,dprocess,i),endo_nodes,p) for i=1:n_m]
        x_ub = [controls_ub(model, node(Point,dprocess,i),endo_nodes,p) for i=1:n_m]
        BIG = 100000
        for i=1:n_m
          for n=1:N
            x_lb[i][n] = max.(x_lb[i][n],-BIG)
            x_ub[i][n] = min.(x_ub[i][n], BIG)
          end
        end
    end

    # create decision rule (which interpolates x0)
    dr = CachedDecisionRule(dprocess, grid, x0)

    # Define controls of tomorrow
    x1 = deepcopy(x0)

    # define states of today
    s = deepcopy(endo_nodes);

    # loop option
    it = 0
    err_0 = NaN
    err = 1.0

    log = TimeIterationLog()
    initialize(log, verbose=verbose)
    append!(log; verbose=verbose, it=0, err=NaN, gain=NaN, epsilon=NaN, time=0.0, nit=NaN)

    ###############################   Iteration loop
    n_h = length(model.symbols[:expectations])

    while it<maxit && err>tol_η

        it += 1

        t1 = time_ns()

        E_f = [zeros(Point{n_h},N) for i=1:number_of_smooth_drs(dprocess)]

        set_values!(dr, x0)
        # Compute expectations function E_f and states of tomorrow

        # S = zeros(size(s))

        for i in 1:size(E_f, 1)
            m = node(Point,dprocess, i)
            for (w, M, j) in get_integration_nodes(Point,dprocess,i)
                # Update the states
                # S[:,:] = Dolo.transition(model, m, s, x0[i], M, p)
                S = Dolo.transition(model, m, s, x0[i], M, p)
                # interpolate controles conditional states of tomorrow
                X = dr(i, j, S)
                # Compute expectations as a weighted average of the exo states w_j
                E_f[i] += w*Dolo.expectation(model, M, S, X, p)

            end
            # compute controles of tomorrow
            x1[i][:] = Dolo.direct_response(model, m, s, E_f[i], p)
        end

        err = 0.0
        for i in 1:size(x1, 1)
            # apply bounds # TODO: improve
            x1[i] = clamp.(x1[i], x_lb[i], x_ub[i])
            # update error
            err = max(err, maxabs(x1[i] - x0[i]))
            # copy controls back into x0
            copyto!(x0[i], x1[i])
        end

        gain = err/err_0
        err_0 = err

        elapsed = time_ns()-t1

        append!(log; verbose=verbose, it=it, err=err, gain=gain, time=elapsed, epsilon=NaN, nit=NaN)

    end

    finalize(log, verbose=verbose)

    converged = err<tol_η
    TimeIterationResult(dr.dr, it, true, dprocess, converged, tol_η, err, log, ti_trace)

end

# get grid for endogenous
function time_iteration_direct(model, dprocess, init_dr::AbstractDecisionRule; grid=Dict(), kwargs...)
    grid = get_grid(model, options=grid)
    return time_iteration_direct(model, dprocess, grid, init_dr;  kwargs...)
end

# get stupid initial rule
function time_iteration_direct(model, dprocess::AbstractDiscretizedProcess; grid=Dict(), kwargs...)

    init_dr = ConstantDecisionRule(model.calibration[:controls])
    return time_iteration_direct(model, dprocess, init_dr;  grid=grid, kwargs...)
end

function time_iteration_direct(model, init_dr::AbstractDecisionRule; grid=Dict(), kwargs...)
    dprocess = discretize( model.exogenous )
    return time_iteration_direct(model, dprocess, init_dr; grid=grid, kwargs...)
end

function time_iteration_direct(model; grid=Dict(), kwargs...)
    dprocess = discretize( model.exogenous )
    init_dr = ConstantDecisionRule(model.calibration[:controls])
    return time_iteration_direct(model, dprocess, init_dr; grid=grid, kwargs...)
end


text:
# ---
# title: GaussianBlur
# cover: gaussianblur.gif
# description: blur the input image using a gaussian kernel
# ---

# [`GaussianBlur`](@ref) can be used to blur the input image using a gaussian
# kernel with a specified kernel size and standard deviation.

using Augmentor
using ImageShow, ImageCore

img_in = testpattern(RGB, ratio=0.5)

mosaicview(
    img_in,
    augment(img_in, GaussianBlur(3)),
    augment(img_in, GaussianBlur(5, 2.5));
    fillvalue=colorant"white", nrow=1, npad=10
)

# ## References

#md # ```@docs
#md # GaussianBlur
#md # ```

## save covers #src
using ImageMagick #src
using FileIO #src
include(joinpath("..", "assets", "utilities.jl")) #src
cover = make_gif(testpattern(RGB, ratio=0.5), GaussianBlur(5), 2) #src
ImageMagick.save("gaussianblur.gif", cover; fps=1) #src


text:
using ForwardDiff
using ForwardDiff: Dual

seed(x::Real, ::Val) = Dual(x, true)

function seed(x, ::Val{N}, offset = 0) where N
  map(x, reshape(1:length(x), size(x))) do x, i
    Dual(x, ntuple(j -> j+offset == i, Val(N)))
  end
end

extract(x::ForwardDiff.Dual) = x.value, [x.partials...]

function extract(xs::AbstractArray{ForwardDiff.Dual{T,V,N}}) where {T,V,N}
  J = similar(xs, V, N, length(xs))
  for i = 1:length(xs), j = 1:N
    J[j, i] = xs[i].partials.values[j]
  end
  return map(x -> x.value, xs), J
end

function forward_jacobian(f, x, ::Val{N}) where N
  y, _J = extract(f(seed(x, Val(N))))
  J = similar(_J, length(x), length(y))
  J[1:N,:] = _J
  offset = 0
  while offset + N < length(x)
    offset += N
    _, _J = extract(f(seed(x, Val(N), offset)))
    range = (1+offset):min(N+offset,length(x))
    J[range,:] = @view _J[range.-offset,:]
  end
  return y, J
end

function forward_jacobian(f, x)
  if length(x) < ForwardDiff.DEFAULT_CHUNK_THRESHOLD
    forward_jacobian(f, x, Val(length(x)))
  else
    forward_jacobian(f, x, Val(ForwardDiff.DEFAULT_CHUNK_THRESHOLD))
  end
end

vec_scalar(x) = vec(x)
vec_scalar(x::Real) = [x]
reshape_scalar(x, y) = reshape(y, size(x))
reshape_scalar(x::Real, y) = y[]

forwarddiff(f, x) = f(x)

@adjoint function forwarddiff(f, x)
  y, J = forward_jacobian(f, x)
  return y, ȳ -> (nothing, reshape_scalar(x, J*vec_scalar(ȳ)))
end


text:
# This script build the API docs and should be run
# only if both Nervoso and Lexicon are installed
using Nervoso, Lexicon

save("temp.md", Nervoso, Config(
   md_subheader=:skip,
))

open("temp.md") do temp
    d = readall(temp)

    # Remove docstring signatures
    d = replace(d, r"\n+\s*`[^`]*`\s*\n+", "\n\n")

    # Change objname from header to code
    header_to_code(s) = string("`", s[6:end-2], "`[")
    d = replace(d, r"#### .* \[", header_to_code)

    # Remove "Nervoso." from readability
    d = replace(d, r"Nervoso\.", "")

    # Remove source links
    d = replace(d, r"\n\*source:\*\n.*\n", "")

    # Add newline after permalinks
    d = replace(d, r"\[¶\](.*)", s -> string(s,"\n"))

    # Change title
    d = replace(d, r"# Nervoso", "# API Reference")

    # Save
    open(joinpath("reference", "api.md"), "w") do f
        write(f, d)
    end

    # Remove temp
    rm("temp.md")
end


text:
"""
    plot_mesh(femesh)

Plot finite element mesh.
"""
function plot_mesh(femesh::FEMesh)
    ncompartment, nboundary = size(femesh.facets)
    scene = nothing
    first = true
    for icmpt = 1:ncompartment, iboundary = 1:nboundary
        facets = femesh.facets[icmpt, iboundary]
        points = femesh.points[icmpt]
        color = points[3, :]
        # color = fill(1.0 * iboundary, length(colors))
        # color = iboundary / nboundary
        if first
            # scene = mesh(points', facets'; color, shading = false)
            scene = poly(points', facets'; color, strokewidth = 1, shading = false)
            first = false
        else
            # mesh!(points', facets'; color, shading = false)
            poly!(points', facets'; color, strokewidth = 1, shading = false)
        end
    end
    scene
end


text:
include("double_integrator.jl")

## Test Euler
@autodiff struct DoubleIntegrator{D} <: DI{D} end
cmodel = DoubleIntegrator{2}()
n,m = RD.dims(cmodel)

# Evaluate Euler step explicitly
x0,u0 = randn(cmodel)
t,h = 0.0, 0.1
xn = zeros(n)
z0 = KnotPoint(x0,u0,t,h)
xdot = RD.dynamics(cmodel, z0)
x2 = x0 + h * xdot

# Create discretized model
model = RD.DiscretizedDynamics{RD.Euler}(cmodel)
@test x2 ≈ RD.discrete_dynamics(model, z0)
RD.discrete_dynamics!(model, xn, z0)
@test x2 ≈ xn

test_model(model)

## Test RK4
x,u = x0, u0
k1 = RD.dynamics(cmodel, x,        u, t      )*h
k2 = RD.dynamics(cmodel, x + k1/2, u, t + h/2)*h
k3 = RD.dynamics(cmodel, x + k2/2, u, t + h/2)*h
k4 = RD.dynamics(cmodel, x + k3,   u, t + h  )*h
x_rk4 = x + (k1 + 2k2 + 2k3 + k4)/6

model = RD.DiscretizedDynamics{RD.RK4}(cmodel)
@test x_rk4 ≈ RD.discrete_dynamics(model, z0)
RD.discrete_dynamics!(model, xn, z0)
@test x_rk4 ≈ xn

test_model(model)
test_error_allocs(model)

## Implicit Midpoint
model = RD.DiscretizedDynamics{RD.ImplicitMidpoint}(cmodel)
z1 = KnotPoint(randn(model)..., t, h)
z2 = KnotPoint(randn(model)..., t+h, h)
x1,u1 = RD.state(z1), RD.control(z1)
x2,u2 = RD.state(z2), RD.control(z2)
y1,y2 = zeros(n), zeros(n)
J1,J2 = zeros(n,n+m), zeros(n,n+m)
J10,J20 = zeros(n,n+m), zeros(n,n+m)

for i = 1:2
    if i == 2
        z1 = KnotPoint{n,m}(Vector(z1.z), t,h)
        z2 = KnotPoint{n,m}(Vector(z2.z), t+h,h)
    end

    e_mid = x1 - x2 + h*RD.dynamics(cmodel, (x1 + x2)/2, u1, 0.0)
    @test e_mid ≈ RD.dynamics_error(model, z2, z1)
    RD.dynamics_error!(model, y2, y1, z2, z1)
    @test e_mid ≈ y2

    RD.dynamics_error_jacobian!(StaticReturn(), RD.UserDefined(), model, J20, J10, y2, y1, z2, z1)
    RD.dynamics_error_jacobian!(StaticReturn(), ForwardAD(), model, J2, J1, y2, y1, z2, z1)
    @test J2 ≈ J20
    @test J1 ≈ J10

    RD.dynamics_error_jacobian!(InPlace(), RD.UserDefined(), model, J20, J10, y2, y1, z2, z1)
    RD.dynamics_error_jacobian!(InPlace(), ForwardAD(), model, J2, J1, y2, y1, z2, z1)
    @test J2 ≈ J20
    @test J1 ≈ J10

    RD.dynamics_error_jacobian!(StaticReturn(), FiniteDifference(), model, J2, J1, y2, y1, z2, z1)
    @test J2 ≈ J20
    @test J1 ≈ J10
    RD.dynamics_error_jacobian!(InPlace(), FiniteDifference(), model, J2, J1, y2, y1, z2, z1)
    @test J2 ≈ J20
    @test J1 ≈ J10
end

test_error_allocs(model)


text:
function est_br(x::AbstractVector, fs=10)
    # Returns breathing rate in bpm
    bounds = Int.((1.2 * fs, 30 * fs))
    s = autocor(x, bounds[1]:bounds[2]-1)
    delay = argmax(s) + bounds[1]
    return 60 * fs / delay
end

function est_br_fft(x::AbstractVector, fs=10)
    # Returns breathing rate in bpm
    n = length(x)
    s = abs2.(rfft(x))
    f = rfftfreq(n, fs)
    return f[argmax(s)] * 60
end

function est_br_fft2(x::AbstractVector, fs=10)
    # using periodogram, not accurate
    s = periodogram(x, fs=fs)
    return freq(s)[argmax(power(s))] * 60
end


text:
using FermiCG
using Printf
using Test
using LinearAlgebra
using Profile
using HDF5
using Random
using PyCall

@testset "CompressedTuckerState" begin
    atoms = []

    r = 1
    a = 1
    push!(atoms,Atom(1,"H", [0, 0*a, 0*r]))
    push!(atoms,Atom(2,"H", [0, 0*a, 1*r]))
    push!(atoms,Atom(3,"H", [0, 1*a, 2*r]))
    push!(atoms,Atom(4,"H", [0, 1*a, 3*r]))
    push!(atoms,Atom(5,"H", [0, 2*a, 4*r]))
    push!(atoms,Atom(6,"H", [0, 2*a, 5*r]))
    push!(atoms,Atom(7,"H", [0, 3*a, 6*r]))
    push!(atoms,Atom(8,"H", [0, 3*a, 7*r]))
    push!(atoms,Atom(9,"H", [0, 4*a, 8*r]))
    push!(atoms,Atom(10,"H",[0, 4*a, 9*r]))
    push!(atoms,Atom(11,"H",[0, 5*a, 10*r]))
    push!(atoms,Atom(12,"H",[0, 5*a, 11*r]))

    clusters    = [(1:4),(5:8),(9:12)]
    init_fspace = [(2,2),(2,2),(2,2)]
    na = 6
    nb = 6

    basis = "sto-3g"
    mol     = Molecule(0,1,atoms,basis)

    # get integrals
    mf = FermiCG.pyscf_do_scf(mol)
    nbas = size(mf.mo_coeff)[1]
    ints = FermiCG.pyscf_build_ints(mol,mf.mo_coeff, zeros(nbas,nbas));
    #e_fci, d1_fci, d2_fci = FermiCG.pyscf_fci(ints, na, nb, conv_tol=1e-10,max_cycle=100, nroots=1)
    e_fci = -18.33022092

    # localize orbitals
    C = mf.mo_coeff
    Cl = FermiCG.localize(mf.mo_coeff,"lowdin",mf)
    FermiCG.pyscf_write_molden(mol,Cl,filename="lowdin.molden")
    S = FermiCG.get_ovlp(mf)
    U =  C' * S * Cl
    println(" Rotate Integrals")
    flush(stdout)
    ints = FermiCG.orbital_rotation(ints,U)
    println(" done.")
    flush(stdout)

    #
    # define clusters
    clusters = [Cluster(i,collect(clusters[i])) for i = 1:length(clusters)]
    display(clusters)

    rdm1 = zeros(size(ints.h1))

    e_cmf, U, Da, Db  = FermiCG.cmf_oo(ints, clusters, init_fspace, rdm1,
                                       max_iter_oo=40, verbose=0, gconv=1e-6, method="bfgs")
    FermiCG.pyscf_write_molden(mol,Cl*U,filename="cmf.molden")
    ints = FermiCG.orbital_rotation(ints,U)

    e_ref = e_cmf - ints.h0

    max_roots = 100
    # build Hamiltonian, cluster_basis and cluster ops
    #display(Da)
    #cluster_bases = FermiCG.compute_cluster_eigenbasis(ints, clusters, verbose=2, max_roots=max_roots)
    cluster_bases = FermiCG.compute_cluster_eigenbasis(ints, clusters, verbose=0, max_roots=max_roots,
                                                       init_fspace=init_fspace, rdm1a=Da, rdm1b=Db)
    clustered_ham = FermiCG.extract_ClusteredTerms(ints, clusters)
    cluster_ops = FermiCG.compute_cluster_ops(cluster_bases, ints);

    p_spaces = Vector{FermiCG.ClusterSubspace}()
    q_spaces = Vector{FermiCG.ClusterSubspace}()

    # define p spaces
    for ci in clusters
        tss = FermiCG.ClusterSubspace(ci)
        tss[init_fspace[ci.idx]] = 1:1
        push!(p_spaces, tss)
    end

    # define q spaces
    for tssp in p_spaces
        tss = FermiCG.get_ortho_compliment(tssp, cluster_bases[tssp.cluster.idx])
        push!(q_spaces, tss)
    end

    println(" ================= Cluster P Spaces ===================")
    display.(p_spaces)
    println(" ================= Cluster Q Spaces ===================")
    display.(q_spaces)

    FermiCG.add_cmf_operators!(cluster_ops, cluster_bases, ints, Da, Db);

    nroots = 1

    #
    # initialize with eye
    ref_vector = FermiCG.TuckerState(clusters, p_spaces, q_spaces, na, nb)
    FermiCG.set_vector!(ref_vector, Matrix(1.0I, length(ref_vector),nroots))

    ref_vec  = FermiCG.CompressedTuckerState(ref_vector, thresh=-1);

    e_var, v_var = FermiCG.block_sparse_tucker(ref_vec, cluster_ops, clustered_ham, nbody=2,
                                                      thresh_foi=1e-7,
                                                      thresh_pt =1e-5,
                                                      thresh_var=1e-4,
                                                      tol_ci=1e-10, tol_tucker=1e-4, do_pt=true)

    @test isapprox(e_var[1], -18.33000292416142, atol=1e-8)

    e_cepa, v_cepa = FermiCG.do_fois_cepa(ref_vec, cluster_ops, clustered_ham, thresh_foi=1e-8, max_iter=50, tol=1e-8)
    @test isapprox(e_cepa[1], -18.330092885975663, atol=1e-8)

    #e_cepa, v_cepa = FermiCG.do_fois_cepa(v_var, cluster_ops, clustered_ham, thresh_foi=1e-10, max_iter=50, tol=1e-8)
    #@test isapprox(e_cepa[1], -18.330092885975663, atol=1e-8)

    e_pt, v_pt = FermiCG.do_fois_pt2(ref_vec, cluster_ops, clustered_ham, thresh_foi=1e-8, max_iter=50, tol=1e-8)
    @test isapprox(e_pt, -18.32863783512617, atol=1e-8)

    e_ci, v_ci = FermiCG.tucker_ci_solve(v_cepa, cluster_ops, clustered_ham)
    @test isapprox(e_ci[1], -18.330044872505518, atol=1e-8)

    #v_pt, e_pt = FermiCG.hylleraas_compressed_mp2(v_cepa, ref_vec, cluster_ops, clustered_ham; tol=1e-8, max_iter=100, H0="Hcmf", verbose=1)
    #@test isapprox(e_pt, -18.32863783512617, atol=1e-8)
end

max_stars_repo_name: nmayhall-vt/FermiCG

max_issues_repo_name: nmayhall-vt/FermiCG

max_forks_repo_name: nmayhall-vt/FermiCG


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     cell_metadata_json: true
#     formats: ipynb,jl:light
#     text_representation:
#       extension: .jl
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.11.3
#   kernelspec:
#     display_name: Julia 1.6.2
#     language: julia
#     name: julia-1.6
# ---

# + [markdown] {"slideshow": {"slide_type": "slide"}}
# # Ocean Transports -- Integration And Mapping
#
# Transports in the climate system are often represented as gridded vector fields (e.g. on a `C-grid`) and integrated across model grid lines. These paths may for examples connect a pair of coordinates, or track latitude circles. For more detail, please refer to [Forget et al, 2015](https://doi.org/10.5194/gmd-8-3071-2015) and appendices.

# + [markdown] {"slideshow": {"slide_type": "slide"}}
# ### Read grid & transports from file
#
# 1. pre-requisites
# 2. read transports

# + {"slideshow": {"slide_type": "subslide"}}
include("helper_functions.jl")

pth=MeshArrays.GRID_LLC90
γ=GridSpec("LatLonCap",pth)
Γ=GridLoad(γ;option="full")
(Tx,Ty,τx,τy,η)=trsp_read("LatLonCap",pth);
# -

# Prepare for visualization:

# +
μ =Γ.hFacC[:,1]
μ[findall(μ.>0.0)].=1.0
μ[findall(μ.==0.0)].=NaN

lon=[i for i=-179.:2.0:179., j=-89.:2.0:89.]
lat=[j for i=-179.:2.0:179., j=-89.:2.0:89.]

#(f,i,j,w)=InterpolationFactors(Γ,vec(lon),vec(lat))
#λ=(lon=lon,lat=lat,f=f,i=i,j=j,w=w);
#df = DataFrame(f=λ.f[:], i=λ.i[:], j=λ.j[:], w=Float32.(λ.w[:]));
#CSV.write("interp_coeffs.csv", df)

df=DataFrame(CSV.File("interp_coeffs.csv"))
λ=(f=reshape(df.f,length(lon[:]),4), i=reshape(df.i,length(lon[:]),4),
    j=reshape(df.j,length(lon[:]),4), w=reshape(df.w,length(lon[:]),4));

# + [markdown] {"slideshow": {"slide_type": "slide"}}
# ### Integrate transport across latitude lines
#
# 1. `LatitudeCircles` computes `grid edge path`s that track latitude circles
# 2. `ThroughFlow` integrates transports accross the specified `grid edge path`s
# 3. Plot integrated meridional transport in `Sverdrup` units (1Sv=10^6 m³/s)

# + {"slideshow": {"slide_type": "subslide"}}
uv=Dict("U"=>Tx,"V"=>Ty,"dimensions"=>["x","y"])
L=-89.0:89.0; LC=LatitudeCircles(L,Γ)

T=Array{Float64,1}(undef,length(LC))
[T[i]=1e-6*ThroughFlow(uv,LC[i],Γ) for i=1:length(LC)]

plot(L,T,xlabel="latitude",ylabel="Northward Transport (in Sv)",label="ECCO estimate")

# + [markdown] {"slideshow": {"slide_type": "slide"}}
# ## Transport Directions
#
# 1. `u,v` are oriented in standard `Eastward,Northward` directions
# 2. `uC,vC` are oriented along the `x,y` directions of each subdomain

# + {"slideshow": {"slide_type": "-"}}
u,v,uC,vC=rotate_uv(uv,Γ);

# + {"slideshow": {"slide_type": "subslide"}, "cell_style": "split"}
heatmap(u,clims=(-20.0,20.0),title="East-ward")
#heatmap(v,clims=(-20.0,20.0),title="North-ward")

# + {"slideshow": {"slide_type": "fragment"}, "cell_style": "split"}
heatmap(uC,clims=(-20.0,20.0),title="x-ward")
#heatmap(vC,clims=(-20.0,20.0),title="y-ward")

# + [markdown] {"slideshow": {"slide_type": "slide"}}
# ### Global Maps
#
# 1. interpolate `u,v` to a `1/2 x 1/2` degree grid for plotting
# 2. map out the Eastward,Northward transport components
# -

uI=reshape(Interpolate(u,λ.f,λ.i,λ.j,λ.w),size(lon))
vI=reshape(Interpolate(v,λ.f,λ.i,λ.j,λ.w),size(lon))
size(lon)

# + {"slideshow": {"slide_type": "subslide"}, "cell_style": "split"}
heatmap(lon[:,1],lat[1,:],
    permutedims(uI),clims=(-20.0,20.0),
    title="Eastward transport (in Sv / cell)")
# + {"slideshow": {"slide_type": "fragment"}, "cell_style": "split"}
heatmap(lon[:,1],lat[1,:],
    permutedims(vI),clims=(-20.0,20.0),
    title="Northward transport (in Sv / cell)")
# -

max_stars_repo_name: gaelforget/Notebooks

max_issues_repo_name: gaelforget/Notebooks

max_forks_repo_name: gaelforget/Notebooks


text:
# Manually added
const PNGCAPI = nothing
const PNG_BYTES_TO_CHECK = 8
const png_FILE_p = Ptr{Cvoid}

# Automatically generated using Clang.jl

# Skipping MacroDefinition: PNGARG ( arglist ) arglist

const PNGCBAPI = PNGCAPI
const PNGAPI = PNGCAPI

# Skipping MacroDefinition: PNG_FUNCTION ( type , name , args , attributes ) attributes type name args
# Skipping MacroDefinition: PNG_EXPORT_TYPE ( type ) PNG_IMPEXP type
# Skipping MacroDefinition: PNG_EXPORTA ( ordinal , type , name , args , attributes ) PNG_FUNCTION ( PNG_EXPORT_TYPE ( type ) , ( PNGAPI name ) , PNGARG ( args ) , PNG_LINKAGE_API attributes )
# Skipping MacroDefinition: PNG_EXPORT ( ordinal , type , name , args ) PNG_EXPORTA ( ordinal , type , name , args , PNG_EMPTY )
# Skipping MacroDefinition: PNG_REMOVED ( ordinal , type , name , args , attributes )
# Skipping MacroDefinition: PNG_CALLBACK ( type , name , args ) type ( PNGCBAPI name ) PNGARG ( args )
# Skipping MacroDefinition: PNG_USE_RESULT __attribute__ ( ( __warn_unused_result__ ) )
# Skipping MacroDefinition: PNG_NORETURN __attribute__ ( ( __noreturn__ ) )
# Skipping MacroDefinition: PNG_ALLOCATED __attribute__ ( ( __malloc__ ) )
# Skipping MacroDefinition: PNG_DEPRECATED __attribute__ ( ( __deprecated__ ) )
# Skipping MacroDefinition: PNG_PRIVATE __attribute__ ( ( __unavailable__ ( "This function is not exported by libpng." ) ) )
# Skipping MacroDefinition: PNG_FP_EXPORT ( ordinal , type , name , args ) PNG_EXPORT ( ordinal , type , name , args ) ;
# Skipping MacroDefinition: PNG_FIXED_EXPORT ( ordinal , type , name , args )

const png_byte = Cuchar
const png_int_16 = Int16
const png_uint_16 = UInt16
const png_int_32 = Cint
const png_uint_32 = UInt32
const png_size_t = Csize_t
const png_ptrdiff_t = Cptrdiff_t
const png_alloc_size_t = Csize_t
const png_fixed_point = png_int_32
const png_voidp = Ptr{Cvoid}
const png_const_voidp = Ptr{Cvoid}
const png_bytep = Ptr{png_byte}
const png_const_bytep = Ptr{png_byte}
const png_uint_32p = Ptr{png_uint_32}
const png_const_uint_32p = Ptr{png_uint_32}
const png_int_32p = Ptr{png_int_32}
const png_const_int_32p = Ptr{png_int_32}
const png_uint_16p = Ptr{png_uint_16}
const png_const_uint_16p = Ptr{png_uint_16}
const png_int_16p = Ptr{png_int_16}
const png_const_int_16p = Ptr{png_int_16}
const png_charp = Cstring
const png_const_charp = Cstring
const png_fixed_point_p = Ptr{png_fixed_point}
const png_const_fixed_point_p = Ptr{png_fixed_point}
const png_size_tp = Ptr{Csize_t}
const png_const_size_tp = Ptr{Csize_t}
const png_doublep = Ptr{Cdouble}
const png_const_doublep = Ptr{Cdouble}
const png_bytepp = Ptr{Ptr{png_byte}}
const png_uint_32pp = Ptr{Ptr{png_uint_32}}
const png_int_32pp = Ptr{Ptr{png_int_32}}
const png_uint_16pp = Ptr{Ptr{png_uint_16}}
const png_int_16pp = Ptr{Ptr{png_int_16}}
const png_const_charpp = Ptr{Cstring}
const png_charpp = Ptr{Cstring}
const png_fixed_point_pp = Ptr{Ptr{png_fixed_point}}
const png_doublepp = Ptr{Ptr{Cdouble}}
const png_charppp = Ptr{Ptr{Cstring}}
const PNG_API_RULE = 0
const PNG_DEFAULT_READ_MACROS = 1
const PNG_GAMMA_THRESHOLD_FIXED = 5000
const PNG_ZBUF_SIZE = 8192
const PNG_IDAT_READ_SIZE = PNG_ZBUF_SIZE
const PNG_INFLATE_BUF_SIZE = 1024
const PNG_MAX_GAMMA_8 = 11
const PNG_QUANTIZE_BLUE_BITS = 5
const PNG_QUANTIZE_GREEN_BITS = 5
const PNG_QUANTIZE_RED_BITS = 5
const PNG_TEXT_Z_DEFAULT_COMPRESSION = -1
const PNG_TEXT_Z_DEFAULT_STRATEGY = 0
const PNG_USER_CHUNK_CACHE_MAX = 1000
const PNG_USER_CHUNK_MALLOC_MAX = 8000000
const PNG_USER_HEIGHT_MAX = 1000000
const PNG_USER_WIDTH_MAX = 1000000
const PNG_ZLIB_VERNUM = 0x1250
const PNG_Z_DEFAULT_COMPRESSION = -1
const PNG_Z_DEFAULT_NOFILTER_STRATEGY = 0
const PNG_Z_DEFAULT_STRATEGY = 1
const PNG_sCAL_PRECISION = 5
const PNG_sRGB_PROFILE_CHECKS = 2
const PNG_LIBPNG_VER_STRING = "1.6.37"
const PNG_HEADER_VERSION_STRING = " libpng version 1.6.37 - April 14, 2019\n"
const PNG_LIBPNG_VER_SONUM = 16
const PNG_LIBPNG_VER_DLLNUM = 16
const PNG_LIBPNG_VER_MAJOR = 1
const PNG_LIBPNG_VER_MINOR = 6
const PNG_LIBPNG_VER_RELEASE = 37
const PNG_LIBPNG_VER_BUILD = 0
const PNG_LIBPNG_BUILD_ALPHA = 1
const PNG_LIBPNG_BUILD_BETA = 2
const PNG_LIBPNG_BUILD_RC = 3
const PNG_LIBPNG_BUILD_STABLE = 4
const PNG_LIBPNG_BUILD_RELEASE_STATUS_MASK = 7
const PNG_LIBPNG_BUILD_PATCH = 8
const PNG_LIBPNG_BUILD_PRIVATE = 16
const PNG_LIBPNG_BUILD_SPECIAL = 32
const PNG_LIBPNG_BUILD_BASE_TYPE = PNG_LIBPNG_BUILD_STABLE
const PNG_LIBPNG_VER = 10637
const PNG_LIBPNG_BUILD_TYPE = PNG_LIBPNG_BUILD_BASE_TYPE

# Skipping MacroDefinition: png_libpng_ver png_get_header_ver ( NULL )

const PNG_TEXT_COMPRESSION_NONE_WR = -3
const PNG_TEXT_COMPRESSION_zTXt_WR = -2
const PNG_TEXT_COMPRESSION_NONE = -1
const PNG_TEXT_COMPRESSION_zTXt = 0
const PNG_ITXT_COMPRESSION_NONE = 1
const PNG_ITXT_COMPRESSION_zTXt = 2
const PNG_TEXT_COMPRESSION_LAST = 3
const PNG_HAVE_IHDR = 0x01
const PNG_HAVE_PLTE = 0x02
const PNG_AFTER_IDAT = 0x08

# Skipping MacroDefinition: PNG_UINT_31_MAX ( ( png_uint_32 ) 0x7fffffffL )
# Skipping MacroDefinition: PNG_UINT_32_MAX ( ( png_uint_32 ) ( - 1 ) )
# Skipping MacroDefinition: PNG_SIZE_MAX ( ( size_t ) ( - 1 ) )

const PNG_FP_1 = 100000
const PNG_FP_HALF = 50000

# Skipping MacroDefinition: PNG_FP_MAX ( ( png_fixed_point ) 0x7fffffffL )

const PNG_FP_MAX =  typemax(png_fixed_point)
const PNG_FP_MIN = typemin(png_fixed_point)
const PNG_COLOR_MASK_PALETTE = 1
const PNG_COLOR_MASK_COLOR = 2
const PNG_COLOR_MASK_ALPHA = 4
const PNG_COLOR_TYPE_GRAY = 0
const PNG_COLOR_TYPE_PALETTE = PNG_COLOR_MASK_COLOR | PNG_COLOR_MASK_PALETTE
const PNG_COLOR_TYPE_RGB = PNG_COLOR_MASK_COLOR
const PNG_COLOR_TYPE_RGB_ALPHA = PNG_COLOR_MASK_COLOR | PNG_COLOR_MASK_ALPHA
const PNG_COLOR_TYPE_GRAY_ALPHA = PNG_COLOR_MASK_ALPHA
const PNG_COLOR_TYPE_RGBA = PNG_COLOR_TYPE_RGB_ALPHA
const PNG_COLOR_TYPE_GA = PNG_COLOR_TYPE_GRAY_ALPHA
const PNG_COMPRESSION_TYPE_BASE = 0
const PNG_COMPRESSION_TYPE_DEFAULT = PNG_COMPRESSION_TYPE_BASE
const PNG_FILTER_TYPE_BASE = 0
const PNG_INTRAPIXEL_DIFFERENCING = 64
const PNG_FILTER_TYPE_DEFAULT = PNG_FILTER_TYPE_BASE
const PNG_INTERLACE_NONE = 0
const PNG_INTERLACE_ADAM7 = 1
const PNG_INTERLACE_LAST = 2
const PNG_OFFSET_PIXEL = 0
const PNG_OFFSET_MICROMETER = 1
const PNG_OFFSET_LAST = 2
const PNG_EQUATION_LINEAR = 0
const PNG_EQUATION_BASE_E = 1
const PNG_EQUATION_ARBITRARY = 2
const PNG_EQUATION_HYPERBOLIC = 3
const PNG_EQUATION_LAST = 4
const PNG_SCALE_UNKNOWN = 0
const PNG_SCALE_METER = 1
const PNG_SCALE_RADIAN = 2
const PNG_SCALE_LAST = 3
const PNG_RESOLUTION_UNKNOWN = 0
const PNG_RESOLUTION_METER = 1
const PNG_RESOLUTION_LAST = 2
const PNG_sRGB_INTENT_PERCEPTUAL = 0
const PNG_sRGB_INTENT_RELATIVE = 1
const PNG_sRGB_INTENT_SATURATION = 2
const PNG_sRGB_INTENT_ABSOLUTE = 3
const PNG_sRGB_INTENT_LAST = 4
const PNG_KEYWORD_MAX_LENGTH = 79
const PNG_MAX_PALETTE_LENGTH = 256
const PNG_INFO_gAMA = UInt32(0x0001)
const PNG_INFO_sBIT = UInt32(0x0002)
const PNG_INFO_cHRM = UInt32(0x0004)
const PNG_INFO_PLTE = UInt32(0x0008)
const PNG_INFO_tRNS = UInt32(0x0010)
const PNG_INFO_bKGD = UInt32(0x0020)
const PNG_INFO_hIST = UInt32(0x0040)
const PNG_INFO_pHYs = UInt32(0x0080)
const PNG_INFO_oFFs = UInt32(0x0100)
const PNG_INFO_tIME = UInt32(0x0200)
const PNG_INFO_pCAL = UInt32(0x0400)
const PNG_INFO_sRGB = UInt32(0x0800)
const PNG_INFO_iCCP = UInt32(0x1000)
const PNG_INFO_sPLT = UInt32(0x2000)
const PNG_INFO_sCAL = UInt32(0x4000)
const PNG_INFO_IDAT = UInt32(0x8000)
const PNG_INFO_eXIf = UInt32(0x00010000)
const PNG_TRANSFORM_IDENTITY = 0x0000
const PNG_TRANSFORM_STRIP_16 = 0x0001
const PNG_TRANSFORM_STRIP_ALPHA = 0x0002
const PNG_TRANSFORM_PACKING = 0x0004
const PNG_TRANSFORM_PACKSWAP = 0x0008
const PNG_TRANSFORM_EXPAND = 0x0010
const PNG_TRANSFORM_INVERT_MONO = 0x0020
const PNG_TRANSFORM_SHIFT = 0x0040
const PNG_TRANSFORM_BGR = 0x0080
const PNG_TRANSFORM_SWAP_ALPHA = 0x0100
const PNG_TRANSFORM_SWAP_ENDIAN = 0x0200
const PNG_TRANSFORM_INVERT_ALPHA = 0x0400
const PNG_TRANSFORM_STRIP_FILLER = 0x0800
const PNG_TRANSFORM_STRIP_FILLER_BEFORE = PNG_TRANSFORM_STRIP_FILLER
const PNG_TRANSFORM_STRIP_FILLER_AFTER = 0x1000
const PNG_TRANSFORM_GRAY_TO_RGB = 0x2000
const PNG_TRANSFORM_EXPAND_16 = 0x4000
const PNG_TRANSFORM_SCALE_16 = 0x8000
const PNG_FLAG_MNG_EMPTY_PLTE = 0x01
const PNG_FLAG_MNG_FILTER_64 = 0x04
const PNG_ALL_MNG_FEATURES = 0x05

# Skipping MacroDefinition: png_check_sig ( sig , n ) ! png_sig_cmp ( ( sig ) , 0 , ( n ) )
# Skipping MacroDefinition: png_jmpbuf ( png_ptr ) ( * png_set_longjmp_fn ( ( png_ptr ) , longjmp , ( sizeof ( jmp_buf ) ) ) )

const PNG_ERROR_ACTION_NONE = 1
const PNG_ERROR_ACTION_WARN = 2
const PNG_ERROR_ACTION_ERROR = 3
const PNG_RGB_TO_GRAY_DEFAULT = -1
const PNG_ALPHA_PNG = 0
const PNG_ALPHA_STANDARD = 1
const PNG_ALPHA_ASSOCIATED = 1
const PNG_ALPHA_PREMULTIPLIED = 1
const PNG_ALPHA_OPTIMIZED = 2
const PNG_ALPHA_BROKEN = 3
const PNG_DEFAULT_sRGB = -1
const PNG_GAMMA_MAC_18 = -2
const PNG_GAMMA_sRGB = 220000
const PNG_GAMMA_LINEAR = PNG_FP_1
const PNG_FILLER_BEFORE = 0
const PNG_FILLER_AFTER = 1
const PNG_BACKGROUND_GAMMA_UNKNOWN = 0
const PNG_BACKGROUND_GAMMA_SCREEN = 1
const PNG_BACKGROUND_GAMMA_FILE = 2
const PNG_BACKGROUND_GAMMA_UNIQUE = 3
const PNG_GAMMA_THRESHOLD = PNG_GAMMA_THRESHOLD_FIXED * 1.0e-5
const PNG_CRC_DEFAULT = 0
const PNG_CRC_ERROR_QUIT = 1
const PNG_CRC_WARN_DISCARD = 2
const PNG_CRC_WARN_USE = 3
const PNG_CRC_QUIET_USE = 4
const PNG_CRC_NO_CHANGE = 5
const PNG_NO_FILTERS = 0x00
const PNG_FILTER_NONE = 0x08
const PNG_FILTER_SUB = 0x10
const PNG_FILTER_UP = 0x20
const PNG_FILTER_AVG = 0x40
const PNG_FILTER_PAETH = 0x80
const PNG_FAST_FILTERS = (PNG_FILTER_NONE | PNG_FILTER_SUB) | PNG_FILTER_UP
const PNG_ALL_FILTERS = (PNG_FAST_FILTERS | PNG_FILTER_AVG) | PNG_FILTER_PAETH
const PNG_FILTER_VALUE_NONE = 0
const PNG_FILTER_VALUE_SUB = 1
const PNG_FILTER_VALUE_UP = 2
const PNG_FILTER_VALUE_AVG = 3
const PNG_FILTER_VALUE_PAETH = 4
const PNG_FILTER_VALUE_LAST = 5
const PNG_FILTER_HEURISTIC_DEFAULT = 0
const PNG_FILTER_HEURISTIC_UNWEIGHTED = 1
const PNG_FILTER_HEURISTIC_WEIGHTED = 2
const PNG_FILTER_HEURISTIC_LAST = 3
const PNG_DESTROY_WILL_FREE_DATA = 1
const PNG_SET_WILL_FREE_DATA = 1
const PNG_USER_WILL_FREE_DATA = 2
const PNG_FREE_HIST = UInt32(0x0008)
const PNG_FREE_ICCP = UInt32(0x0010)
const PNG_FREE_SPLT = UInt32(0x0020)
const PNG_FREE_ROWS = UInt32(0x0040)
const PNG_FREE_PCAL = UInt32(0x0080)
const PNG_FREE_SCAL = UInt32(0x0100)
const PNG_FREE_UNKN = UInt32(0x0200)
const PNG_FREE_PLTE = UInt32(0x1000)
const PNG_FREE_TRNS = UInt32(0x2000)
const PNG_FREE_TEXT = UInt32(0x4000)
const PNG_FREE_EXIF = UInt32(0x8000)
const PNG_FREE_ALL = UInt32(0xffff)
const PNG_FREE_MUL = UInt32(0x4220)
const PNG_HANDLE_CHUNK_AS_DEFAULT = 0
const PNG_HANDLE_CHUNK_NEVER = 1
const PNG_HANDLE_CHUNK_IF_SAFE = 2
const PNG_HANDLE_CHUNK_ALWAYS = 3
const PNG_HANDLE_CHUNK_LAST = 4
const PNG_IO_NONE = 0x0000
const PNG_IO_READING = 0x0001
const PNG_IO_WRITING = 0x0002
const PNG_IO_SIGNATURE = 0x0010
const PNG_IO_CHUNK_HDR = 0x0020
const PNG_IO_CHUNK_DATA = 0x0040
const PNG_IO_CHUNK_CRC = 0x0080
const PNG_IO_MASK_OP = Float32(0x0000)
const PNG_IO_MASK_LOC = 0x00f0
const PNG_INTERLACE_ADAM7_PASSES = 7

# Skipping MacroDefinition: PNG_PASS_START_ROW ( pass ) ( ( ( 1 & ~ ( pass ) ) << ( 3 - ( ( pass ) >> 1 ) ) ) & 7 )
# Skipping MacroDefinition: PNG_PASS_START_COL ( pass ) ( ( ( 1 & ( pass ) ) << ( 3 - ( ( ( pass ) + 1 ) >> 1 ) ) ) & 7 )
# Skipping MacroDefinition: PNG_PASS_ROW_OFFSET ( pass ) ( ( pass ) > 2 ? ( 8 >> ( ( ( pass ) - 1 ) >> 1 ) ) : 8 )
# Skipping MacroDefinition: PNG_PASS_COL_OFFSET ( pass ) ( 1 << ( ( 7 - ( pass ) ) >> 1 ) )
# Skipping MacroDefinition: PNG_PASS_ROW_SHIFT ( pass ) ( ( pass ) > 2 ? ( 8 - ( pass ) ) >> 1 : 3 )
# Skipping MacroDefinition: PNG_PASS_COL_SHIFT ( pass ) ( ( pass ) > 1 ? ( 7 - ( pass ) ) >> 1 : 3 )
# Skipping MacroDefinition: PNG_PASS_ROWS ( height , pass ) ( ( ( height ) + ( ( ( 1 << PNG_PASS_ROW_SHIFT ( pass ) ) - 1 ) - PNG_PASS_START_ROW ( pass ) ) ) >> PNG_PASS_ROW_SHIFT ( pass ) )
# Skipping MacroDefinition: PNG_PASS_COLS ( width , pass ) ( ( ( width ) + ( ( ( 1 << PNG_PASS_COL_SHIFT ( pass ) ) - 1 ) - PNG_PASS_START_COL ( pass ) ) ) >> PNG_PASS_COL_SHIFT ( pass ) )
# Skipping MacroDefinition: PNG_ROW_FROM_PASS_ROW ( y_in , pass ) ( ( ( y_in ) << PNG_PASS_ROW_SHIFT ( pass ) ) + PNG_PASS_START_ROW ( pass ) )
# Skipping MacroDefinition: PNG_COL_FROM_PASS_COL ( x_in , pass ) ( ( ( x_in ) << PNG_PASS_COL_SHIFT ( pass ) ) + PNG_PASS_START_COL ( pass ) )
# Skipping MacroDefinition: PNG_PASS_MASK ( pass , off ) ( ( ( 0x110145AF >> ( ( ( 7 - ( off ) ) - ( pass ) ) << 2 ) ) & 0xF ) | ( ( 0x01145AF0 >> ( ( ( 7 - ( off ) ) - ( pass ) ) << 2 ) ) & 0xF0 ) )
# Skipping MacroDefinition: PNG_ROW_IN_INTERLACE_PASS ( y , pass ) ( ( PNG_PASS_MASK ( pass , 0 ) >> ( ( y ) & 7 ) ) & 1 )
# Skipping MacroDefinition: PNG_COL_IN_INTERLACE_PASS ( x , pass ) ( ( PNG_PASS_MASK ( pass , 1 ) >> ( ( x ) & 7 ) ) & 1 )
# Skipping MacroDefinition: png_composite ( composite , fg , alpha , bg ) { png_uint_16 temp = ( png_uint_16 ) ( ( png_uint_16 ) ( fg ) * ( png_uint_16 ) ( alpha ) + ( png_uint_16 ) ( bg ) * ( png_uint_16 ) ( 255 - ( png_uint_16 ) ( alpha ) ) + 128 ) ; ( composite ) = ( png_byte ) ( ( ( temp + ( temp >> 8 ) ) >> 8 ) & 0xff ) ; }
# Skipping MacroDefinition: png_composite_16 ( composite , fg , alpha , bg ) { png_uint_32 temp = ( png_uint_32 ) ( ( png_uint_32 ) ( fg ) * ( png_uint_32 ) ( alpha ) + ( png_uint_32 ) ( bg ) * ( 65535 - ( png_uint_32 ) ( alpha ) ) + 32768 ) ; ( composite ) = ( png_uint_16 ) ( 0xffff & ( ( temp + ( temp >> 16 ) ) >> 16 ) ) ; }
# Skipping MacroDefinition: PNG_get_uint_32 ( buf ) ( ( ( png_uint_32 ) ( * ( buf ) ) << 24 ) + ( ( png_uint_32 ) ( * ( ( buf ) + 1 ) ) << 16 ) + ( ( png_uint_32 ) ( * ( ( buf ) + 2 ) ) << 8 ) + ( ( png_uint_32 ) ( * ( ( buf ) + 3 ) ) ) )
# Skipping MacroDefinition: PNG_get_uint_16 ( buf ) ( ( png_uint_16 ) ( ( ( unsigned int ) ( * ( buf ) ) << 8 ) + ( ( unsigned int ) ( * ( ( buf ) + 1 ) ) ) ) )
# Skipping MacroDefinition: PNG_get_int_32 ( buf ) ( ( png_int_32 ) ( ( * ( buf ) & 0x80 ) ? - ( ( png_int_32 ) ( ( ( png_get_uint_32 ( buf ) ^ 0xffffffffU ) + 1U ) & 0x7fffffffU ) ) : ( png_int_32 ) png_get_uint_32 ( buf ) ) )
# Skipping MacroDefinition: png_get_uint_32 ( buf ) PNG_get_uint_32 ( buf )
# Skipping MacroDefinition: png_get_uint_16 ( buf ) PNG_get_uint_16 ( buf )
# Skipping MacroDefinition: png_get_int_32 ( buf ) PNG_get_int_32 ( buf )

const PNG_IMAGE_VERSION = 1
const PNG_IMAGE_WARNING = 1
const PNG_IMAGE_ERROR = 2

# Skipping MacroDefinition: PNG_IMAGE_FAILED ( png_cntrl ) ( ( ( ( png_cntrl ) . warning_or_error ) & 0x03 ) > 1 )

const PNG_FORMAT_FLAG_ALPHA = UInt32(0x01)
const PNG_FORMAT_FLAG_COLOR = UInt32(0x02)
const PNG_FORMAT_FLAG_LINEAR = UInt32(0x04)
const PNG_FORMAT_FLAG_COLORMAP = UInt32(0x08)
const PNG_FORMAT_FLAG_BGR = UInt32(0x10)
const PNG_FORMAT_FLAG_AFIRST = UInt32(0x20)
const PNG_FORMAT_FLAG_ASSOCIATED_ALPHA = UInt32(0x40)
const PNG_FORMAT_GRAY = 0
const PNG_FORMAT_GA = PNG_FORMAT_FLAG_ALPHA
const PNG_FORMAT_AG = PNG_FORMAT_GA | PNG_FORMAT_FLAG_AFIRST
const PNG_FORMAT_RGB = PNG_FORMAT_FLAG_COLOR
const PNG_FORMAT_BGR = PNG_FORMAT_FLAG_COLOR | PNG_FORMAT_FLAG_BGR
const PNG_FORMAT_RGBA = PNG_FORMAT_RGB | PNG_FORMAT_FLAG_ALPHA
const PNG_FORMAT_ARGB = PNG_FORMAT_RGBA | PNG_FORMAT_FLAG_AFIRST
const PNG_FORMAT_BGRA = PNG_FORMAT_BGR | PNG_FORMAT_FLAG_ALPHA
const PNG_FORMAT_ABGR = PNG_FORMAT_BGRA | PNG_FORMAT_FLAG_AFIRST
const PNG_FORMAT_LINEAR_Y = PNG_FORMAT_FLAG_LINEAR
const PNG_FORMAT_LINEAR_Y_ALPHA = PNG_FORMAT_FLAG_LINEAR | PNG_FORMAT_FLAG_ALPHA
const PNG_FORMAT_LINEAR_RGB = PNG_FORMAT_FLAG_LINEAR | PNG_FORMAT_FLAG_COLOR
const PNG_FORMAT_LINEAR_RGB_ALPHA = (PNG_FORMAT_FLAG_LINEAR | PNG_FORMAT_FLAG_COLOR) | PNG_FORMAT_FLAG_ALPHA
const PNG_FORMAT_RGB_COLORMAP = PNG_FORMAT_RGB | PNG_FORMAT_FLAG_COLORMAP
const PNG_FORMAT_BGR_COLORMAP = PNG_FORMAT_BGR | PNG_FORMAT_FLAG_COLORMAP
const PNG_FORMAT_RGBA_COLORMAP = PNG_FORMAT_RGBA | PNG_FORMAT_FLAG_COLORMAP
const PNG_FORMAT_ARGB_COLORMAP = PNG_FORMAT_ARGB | PNG_FORMAT_FLAG_COLORMAP
const PNG_FORMAT_BGRA_COLORMAP = PNG_FORMAT_BGRA | PNG_FORMAT_FLAG_COLORMAP
const PNG_FORMAT_ABGR_COLORMAP = PNG_FORMAT_ABGR | PNG_FORMAT_FLAG_COLORMAP

# Skipping MacroDefinition: PNG_IMAGE_SAMPLE_CHANNELS ( fmt ) ( ( ( fmt ) & ( PNG_FORMAT_FLAG_COLOR | PNG_FORMAT_FLAG_ALPHA ) ) + 1 )
# Skipping MacroDefinition: PNG_IMAGE_SAMPLE_COMPONENT_SIZE ( fmt ) ( ( ( ( fmt ) & PNG_FORMAT_FLAG_LINEAR ) >> 2 ) + 1 )
# Skipping MacroDefinition: PNG_IMAGE_SAMPLE_SIZE ( fmt ) ( PNG_IMAGE_SAMPLE_CHANNELS ( fmt ) * PNG_IMAGE_SAMPLE_COMPONENT_SIZE ( fmt ) )
# Skipping MacroDefinition: PNG_IMAGE_MAXIMUM_COLORMAP_COMPONENTS ( fmt ) ( PNG_IMAGE_SAMPLE_CHANNELS ( fmt ) * 256 )
# Skipping MacroDefinition: PNG_IMAGE_PIXEL_ ( test , fmt ) ( ( ( fmt ) & PNG_FORMAT_FLAG_COLORMAP ) ? 1 : test ( fmt ) )
# Skipping MacroDefinition: PNG_IMAGE_PIXEL_CHANNELS ( fmt ) PNG_IMAGE_PIXEL_ ( PNG_IMAGE_SAMPLE_CHANNELS , fmt )
# Skipping MacroDefinition: PNG_IMAGE_PIXEL_COMPONENT_SIZE ( fmt ) PNG_IMAGE_PIXEL_ ( PNG_IMAGE_SAMPLE_COMPONENT_SIZE , fmt )
# Skipping MacroDefinition: PNG_IMAGE_PIXEL_SIZE ( fmt ) PNG_IMAGE_PIXEL_ ( PNG_IMAGE_SAMPLE_SIZE , fmt )
# Skipping MacroDefinition: PNG_IMAGE_ROW_STRIDE ( image ) ( PNG_IMAGE_PIXEL_CHANNELS ( ( image ) . format ) * ( image ) . width )
# Skipping MacroDefinition: PNG_IMAGE_BUFFER_SIZE ( image , row_stride ) ( PNG_IMAGE_PIXEL_COMPONENT_SIZE ( ( image ) . format ) * ( image ) . height * ( row_stride ) )
# Skipping MacroDefinition: PNG_IMAGE_SIZE ( image ) PNG_IMAGE_BUFFER_SIZE ( image , PNG_IMAGE_ROW_STRIDE ( image ) )
# Skipping MacroDefinition: PNG_IMAGE_COLORMAP_SIZE ( image ) ( PNG_IMAGE_SAMPLE_SIZE ( ( image ) . format ) * ( image ) . colormap_entries )

const PNG_IMAGE_FLAG_COLORSPACE_NOT_sRGB = 0x01
const PNG_IMAGE_FLAG_FAST = 0x02
const PNG_IMAGE_FLAG_16BIT_sRGB = 0x04

# Skipping MacroDefinition: png_image_write_get_memory_size ( image , size , convert_to_8_bit , buffer , row_stride , colormap ) png_image_write_to_memory ( & ( image ) , 0 , & ( size ) , convert_to_8_bit , buffer , row_stride , colormap )
# Skipping MacroDefinition: PNG_IMAGE_DATA_SIZE ( image ) ( PNG_IMAGE_SIZE ( image ) + ( image ) . height )
# Skipping MacroDefinition: PNG_ZLIB_MAX_SIZE ( b ) ( ( b ) + ( ( ( b ) + 7U ) >> 3 ) + ( ( ( b ) + 63U ) >> 6 ) + 11U )
# Skipping MacroDefinition: PNG_IMAGE_COMPRESSED_SIZE_MAX ( image ) PNG_ZLIB_MAX_SIZE ( ( png_alloc_size_t ) PNG_IMAGE_DATA_SIZE ( image ) )
# Skipping MacroDefinition: PNG_IMAGE_PNG_SIZE_MAX_ ( image , image_size ) ( ( 8U /*sig*/ + 25U /*IHDR*/ + 16U /*gAMA*/ + 44U /*cHRM*/ + 12U /*IEND*/ + ( ( ( image ) . format & PNG_FORMAT_FLAG_COLORMAP ) ? /*colormap: PLTE, tRNS*/ 12U + 3U * ( image ) . colormap_entries /*PLTE data*/ + ( ( ( image ) . format & PNG_FORMAT_FLAG_ALPHA ) ? 12U /*tRNS*/ + ( image ) . colormap_entries : 0U ) : 0U ) + 12U ) + ( 12U * ( ( image_size ) / PNG_ZBUF_SIZE ) ) /*IDAT*/ + ( image_size ) )
# Skipping MacroDefinition: PNG_IMAGE_PNG_SIZE_MAX ( image ) PNG_IMAGE_PNG_SIZE_MAX_ ( image , PNG_IMAGE_COMPRESSED_SIZE_MAX ( image ) )

const PNG_MAXIMUM_INFLATE_WINDOW = 2
const PNG_SKIP_sRGB_CHECK_PROFILE = 4
const PNG_IGNORE_ADLER32 = 8
const PNG_OPTION_NEXT = 12
const PNG_OPTION_UNSET = 0
const PNG_OPTION_INVALID = 1
const PNG_OPTION_OFF = 2
const PNG_OPTION_ON = 3
const png_libpng_version_1_6_37 = Cstring
const png_struct_def = Cvoid
const png_struct = png_struct_def
const png_const_structp = Ptr{png_struct}
const png_structp = Ptr{png_struct}
const png_structpp = Ptr{Ptr{png_struct}}
const png_info_def = Cvoid
const png_info = png_info_def
const png_infop = Ptr{png_info}
const png_const_infop = Ptr{png_info}
const png_infopp = Ptr{Ptr{png_info}}
const png_structrp = Ptr{png_struct}
const png_const_structrp = Ptr{png_struct}
const png_inforp = Ptr{png_info}
const png_const_inforp = Ptr{png_info}

struct png_color_struct
    red::png_byte
    green::png_byte
    blue::png_byte
end

const png_color = png_color_struct
const png_colorp = Ptr{png_color}
const png_const_colorp = Ptr{png_color}
const png_colorpp = Ptr{Ptr{png_color}}

struct png_color_16_struct
    index::png_byte
    red::png_uint_16
    green::png_uint_16
    blue::png_uint_16
    gray::png_uint_16
end

const png_color_16 = png_color_16_struct
const png_color_16p = Ptr{png_color_16}
const png_const_color_16p = Ptr{png_color_16}
const png_color_16pp = Ptr{Ptr{png_color_16}}

struct png_color_8_struct
    red::png_byte
    green::png_byte
    blue::png_byte
    gray::png_byte
    alpha::png_byte
end

const png_color_8 = png_color_8_struct
const png_color_8p = Ptr{png_color_8}
const png_const_color_8p = Ptr{png_color_8}
const png_color_8pp = Ptr{Ptr{png_color_8}}

struct png_sPLT_entry_struct
    red::png_uint_16
    green::png_uint_16
    blue::png_uint_16
    alpha::png_uint_16
    frequency::png_uint_16
end

const png_sPLT_entry = png_sPLT_entry_struct
const png_sPLT_entryp = Ptr{png_sPLT_entry}
const png_const_sPLT_entryp = Ptr{png_sPLT_entry}
const png_sPLT_entrypp = Ptr{Ptr{png_sPLT_entry}}

struct png_sPLT_struct
    name::png_charp
    depth::png_byte
    entries::png_sPLT_entryp
    nentries::png_int_32
end

const png_sPLT_t = png_sPLT_struct
const png_sPLT_tp = Ptr{png_sPLT_t}
const png_const_sPLT_tp = Ptr{png_sPLT_t}
const png_sPLT_tpp = Ptr{Ptr{png_sPLT_t}}

struct png_text_struct
    compression::Cint
    key::png_charp
    text::png_charp
    text_length::Csize_t
    itxt_length::Csize_t
    lang::png_charp
    lang_key::png_charp
end

const png_text = png_text_struct
const png_textp = Ptr{png_text}
const png_const_textp = Ptr{png_text}
const png_textpp = Ptr{Ptr{png_text}}

struct png_time_struct
    year::png_uint_16
    month::png_byte
    day::png_byte
    hour::png_byte
    minute::png_byte
    second::png_byte
end

const png_time = png_time_struct
const png_timep = Ptr{png_time}
const png_const_timep = Ptr{png_time}
const png_timepp = Ptr{Ptr{png_time}}

struct png_unknown_chunk_t
    name::NTuple{5, png_byte}
    data::Ptr{png_byte}
    size::Csize_t
    location::png_byte
end

const png_unknown_chunk = png_unknown_chunk_t
const png_unknown_chunkp = Ptr{png_unknown_chunk}
const png_const_unknown_chunkp = Ptr{png_unknown_chunk}
const png_unknown_chunkpp = Ptr{Ptr{png_unknown_chunk}}

struct png_row_info_struct
    width::png_uint_32
    rowbytes::Csize_t
    color_type::png_byte
    bit_depth::png_byte
    channels::png_byte
    pixel_depth::png_byte
end

const png_row_info = png_row_info_struct
const png_row_infop = Ptr{png_row_info}
const png_row_infopp = Ptr{Ptr{png_row_info}}
const png_error_ptr = Ptr{Cvoid}
const png_rw_ptr = Ptr{Cvoid}
const png_flush_ptr = Ptr{Cvoid}
const png_read_status_ptr = Ptr{Cvoid}
const png_write_status_ptr = Ptr{Cvoid}
const png_progressive_info_ptr = Ptr{Cvoid}
const png_progressive_end_ptr = Ptr{Cvoid}
const png_progressive_row_ptr = Ptr{Cvoid}
const png_user_transform_ptr = Ptr{Cvoid}
const png_user_chunk_ptr = Ptr{Cvoid}
const png_longjmp_ptr = Ptr{Cvoid}
const png_malloc_ptr = Ptr{Cvoid}
const png_free_ptr = Ptr{Cvoid}
const png_control = Cvoid
const png_controlp = Ptr{png_control}

struct png_image
    opaque::png_controlp
    version::png_uint_32
    width::png_uint_32
    height::png_uint_32
    format::png_uint_32
    flags::png_uint_32
    colormap_entries::png_uint_32
    warning_or_error::png_uint_32
    message::NTuple{64, UInt8}
end

const png_imagep = Ptr{png_image}


text:
@testset "Stemmer" begin

    algs = stemmer_types()
    @test !isempty(algs)

    for alg in algs
        stmr = Stemmer(alg)
        TextAnalysis.release(stmr)
    end

    test_cases = Dict{String, Any}(
        "english" => Dict{AbstractString, AbstractString}(
            "working" => "work",
            "worker" => "worker",
            "aβc" => "aβc",
            "a∀c" => "a∀c"
        )
    )

    for (alg, test_words) in test_cases
        stmr = Stemmer(alg)
        for (n,v) in test_words
            @test v == stem(stmr, n)
        end
    end
end


text:
# Test script for scalings.
include("Utility.jl")
using Plots

results = postprocess("results/cray_BATL_update/")

compiler= (results[:Compiler].=="cray")

results_1 = sort(results[compiler .& (results[:Thread].==1),:],(:Proc))
results_2 = sort(results[compiler .& (results[:Thread].==2),:],(:Proc))
results_4 = sort(results[compiler .& (results[:Thread].==4),:],(:Proc))
results_8 = sort(results[compiler .& (results[:Thread].==8),:],(:Proc))
results_16 = sort(results[compiler .& (results[:Thread].==16),:],(:Proc))
results_32 = sort(results[compiler .& (results[:Thread].==32),:],(:Proc))

runtime = fill(NaN,13,6)
runtime[1:size(results_1)[1],1] = results_1[:Time]
runtime[1:size(results_2)[1],2] = results_2[:Time]
runtime[1:size(results_4)[1],3] = results_4[:Time]
runtime[1:size(results_8)[1],4] = results_8[:Time]
runtime[1:size(results_16)[1],5] = results_16[:Time]
runtime[1:size(results_32)[1],6] = results_32[:Time]

# Calculate the number of cell updates per second for 2nd order scheme
cell_updates = sort(8^3 .* 256 .* results_32[:Thread] .* results_32[:Proc] .* 2)

speed = hcat([cell_updates ./ runtime[:,i] for i in 1:6]...)

proc = results_32[:Proc] .* results_32[:Thread]

pyplot()
mylabel = reshape(["1 thread", "2 thread", "4 thread", "8 thread", "16 threads",
   "32 threads", "", "", "", "", "", ""],1,12)

#=
plt1 = plot(proc,speed, seriestype=[:line, :scatter], markersize=1, lw=1,
     xscale=:log10, yscale=:log10, legend=:best,
     xlim=[10^1, 10^5.5], ylim=[10^5.5, 10^9.5],
     alpha=0.5,
     markercolor=reshape([(1:13)...],1,13),
     linecolor=reshape([(1:13)...],1,13),
     label=mylabel,
     xlabel="number of cores",
     ylabel="cell updates per second",
     title="BATSRUS Hybrid Weak Scaling",
     xtickfontsize=16, ytickfontsize=16,
     legendfontsize=16)
=#

plt1 = plot(proc,speed, seriestype=[:line, :scatter], markersize=1, lw=1,
          legend=:best,
          markercolor=reshape([(1:13)...],1,13),
          linecolor=reshape([(1:13)...],1,13),
          label=mylabel,
          xlabel="number of cores",
          ylabel="cell updates per second",
          title="BATSRUS Hybrid Weak Scaling",
          xtickfontsize=16, ytickfontsize=16,
          legendfontsize=16)

speed_ref = speed[1] ./ proc[1] .* proc
plot!(proc,speed_ref, linestyle=:dashdot, linecolor=:black,
     label="linear ref.")

results[:core] = results[:Thread] .* results[:Proc]

using StatsPlots
mylabel = reshape(["2^5 cores","2^6 cores","2^7 cores","2^8 cores",
   "2^9 cores","2^10 cores","2^11 cores","2^12 cores","2^13 cores",
   "2^14 cores","2^15 cores","2^16 cores","2^17 cores",
   "","","","","","","","","","","","",""],1,26)

plt2 = plot(reuse = false)
@df results plot(:Thread, :Time, xticks=[2^i for i in 0:5],
   group = (:core),
   label=mylabel,
   seriestype=[:line, :scatter],
   markersize=0.6 .* log2.(:Proc),
   markercolor=reshape([(1:13)...],1,13),
   linecolor=reshape([(1:13)...],1,13),
   xlabel="thread",
   ylabel="execution time [s]",
   title="Timing",
   xtickfontsize=16, ytickfontsize=16,
   legendfontsize=12)