text:
include("LESmodels/les_real_calc.jl")

@par function advance_in_time!(s::A) where {A<:@par(AbstractSimulation)}
    calculate_rhs!(s)

    if hasforcing(A)
        s.forcing(s)
        is_spec_time(s) && write_spec_forcing(s)
    end

    mod(s.iteration[],s.dtstats) == 0 && writestats(s)
    time_step!(s)

    s.iteration[] += 1
    s.time[] += get_dt(s)
    return nothing
end

@par function calculate_rhs!(s::A) where {A<:@par(AbstractSimulation)}
    # Set necessary fields to fourier
    mysetfourier!(s.rhs)
    hasdensity(A) && mysetfourier!(s.densitystratification.rhs)
    hasdensityles(A) && mysetfourier!(s.densitystratification.flux)
    haspassivescalar(A) && mysetfourier!(s.passivescalar.rhs)
    haspassivescalarles(A) && mysetfourier!(s.passivescalar.flux)
    hasles(A) && mysetfourier!(s.lesmodel.tau)
    is_dynamic_les(A) && (mysetfourier!(s.lesmodel.û); mysetfourier!(s.lesmodel.M))
    is_dynP_les(A) && (mysetfourier!(s.lesmodel.ŵ))

    #actual calculation
    fourierspacep1!(s)
    realspace!(s)

    has_variable_timestep(s) && set_dt!(s)

    is_spec_time(s) && write_spec(s)
    if is_stats_time(s)
        if !STAT_FIL
            s.nlstats[] = non_linear_stats(s.reductionh,s.reductionv,s.u,s.rhs,s.kspec)
            s.pressstats[] = pressure_stats(s.reductionh,s.reductionv,s)
        else
            s.nlstats[] = non_linear_stats(s.reductionh,s.reductionv,s.reductionh_fil,s.reductionv_fil,s.u,s.rhs,s.kspec)
            s.pressstats[] = pressure_stats(s.reductionh,s.reductionv,s.reductionh_fil,s.reductionv_fil,s)
        end
    end

    fourierspacep2!(s)

    return nothing
end

@par function fourierspacep1!(s::A) where {A<:@par(AbstractSimulation)}
    @mthreads for ind in ZYRANGE
        jj,kk = Tuple(ind)
        k = ZRANGE[kk]
        j = YRANGE[jj]
        fourierspacep1!(j,k,s)
    end
    return nothing
end

@par function fourierspacep1!(j,k,s::A) where {A<:@par(AbstractSimulation)}
    rhs = s.rhs.c
    u = s.u.c
    if hasles(A)
        τ = s.lesmodel.tau.c
        if is_dynP_les(A)
            ŵ = s.lesmodel.ŵ.c
        end
        if is_dynamic_les(A)
            û = s.lesmodel.û.c
            Δ̂² = s.lesmodel.Δ̂²
            M = s.lesmodel.M.c
        end
    end

    if hasdensity(A) && hasdensityles(A)
        ρ = s.densitystratification.ρ
        f = s.densitystratification.flux
    end

    if haspassivescalar(A) && haspassivescalarles(A)
        φ = s.passivescalar.φ
        fφ = s.passivescalar.flux
    end

        @msimd for i in XRANGE
            v = u[i,j,k]
            ∇ = im*K[i,j,k]

            if need_vorticity(A)
                w = ∇ × v
                rhs[i,j,k] = w
            end

            if hasles(A)
                S = symouter(∇,v)
                τ[i,j,k] = S
                if is_dynamic_les(A)
                    G = Gaussfilter(Δ̂²,i,j,k)
                    M[i,j,k] = G*S
                    û[i,j,k] = G*v
                    if is_dynP_les(A)
                        ŵ[i,j,k] = G*w
                    end
                end
            end

            if hasdensity(A) && hasdensityles(A)
                f[i,j,k] = ∇*ρ[i,j,k]
            end

            if haspassivescalar(A) && haspassivescalarles(A)
                fφ[i,j,k] = ∇*φ[i,j,k]
            end

        end
end

@par function realspace!(s::A) where {A<:@par(AbstractSimulation)}
    real!(s.u)
    if need_vorticity(A)
        real!(s.rhs)
    end
    if !is_vorticityEquation(A)
        setreal!(s.equation.uu)
    end

    if is_spec_time(s) || is_stats_time(s)
        setreal!(s.kspec)
        calc_k!(s.kspec.rr,s.u.rr)
        myfourier!(s.kspec)
    end

    haspassivescalar(A) && real!(s.passivescalar.φ)
    hasdensity(A) && real!(s.densitystratification.ρ)

    hasles(s) && real!(s.lesmodel.tau)

    if hasdensity(A)
        if hasdensityles(A)
            real!(s.densitystratification.flux)
        else
            setreal!(s.densitystratification.flux)
        end
    end

    if haspassivescalar(A)
        if haspassivescalarles(A)
            real!(s.passivescalar.flux)
        else
            setreal!(s.passivescalar.flux)
        end
    end

    realspace_LES_calculation!(s)
    is_output_time(s) && writeoutput(s)

    realspacecalculation!(s)

    if has_variable_timestep(s)
        find_max(s.reductionh,s.u.r)
        if hasdensity(s)
            find_max(s.densitystratification.reduction,s.densitystratification.ρ.r)
        end
    end

    #is_output_time(s) && writeoutput(s)

    if is_vorticityEquation(A)
        myfourier!(s.rhs)
        # fix for erros on u×ω calculation
        s.rhs.c[1] = zero(Vec{ComplexF64})
        #dealias!(s.rhs)
    else
        myfourier!(s.equation.uu)
        #dealias!(s.equation.uu)
        mysetfourier!(s.rhs)
        div!(s.rhs,s.equation.uu)
    end

    #if is_output_time(s)
    #    real!(s.rhs)
    #    init = s.iteration[]
    #    write("nn1.$init",s.rhs.rr.x)
    #    write("nn2.$init",s.rhs.rr.y)
    #    write("nn3.$init",s.rhs.rr.z)
    #    myfourier!(s.rhs)
    #end

    myfourier!(s.u)
    pressure_projection!(s.u,s)

    if haspassivescalar(A)
        myfourier!(s.passivescalar.flux)
        #dealias!(s.passivescalar.flux)
        myfourier!(s.passivescalar.φ)
    end
    if hasdensity(A)
        myfourier!(s.densitystratification.flux)
        #dealias!(s.densitystratification.flux)
        myfourier!(s.densitystratification.ρ)
    end

    if hasles(s)
        if !(!is_SandP(A) && is_FakeSmagorinsky(A))
            myfourier!(s.lesmodel.tau)
            #dealias!(s.lesmodel.tau)
        end
    end

    return nothing
end

function calc_k!(k,u)
    @mthreads for j in TRANGE
        calc_k!(k,u,j)
    end
end

@inline function calc_k!(k,u,j)
    @inbounds for i in REAL_RANGES[j]
        k[i] = 0.5*(u[i]⋅u[i])
    end
end

@par function realspacecalculation!(s::A) where {A<:@par(AbstractSimulation)}
    #@assert !(hasdensity(A) & haspassivescalar(A))
    @mthreads for j in TRANGE
        realspacecalculation!(s,j)
    end
    return nothing
end

@par function realspacecalculation!(s::A,j::Integer) where {A<:@par(AbstractSimulation)}

    u = s.u.rr
    rhs = s.rhs.rr

    if !is_vorticityEquation(A)
        uu = s.equation.uu.rr
    end

    if haspassivescalar(A)
        φ = s.passivescalar.φ.field.data
        fφ = s.passivescalar.flux.rr
    end

    if hasdensity(A)
        ρ = s.densitystratification.ρ.field.data
        f = s.densitystratification.flux.rr
    end

    @inbounds @msimd for i in REAL_RANGES[j]

        v = u[i]
        if is_vorticityEquation(A)
            w = rhs[i]
            rhs[i] = v × w
        else
            uu[i] = traceless(symouter(v,-v))
        end

        if hasdensity(A)
            rhsden = (-ρ[i]) * v
            hasdensityles(A) && (rhsden += f[i])
            f[i] = rhsden
        end

        if haspassivescalar(A)
            rhsp = (-φ[i]) * v
            haspassivescalarles(A) && (rhsp += fφ[i])
            fφ[i] = rhsp
        end

    end

    return nothing
end

@par function fourierspacep2_velocity!(s::A) where {A<:@par(AbstractSimulation)}
    @mthreads for ind in ZYRANGE
        jj,kk = Tuple(ind)
        k = ZRANGE[kk]
        j = YRANGE[jj]
        fourierspacep2_velocity!(j,k,s)
    end
    return nothing
end

@inline @par function fourierspacep2_velocity!(j,k,s::A) where {A<:@par(AbstractSimulation)}
    u = s.u.c
    rhsv = s.rhs.c

    if !is_implicit(A)
        mν = -ν
        if hashyperviscosity(A)
            mνh = -nuh(A)
            M = get_hyperviscosity_exponent(A)
        end
    end

    if hasles(A)
        τ = s.lesmodel.tau.c
    end

    if hasdensity(A)
        ρ = s.densitystratification.ρ.field
        g = gravity(s.densitystratification)
    end

        @inbounds @msimd for i in XRANGE

            kh = K[i,j,k]
            K2 = kh⋅kh
            v = u[i,j,k]
            rhs = rhsv[i,j,k]

            if !is_implicit(A)
                if hashyperviscosity(A)
                    rhs += (K2*mν + (K2^M)*mνh)*v
                else
                    rhs += (K2*mν)*v
                end
            end

            if hasles(A)
                if !(!is_SandP(A) && is_FakeSmagorinsky(A))
                    rhs += (im*kh)⋅τ[i,j,k]
                end
            end

            if hasdensity(A)
                rhs += ρ[i,j,k]*g
            end

            p1 = ifelse(k==j==i==1,zero(ComplexF64),-(kh⋅rhs)/K2)
            rhsv[i,j,k] = p1*kh + rhs
        end
end

@par function fourierspacep2!(s::A) where {A<:@par(AbstractSimulation)}
    #is_implicit(typeof(s.timestep)) || add_viscosity!(s.rhs,s.u,ν,s)
    #hasles(s) && add_residual_tensor!(s.rhs,s.lesmodel.tau,s)
    #if hasdensity(A)
        #Gdirec = graddir(s.densitystratification)
        #gdir = Gdirec === :x ? s.rhs.cx : Gdirec === :y ? s.rhs.cy : s.rhs.cz
        #addgravity!(gdir, complex(s.densitystratification.ρ), -gravity(s.densitystratification), s)
    #end
    #pressure_projection!(s.rhs.c.x,s.rhs.c.y,s.rhs.c.z,s)
    fourierspacep2_velocity!(s)
    if haspassivescalar(A)
        gdir = graddir(s.passivescalar) === :x ? s.u.c.x : graddir(s.passivescalar) === :y ? s.u.c.y : s.u.c.z
        div!(complex(s.passivescalar.rhs), s.passivescalar.flux.c.x, s.passivescalar.flux.c.y, s.passivescalar.flux.c.z, gdir, -meangradient(s.passivescalar), s)
        is_implicit(typeof(s.passivescalar.timestep)) || add_scalar_difusion!(complex(s.passivescalar.rhs),complex(s.passivescalar.φ),diffusivity(s.passivescalar),s)
    end
    if hasdensity(A)
        gdir = graddir(s.densitystratification) === :x ? s.u.c.x : graddir(s.densitystratification) === :y ? s.u.c.y : s.u.c.z
        div!(complex(s.densitystratification.rhs), s.densitystratification.flux.c.x, s.densitystratification.flux.c.y, s.densitystratification.flux.c.z, gdir, -meangradient(s.densitystratification), s)
        is_implicit(typeof(s.densitystratification.timestep)) || add_scalar_difusion!(complex(s.densitystratification.rhs),complex(s.densitystratification.ρ),diffusivity(s.densitystratification),s)
    end
    return nothing
end

function div!(out,input)
    @mthreads for ind in ZYRANGE
        jj,kk = Tuple(ind)
        k = ZRANGE[kk]
        j = YRANGE[jj]
            @inbounds @msimd for i in XRANGE
                out[i,j,k] = (im*K[i,j,k]) ⋅ input[i,j,k]
            end
    end
end

#@par function add_residual_tensor!(rhs::VectorField,τ::SymmetricTracelessTensor,s::@par(AbstractSimulation))
#    @mthreads for k in 1:Nz
#        add_residual_tensor!(rhs,τ,k,s)
#    end
#end
#
#@inline @par function add_residual_tensor!(rhs::VectorField,tau::SymmetricTracelessTensor,k::Int,s::@par(AbstractSimulation))
#    rx = rhs.cx
#    ry = rhs.cy
#    rz = rhs.cz
#    txx = tau.cxx
#    txy = tau.cxy
#    txz = tau.cxz
#    tyy = tau.cyy
#    tyz = tau.cyz
#    @inbounds for j in 1:Ny
#        @msimd for i in 1:Nx
#            rx[i,j,k] += im*(kx[i]*txx[i,j,k] + ky[j]*txy[i,j,k] + kz[k]*txz[i,j,k])
#            ry[i,j,k] += im*(kx[i]*txy[i,j,k] + ky[j]*tyy[i,j,k] + kz[k]*tyz[i,j,k])
#            rz[i,j,k] += im*(kx[i]*txz[i,j,k] + ky[j]*tyz[i,j,k] + kz[k]*(-txx[i,j,k]-tyy[i,j,k]))
#        end
#    end
#end

#function add_viscosity!(rhs::VectorField,u::VectorField,ν::Real,s::AbstractSimulation)
#    if hashyperviscosity(s)
#        _add_hviscosity!(rhs.c.x,u.c.x,s)
#        _add_hviscosity!(rhs.c.y,u.c.y,s)
#        _add_hviscosity!(rhs.c.z,u.c.z,s)
#    else
#        _add_viscosity!(rhs.c.x,u.c.x,-ν,s)
#        _add_viscosity!(rhs.c.y,u.c.y,-ν,s)
#        _add_viscosity!(rhs.c.z,u.c.z,-ν,s)
#    end
#end
#
#@par function _add_viscosity!(rhs::AbstractArray,u::AbstractArray,mν::Real,s::@par(AbstractSimulation))
#    @mthreads for k in ZRANGE
#        for j in YRANGE
#            @inbounds @msimd for i in XRANGE
#                #rhs[i,j,k] = (kx[i]*kx[i] + ky[j]*ky[j] + kz[k]*kz[k])*mŒΩ*u[i,j,k] + rhs[i,j,k]
#                rhs[i,j,k] = muladd(muladd(KX[i], KX[i], muladd(KY[j], KY[j], KZ[k]*KZ[k])), mν*u[i,j,k], rhs[i,j,k])
#            end
#        end
#    end
#end
#
#@par function _add_hviscosity!(rhs::AbstractArray,u::AbstractArray,s::@par(AbstractSimulation))
#    mν = -nu(s)
#    mνh = -nuh(s)
#    M = get_hyperviscosity_exponent(s)
#    @mthreads for k in ZRANGE
#        for j in YRANGE
#            @inbounds @msimd for i in XRANGE
#                modk2 = muladd(KX[i], KX[i], muladd(KY[j], KY[j], KZ[k]*KZ[k]))
#                rhs[i,j,k] = muladd(muladd(modk2, mν, modk2^M * mνh), u[i,j,k], rhs[i,j,k])
#            end
#        end
#    end
#end
#
#function add_scalar_difusion!(rhs::AbstractArray,u::AbstractArray,ν::Real,s::AbstractSimulation)
#    _add_viscosity!(rhs,u,-ν,s)
#end

@par function pressure_projection!(rhs,s::@par(AbstractSimulation))
    @mthreads for ind in ZYRANGE
        jj,kk = Tuple(ind)
        k = ZRANGE[kk]
        j = YRANGE[jj]
            @inbounds @msimd for i in XRANGE
                rhsh = rhs[i,j,k]
                kh = K[i,j,k]
                K2 = kh⋅kh
                p1 = ifelse(k==j==i==1,zero(ComplexF64),-(kh⋅rhsh)/K2)
                rhs[i,j,k] = p1*kh + rhsh
            end
    end
end

#@par function addgravity!(rhs,ρ,g::Real,s::@par(AbstractSimulation))
    #@mthreads for k in ZRANGE
        #for j in YRANGE
            #@inbounds @msimd for i in XRANGE
                #rhs[i,j,k] = muladd(ρ[i,j,k],g,rhs[i,j,k])
            #end
        #end
    #end
#end

@par function time_step!(s::A) where {A<:@par(AbstractSimulation)}
    s.timestep(s.u,s.rhs,s)
    haspassivescalar(s) && s.passivescalar.timestep(s.passivescalar.φ,s.passivescalar.rhs,s)
    hasdensity(s) && s.densitystratification.timestep(s.densitystratification.ρ,s.densitystratification.rhs,s)
end


text:
function plotssavefig(plts, figname="visclaw.svg"; num_start::Integer=1, kwargs...)
    dn = dirname(figname)
    bn = basename(figname)

    for i=1:length(plts)
        numstr = @sprintf("%03d",(i-1)+num_start)
        bnnum = occursin(".", bn) ? replace(bn, "." => "-"*numstr*".") : bn*"-"*numstr
        Plots.savefig(plts[i], joinpath(dn,bnnum); kwargs...)
    end

end

function plotsgif(plts, gifname::AbstractString="visclaw.gif"; kwargs...)
    anim = Plots.Animation()
    map(p->Plots.frame(anim, p), plts)
    Plots.gif(anim, gifname; kwargs...)
end

ext: jl


text:
"""
    Histogram

Histogram represents the data as a vector where each coordinate corresponds to
one element of the data universe. This is the default representation for MWEM.
"""
mutable struct Histogram <: Data
    weights::Array{Float64, 1}
    num_samples::Int64
end

function Histogram(weights::Array{Float64, 1})
    Histogram(weights, 0)
end

struct HistogramQuery <: Query
    weights::Array{Float64, 1}
end

struct HistogramQueries <: Queries
    queries::Array{Float64, 2}
end

function get(queries::HistogramQueries, i::QueryIndex)
    HistogramQuery(queries.queries[:, i])
end

function evaluate(query::HistogramQuery, h::Histogram)
    dot(query.weights, h.weights)
end

function evaluate(queries::HistogramQueries, h::Histogram)
    queries.queries' * vec
end

function normalize!(h::Histogram)
    h.weights /= sum(h.weights)
    h
end

# if error is > 0, increase weight on positive elements of queries[query] and
# decrease weight on negative elements. Magnitude of error determines step size.
function update!(q::HistogramQuery, h::Histogram, error::Float64)

    Threads.@threads for j = 1:length(h.weights)# @threads can help to speed up
        @inbounds h.weights[j] *= exp(error * q.weights[j] / 2.0)
    end
    normalize!(h)
end

function initialize(queries::Queries, data::Histogram, ps::MWParameters)
    data = normalize!(data)
    histogram_length = length(data.weights)
    num_samples = data.num_samples
    if ps.noisy_init
        # Noisy init incurs an additional `epsilon` privacy cost
        weights = zeros(Float64, histogram_length)
        noise = rand(Laplace(0.0, 1/(num_samples*ps.epsilon*ps.init_budget)), histogram_length)
        @simd for i = 1:histogram_length
             @inbounds weights[i] =
                 max(data.weights[i] + noise[i] - 1.0/(e*num_samples*ps.epsilon), 0.0)
        end
        weights /= sum(weights)
        synthetic = Histogram(0.5 * weights + 0.5/histogram_length)
        ps.epsilon = (1-ps.init_budget)*ps.epsilon
    else
        weights = ones(Float64, histogram_length)
        synthetic = Histogram(weights/histogram_length)
    end
    real_answers = evaluate(queries, data)
    scale = 2.0/(ps.epsilon*num_samples)
    MWState(data, synthetic, queries, real_answers, Dict{Int, Float64}(), scale)
end

function initialize(queries::HistogramQueries, data::Tabular, ps::MWParameters)
    initialize(queries, Histogram(data), ps)
end

"""
    Histogram(table)

Create histogram representation from tabular data.
"""
function Histogram(table::Tabular)
    d, n = size(table.data)
    histogram = zeros(Float64,2^d)
    for i = 1:n
        num = 0
        x = vec(table.data[:, i])
        for j = 1:d
            num += convert(Int64, x[d-j+1]) * 2^(j-1)
        end
        histogram[num+1] += 1.0
    end
    normalize!(Histogram(histogram, n))
end

"""
    Tabular(histogram, n)

Create tabular data by sampling n times weighted by histogram.
"""
function Tabular(histogram::Histogram, n::Int)
    N = length(histogram.weights)
    d = convert(Int64, log(2, N))
    idx = wsample(collect(0:N-1), histogram.weights, n)
    data_matrix = zeros(d, n)
    for i = 1:n
        #println(idx[i], 2, d)
        data_matrix[:, i] = reverse(digits(idx[i], base = 2, pad = d))
    end
    Tabular(data_matrix)
end


text:
module TestGetDataframe

using Mimi
using Test

import Mimi:
    reset_compdefs, _load_dataframe

reset_compdefs()

#------------------------------------------------------------------------------
#   1. Test with 1 dimension
#------------------------------------------------------------------------------

model1 = Model()

@defcomp testcomp1 begin
    var1 = Variable(index=[time])
    par1 = Parameter(index=[time])
    par_scalar = Parameter()

    function run_timestep(p, v, d, t)
        v.var1[t] = p.par1[t]
    end
end

@defcomp testcomp2 begin
    var2 = Variable(index=[time])
    par2 = Parameter(index=[time])

    function run_timestep(p, v, d, t)
        v.var2[t] = p.par2[t]
    end
end

years = collect(2015:5:2110)
late_first = 2030
early_last = 2100

set_dimension!(model1, :time, years)
add_comp!(model1, testcomp1)
set_param!(model1, :testcomp1, :par1, years)
set_param!(model1, :testcomp1, :par_scalar, 5.)
add_comp!(model1, testcomp2; first = late_first, last = early_last)
set_param!(model1, :testcomp2, :par2, late_first:5:early_last)

# Test running before model built
@test_throws ErrorException df = getdataframe(model1, :testcomp1, :var1)

# Now run model
run(model1)

# Test trying to load a dataframe for a scalar parameter
@test_throws ErrorException getdataframe(model1, :testcomp1, :par_scalar)

# Test trying to getdataframe from a variable that does not exist in the component
@test_throws ErrorException getdataframe(model1, :testcomp1, :var2)

# Regular getdataframe
df = getdataframe(model1, :testcomp1=>:var1, :testcomp1=>:par1, :testcomp2=>:var2, :testcomp2=>:par2)
dim = Mimi.dimension(model1, :time)
@test df[:var1] == df[:par1] == years
@test all(ismissing, df[:var2][1 : dim[late_first]-1])
@test all(ismissing, df[:par2][1 : dim[late_first]-1])
@test df[:var2][dim[late_first] : dim[early_last]] == df[:par2][dim[late_first] : dim[early_last]] == late_first:5:early_last
@test all(ismissing, df[:var2][dim[years[end]] : dim[early_last]])
@test all(ismissing, df[:par2][dim[years[end]] : dim[early_last]])

# Test trying to load an item into an existing dataframe where that item key already exists
@test_throws ErrorException _load_dataframe(model1, :testcomp1, :var1, df)

#------------------------------------------------------------------------------
#   2. Test with > 2 dimensions
#------------------------------------------------------------------------------

@defcomp testcomp3 begin
    par3 = Parameter(index=[time, regions, rates])
    var3 = Variable(index=[time])
end

regions = [:reg1, :reg2]
rates   = [0.025, 0.05]

# A. Simple case where component has same time length as model

model2 = Model()

set_dimension!(model2, :time, years)
set_dimension!(model2, :regions, regions)
set_dimension!(model2, :rates, rates)

data = Array{Int}(undef, length(years), length(regions), length(rates))
data[:] = 1:(length(years) * length(regions) * length(rates))

add_comp!(model2, testcomp3)
set_param!(model2, :testcomp3, :par3, data)

run(model2)

df2 = getdataframe(model2, :testcomp3, :par3)
@test size(df2) == (length(data), 4)

# Test trying to combine two items with different dimensions into one dataframe
@test_throws ErrorException getdataframe(model2, Pair(:testcomp3, :par3), Pair(:testcomp3, :var3))

# B. Test with shorter time than model

model3 = Model()
set_dimension!(model3, :time, years)
set_dimension!(model3, :regions, regions)
set_dimension!(model3, :rates, rates)
add_comp!(model3, testcomp3; first = late_first, last = early_last)

par3 = Array{Float64}(undef, length(late_first:5:early_last), length(regions), length(rates))
par3[:] = 1:(length(late_first:5:early_last) * length(regions) * length(rates))
set_param!(model3, :testcomp3, :par3, par3)
run(model3)

df3 = getdataframe(model3, :testcomp3=>:par3)
@test size(df3) == (length(rates)*length(regions)*length(years), 4)

# Test that times outside the component's time span are padded with `missing` values
@test all(ismissing, df3[:par3][1 : (length(rates)*length(regions)*(dim[late_first]-1))])
@test all(ismissing, df3[:par3][end - (length(rates)*length(regions)*(dim[end]-dim[early_last]))+1: end])

end #module


text:
using DiskFrame
using Test

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


text:
abstract type PotentialMatrix{T} <: NoAccessArray{T, 2} end

@inline LinearAlgebra.mul!(
    Y::AbstractArray{T, 1},
    A::PotentialMatrix{T},
    x::AbstractArray{T, 1}
) where T = Y .= A * x

@inline _kcfg(A::PotentialMatrix{T}) where T = _kcfg(size(A, 1))

struct LaplacePotentialMatrix{T, L <: PotentialType} <: PotentialMatrix{T}
    Ξ       ::PositionVector{T}
    elements::TriangleVector{T}
end

@inline LaplacePotentialMatrix{L}(
    Ξ       ::PositionVector{T},
    elements::TriangleVector{T}
) where {T, L <: PotentialType} = LaplacePotentialMatrix{T, L}(Ξ, elements)

@inline Base.size(A::LaplacePotentialMatrix{T}) where T = (length(A.Ξ), length(A.elements))

@inline Base.:*(
    A::LaplacePotentialMatrix{T, SingleLayer},
    x::AbstractArray{T, 1}
) where T = _mul(A, x, _laplace_single_mul_kernel!)

@inline Base.:*(
    A::LaplacePotentialMatrix{T, DoubleLayer},
    x::AbstractArray{T, 1}
) where T = _mul(A, x, _laplace_double_mul_kernel!)

@inline LinearAlgebra.diag(
    A::LaplacePotentialMatrix{T, SingleLayer},
    k::Int = 0
) where T = _diag(A, k, _laplace_single_diag_kernel!)

@inline LinearAlgebra.diag(
    A::LaplacePotentialMatrix{T, DoubleLayer},
    k::Int = 0
) where T = _diag(A, k, _laplace_double_diag_kernel!)

function _diag(A::LaplacePotentialMatrix{T}, k::Int, pot::F) where {T, F <: Function}
    @assert size(A, 1) == size(A, 2) "diag requires a square matrix"
    k != 0 && error("diag not defined for k != 0 on ", typeof(A))
    dst = CuArray{T}(undef, size(A, 1))
    cfg = _kcfg(A)
    @cuda threads=cfg.threads blocks=cfg.blocks pot(dst, A.Ξ, A.elements)
    Array(dst)
end

function _mul(A::LaplacePotentialMatrix{T}, x::CuArray{T, 1}, pot::F) where {T, F <: Function}
    dst = CuArray{T}(undef, size(A, 1))
    cfg = _kcfg(A)
    @cuda threads=cfg.threads blocks=cfg.blocks pot(dst, A.Ξ, A.elements, x)
    dst
end

@inline _mul(
    A  ::LaplacePotentialMatrix{T},
    x  ::AbstractArray{T, 1},
    pot::F
) where {T, F <: Function} = Array(_mul(A, CuArray(x), pot))

struct ReYukawaPotentialMatrix{T, L <: PotentialType} <: PotentialMatrix{T}
    Ξ       ::PositionVector{T}
    elements::TriangleVector{T}
    yuk     ::T
end

@inline ReYukawaPotentialMatrix{L}(
    Ξ       ::PositionVector{T},
    elements::TriangleVector{T},
    yuk     ::T
) where {T, L <: PotentialType} = ReYukawaPotentialMatrix{T, L}(Ξ, elements, yuk)

@inline Base.size(A::ReYukawaPotentialMatrix{T}) where T = (length(A.Ξ), length(A.elements))

@inline Base.:*(
    A::ReYukawaPotentialMatrix{T, SingleLayer},
    x::AbstractArray{T, 1}
) where T = _mul(A, x, _reyukawa_single_mul_kernel!)

@inline Base.:*(
    A::ReYukawaPotentialMatrix{T, DoubleLayer},
    x::AbstractArray{T, 1}
) where T = _mul(A, x, _reyukawa_double_mul_kernel!)

@inline LinearAlgebra.diag(
    A::ReYukawaPotentialMatrix{T, SingleLayer},
    k::Int = 0
) where T = _diag(A, k, _reyukawa_single_diag_kernel!)

@inline LinearAlgebra.diag(
    A::ReYukawaPotentialMatrix{T, DoubleLayer},
    k::Int = 0
) where T = _diag(A, k, _reyukawa_double_diag_kernel!)

function _diag(A::ReYukawaPotentialMatrix{T}, k::Int, pot::F) where {T, F <: Function}
    @assert size(A, 1) == size(A, 2) "diag requires a square matrix"
    k != 0 && error("diag not defined for k != 0 on ", typeof(A))
    dst = CuArray{T}(undef, size(A, 1))
    cfg = _kcfg(A)
    @cuda threads=cfg.threads blocks=cfg.blocks pot(dst, A.Ξ, A.elements, A.yuk)
    Array(dst)
end

function _mul(
    A  ::ReYukawaPotentialMatrix{T},
    x  ::CuArray{T, 1},
    pot::F
) where {T, F <: Function}
    dst = CuArray{T}(undef, size(A, 1))
    cfg = _kcfg(A)
    @cuda threads=cfg.threads blocks=cfg.blocks pot(dst, A.Ξ, A.elements, x, A.yuk)
    dst
end

@inline _mul(
    A  ::ReYukawaPotentialMatrix{T},
    x  ::AbstractArray{T, 1},
    pot::F
) where {T, F <: Function} = Array(_mul(A, CuArray(x), pot))


text:
struct Time
    total_time_ns::TimeNS
    ncalls::Int
    tag::Union{Symbol,Nothing}
    file::Union{Symbol,Nothing}
    line::Int
    fullname::Tuple{Vararg{Symbol}}
    pkgid::Union{UUID,Nothing}
end

time_table() = append_times!(Time[])

function append_times!(table)
    for thread_records in __RECORDS
        for record in thread_records
            ncalls = sum(getindex, record.total_calls)
            ncalls == 0 && continue
            total_time_ns = sum(getindex, record.total_times)
            row = Time(
                total_time_ns,
                ncalls,
                record.location.tag,
                record.location.file,
                record.location.line,
                record.location.fullname,
                record.location.pkgid,
            )
            push!(table, row)
        end
    end
    return table
end


text:
function get_neighbours(layout, ix, jy)
    nr_rows, nr_cols = size(layout)
    neighbours = Char[]
    for dx in -1:1, dy in -1:1
        if dx != 0 || dy != 0
            x = ix + dx
            y = jy + dy
            if 1 <= x <= nr_rows && 1 <= y <= nr_cols
                push!(neighbours, layout[x, y])
            end
        end
    end
    neighbours
end

function single_step(layout)
    updated = Array{Char, 2}(undef, size(layout))
    nr_changed = 0
    for jy in 1:size(layout, 2)
        for ix in 1:size(layout, 1)
            item = layout[ix, jy]
            if item == '.'
                updated[ix, jy] = '.'
            else
                neighbours = get_neighbours(layout, ix, jy)
                if item == 'L' && all(neighbour != '#' for neighbour in neighbours)
                    updated[ix, jy] = '#'
                    nr_changed += 1
                elseif item == '#' && count(n -> n == '#', neighbours) >= 4
                    updated[ix, jy] = 'L'
                    nr_changed += 1
                else
                    updated[ix, jy] = item
                end
            end
        end
    end
    updated, nr_changed
end

function seating_system(seat_layout)
    while true
        seat_layout, nr_changed_seats = single_step(seat_layout)
        if nr_changed_seats == 0
            return count(x -> x == '#', seat_layout)
        end
    end
end

function main()
    lines = readlines("../resources/input_11.txt")
    seat_layout = permutedims(hcat(map(collect, lines)...))
    result = seating_system(seat_layout)
    println(result)
end

main()


text:
abstract type Healing end
abstract type NoHealing <: Healing end
abstract type AnnealingOnly <: Healing end
abstract type IsotropicDomainOnly <: Healing end
abstract type AnnealingIsotropicDomain <: Healing end
abstract type FabricDestructionOnly <: Healing end
abstract type AnnealingFabricDestruction <: Healing end

struct FiniteStrainEllipsoid{T}
    x1::T
    x2::T
    y1::T
    y2::T
    a1::T
    a2::T
end

# Define isotropic domain
struct IsotropicDomain{T}
    r::T  # depth above which Ω is isotropic
end

# Rate of annealing
struct Annealing{T}
    rate::T # rate of annealing
end

# Destruction of fabric
struct FabricDestruction{T}
    ϵ::T # ratio a1/a2 at which aggregate goes back to its anisotropic form
end

struct FiniteStrain{T} <: Healing
    fse::Matrix{FiniteStrainEllipsoid{Float64}} # matrix containing the FSE at every ip
    annealing::Annealing # rate of annealing
    destruction::FabricDestruction # ratio a1/a2 at which aggregate goes back to its anisotropic form
    isotropic_domain::IsotropicDomain # depth above which Ω is isotropic

    function FiniteStrain(nel; nip = 6, ϵ = 1e3, annealing_rate = 0, r_iso = 0)
        # instantiate isotropic finite strain ellipsoid
        fse = [FiniteStrainEllipsoid(1.0, 0.0, 0.0, 1.0, 1.0, 1.0) for _ in 1:nel, _ in 1:nip]
        # define active healing mechanism(s)
        if (annealing_rate != 0) && (ϵ != 0)
            type = AnnealingFabricDestruction

        elseif (annealing_rate != 0) && (ϵ == 0)
            type = AnnealingOnly

        elseif (annealing_rate == 0) && (ϵ != 0)
            type = FabricDestructionOnly

        else
            type = NoHealing

        end
        # create object
        new{type}(
            fse,
            Annealing(annealing_rate),
            FabricDestruction(ϵ),
            IsotropicDomain(r_iso)
        )

    end
end

*(a::Number, annealing::Annealing) = a*annealing.rate
*(annealing::Annealing, a::Number) = a*annealing

rebuild_FSE(vx1, vx2, vy1, vy2, a1, a2) = [FiniteStrainEllipsoid(vx1[i], vx2[i], vy1[i], vy2[i], a1[i], a2[i]) for i in CartesianIndices(a1)]

function isotropic_lithosphere!(F, idx)
    Is = @SMatrix [1.0 0.0; 0.0 1.0]
    Threads.@threads for i in idx
        @inbounds F[i] = Is
    end
end

function healing(F, FSE)
    Is = @SMatrix [1.0 0.0; 0.0 1.0]
    Threads.@threads for i in CartesianIndices(F)
        @inbounds if FSE[i].a1./FSE[i].a2 > 1e2
            F[i] = Is
        end
    end
    return F
end

function getFSE(F, FSE::FiniteStrain{NoHealing})
    # F can grow A LOT in long computations, eventually overflowing at ~1e309
    # Thus we need to normalize F from time to time. We normalize F ∈ Ω w.r.t.
    # the same number (max(F[1])), otherwise a heterogeneous normalization
    # will screw up the particles interpolation
    normalize_F!(F)
    Threads.@threads for iel in CartesianIndices(F)
        @inbounds FSE.fse[iel] = _FSE(F[iel])
    end
    FSE, F
end

function _FSE(Fi)

    # Compute FSE
    # Fi = normalize_F(F) # normalize F if its too big, otherwise F*F' will be Inf
    L = Fi * Fi' # left-strecth tensor
    # eigval, evect = eigen(Fi * Fi') # get length of FSE semi-axes and orientation
    # imax, imin = eigval_order(eigval) # get the right order of the semi-axis length
    eigval, evect, imax, imin = _eigen(L) # get length of FSE semi-axes and orientation
                                          # and the right order of eigvals and eigvects

    # Fill FSE
    FiniteStrainEllipsoid(
        evect[1,imax], # vx1
        evect[1,imin], # vx2
        evect[2,imax], # vy1
        evect[2,imin], # vy2
        √(abs(eigval[imax])), # a1
        √(abs(eigval[imin])), # a2
    )

end

function normalize_F!(F; ϵ = 1e30)
    normalizer = abs(maximum(F[1])) # make sure we take the absolute value, don't want to change signs
                                    # need to check only one component, as the rest will havea similar exponent
    if normalizer > ϵ
        normalizer = 1/normalizer
        Threads.@threads for i in eachindex(F)
            @inbounds @fastmath F[i] *= normalizer
        end
    end
end

function normalize_F(F; ϵ = 1e30)
    Fmax = abs(maximum(F)) # take absolute value because we do not want to change the sign of Fij
    if Fmax > ϵ
        return F./Fmax
    else
        return F
    end
end

eigval_order(eigval) = ifelse(
    @inbounds(eigval[2] > eigval[1]),
    (2, 1),
    (1, 2)
)

function getFSE(F, FSE::FiniteStrain{FabricDestructionOnly})
    Threads.@threads for iel in eachindex(F)
        healing_FSE!(FSE, F, iel)
    end
    return FSE, F
end

@inbounds function healing_FSE!(FSE, F, iel)

    # Compute FSE
    Fi = normalize_F(F[iel]) # normalize F if its too big, otherwise F*F' will be Inf
    eigval, evect = eigen(Fi * Fi') # get length of FSE semi-axes and orientation
    imax, imin = eigval_order(eigval) # get the right order of the semi-axis length

    a1 = √(abs(eigval[imax]))
    a2 = √(abs(eigval[imin]))
    @inbounds if a1/a2 < FSE.destruction.ϵ
        # Fill FSE
        FSE.fse[iel] = FiniteStrainEllipsoid(
            evect[1,imax]::Float64, # vx1
            evect[1,imin]::Float64, # vx2
            evect[2,imax]::Float64, # vy1
            evect[2,imin]::Float64, # vy2
            a1, # a1
            a2, # a2
        )

    else
        # destruction of the fabric -> isotropic aggregate
        F[iel] = @SMatrix [1.0 0.0; 0.0 1.0]
        # Fill FSE
        FSE.fse[iel] = FiniteStrainEllipsoid(
            1.0, # vx1
            0.0, # vx2
            0.0, # vy1
            1.0, # vy2
            1.0, # a1
            1.0, # a2
        )

    end

end

function getFSE(F, FSE::FiniteStrain{AnnealingOnly})
    annealing = FSE.annealing.rate
    Threads.@threads for iel in eachindex(F)
        @inbounds FSE.fse[iel] = _FSE_annealing(F[iel], annealing)
    end
    return FSE, F
end

@inbounds function _FSE_annealing(F, s)

    # Compute FSE
    Fi = normalize_F(F) # normalize F if its too big, otherwise F*F' will be Inf
    L = Fi * Fi' # left-strecth tensor
    # eigval, evect = eigen(Fi * Fi') # get length of FSE semi-axes and orientation
    # imax, imin = eigval_order(eigval) # get the right order of the semi-axis length
    eigval, evect, imax, imin = _eigen(L) # get length of FSE semi-axes and orientation
                                          # and the right order of eigvals and eigvects

    # principal semi-axes
    a1, a2 = √(abs(eigval[imax])),  √(abs(eigval[imin]))
    # eigenvectors
    λ11, λ12, λ21, λ22 = evect[1,imax], evect[2,imax], evect[1,imin], evect[2,imin]
    # unstretch semi-axes due to annealing
    a1u, a2u = unstretch_axes(a1, a2, s)
    # Fu = recover_F(Fi, L, λ11, λ12, λ21, λ22, a1u, a2u)
    # Fu = recover_F(Fi, L, λ11, λ12, λ21, λ22, a1, a2)

    # force isotropy in case we unstretched too much
    a1u, a2u = force_isotropy(a1, a2)

    # Fill FSE
    FSE = FiniteStrainEllipsoid(
        λ11, # vx1
        λ21, # vx2
        λ12, # vy1
        λ22, # vy2
        a1u,
        a2u,
    )

    return FSE

end

function force_isotropy(a1::T, a2::T) where T
    if (a1 ≤ 1) || (a2 ≥ 1)
        a1 = one(T)
        a2 = zero(T)
    end
    return a1, a2
end

function getFSE(F, FSE::FiniteStrain{AnnealingFabricDestruction})
    # F can grow A LOT in long computations, eventually overflowing at ~1e309
    # Thus we need to normalize F from time to time. We normalize F ∈ Ω w.r.t.
    # the same number (max(F[1])), otherwise a heterogeneous normalization
    # will screw up the particles interpolation
    normalize_F!(F)
    Threads.@threads for iel in CartesianIndices(F)
        @inbounds FSE.fse[iel], F[iel] = _FSE(F[iel], FSE.destruction, FSE.annealing)
    end
    return FSE, F
end

function _FSE(F, destruction::FabricDestruction, annealing::Annealing)

    ϵ, s = destruction.ϵ, annealing.rate

    # Compute FSE
    Fi = normalize_F(F) # normalize F if its too big, otherwise F*F' will be Inf
    L = Fi * Fi' # left-strecth tensor
    # eigval, evect = eigen(Fi * Fi') # get length of FSE semi-axes and orientation
    # imax, imin = eigval_order(eigval) # get the right order of the semi-axis length
    eigval, evect, imax, imin = _eigen(L) # get length of FSE semi-axes and orientation
                                          # and the right order of eigvals and eigvects

    # principal semi-axes
    a1, a2 = √(abs(eigval[imax])),  √(abs(eigval[imin]))
    # eigenvectors
    λ11, λ12, λ21, λ22 = evect[1,imax], evect[2,imax], evect[1,imin], evect[2,imin]
    # unstretch semi-axes due to annealing
    a1u, a2u = unstretch_axes(a1, a2, s)
    # force isotropy in case we unstretched too much
    a1u, a2u = force_isotropy(a1, a2)

    if a1u/a2u < ϵ # check whether fabric is destroyed or not
        # Fill FSE
        FSE = FiniteStrainEllipsoid(
            λ11, # vx1
            λ21, # vx2
            λ12, # vy1
            λ22, # vy2
            a1u,
            a2u,
        )

    else
        # destruction of the fabric -> isotropic aggregate
        F = @SMatrix [1.0 0.0; 0.0 1.0]
        # Fill FSE
        FSE = FiniteStrainEllipsoid(
            1.0, # vx1
            0.0, # vx2
            0.0, # vy1
            1.0, # vy2
            1.0, # a1
            1.0, # a2
        )

    end

    return FSE, F

end

function _eigen(A)
    eigval::SVector{2, Float64}, evect::SMatrix{2, 2, Float64} = eigen(A) # get length of FSE semi-axes and orientation
    # eigval, evect = eigen(A) # get length of FSE semi-axes and orientation
    imax, imin = eigval_order(eigval) # get the right order of the semi-axis length
    return eigval, evect, imax, imin
end

#=
Unstrecth the FSE semi-axes by a factor of s.
The area of the ellipse must remain the same
after unstretching
=#
function unstretch_axes(a1, a2, s)
    # a2, a1 are the old semiaxes; s is the unstretching factor
    r0 = a2/a1 # old ratio
    r = r0 + (1-r0)*s # new aspect ratio
    A = π*a2*a1 # Area of the FSE
    # find new unstrecthed axes such that the
    # aspect ratio is r and A remains untouched
    a1u = √(A/(π*r))
    a2u = r*a1u

    return a1u, a2u
end

#=
Recover the new F tensor after unstretching the FSE semi-axes.
We know that F = LR, where L is the left-stretch tensor and R⁻¹=Rᵀ
is the orthogonal rotation tensor. First we recover the new L using
the old eigenvalues (we do not want to rotate the FSE, only unstretch
it) and new eigenvalues:
    L = PDP⁻¹
where P = [λ₁ λ₃ λ₃] and D = I⋅[a₁ a₃ a₃]ᵀ. R is obtained using the
old F and L0 tensors:
    R = F\L0
and finally we obtain the new deformation gradient tensor:
    Fn = LR
=#
function recover_F(F, L0, λ11, λ12, λ21, λ22, a1, a2)
    P = @SMatrix [
        λ11 λ21
        λ12 λ22
    ]
    D = @SMatrix [
        a1*a1  0.0
        0.0    a2*a2
    ]
    # recover left-stretch tensor
    L = P*D\P
    # rotation tensor R = F\(F*Fᵀ) = F\l0
    R = F\L0
    # unstretched tensor of deformation
    Fu = L*R
end

function volume_integral(V, EL2NOD, θ, r)
    ## MODEL PARAMETERS
    nip = 6
    nnodel = size(EL2NOD,1)
    nnodel = 6
    nel = size(EL2NOD,2)

    ## INTEGRATION POINTS & DERIVATIVES wrt LOCAL COORDINATES
    ni, nn, nn3 = Val(nip), Val(nnodel), Val(3)
    N, ∇N,_ , w_ip = _get_SF(ni,nn)
    N3, b, dN3ds,c = _get_SF(ni,nn3)

    ## LOOP OVER ELEMENETS
    integral = zeros(nel)
    Threads.@threads for iel in 1:nel
        #= NOTE: For triangular elements with non-curved edges the Jacobian is
              the same for all integration points (i.e. calculated once
              before the integration loop). Further, linear 3-node shape are
              sufficient to calculate the Jacobian. =#

        # Polar coordinates of element nodes
        θ_el = @SVector [θ[i,iel] for i in 1:3]
        r_el = @SVector [r[i,iel] for i in 1:3]
        iel_coords = SMatrix{3,2}([θ_el r_el])

        # Jacobian n. 1 (p:=polar, l:=local): reference element --> current element
        J_pl = dN3ds*iel_coords
        detJ_pl = det(J_pl) # fast |-> unrolled

        # Jacobian ∂ξ∂θ to transform local (ξ, η) into global (θ,r) derivatives
        #     ∂ξ∂θ = [ R_31    -R_21
        #             -Th_31    Th_21] / detJa_PL
        # ∂ξ∂θ = @SMatrix [(r_el[3]-r_el[1])/detJ_pl -(r_el[2]-r_el[1])/detJ_pl
        #                 -(θ_el[3]-θ_el[1])/detJ_pl  (θ_el[2]-θ_el[1])/detJ_pl]

        idx = view(EL2NOD, 1:6, iel)
        V_el = SVector{6}(view(V, idx))

        ## INTEGRATION LOOP
        @inbounds for ip in 1:nip
            # Polar coordinates of the integration points
            r_ip = mydot(r_el,N3[ip])

            # Integration weight
            ω = r_ip*detJ_pl*w_ip[ip]

            # field at integration point
            V_ip = mydot(V_el,N[ip])

            # Integral reduction
            integral[iel] += V_ip*ω
        end

    end # end block loop

    return sum(integral)

end # END OF ASSEMBLY FUNCTION

function volume_integral_cartesian(V, gr)
    ## MODEL PARAMETERS
    EL2NOD = gr.e2n
    x = gr.x[gr.e2n]
    z = gr.z[gr.e2n]
    nip = 6
    nnodel = size(EL2NOD,1)
    nnodel = 6
    nel = size(EL2NOD,2)

    ## INTEGRATION POINTS & DERIVATIVES wrt LOCAL COORDINATES
    ni, nn, nn3 = Val(nip), Val(nnodel), Val(3)
    N, ∇N,_ , w_ip = _get_SF(ni,nn)
    _, b, dN3ds,c = _get_SF(ni,nn3)

    ## LOOP OVER ELEMENETS
    integral = zeros(nel)
    Threads.@threads for iel in 1:nel
        #= NOTE: For triangular elements with non-curved edges the Jacobian is
              the same for all integration points (i.e. calculated once
              before the integration loop). Further, linear 3-node shape are
              sufficient to calculate the Jacobian. =#

        # Polar coordinates of element nodes
        x_el = @SVector [x[i,iel] for i in 1:3]
        z_el = @SVector [z[i,iel] for i in 1:3]
        iel_coords = SMatrix{3,2}([x_el z_el])

        # Jacobian n. 1 (p:=polar, l:=local): reference element --> current element
        J_pl = dN3ds*iel_coords
        detJ_pl = det(J_pl) # fast |-> unrolled

        idx = view(EL2NOD, 1:6, iel)
        V_el = SVector{6}(view(V, idx))

        ## INTEGRATION LOOP
        @inbounds for ip in 1:nip
            # Integration weight
            ω = detJ_pl*w_ip[ip]

            # field at integration point
            V_ip = mydot(V_el, N[ip])
            # Integral reduction
            integral[iel] += V_ip*ω
        end

    end # end block loop

    return sum(integral)

end # END OF ASSEMBLY FUNCTION

max_stars_repo_name: albert-de-montserrat/Persephone

max_issues_repo_name: albert-de-montserrat/Persephone

max_forks_repo_name: albert-de-montserrat/Persephone


text:
const WORKLOAD_MAP = Dict(
    "test_vgg" => AutoTM.Experiments.test_vgg(),
    "large_vgg" => AutoTM.Experiments.large_vgg(),
    "large_inception" => AutoTM.Experiments.large_inception(),
    "large_resnet" => AutoTM.Experiments.large_resnet(),
    "large_densenet" => AutoTM.Experiments.large_densenet(),
)

workloads() = keys(WOAKLOAD_MAP)

const OPTIMIZER_MAP = Dict(
    "1lm" => AutoTM.Optimizer.Synchronous(185_000_000_000),
    "2lm" => AutoTM.Optimizer.Optimizer2LM(),
)

struct AutoTMParams
    workload::String
    cache::Any
    # 1LM or 2LM
    mode::String
    # Bool flag to use scratchpad
    use_2lm_scratchpad::Bool
    pipe::Any
end

getworkload(A::AutoTMParams) = WORKLOAD_MAP[A.workload]
getoptimizer(A::AutoTMParams) = OPTIMIZER_MAP[A.mode]

maybeunwrap(x) = x
maybeunwrap(x::Tuple) = first(x)

# Hoist into a function so GC works correctly
function run(backend, params::AutoTMParams, measurements)
    # Get workload and optimizer.
    f = getworkload(params)
    opt = getoptimizer(params)

    # Build up a list of optional keyword arguments.
    kw = []
    if params.use_2lm_scratchpad
        push!(kw, :use_scratchpad => true)
    end

    # Instantiate the ngraph function to run - go through `factory` so the node affinity
    # optimization pass runs

    fex = AutoTM.Optimizer.factory(
        backend,
        f,
        opt;
        cache = cache,
        kw...
    )

    # Sometimes, we might get a tuple back from factory - unpack it so we just get the
    # FlusExecutable.
    fex = maybeunwrap(fex)

    # Run the function once to warm it up.
    @time fex()

    # Reuse the initial compilation for all counter sets
    for counter_set in parsed_args["counter_type"]
        println("Counter Set: $counter_set")
        filepath = make_filename(
            workload,
            parsed_args["mode"],
            counter_set,
            parsed_args["use_2lm_scratchpad"]
        )

        # Configure `counters.jl`
        println(pipe, "filepath $filepath")
        println(pipe, "counters $(counter_set)")

        # Invoke the sampling subprocess
        println("Invoking Sampler")
        println(pipe, "start")

        # Sleep for 5 seconds to give the subprocess time to startup and invoke the
        # FluxExecutable again.
        sleep(5)
        @time fex()
        sleep(5)

        # Signal worker process that we're done by writing to the NamedPipe
        println(pipe, "stop")
        println("Successfully Ran Worker Process")
    end
    return nothing
end

max_stars_repo_name: darchr/AutoTM

max_issues_repo_name: darchr/AutoTM

max_forks_repo_name: darchr/AutoTM


text:
call(f, xs...) = f(xs...)

function optimiser(ps, fs...)
  ps = [Param(p) for p in ps]
  fs = map(ps) do p
    os = map(f -> f(p), fs)
    () -> foreach(call, os)
  end
  () -> foreach(call, fs)
end

"""
    SGD(params, η = 1; decay = 0)

Classic gradient descent optimiser. For each parameter `p` and its
gradient `δp`, this runs `p -= η*δp`.

Supports decayed learning rate decay if the `decay` argument is provided.
"""
SGD(ps, η = 1; decay = 0) =
  optimiser(ps, p -> invdecay(p, decay), p -> descent(p, η))

"""
    Momentum(params, ρ, decay = 0)

SGD with momentum `ρ` and optional learning rate decay.
"""
Momentum(ps, ρ; decay = 0) =
  optimiser(ps, p -> momentum(p, ρ), p -> invdecay(p, decay), p -> descent(p, 1))

"""
    Nesterov(params, ρ, decay = 0)

SGD with Nesterov momentum `ρ` and optional learning rate decay.
"""
Nesterov(ps, ρ; decay = 0) =
  optimiser(ps, p -> nesterov(p, ρ), p -> invdecay(p, decay), p -> descent(p, 1))

"""
    RMSProp(params; η = 0.001, ρ = 0.9, ϵ = 1e-8, decay = 0)

[RMSProp](http://www.cs.toronto.edu/~tijmen/csc321/slides/lecture_slides_lec6.pdf)
optimiser. Parameters other than learning rate don't need tuning. Often a good
choice for recurrent networks.
"""
RMSProp(ps, η = 0.001; ρ = 0.9, ϵ = 1e-8, decay = 0) =
  optimiser(ps, p -> rmsprop(p; η = η, ρ = ρ, ϵ = ϵ), p -> invdecay(p, decay), p -> descent(p, 1))

"""
    ADAM(params; η = 0.001, β1 = 0.9, β2 = 0.999, ϵ = 1e-08, decay = 0)

[ADAM](https://arxiv.org/abs/1412.6980v8) optimiser.
"""
ADAM(ps, η = 0.001; β1 = 0.9, β2 = 0.999, ϵ = 1e-08, decay = 0) =
  optimiser(ps, p -> adam(p; η = η, β1 = β1, β2 = β2, ϵ = ϵ), p -> invdecay(p, decay), p -> descent(p, 1))

"""
    ADAGrad(params; η = 0.01, ϵ = 1e-8, decay = 0)

[ADAGrad](http://www.jmlr.org/papers/volume12/duchi11a/duchi11a.pdf) optimiser.
Parameters don't need tuning.
"""
ADAGrad(ps; η = 0.01, ϵ = 1e-8, decay = 0) =
  optimiser(ps, p -> adagrad(p; η = η, ϵ = ϵ), p -> invdecay(p, decay), p -> descent(p, 1))

"""
    ADADelta(params; η = 0.01, ρ = 0.95, ϵ = 1e-8, decay = 0)

[ADADelta](http://arxiv.org/abs/1212.5701) optimiser. Parameters don't need
tuning.
"""
ADADelta(ps; η = 0.01, ρ = 0.95, ϵ = 1e-8, decay = 0) =
  optimiser(ps, p -> adadelta(p; ρ = ρ, ϵ = ϵ), p -> invdecay(p, decay), p -> descent(p, 1))


text:
@enum NODETYPE nature=0 player=1 terminal=2
type GameTree{IntType<:Integer,FloatType<:Real}
    numnodes::IntType
    root::IntType
    numchancehistories::IntType
    numterminalhistories::IntType

    # individual node information
    nodetype::Array{NODETYPE}
    playerid::Array{IntType} # 0 is nature, players are 1,2,...
    nodeutility::Array{IntType}
    nodeparentnode::Array{IntType}
    nodeinfoset::Array{IntType}
    nodenatureprobability::Array{FloatType}
    nodenames::Array{String}
    nodechildren::Array{IntType,2}
end


text:
module Benitor

function fibonacci(n::Int)
    if n <= 1
        return n
    end
    return fibonacci(n-1) + fibonacci(n-2)
end

end # module


text:
function force(i::Int, j::Int)
  if i == j
    return 0
  end

  if formulation == "VanillaIce"
    force_hydro =
      -Z(i) * Z(j) / y(i) / y(j) * (sph[i].P + sph[j].P) * gradW_sym(i, j)
  elseif formulation == "Lagrangian"
    force_hydro =
      -Z(i) * Z(j) * (
        grad_h(i) * sph[i].P / y(i)^2 * gradW(i, j) +
        grad_h(j) * sph[j].P / y(j)^2 * gradW(i, j, j)
      )
  else
    @assert 0
  end
  force_visc = -sph[i].m * sph[j].m * PI_visc(i, j) * gradW_sym(i, j)
  return force_hydro + force_visc
end

function dU(i::Int, j::Int)
  if i == j
    return 0
  end

  if formulation == "VanillaIce"
    dU_hydro =
      sph[i].P * Z(i) * Z(j) / y(i) / y(j) * vij(i, j) * gradW_sym(i, j)
  elseif formulation == "Lagrangian"
    dU_hydro =
      grad_h(i) * sph[i].P * Z(i) * Z(j) / y(i)^2 * vij(i, j) * gradW(i, j)
  else
    @assert 0
  end
  dU_visc =
    0.5 * sph[i].m * sph[j].m * PI_visc(i, j) * vij(i, j) * gradW_sym(i, j)
  return dU_hydro + dU_visc
end

function grad_h(i)
  invf = 1 + sph[i].hsml / y(i) * sph[i].dydh
  return 1/invf
end

function PI_visc(i, j)
  if eij(i, j) * vij(i, j) > 0
    return 0
  else
    if !time_dependent_viscosity
      alpha = 1
      beta = 2
    else
      alpha = (sph[i].alpha + sph[j].alpha) / 2
      beta = 2 * alpha
    end
    epsilon = 0.01
    hbar = (sph[i].hsml + sph[j].hsml) / 4
    cbar = (cs(i) + cs(j)) / 2
    rhobar = (sph[i].rho + sph[j].rho) / 2
    mu = (hbar * dist(i, j) * eij(i, j) * vij(i, j)) / (dist(i, j)^2 + epsilon * hbar^2)
    PI_visc = (-alpha * cbar * mu + beta * mu^2) / rhobar
    return PI_visc
  end
end

function divergence_v(i::Int, j::Int)
  if i == j
    return 0
  else
    return -Z(j) / y(j) * vij(i, j) * eij(i, j) * gradW_sym(i, j)
  end
end

function dalpha_visc(i, divv)
  S = max(-divv * (alpha_max - sph[i].alpha), 0)
  tau = sph[i].hsml / 0.2 / cs(i)
  dalpha_visc = -(sph[i].alpha - alpha_min) / tau + S
  return dalpha_visc
end


text:
"""
    CertainValue

A simple wrapper type for values with no uncertainty (i.e. represented by a scalar).

## Examples

The two following ways of constructing values without uncertainty are equivalent.

```julia
u1, u2 = CertainValue(2.2), CertainValue(6)
w1, w2 = UncertainValue(2.2), UncertainValue(6)
```
"""
struct CertainValue{T} <: AbstractUncertainValue
    value::T
end

Broadcast.broadcastable(x::CertainValue) = Ref(x.value)

function summarise(uval::CertainValue)
    _type = typeof(uval)
    val = uval.value
    "$_type($val)"
end
Base.show(io::IO, uval::CertainValue) = print(io, summarise(uval))

eltype(v::CertainValue{T}) where {T} = T

Base.size(x::CertainValue) = ()
Base.size(x::CertainValue,d) = convert(Int,d)<1 ? throw(BoundsError()) : 1
Base.axes(x::CertainValue) = ()
Base.axes(x::CertainValue,d) = convert(Int,d)<1 ? throw(BoundsError()) : Base.OneTo(1)
Base.ndims(x::CertainValue) = 0
Base.ndims(::Type{<:CertainValue}) = 0
Base.length(x::CertainValue) = 1
Base.firstindex(x::CertainValue) = 1
Base.lastindex(x::CertainValue) = 1
Base.IteratorSize(::Type{<:CertainValue}) = Base.HasShape{0}()
Base.keys(::CertainValue) = Base.OneTo(1)
Base.getindex(x::CertainValue) = x

function Base.getindex(x::CertainValue, i::Integer)
    Base.@_inline_meta
    @boundscheck i == 1 || throw(BoundsError())
    x
end
function Base.getindex(x::CertainValue, I::Integer...)
    Base.@_inline_meta
    @boundscheck all([i == 1 for i in I]) || throw(BoundsError())
    x
end

Base.first(x::CertainValue) = x
Base.last(x::CertainValue) = x
Base.copy(x::CertainValue) = x

Base.minimum(v::CertainValue) = v.value
Base.maximum(v::CertainValue) = v.value
Base.isnan(x::CertainValue) = Base.isnan(x.value)
Base.abs2(x::CertainValue) = Base.abs2(x.value)

StatsBase.mean(v::CertainValue) = v.value
StatsBase.median(v::CertainValue) = v.value
StatsBase.middle(v::CertainValue) = v.value
StatsBase.quantile(v::CertainValue, q) = v.value
StatsBase.quantile(v::CertainValue, q, n::Int) = v.value
StatsBase.std(v::CertainValue{T}) where {T} = zero(T)

Base.rand(v::CertainValue) = v.value
Base.rand(v::CertainValue{T}, n::Int) where T = repeat([v.value], n)

Base.float(v::CertainValue) = float(v.value)

function Base.:<(x::CertainValue{T1}, y::CertainValue{T2}) where {
        T1 <: Real, T2 <: Real}
    x.value < y.value
end

function IntervalArithmetic.interval(x::CertainValue{T1}, y::CertainValue{T2}) where {
        T1 <: Real, T2 <: Real}
    interval(x.value, y.value)
end

export
CertainValue,
UncertainValue


text:
"""
    GraphEdge

# Examples
```julia
function draw(;position_1, position_2, line_width=5)
    setline(line_width)
    line(position_1, position_2, action=:stroke)
    return O
end
function e(g, node1, node2, attr)
    return g[node1][node2]
end

# g mimics an adjacency list with edge weights
g = [Dict(2=>5),
     Dict(1=>3)]
ga = GraphAnimation(g, true, 100, 100, O; get_edge_attribute=e)
node1 = Object(1:30, GraphNode(ga, 1; animate_on=:scale))
node2 = Object(10:40, GraphNode(ga, 2; animate_on=:scale))
edge = Object(30:60, GraphEdge(ga, 1, 2, draw; animate_on=:scale, property_style_map=Dict(""=>:line_width)))
render(video; pathname="graph_node.gif")
```
"""
struct GraphEdge <: AbstractGraphEdge
    from_vertex::AbstractGraphVertex
    to_vertex::AbstractGraphVertex
    object::AbstractObject
    animate_on::Symbol
    style_property_map::Dict{Symbol, Any}
    opts::Dict{Symbol, Any}
end

Base.convert(::Type{Pair}, x::AbstractGraphEdge) = Pair(x.from_vertex.vertex_id, x.to_vertex.vertex_id)

GRAPH_EDGES = Vector{AbstractGraphEdge}()
CURRENT_GRAPH_EDGE = Array{AbstractGraphEdge, 1}()

# """
#     GraphEdge(graph::AbstractObject, from_node::Int, to_node::Int, draw::Function; <keyword arguments>)

# Create a graph edge with additional drawing function and options.

# # Arguments
# - `graph::AbstractObject`: The graph created using [`Graph`](@ref) to which this edge should be added to.
# - `from_node::Int`
# - `to_node::Int`
# - `draw::Function`: The drawing function used to draw the edge.
#     - Implementing the drawing function in a special way to expose the drawing parameters helps in better animation.

# # Keywords
# - `animate_on::Symbol`: Control the animation effect on the edges using pre-defined drawing parameters
#     - For edges, it can be
#         - `:opacity`
#         - `:line_width`
#         - `:length`
#     - For known graph types, it can additionally be
#         - `:color`
#         - `:weights`: only possible for weighted graphs i.e. `SimpleWeightedGraphs`.
# - `properties_to_style_map::Dict{Any,Symbol}`: A mapping to of how edge attributes map to edge drawing styles.
# """
GraphEdge(from_vertex, to_vertex, args...; kwargs...) =
    GraphEdge(CURRENT_GRAPH[1], from_vertex, to_vertex, args...; kwargs...)

function GraphEdge(
    jg::JGraph,
    from_vertex::AbstractGraphVertex,
    to_vertex::AbstractGraphVertex,
    args...;
    kwargs...
)
    return GraphEdge(jg, from_vertex.vertex_id, to_vertex.vertex_id, args...; kwargs...)
end

function GraphEdge(
    jg::JGraph,
    from_vertex::Integer,
    to_vertex::Integer,
    frames;
    kwargs...
)
    check_frames_within(jg.object.frames, frames, from_vertex => to_vertex)
    object = Object(frames, get_draw(:edge); kwargs...)
    return GraphEdge(jg, from_vertex, to_vertex, object; kwargs...)
end

function GraphEdge(
    jg::JGraph,
    from_vertex::Integer,
    to_vertex::Integer,
    object::AbstractObject;
    animate_on::Symbol = :opacity,
    style_property_map::Dict{Any,Symbol} = Dict{Any,Symbol}()
)
    if jg.mode == :static
        if has_edge(jg.graph.adjacency_graph, from_vertex, to_vertex)
            @warn "Edge $(from_vertex)=>$(to_vertex) is already created on canvas. Recreating it will leave orphan edges in the animation. To undo, call `rem_edge!`"
        end
        object.opts[:_graph_idx] = jg.object.opts[:_graph_idx]
        object.opts[:_edge_idx] = length(GRAPH_EDGES) + 1
        opts = Dict{Symbol, Any}()
        opts[:styles] = OrderedDict{Symbol, Function}()
        object.opts[:_style_opts_cache] = Dict{Symbol, Any}()
        add_edge!(jg.graph.adjacency_graph, from_vertex, to_vertex)
        set_prop!(jg.graph.adjacency_graph, from_vertex, to_vertex, length(jg.ordering)+1)
        if !is_directed(jg.graph.adjacency_graph)
            set_prop!(jg.graph.adjacency_graph, to_vertex, from_vertex, length(jg.ordering)+1)
        end
        if from_vertex == to_vertex
            object.opts[:self_loop] = true
        elseif is_directed(jg.graph.adjacency_graph) && has_edge(jg.graph.adjacency_graph, to_vertex, from_vertex)
            object.opts[:loop] = true
            other_edge = jg.ordering[get_prop(jg.graph.adjacency_graph, to_vertex, from_vertex)]
            other_edge.object.opts[:loop] = true
        end
        from_vertex_object = jg.ordering[get_prop(jg.graph.adjacency_graph, from_vertex)]
        to_vertex_object = jg.ordering[get_prop(jg.graph.adjacency_graph, to_vertex)]
        edge = GraphEdge(from_vertex_object, to_vertex_object, object, animate_on, style_property_map, opts)
        push!(GRAPH_EDGES, edge)
        push!(jg.ordering, edge)
        if isempty(CURRENT_GRAPH_EDGE)
            push!(CURRENT_GRAPH_EDGE, edge)
        else
            CURRENT_GRAPH_EDGE[1] = edge
        end
        return edge
    elseif jg.mode == :dynamic
    end
end


text:
@deprecate point_in_domain(basis::Dictionary, scalar) affine_point_in_domain(basis, scalar) false

# Interpolate linearly between the extrema of the support, using the value 0 <= scalar <= 1
affine_point_in_domain(basis::Dictionary, scalar) = affine_point_in_domain(support(basis), scalar)

function affine_point_in_domain(domain::Domain, scalar)
    # Try to find an interval within the support of the basis
    a = infimum(domain)
    b = supremum(domain)

    if isinf(norm(a)) || isinf(norm(b))
        va = MVector(a)
        vb = MVector(b)
        # Avoid infinities for some bases on the real line
        for i in 1:length(va)
            if isinf(va[i])
                va[i] = -1
            end
            if isinf(vb[i])
                vb[i] = 1
            end
        end
        a = SVector(va)
        b = SVector(vb)
    end
    x = (1-scalar) * a + scalar * b
end

function affine_point_in_domain(domain::Domain{<:Number}, scalar)
    # Try to find an interval within the support of the basis
    a = infimum(domain)
    b = supremum(domain)

    T = eltype(domain)
    # Avoid infinities for some bases on the real line
    if isinf(a)
        a = -T(1)
    end
    if isinf(b)
        b = T(1)
    end
    x = (1-scalar) * a + scalar * b
end

function affine_point_in_domain(domain::MappedDomain, scalar)
    m = forward_map(domain)
    m(affine_point_in_domain(superdomain(domain), scalar))
end

Base.rationalize(x::SVector{N,Float64}) where {N} = SVector{N,Rational{Int}}([rationalize(x_i) for x_i in x])
Base.rationalize(x::SVector{N,BigFloat}) where {N} = SVector{N,Rational{BigInt}}([rationalize(x_i) for x_i in x])

function Delimit(s::AbstractString)
    println()
    println("############")
    println("# ",s)
    println("############")
end

function delimit(s::AbstractString)
    println()
    println("## ",s)
end

# Abuse affine_point_in_domain with a scalar greater than one in order to get
# a point outside the domain.
point_outside_domain(basis::Dictionary) = affine_point_in_domain(basis, convert(prectype(basis), 1.1))

point_outside_domain(basis::Laguerre) = -one(domaintype(basis))
point_outside_domain(basis::Hermite) = one(domaintype(basis))+im

function random_point_in_domain(basis::Dictionary)
    T = prectype(basis)
    w = one(T) * rand()
    p = affine_point_in_domain(basis, w)
    in_support(basis, p) ? p : random_point_in_domain(basis)
end

function fixed_point_in_domain(basis::Dictionary)
    T = prectype(basis)
    w = 1/sqrt(T(2))
    affine_point_in_domain(basis, w)
end

random_index(dict::Dictionary) = 1 + floor(Int, rand()*length(dict))

widen_type(::Type{T}) where {T <: Number} = widen(T)
widen_type(::Type{Tuple{A,B}}) where {A,B} = Tuple{widen(A),widen(B)}
widen_type(::Type{Tuple{A,B,C}}) where {A,B,C} = Tuple{widen(A),widen(B),widen(C)}

test_tolerance(::Type{T}) where {T <: AbstractFloat} = sqrt(eps(T))
test_tolerance(::Type{T}) where {T <: Number} = test_tolerance(float(T))
test_tolerance(::Type{Complex{T}}) where {T <: Number} = test_tolerance(T)
test_tolerance(::Type{T}) where {T} = test_tolerance(prectype(T))


text:
import Base.(^)
function (^){N <: Number}(x::Var, n::N)
    y = map(v -> v^n, x.data)
    function df{T}(gy::UniArray{T})
        if hasgrad(x)
            x.data += n * map(v -> v^(n-1), x.data) * gy
        end
    end
    Var(y, [x], ^, df)
end

max_stars_repo_name: hshindo/JukaiNLP

max_issues_repo_name: hshindo/JukaiNLP

max_forks_repo_name: hshindo/JukaiNLP


text:
include("../model/VelocityControlledDifferentialDrive.jl")
using Base.Test

let state = initial_state(1, 0, π/2)
	cycle_time = 0.1
	angular_velocity = 5.0
	expected_position_x = 1.0
	expected_position_y = 0.0165
	ε_position = 0.001
	expected_orientation = 1.5707963
	ε_orientation = 0.01
	control = [angular_velocity angular_velocity]

	state = new_state(state, control, cycle_time)

	@test state[1] ≈ expected_position_x atol = ε_position
	@test state[2] ≈ expected_position_y atol = ε_position
	@test state[3] ≈ expected_orientation atol = ε_orientation
end

let state = initial_state(1, 0, π/2)
	cycle_time = 0.1
	angular_velocity = 5.0
	expected_position_y = 0.0165
	expected_position_x = 1.0
	ε_position = 0.001
	expected_orientation = 1.5707963
	ε_orientation = 0.01
	control = [angular_velocity angular_velocity]

	state = new_state(state, control, cycle_time)

	println()
	println("Constant linear move results")
	println(state)
	@test state[1] ≈ expected_position_x atol = ε_position
	@test state[2] ≈ expected_position_y atol = ε_position
	@test state[3] ≈ expected_orientation atol = ε_orientation
	println("Constant linear move pass")
end

let state = initial_state(1, 0, π/2)
	cycle_time = 0.1
	angular_velocity_left = 5.5
	angular_velocity_right = 6.5
	expected_position_x = 0.99974871
	expected_position_y = 0.0197979
	ε_position = 0.001
	expected_orientation = 1.596181
	ε_orientation = 0.01
	control = [angular_velocity_left angular_velocity_right]

	state = new_state(state, control, cycle_time)

	println()
	println("Constant circular move long cycle results")
	println("Expected:   ", [expected_position_x expected_position_y expected_orientation])
	println("Calculated: ", state)
	println("Error[mm]:  ", sqrt((expected_position_x - state[1])^2 + (expected_position_y - state[2])^2) * 1000)
	@test state[1] ≈ expected_position_x atol = ε_position
	@test state[2] ≈ expected_position_y atol = ε_position
	@test state[3] ≈ expected_orientation atol = ε_orientation
	println("Constant circular move long cycle pass")
end

let state = initial_state(1, 0, π/2)
	cycle_time = 0.1
	angular_velocity_left = 4.5
	angular_velocity_right = 6.5
	expected_position_x = 0.99953937
	expected_position_y = 0.0181422
	ε_position = 0.001
	expected_orientation = 1.621566
	ε_orientation = 0.01
	control = [angular_velocity_left angular_velocity_right]

	state = new_state(state, control, cycle_time)

	println()
	println("Constant circular move long cycle results")
	println("Expected:   ", [expected_position_x expected_position_y expected_orientation])
	println("Calculated: ", state)
	println("Error[mm]:  ", sqrt((expected_position_x - state[1])^2 + (expected_position_y - state[2])^2) * 1000)
	@test state[1] ≈ expected_position_x atol = ε_position
	@test state[2] ≈ expected_position_y atol = ε_position
	@test state[3] ≈ expected_orientation atol = ε_orientation
	println("Constant circular move long cycle high curvature pass")
end

let state = initial_state(1, 0, π/2)
	cycle_time = 0.1
	angular_velocity_left = 5.5
	angular_velocity_right = 6.5
	expected_position_x = 0.944136
	expected_position_y = 0.289875
	ε_position = 0.001
	expected_orientation = 1.95157
	ε_orientation = 0.01
	control = [angular_velocity_left angular_velocity_right]

	for i in 1:15
		state = new_state(state, control, cycle_time)
	end

	println()
	println("Constant circular move multiple update long cycle results")
	println("Expected:   ", [expected_position_x expected_position_y expected_orientation])
	println("Calculated: ", state)
	println("Error[mm]:  ", sqrt((expected_position_x - state[1])^2 + (expected_position_y - state[2])^2) * 1000)
	@test state[1] ≈ expected_position_x atol = ε_position
	@test state[2] ≈ expected_position_y atol = ε_position
	@test state[3] ≈ expected_orientation atol = ε_orientation
	println("Constant circular move multiple update long cycle pass")
end

let state = initial_state(1, 0, π/2)
	cycle_time = 0.05
	angular_velocity_left = 5.5
	angular_velocity_right = 6.5
	expected_position_x = 0.999937174
	expected_position_y = 0.00989973
	ε_position = 0.001
	expected_orientation = 1.583489
	ε_orientation = 0.01
	control = [angular_velocity_left angular_velocity_right]

	state = new_state(state, control, cycle_time)

	println()
	println("Constant circular move short cycle results")
	println("Expected:   ", [expected_position_x expected_position_y expected_orientation])
	println("Calculated: ", state)
	println("Error[mm]:  ", sqrt((expected_position_x - state[1])^2 + (expected_position_y - state[2])^2) * 1000)
	@test state[1] ≈ expected_position_x atol = ε_position
	@test state[2] ≈ expected_position_y atol = ε_position
	@test state[3] ≈ expected_orientation atol = ε_orientation
	println("Constant circular move short cycle pass")
end

let state = initial_state(1, 0, π/2)
	cycle_time = 0.05
	angular_velocity_left = 5.5
	angular_velocity_right = 6.5
	expected_position_x = 0.944136
	expected_position_y = 0.289875
	ε_position = 0.001
	expected_orientation = 1.95157
	ε_orientation = 0.01
	control = [angular_velocity_left angular_velocity_right]

	for i in 1:30
		state = new_state(state, control, cycle_time)
	end

	println()
	println("Constant circular move multiple update short cycle results")
	println("Expected:   ", [expected_position_x expected_position_y expected_orientation])
	println("Calculated: ", state)
	println("Error[mm]:  ", sqrt((expected_position_x - state[1])^2 + (expected_position_y - state[2])^2) * 1000)
	@test state[1] ≈ expected_position_x atol = ε_position
	@test state[2] ≈ expected_position_y atol = ε_position
	@test state[3] ≈ expected_orientation atol = ε_orientation
	println("Constant circular move multiple update short cycle pass")
end

max_stars_repo_name: grgomrton/robot-model

max_issues_repo_name: grgomrton/robot-model

max_forks_repo_name: grgomrton/robot-model


text:
module AnalyzeRegistry

# Standard libraries
using Pkg, TOML, UUIDs
# Third-party packages
using FLoops # for the `@floops` macro
using MicroCollections # for `EmptyVector` and `SingletonVector`
using BangBang # for `append!!`
using LicenseCheck # for `find_license` and `is_osi_approved`
using JSON3 # for interfacing with `tokei` to count lines of code
using Tokei_jll # count lines of code

export general_registry, find_packages, analyze, analyze_from_registry, analyze_from_registry!

include("count_loc.jl")
const LicenseTableEltype=@NamedTuple{license_filename::String, licenses_found::Vector{String}, license_file_percent_covered::Float64}

struct Package
    name::String # name of the package
    uuid::UUID # uuid of the package
    repo::String # URL of the repository
    subdir::String # subdirectory of the package in the repo
    reachable::Bool # can the repository be cloned?
    docs::Bool # does it have documentation?
    runtests::Bool # does it have the test/runtests.jl file?
    github_actions::Bool # does it use GitHub Actions?
    travis::Bool # does it use Travis CI?
    appveyor::Bool # does it use AppVeyor?
    cirrus::Bool # does it use Cirrus CI?
    circle::Bool # does it use Circle CI?
    drone::Bool # does it use Drone CI?
    buildkite::Bool # does it use Buildkite?
    azure_pipelines::Bool # does it use Azure Pipelines?
    gitlab_pipeline::Bool # does it use Gitlab Pipeline?
    license_files::Vector{LicenseTableEltype} # a table of all possible license files
    licenses_in_project::Vector{String} # any licenses in the `license` key of the Project.toml
    lines_of_code::Vector{LoCTableEltype}
end
function Package(name, uuid, repo;
                 subdir="",
                 reachable=false,
                 docs=false,
                 runtests=false,
                 github_actions=false,
                 travis=false,
                 appveyor=false,
                 cirrus=false,
                 circle=false,
                 drone=false,
                 buildkite=false,
                 azure_pipelines=false,
                 gitlab_pipeline=false,
                 license_files=LicenseTableEltype[],
                 licenses_in_project=String[],
                 lines_of_code=Vector{LoCTableEltype}(),
                 )
    return Package(name, uuid, repo, subdir, reachable, docs, runtests, github_actions, travis,
                   appveyor, cirrus, circle, drone, buildkite, azure_pipelines, gitlab_pipeline,
                   license_files, licenses_in_project, lines_of_code)
end

# define `isequal`, `==`, and `hash` just in terms of the fields
for f in (:isequal, :(==))
    @eval begin
        function Base.$f(A::Package, B::Package)
            for i = 1:fieldcount(Package)
                $f(getfield(A, i), getfield(B, i)) || return false
            end
            true
        end
    end
end

Base.hash(A::Package, h::UInt) = hash(:Package, hash(ntuple(i -> getfield(A, i), fieldcount(Package)), h))

function Base.show(io::IO, p::Package)
    body = """
        Package $(p.name):
          * repo: $(p.repo)
          * uuid: $(p.uuid)
          * is reachable: $(p.reachable)
        """
    if p.reachable
        if !isempty(p.lines_of_code)
            body *= """
                  * lines of Julia code in `src`: $(count_julia_loc(p.lines_of_code, "src"))
                  * lines of Julia code in `test`: $(count_julia_loc(p.lines_of_code, "test"))
                """
        end
        if isempty(p.license_files)
            body *= "  * no license found\n"
        else
            lic = p.license_files[1]
            license_string = join(lic.licenses_found, ", ")
            body *= "  * has license(s) in file: $(license_string)\n"
            body *= "    * filename: $(lic.license_filename)\n"
            body *= "    * OSI approved: $(all(is_osi_approved, lic.licenses_found))\n"
        end
        if !isempty(p.licenses_in_project)
            lic_project = join(p.licenses_in_project, ", ")
            body *= "  * has license(s) in Project.toml: $(lic_project)\n"
            body *= "    * OSI approved: $(all(is_osi_approved, p.licenses_in_project))\n"
        end
        body *= """
              * has documentation: $(p.docs)
              * has tests: $(p.runtests)
            """
        ci_services = (p.github_actions => "GitHub Actions",
                       p.travis => "Travis",
                       p.appveyor => "AppVeyor",
                       p.cirrus => "Cirrus",
                       p.circle => "Circle",
                       p.drone => "Drone CI",
                       p.buildkite => "Buildkite",
                       p.azure_pipelines => "Azure Pipelines",
                       p.gitlab_pipeline => "GitLab Pipeline",
                       )
        if any(first.(ci_services))
            body *= "  * has continuous integration: true\n"
            for (k, v) in ci_services
                if k
                    body *= "    * $(v)\n"
                end
            end
        else
            body *= "  * has continuous integration: false\n"
        end
    end
    print(io, strip(body))
end

"""
    general_registry() -> String

Guess the path of the General registry.
"""
general_registry() =
    first([joinpath(d, "registries", "General") for d in Pkg.depots() if isfile(joinpath(d, "registries", "General", "Registry.toml"))])

"""
    find_packages(; registry = general_registry()) -> Vector{String}
    find_packages(names::AbstractString...; registry = general_registry()) -> Vector{String}
    find_packages(names; registry = general_registry()) -> Vector{String}

Find all packages in the given registry (specified by the `registry` keyword argument),
the General registry by default. Return a vector with the paths to the directories
of each package in the registry.

Pass a list of package `names` as the first argument to return the paths corresponding to those packages,
or individual package names as separate arguments.
"""
find_packages

find_packages(names::AbstractString...; registry = general_registry()) =  find_packages(names; registry=registry)

function find_packages(names; registry = general_registry())
    if names !== nothing
        paths = String[]
        for name in names
            path = joinpath(registry, string(uppercase(first(name))), name)
            if isdir(path)
               push!(paths, path)
            else
                @error("Could not find package in registry!", name, path)
            end
        end
        return paths
    end
end

function find_packages(; registry = general_registry())
    # Get the list of packages in the registry by parsing the `Registry.toml`
    # file in the given directory.
    packages = TOML.parsefile(joinpath(registry, "Registry.toml"))["packages"]
    # Get the directories of all packages.  Filter out JLL packages: they are
    # automatically generated and we know that they don't have testing nor
    # documentation.
    return [joinpath(registry, splitpath(p["path"])...) for (_, p) in packages if !endswith(p["name"], "_jll")]
end

"""
    analyze_from_registry!(root, dir::AbstractString) -> Package

Analyze the package whose entry in the registry is in the `dir` directory,
cloning the package code to `joinpath(root, uuid)` where `uuid` is the UUID
of the package, if such a directory does not already exist.

"""
function analyze_from_registry!(root, dir::AbstractString)
    # Parse the `Package.toml` file in the given directory.
    toml = TOML.parsefile(joinpath(dir, "Package.toml"))
    name = toml["name"]::String
    uuid_string = toml["uuid"]::String
    uuid = UUID(uuid_string)
    repo = toml["repo"]::String
    subdir = get(toml, "subdir", "")::String

    dest = joinpath(root, uuid_string)

    isdir(dest) && return analyze(dest; repo, subdir)

    reachable = try
        # Clone only latest commit on the default branch.  Note: some
        # repositories aren't reachable because the author made them private
        # or deleted them.  In these cases git would ask for username and
        # password, provide it with fake values just to move on:
        # https://stackoverflow.com/a/65705346/2442087
        run(pipeline(`git clone -q --depth 1 --config credential.helper='!f() { echo -e "username=git\npassword="; }; f' $(repo) $(dest)`; stderr=devnull))
        true
    catch
        # The repository may be unreachable
        false
    end
    return reachable ? analyze(dest; repo, reachable, subdir) : Package(name, uuid, repo; subdir)
end

"""
    analyze_from_registry!(root, packages::AbstractVector{<:AbstractString}) -> Vector{Package}

Analyze all packages in the iterable `packages`, using threads, cloning them to `root`
if a directory with their `uuid` does not already exist.  Returns a
`Vector{Package}`.
"""
function analyze_from_registry!(root, packages::AbstractVector{<:AbstractString})
    @floop for p in packages
        ps = SingletonVector((analyze_from_registry!(root, p),))
        @reduce(result = append!!(EmptyVector(), ps))
    end
    result
end

"""
    analyze_from_registry(dir::AbstractString) -> Package
    analyze_from_registry(packages::AbstractVector{<:AbstractString}) -> Vector{Package}

Analyzes a package or list of packages using the information in their directory
in a registry by creating a temporary directory and calling `analyze_from_registry!`,
cleaning up the temporary directory afterwards.

## Example
```julia
julia> analyze_from_registry(joinpath(general_registry(), "B", "BinaryBuilder"))
Package BinaryBuilder:
  * repo: https://github.com/JuliaPackaging/BinaryBuilder.jl.git
  * uuid: 12aac903-9f7c-5d81-afc2-d9565ea332ae
  * is reachable: true
  * lines of Julia code in `src`: 4733
  * lines of Julia code in `test`: 1520
  * has license(s) in file: MIT
    * filename: LICENSE.md
    * OSI approved: true
  * has documentation: true
  * has tests: true
  * has continuous integration: true
    * GitHub Actions
    * Azure Pipelines

```
"""
function analyze_from_registry(p)
    mktempdir() do root
        analyze_from_registry!(root, p)
    end
end

function parse_project(dir)
    bad_project = (; name = "Invalid Project.toml", uuid = UUID(UInt128(0)), licenses_in_project=String[])
    project_path = joinpath(dir, "Project.toml")
    if !isfile(project_path)
        project_path = joinpath(dir, "JuliaProject.toml")
    end
    isfile(project_path) || return bad_project
    project = TOML.tryparsefile(project_path)
    project isa TOML.ParserError && return bad_project
    haskey(project, "name") && haskey(project, "uuid") || return bad_project
    uuid = tryparse(UUID, project["uuid"]::String)
    uuid === nothing && return bad_project
    licenses_in_project = get(project, "license", String[])
    if licenses_in_project isa String
        licenses_in_project = [licenses_in_project]
    end
    return (; name = project["name"]::String, uuid, licenses_in_project)
end

"""
    analyze(dir::AbstractString; repo = "", reachable=true, name=nothing, uuid=nothing)

Analyze the package whose source code is located at `dir`. Optionally `repo`
and `reachable` a boolean indicating whether or not the package is reachable online, since
these can't be inferred from the source code. If `name` or `uuid` are `nothing`, the
directories `Project.toml` is parsed to infer the package's name and UUID.

## Example
```julia
julia> analyze(pkgdir(AnalyzeRegistry))
Package AnalyzeRegistry:
  * repo:
  * uuid: e713c705-17e4-4cec-abe0-95bf5bf3e10c
  * is reachable: true
  * lines of Julia code in `src`: 327
  * lines of Julia code in `test`: 58
  * has license(s) in file: MIT
    * filename: LICENSE
    * OSI approved: true
  * has documentation: false
  * has tests: true
  * has continuous integration: true
    * GitHub Actions

```
"""
function analyze(dir::AbstractString; repo = "", reachable=true, subdir="")
    # we will look for docs, tests, license, and count lines of code
    # in the `pkgdir`; we will look for CI in the `dir`.
    pkgdir = joinpath(dir, subdir)
    name, uuid, licenses_in_project = parse_project(pkgdir)
    docs = isfile(joinpath(pkgdir, "docs", "make.jl")) || isfile(joinpath(pkgdir, "doc", "make.jl"))
    runtests = isfile(joinpath(pkgdir, "test", "runtests.jl"))
    travis = isfile(joinpath(dir, ".travis.yml"))
    appveyor = isfile(joinpath(dir, "appveyor.yml"))
    cirrus = isfile(joinpath(dir, ".cirrus.yml"))
    circle = isfile(joinpath(dir, ".circleci", "config.yml"))
    drone = isfile(joinpath(dir, ".drone.yml"))
    azure_pipelines = isfile(joinpath(dir, "azure-pipelines.yml"))
    buildkite = isfile(joinpath(dir, ".buildkite", "pipeline.yml"))
    gitlab_pipeline = isfile(joinpath(dir, ".gitlab-ci.yml"))
    github_workflows = joinpath(dir, ".github", "workflows")
    if isdir(github_workflows)
        # Find all workflows
        files = readdir(github_workflows)
        # Exclude TagBot and CompatHelper
        filter(f -> lowercase(f) ∉ ("compathelper.yml", "tagbot.yml"), files)
        # Assume all other files are GitHub Actions for CI.  May not
        # _always_ be the case, but it's a good first-order approximation.
        github_actions = length(files) > 0
    else
        github_actions = false
    end
    license_files = find_licenses(dir)
    lines_of_code = count_loc(dir)
    Package(name, uuid, repo; subdir, reachable, docs, runtests, travis, appveyor, cirrus,
            circle, drone, buildkite, azure_pipelines, gitlab_pipeline, github_actions,
            license_files, licenses_in_project, lines_of_code)
end

end # module


text:
module AstroImages

using FITSIO
using FileIO
# Rather than pulling in all of Images.jl, just grab the packages
# we need to extend to our basic functionality.
# We also need ImageShow so that user's images appear automatically.
using ImageBase, ImageShow#, ImageAxes

using WCS
using Statistics
using MappedArrays
using ColorSchemes
using DimensionalData
using Tables
using RecipesBase
using AstroAngles
using Printf
using PlotUtils: PlotUtils
using PlotUtils: optimize_ticks, AbstractColorList
using UUIDs # can remove once reigstered with FileIO

export load,
    save,
    AstroImage,
    AstroImageVec,
    AstroImageMat,
    WCSGrid,
    composecolors,
    Zscale,
    Percent,
    logstretch,
    powstretch,
    sqrtstretch,
    squarestretch,
    asinhstretch,
    sinhstretch,
    powerdiststretch,
    imview,
    render, # deprecated
    header,
    copyheader,
    shareheader,
    wcs,
    Comment,
    History,
    # Dimensions
    Centered,
    Spec,
    Pol,
    Ti,
    X, Y, Z, Dim,
    At, Near, ..,
    dims, refdims,
    recenter,
    pix_to_world,
    pix_to_world!,
    world_to_pix,
    world_to_pix!

# Images.jl expects data to be either a float or a fixed-point number.  Here we define some
# utilities to convert all data types supported by FITS format to float or fixed-point:
#
#   * Float numbers are left as they are
#   * Unsigned integers are mapped to [0, 1] with Normed type
#   * Signed integers are mapped to unsigned integers and then to Normed
_float(x::AbstractFloat) = x
for n in (8, 16, 32, 64)
    SIT = Symbol("Int", n) # signed integer type
    UIT = Symbol("UInt", n) # unsigned integer type
    NIT = Symbol("N0f", n) # fixed-point type for unsigned float
    @eval maxint = $UIT(big(2) ^ ($n - 1)) #
    @eval begin
        _float(x::$UIT) = reinterpret($NIT, x)
        _float(x::$SIT) = _float(xor(reinterpret($UIT, x), $maxint))
    end
end

"""
Provides access to a FITS image along with its accompanying
header and WCS information, if applicable.
"""
struct AstroImage{T,N,D<:Tuple,R<:Tuple,A<:AbstractArray{T,N},W<:Tuple} <: AbstractDimArray{T,N,D,A}
    # Parent array we are wrapping
    data::A
    # Fields for DimensionalData
    dims::D
    refdims::R
    # FITS Heads beloning to this image, if any
    header::FITSHeader
    # cached WCSTransform objects for this data.
    wcs::Vector{WCSTransform}
    # A flag that is set when a user modifies a WCS header.
    # The next access to the wcs object will regenerate from
    # the new header on demand.
    wcs_stale::Base.RefValue{Bool}
    # Correspondance between dims & refdims -> WCS Axis numbers
    wcsdims::W
end
# Provide type aliases for 1D and 2D versions of our data structure.
const AstroImageVec{T,D} = AstroImage{T,1} where {T}
const AstroImageMat{T,D} = AstroImage{T,2} where {T}

"""
    Centered()

Pass centered as a dimesion range to automatically center a dimension
along that axis.

Example:
```julia
cube = load("abc.fits", (X=Centered(), Y=Centered(), Pol=[:I, :Q, :U]))
```

In that case, cube will have dimsions with the centre of the image at 0
in both the X and Y axes.
"""
struct Centered end

# Default dimension names if none are provided
const dimnames = (
    X, Y, Z,
    (Dim{i} for i in 4:10)...
)

const Spec = Dim{:Spec}
const Pol = Dim{:Pol}

"""
    wcsax(img, dim)

Return the WCS axis number associated with a dimension, even if the image
has been slices or otherwise transformed.

Internally, the order is stored in the field `wcsdims`.
"""
function wcsax(img::AstroImage, dim)
    return findfirst(di->name(di)==name(dim), img.wcsdims)
end

# Accessors
"""
    header(img::AstroImage)

Return the underlying FITSIO.FITSHeader object wrapped by an AstroImage.
Note that this object has less flexible getindex and setindex methods.
Indexing by symbol, Comment, History, etc are not supported.
"""
header(img::AstroImage) = getfield(img, :header)
"""
    header(array::AbstractArray)

Returns an empty FITSIO.FITSHeader object when called with a non-AstroImage
abstract array.
"""
header(::AbstractArray) = emptyheader()

"""
    wcs(img)

Computes and returns a list of World Coordinate System WCSTransform objects from WCS.jl.
The resultss are cached after the first call, so subsequent calls are fast.
Modifying a WCS header invalidates this cache automatically, so users should call `wcs(...)`
each time rather than keeping the WCSTransform object around.
"""
function wcs(img::AstroImage)
    if getfield(img, :wcs_stale)[]
        empty!(getfield(img, :wcs))
        append!(getfield(img, :wcs), wcsfromheader(img))
        getfield(img, :wcs_stale)[] = false
    end
    return getfield(img, :wcs)
end
"""
    wcs(img, index)

Computes and returns a World Coordinate System WCSTransform objects from WCS.jl by index.
This is to support files with multiple WCS transforms specified.
`wcs(img,1)` is useful for selecting selecting the first WCSTranform object.
The resultss are cached after the first call, so subsequent calls are fast.
Modifying a WCS header invalidates this cache automatically, so users should call `wcs(...)`
each time rather than keeping the WCSTransform object around.
"""
wcs(img, ind) = wcs(img)[ind]
"""
wcs(array)

Returns a list with a single basic WCSTransform object when called with a non-AstroImage
abstract array.
"""
wcs(arr::AbstractArray) = [emptywcs(arr)]

# Implement DimensionalData interface
Base.parent(img::AstroImage) = getfield(img, :data)
DimensionalData.dims(A::AstroImage) = getfield(A, :dims)
DimensionalData.refdims(A::AstroImage) = getfield(A, :refdims)
DimensionalData.data(A::AstroImage) = getfield(A, :data)
DimensionalData.name(::AstroImage) = DimensionalData.NoName()
DimensionalData.metadata(::AstroImage) = DimensionalData.Dimensions.LookupArrays.NoMetadata()

@inline function DimensionalData.rebuild(
    img::AstroImage,
    data,
    # Fields for DimensionalData
    dims::Tuple=DimensionalData.dims(img),
    refdims::Tuple=DimensionalData.refdims(img),
    name::Union{Symbol,DimensionalData.AbstractName,Nothing}=nothing,
    metadata::Union{DimensionalData.LookupArrays.AbstractMetadata,Nothing}=nothing,
    # FITS Header beloning to this image, if any
    header::FITSHeader=deepcopy(header(img)),
    # A cached WCSTransform object for this data
    wcs::Vector{WCSTransform}=getfield(img, :wcs),
    wcs_stale::Bool=getfield(img, :wcs_stale)[],
    wcsdims::Tuple=(dims...,refdims...),
)
    return AstroImage(data, dims, refdims, header, wcs, Ref(wcs_stale), wcsdims)
end
# Keyword argument version.
# We have to define this since our struct contains additional field names.
@inline function DimensionalData.rebuild(
    img::AstroImage;
    data,
    # Fields for DimensionalData
    dims::Tuple=DimensionalData.dims(img),
    refdims::Tuple=DimensionalData.refdims(img),
    name::Union{Symbol,DimensionalData.AbstractName,Nothing}=nothing,
    metadata::Union{DimensionalData.LookupArrays.AbstractMetadata,Nothing}=nothing,
    # FITS Header beloning to this image, if any
    header::FITSHeader=deepcopy(header(img)),
    # A cached WCSTransform object for this data
    wcs::Vector{WCSTransform}=getfield(img, :wcs),
    wcs_stale::Bool=getfield(img, :wcs_stale)[],
    wcsdims::Tuple=(dims...,refdims...),
)
    return AstroImage(data, dims, refdims, header, wcs, Ref(wcs_stale), wcsdims)
end
@inline DimensionalData.rebuildsliced(
    f::Function,
    img::AstroImage,
    data,
    I,
    header=deepcopy(header(img)),
    wcs=getfield(img, :wcs),
    wcs_stale=getfield(img, :wcs_stale)[],
    wcsdims=getfield(img, :wcsdims),
) = rebuild(img, data, DimensionalData.slicedims(f, img, I)..., nothing, nothing, header, wcs, wcs_stale, wcsdims)

# Return result wrapped in AstroImage
# For these functions that return lazy wrappers, we want to share header
for f in [
    :(Base.adjoint),
    :(Base.transpose),
    :(Base.view)
]
    # TODO: these functions are copying headers
    @eval ($f)(img::AstroImage) = shareheader(img, $f(parent(img)))
end

"""
    AstroImage(img::AstroImage)

Returns its argument. Useful to ensure an argument is converted to an
AstroImage before proceeding.
"""
AstroImage(img::AstroImage) = img

"""
    AstroImage(data::AbstractArray, [header::FITSHeader,] [wcs::WCSTransform,])

Create an AstroImage from an array, and optionally header or header and a
WCSTransform.
"""
function AstroImage(
    data::AbstractArray{T,N},
    dims::Union{Tuple,NamedTuple}=(),
    refdims::Union{Tuple,NamedTuple}=(),
    header::FITSHeader=emptyheader(),
    wcs::Union{WCSTransform,Nothing}=nothing;
    wcsdims=nothing
) where {T, N}
    wcs_stale = isnothing(wcs)
    if isnothing(wcs)
        wcs = [emptywcs(data)]
    end
    # If the user passes in a WCSTransform of their own, we use it and mark
    # wcs_stale=false. It will be kept unless they manually change a WCS header.
    # If they don't pass anything, we start with empty WCS information regardless
    # of what's in the header but we mark it as stale.
    # If/when the WCS info is accessed via `wcs(img)` it will be computed and cached.
    # This avoids those computations if the WCS transform is not needed.
    # It also allows us to create images with invalid WCS header,
    # only erroring when/if they are used.

    # Fields for DimensionalData.
    # TODO: cleanup logic
    if dims == ()
        # if wcsdims
            # ourdims = Tuple(Wcs{i} for i in 1:ndims(data))
        # else
            ourdims = dimnames[1:ndims(data)]
        # end
        dims = map(ourdims, axes(data)) do dim, ax
            dim(ax)
        end
    end
    # Replace any occurences of Centered() with an automatic range
    # from the data.
    dimvals = map(dims, axes(data)) do dim, ax
        if dim isa Centered
            ax .- mean(ax)
        else
            dim
        end
    end
    if dims isa NamedTuple
        dims = NamedTuple{keys(dims)}(dimvals)
    elseif !(dims isa NTuple{N,Dimensions.Dimension} where N) &&
        !(all(d-> d isa Union{UnionAll,DataType} && d <: Dimensions.Dimension, dims))
        k = name.(dimnames[1:ndims(data)])
        dims = NamedTuple{k}(dimvals)
    end
    dims = DimensionalData.format(dims, data)
    if length(dims) != ndims(data)
        error("Number of dims does not match the shape of the data")
    end

    if isnothing(wcsdims)
        wcsdims = (dims...,refdims...)
    end

    return AstroImage(data, dims, refdims, header, wcs, Ref(wcs_stale), wcsdims)
end
function AstroImage(
    darr::AbstractDimArray,
    header::FITSHeader=emptyheader(),
    wcs::Union{Vector{WCSTransform},Nothing}=nothing;
)
    wcs_stale = isnothing(wcs)
    if isnothing(wcs)
        wcs = [emptywcs(darr)]
    end
    wcsdims = (dims(darr)..., refdims(darr)...)
    return AstroImage(parent(darr), dims(darr), refdims(darr), header, wcs, Ref(wcs_stale), wcsdims)
end
AstroImage(
    data::AbstractArray,
    dims::Union{Tuple,NamedTuple},
    header::FITSHeader,
    wcs::Union{Vector{WCSTransform},Nothing}=nothing;
) = AstroImage(data, dims, (), header, wcs)
AstroImage(
    data::AbstractArray,
    header::FITSHeader,
    wcs::Union{Vector{WCSTransform},Nothing}=nothing;
) = AstroImage(data, (), (), header, wcs)

# TODO: ensure this gets WCS dims.
AstroImage(data::AbstractArray, wcs::Vector{WCSTransform}) = AstroImage(data, emptyheader(), wcs)

"""
Index for accessing a comment associated with a header keyword
or COMMENT entry.

Example:
```
img = AstroImage(randn(10,10))
img["ABC"] = 1
img["ABC", Comment] = "A comment describing this key"

push!(img, Comment, "The purpose of this file is to demonstrate comments")
img[Comment] # ["The purpose of this file is to demonstrate comments")]
```
"""
struct Comment end

"""
Allows accessing and setting HISTORY header entries

img = AstroImage(randn(10,10))
push!(img, History, "2023-04-19: Added history entry.")
img[History] # ["2023-04-19: Added history entry."]
"""
struct History end

# We might want getproperty for header access in future.
# function Base.getproperty(img::AstroImage, ::Symbol)
#     error("getproperty reserved for future use.")
# end

# Getting and setting comments
Base.getindex(img::AstroImage, inds::AbstractString...) = getindex(header(img), inds...) # accesing header using strings
function Base.setindex!(img::AstroImage, v, ind::AbstractString)  # modifying header using a string
    setindex!(header(img), v, ind)
    # Mark the WCS object as being out of date if this was a WCS header keyword
    if ind ∈ WCS_HEADERS
        getfield(img, :wcs_stale)[] = true
    end
end
Base.getindex(img::AstroImage, inds::Symbol...) = getindex(img, string.(inds)...) # accessing header using symbol
Base.setindex!(img::AstroImage, v, ind::Symbol) = setindex!(img, v, string(ind))
Base.getindex(img::AstroImage, ind::AbstractString, ::Type{Comment}) = get_comment(header(img), ind) # accesing header comment using strings
Base.setindex!(img::AstroImage, v,  ind::AbstractString, ::Type{Comment}) = set_comment!(header(img), ind, v) # modifying header comment using strings
Base.getindex(img::AstroImage, ind::Symbol, ::Type{Comment}) = get_comment(header(img), string(ind)) # accessing header comment using symbol
Base.setindex!(img::AstroImage,  v, ind::Symbol, ::Type{Comment}) = set_comment!(header(img), string(ind), v) # modifying header comment using Symbol

# Support for special HISTORY and COMMENT entries
function Base.getindex(img::AstroImage, ::Type{History})
    hdr = header(img)
    ii = findall(==("HISTORY"), hdr.keys)
    return view(hdr.comments, ii)
end
function Base.getindex(img::AstroImage, ::Type{Comment})
    hdr = header(img)
    ii = findall(==("COMMENT"), hdr.keys)
    return view(hdr.comments, ii)
end
# Adding new comment and history entries
function Base.push!(img::AstroImage, ::Type{Comment}, history::AbstractString)
    hdr = header(img)
    push!(hdr.keys, "COMMENT")
    push!(hdr.values, nothing)
    push!(hdr.comments, history)
end
function Base.push!(img::AstroImage, ::Type{History}, history::AbstractString)
    hdr = header(img)
    push!(hdr.keys, "HISTORY")
    push!(hdr.values, nothing)
    push!(hdr.comments, history)
end

"""
    copyheader(img::AstroImage, data) -> imgnew
Create a new image copying the header of `img` but
using the data of the AbstractArray `data`. Note that changing the
header of `imgnew` does not affect the header of `img`.
See also: [`shareheader`](@ref).
"""
copyheader(img::AstroImage, data::AbstractArray) =
    AstroImage(data, dims(img), refdims(img), deepcopy(header(img)), copy(getfield(img, :wcs)), Ref(getfield(img, :wcs_stale)[]), getfield(img,:wcsdims))

"""
    shareheader(img::AstroImage, data) -> imgnew
Create a new image reusing the header dictionary of `img` but
using the data of the AbstractArray `data`. The two images have
synchronized header; modifying one also affects the other.
See also: [`copyheader`](@ref).
"""
shareheader(img::AstroImage, data::AbstractArray) = AstroImage(data, dims(img), refdims(img), header(img), getfield(img, :wcs), Ref(getfield(img, :wcs_stale)[]), getfield(img,:wcsdims))
# Share header if an AstroImage, do nothing if AbstractArray
maybe_shareheader(img::AstroImage, data) = shareheader(img, data)
maybe_shareheader(::AbstractArray, data) = data
maybe_copyheader(img::AstroImage, data) = copyheader(img, data)
maybe_copyheader(::AbstractArray, data) = data

Base.promote_rule(::Type{AstroImage{T}}, ::Type{AstroImage{V}}) where {T,V} = AstroImage{promote_type{T,V}}

Base.copy(img::AstroImage) = rebuild(img, copy(parent(img)))
Base.convert(::Type{AstroImage}, A::AstroImage) = A
Base.convert(::Type{AstroImage}, A::AbstractArray) = AstroImage(A)
Base.convert(::Type{AstroImage{T}}, A::AstroImage{T}) where {T} = A
Base.convert(::Type{AstroImage{T}}, A::AstroImage) where {T} = shareheader(A, convert(AbstractArray{T}, parent(A)))
Base.convert(::Type{AstroImage{T}}, A::AbstractArray{T}) where {T} = AstroImage(A)
Base.convert(::Type{AstroImage{T}}, A::AbstractArray) where {T} = AstroImage(convert(AbstractArray{T}, A))
Base.convert(::Type{AstroImage{T,N,D,R,AT}}, A::AbstractArray{T,N}) where {T,N,D,R,AT} = AstroImage(convert(AbstractArray{T}, A))

# TODO: share headers in View. Needs support from DimensionalData.

"""
    emptyheader()

Convenience function to create a FITSHeader with no keywords set.
"""
emptyheader() = FITSHeader(String[],[],String[])

"""
	recenter(img::AstroImage, newcentx, newcenty, ...)

Adjust the dimensions of an AstroImage so that they are centered on the pixel locations given by `newcentx`, .. etc.

This does not affect the underlying array, it just updates the dimensions associated with it by the AstroImage.

Example:
```julia
a = AstroImage(randn(11,11))
a[1,1] # Bottom left
a[At(1),At(1)] # Bottom left
r = recenter(a, 6, 6)
r[1,1] # Still bottom left
r[At(1),At(1)] # Center pixel
```
"""
function recenter(img::AstroImage, centers::Number...)
	newdims = map(dims(img), axes(img), centers) do d, a, c
		return AstroImages.name(d) => a .- c
	end
	newdimsformatted = AstroImages.DimensionalData.format(NamedTuple(newdims), parent(img))
	l = length(newdimsformatted)
	if l < ndims(img)
		newdimsformatted = (newdimsformatted..., dims(img)[l+1:end]...)
	end
	AstroImages.rebuild(img, parent(img), newdimsformatted)
end

include("wcs.jl")
include("io.jl")
include("imview.jl")
include("showmime.jl")
include("plot-recipes.jl")

include("contrib/images.jl")
include("contrib/abstract-ffts.jl")
# include("contrib/reproject.jl")

include("ccd2rgb.jl")
include("precompile.jl")

function __init__()

    # You can only `imview` 2D slices. Add an error hint if the user
    # tries to display a cube.
    if isdefined(Base.Experimental, :register_error_hint)
        Base.Experimental.register_error_hint(MethodError) do io, exc, argtypes, kwargs
            if exc.f == imview && first(argtypes) <: AbstractArray && ndims(first(argtypes)) != 2
                print(io, "\nThe `imview` function only supports 2D images. If you have a cube, try viewing one slice at a time: imview(cube[:,:,1])\n")
            end
        end
    end

    # TODO: This should be registered correctly with FileIO
    del_format(format"FITS")
    add_format(format"FITS",
        # See https://www.loc.gov/preservation/digital/formats/fdd/fdd000317.shtml#sign
        [0x53,0x49,0x4d,0x50,0x4c,0x45,0x20,0x20,0x3d,0x20,0x20,0x20,0x20,0x20,0x20,0x20,0x20,0x20,0x20,0x20,0x20,0x20,0x20,0x20,0x20,0x20,0x20,0x20,0x20,0x54],
        [".fit", ".fits", ".fts", ".FIT", ".FITS", ".FTS", ".fit",],
        [:FITSIO => UUID("525bcba6-941b-5504-bd06-fd0dc1a4d2eb")],
        [:AstroImages => UUID("fe3fc30c-9b16-11e9-1c73-17dabf39f4ad")]
    )

end

end # module