text:
using Gnuplot, Random
Random.seed!(123)
test = rand(50,50)
@gp "set auto fix" "set size square" :-
@gp :- test "w image notit"  """set cblabel "CBTitle \\n (my unit)" """ :-
@gp :- "set cblabel  offset -6.5, 10 font ',8' textcolor lt 3 rotate by 0" :-
@gp :- palette(:plasma) :-
@gp :- "set tmargin at screen 0.90"
save(term = "pngcairo size 400,400", output = "colorbar_labels_single.png")

max_stars_repo_name: lazarusA/gnuplot-examples

max_issues_repo_name: lazarusA/gnuplot-examples

max_forks_repo_name: lazarusA/gnuplot-examples


text:
using PowerLaws
using Base.Test

# tests are based moby dick word distribution
# and other data sets from http://tuvalu.santafe.edu/~aaronc/powerlaws/data.htm
# results are compared to result from R poweRlaw package
# test are done with tolerance to 1e-8
tolerance = 1e-8

moby_data = vec(readdlm(string(Pkg.dir("PowerLaws"),"/data/moby_dick.txt"),' ',Int))
cities = vec(readdlm(string(Pkg.dir("PowerLaws"),"/data/cities.txt"),' ',Int))
electrical_blackout = vec(readdlm(string(Pkg.dir("PowerLaws"),"/data/electrical_blackouts_US.txt"),' ',Int))
population = vec(readdlm(string(Pkg.dir("PowerLaws"),"/data/population.txt"),' ','\n'))

println("test: discrete powerlaw")
est = estimate_xmin(moby_data, dis_powerlaw)
@test_approx_eq_eps est[1].α 1.9527275186009831 tolerance
@test est[1].θ == 7.0
@test_approx_eq_eps est[2] 0.008252950457938835 tolerance

est = estimate_xmin(moby_data, dis_powerlaw, xmins = [2,3,4,10,20])
@test_approx_eq_eps est[1].α 1.9550379794745847 tolerance
@test est[1].θ == 10.0
@test_approx_eq_eps est[2] 0.011867106465087929 tolerance

est1 = estimate_xmin(cities, dis_powerlaw)
@test_approx_eq_eps est1[1].α 1.614392642791718 tolerance
@test est1[1].θ == 1021.0
@test_approx_eq_eps est1[2] 0.06088582981503132 tolerance

println("test: discrete powerlaw - bootstrap")
est1 = estimate_xmin(electrical_blackout, dis_powerlaw)
@test_approx_eq_eps est1[1].α 1.2201523540878356 tolerance
@test est1[1].θ == 1000.0
@test_approx_eq_eps est1[2] 0.36278306128485405 tolerance

bootstr = bootstrap(moby_data,dis_powerlaw,no_of_sims = 15)
@test length(bootstr) == 15

bootstr = bootstrap(moby_data,est[1],no_of_sims = 15)
@test length(bootstr) == 15

bootstr = bootstrap(moby_data,dis_powerlaw,3,no_of_sims = 15)
@test length(bootstr) == 15

bootstr = bootstrap(moby_data,est[1],2,no_of_sims = 15)
@test length(bootstr) == 15

println("test: continuous powerlaw")
est = estimate_xmin(moby_data, con_powerlaw)
@test_approx_eq_eps est[1].α 1.9335240463440206 tolerance
@test est[1].θ == 26.0
@test_approx_eq_eps est[2] 0.03204776085711486 tolerance

est = estimate_xmin(moby_data, con_powerlaw, xmins = [2,3,4,10,20])
@test_approx_eq_eps est[1].α 1.9514282451732778 tolerance
@test est[1].θ == 20.0
@test_approx_eq_eps est[2] 0.04416094210009813 tolerance

est = estimate_xmin(collect(1:100),con_powerlaw)
@test_approx_eq_eps est[1].α 9.178824791215408 tolerance
@test est[1].θ == 79.0
@test_approx_eq_eps est[2] 0.1824330509645657 tolerance

est1 = estimate_xmin(population,con_powerlaw)
@test_approx_eq_eps est1[1].α 2.0337071519535046 tolerance
@test est1[1].θ == 96479.0
@test_approx_eq_eps est1[2] 0.0 tolerance

println("test: continuous powerlaw - bootstrap")
bootstr = bootstrap(moby_data,con_powerlaw,no_of_sims = 15)
@test length(bootstr) == 15

bootstr = bootstrap(moby_data,est[1],no_of_sims = 15)
@test length(bootstr) == 15

bootstr = bootstrap(moby_data,con_powerlaw,3,no_of_sims = 15)
@test length(bootstr) == 15

bootstr = bootstrap(moby_data,est[1],2,no_of_sims = 15)
@test length(bootstr) == 15

println("test: compare distribution")
data = collect(1:100)
d1 = fit(dis_powerlaw,data)
d2 = fit(Poisson,data)
ll_hoods_r = logpdf(d1,data) - logpdf(d2,data)
cmpd = compare_distributions(d1,d2,data)
@test typeof(cmpd) == compare_distributions
@test cmpd.data == data
@test cmpd.log_likehoods_ratio == ll_hoods_r
@test cmpd.xmin == 1
@test cmpd.sig_level == 0.05
@test_approx_eq_eps cmpd.V_test_stat 5.7463540130003254 tolerance
@test_approx_eq_eps cmpd.V_p_val 0.9999999954405856 tolerance
@test cmpd.V_preff_distr == 1
@test cmpd.C_b == 62
@test_approx_eq_eps cmpd.C_p_val 0.020978735677851468 tolerance
@test cmpd.C_preff_distr == 1

cmpd = compare_distributions(d1,Poisson,data)
@test typeof(cmpd) == compare_distributions
@test cmpd.data == data
@test cmpd.log_likehoods_ratio == ll_hoods_r
@test cmpd.xmin == 1
@test cmpd.sig_level == 0.05
@test_approx_eq_eps cmpd.V_test_stat 5.7463540130003254 tolerance
@test_approx_eq_eps cmpd.V_p_val 0.9999999954405856 tolerance
@test cmpd.V_preff_distr == 1
@test cmpd.C_b == 62
@test_approx_eq_eps cmpd.C_p_val 0.020978735677851468 tolerance
@test cmpd.C_preff_distr == 1

cmpd = compare_distributions(dis_powerlaw,Poisson,data)
@test typeof(cmpd) == compare_distributions
@test cmpd.data == data
@test cmpd.log_likehoods_ratio == ll_hoods_r
@test cmpd.xmin == 1
@test cmpd.sig_level == 0.05
@test_approx_eq_eps cmpd.V_test_stat 5.7463540130003254 tolerance
@test_approx_eq_eps cmpd.V_p_val 0.9999999954405856 tolerance
@test cmpd.V_preff_distr == 1
@test cmpd.C_b == 62
@test_approx_eq_eps cmpd.C_p_val 0.020978735677851468 tolerance
@test cmpd.C_preff_distr == 1

cmpd = compare_distributions(dis_powerlaw,con_powerlaw,data)
@test cmpd == Union{}

cmpd = compare_distributions(d1,con_powerlaw,data)
@test cmpd == Union{}

try
    cmpd = compare_distributions(d1,con_powerlaw(),data)
catch e
    @test true == isa(e,MethodError)
end

moby_data = sort(moby_data)
d1 = estimate_xmin(moby_data,con_powerlaw)[1]
d2 = fit(Exponential,moby_data[18072:end])
cmpd = compare_distributions(d1,d2,moby_data,26)
@test cmpd.xmin == 26
@test cmpd.sig_level == 0.05
@test_approx_eq_eps cmpd.V_test_stat 8.710089666437788 tolerance
@test_approx_eq_eps cmpd.V_p_val 1.0 tolerance
@test cmpd.V_preff_distr == 1
@test cmpd.C_b == 566
@test_approx_eq_eps cmpd.C_p_val 0.0 tolerance
@test cmpd.C_preff_distr == 1

d1 = estimate_xmin(moby_data,dis_powerlaw)[1]
d2 = fit(Poisson,moby_data[15898:end])
cmpd = compare_distributions(d1,d2,moby_data,7.0)
@test cmpd.xmin == 7.0
@test cmpd.sig_level == 0.05
@test_approx_eq_eps cmpd.V_test_stat 4.448486373104684 tolerance
@test_approx_eq_eps cmpd.V_p_val 0.9999956761232807 tolerance
@test cmpd.V_preff_distr == 1
@test cmpd.C_b == 2757
@test_approx_eq_eps cmpd.C_p_val 0.0 tolerance
@test cmpd.C_preff_distr == 1

ext: jl


text:
using Simplices.Delaunay: delaunayn
import ..Embeddings: AbstractEmbedding, Embedding, LinearlyInvariantEmbedding
"""
A triangulation of a cloud of embedded points into disjoint simplices.

```julia
# Triangulate a set of random points in 3D space.
t = triangulate(rand(20, 3))

# Refine triangulation until all simplices are below the mean radius of the original
# triangulation.
target_radius = mean(t.radii)
refine_variable_k!(t, target_radius)
```

"""
abstract type AbstractTriangulation <: Partition end

#embedding(t::AbstractTriangulation) = t.embedding
#dimension(t::AbstractTriangulation) = dimension(t.embedding)
#npoints(t::AbstractTriangulation) = npoints(t.embedding)
#radii(t::AbstractTriangulation) = t.radii
#volumes(t::AbstractTriangulation) = t.volumes
#orientations(t::AbstractTriangulation) = t.orientations

"""
    `Triangulation <: AbstractTriangulation`

A triangulation type. Has the following fields:
1. `embedding::Embedding`.
2. `points::Array{Float64, 2}`. The vertices of the triangulation (a subset of the points
    of the embedding).
3. `impoints::Array{Float64, 2}`. The image vertices of the triangulation.
4. `simplex_inds::Array{Int, 2}`. # Array of indices referencing the vertices furnishing
    each simplex, expressed both in terms of the original points and their images under
    the linear forward map.
5. `centroids::Array{Float64, 2}`. Centroids of the simplices.
6. `radii::Vector{Float64}`. Radii of the simplices.
7. `centroids_im::Array{Float64, 2}`. Centroids of the image simplices.
8. `radii_im::Vector{Float64}`. Radii of the image simplices.
9. `orientations::Vector{Float64}`. Orientations of the simplices.
10. `orientations_im::Vector{Float64}`. Orientations of the image simplices.
11. `volumes::Vector{Float64}`. Volumes of the simplices.
12. `volumes_im::Vector{Float64}`. Volumes of the image simplices.
"""
mutable struct Triangulation <: AbstractTriangulation
    embedding::AbstractEmbedding
    points::Array{Float64, 2}
    impoints::Array{Float64, 2}
    simplex_inds::Array{Int, 2}
    centroids::Array{Float64, 2}
    centroids_im::Array{Float64, 2}
    radii::Vector{Float64}
    radii_im::Vector{Float64}
    orientations::Vector{Float64}
    orientations_im::Vector{Float64}
    volumes::Vector{Float64}
    volumes_im::Vector{Float64}
end

export dimension, npoints, radii, volumes, orientations

"""
    LinearlyInvariantTriangulation <: Triangulation`
A triangulation for which we have made sure the point corresponding to the last time
    index falls within the convex hull of the other points.
"""
mutable struct LinearlyInvariantTriangulation <: AbstractTriangulation
    embedding::AbstractEmbedding
    points::Array{Float64, 2}
    impoints::Array{Float64, 2}
    simplex_inds::Array{Int, 2}
    centroids::Array{Float64, 2}
    centroids_im::Array{Float64, 2}
    radii::Vector{Float64}
    radii_im::Vector{Float64}
    orientations::Vector{Float64}
    orientations_im::Vector{Float64}
    volumes::Vector{Float64}
    volumes_im::Vector{Float64}
end

"""
    triangulate(points::Array{Float64, 2})

Triangulate a set of vertices in N dimensions. `points` is an array of vertices,
where each row of the array is a point.
"""
function delaunay_triang(points::AbstractArray{Float64, 2})
    indices = delaunayn(points)
    return indices
end

function triangulate(E::LinearlyInvariantEmbedding)
    points = transpose(E.points[:, 1:end-1])
    simplex_inds = delaunay_triang(points)
    impoints = transpose(E.points[:, 2:end])
    c, r = centroids_radii2(points, simplex_inds)
    c_im, r_im = centroids_radii2(impoints, simplex_inds)
    vol = simplex_volumes(points, simplex_inds)
    vol_im = simplex_volumes(impoints, simplex_inds)
    o = orientations(points, simplex_inds)
    o_im = orientations(impoints, simplex_inds)

    LinearlyInvariantTriangulation(E,
        points,
        impoints,
        simplex_inds,
        c,
        c_im,
        r,
        r_im,
        o,
        o_im,
        vol,
        vol_im)
end

function triangulate(E::AbstractEmbedding)
    points = transpose(E.points[:, 1:end-1])
    simplex_inds = delaunay_triang(points)
    impoints = transpose(E.points[:, 2:end])
    c, r = centroids_radii2(points, simplex_inds)
    c_im, r_im = centroids_radii2(impoints, simplex_inds)
    vol = simplex_volumes(points, simplex_inds)
    vol_im = simplex_volumes(impoints, simplex_inds)
    o = orientations(points, simplex_inds)
    o_im = orientations(impoints, simplex_inds)

    Triangulation(E,
        points,
        impoints,
        simplex_inds,
        c,
        c_im,
        r,
        r_im,
        o,
        o_im,
        vol,
        vol_im)
end

triangulate(pts::AbstractArray{Float64, 2}) = triangulate(embed(pts))

function Base.summary(t::AbstractTriangulation)
    npts = size(t.embedding.points, 2)
    nsimplices = size(t.simplex_inds, 1)
    dim = size(t.embedding.points, 1)
    embeddingtype_tri = typeof(t)
    embeddingtype_emb = typeof(t.embedding)
    return """$dim-dimensional $(embeddingtype_tri) with $nsimplices simplices constructed
            from a $npts-pt $embeddingtype_emb"""
end

function matstring(t::AbstractTriangulation)
    fields = fieldnames(typeof(t))
    fields_str = String.(fields)
    maxlength = maximum([length(str) for str in fields_str]) + 2
    fields_str = [fields_str[i] *
                repeat(" ", maxlength - length(fields_str[i])) for i = 1:length(fields_str)]

    summaries = [join(":"*String(fields_str[i])*summary(getfield(t, fields[i]))*"\n") for i = 1:length(fields_str)] |> join
    infoline = "The following fields are available:\n"

    return summary(t)#*"\n\n"*infoline*summaries
end

Base.show(io::IO, t::AbstractTriangulation) = println(io, matstring(t))

"""
Find the indices of the simplices in the original triangulation that potentially
intersect with the image simplex with index `image_i`.
"""
function maybeintersecting_simplices(t::AbstractTriangulation, image_i::Int)
    inds_potential_simplices = Int[]

    n_simplices = length(t.radii)

    @inbounds for i = 1:n_simplices
        dist_difference = (transpose(t.centroids_im[image_i] - t.centroids[i]) *
                            (t.centroids_im[image_i] - t.centroids[i]) -
                                (t.radii_im[image_i] + t.radii[i])^2)[1]
        if dist_difference < 0
            push!(inds_potential_simplices, i)
        end
    end
    return inds_potential_simplices
end

"""
Find the indices of the image simplices in `t` that potentially intersect with
the original simplex with index `orig_i`.
"""
function maybeintersecting_imsimplices(t::AbstractTriangulation, orig_i::Int)
    inds_potential_simplices = Int[]

    n_simplices = length(t.radii)

    @inbounds for i = 1:n_simplices
        dist_difference = (transpose(t.centroids[orig_i] - t.centroids_im[i]) *
                            (t.centroids[orig_i] - t.centroids_im[i]) -
                                (t.radii[orig_i] + t.radii_im[i])^2)[1]
        if dist_difference < 0
            push!(inds_potential_simplices, i)
        end
    end
    return inds_potential_simplices
end

"""
Draw point representatives from the simplices of a triangulation `t`.

Precedure:
1) Generate one point per simplex.
2) Points are generated from the interior or on the boundary of each simplex.
3) Points are drawn according to a uniform distribution.
"""
function point_representatives(t::AbstractTriangulation)
    dim = size(t.points, 2)
    n_simplices = size(t.simplex_inds, 1)

    # Pre-allocate array to hold the points
    point_representatives = zeros(Float64, n_simplices, dim)

    # Loop over the rows of the simplex_inds array to access all the simplices.
    for i = 1:n_simplices
        simplex = t.points[t.simplex_inds[i, :], :]
        point_representatives[i, :] = Delaunay.childpoint(simplex)
    end

    return point_representatives
end


text:
"""
    newton(f::Function, x0::Real; niters=10^2, tol=1e-8)

Compute roots of scalar `f` by Newton's method

"""

function newton(f::Function, x0::Real; niters=10^2, tol=1e-8)
    x = xnew = x0
    for it in 1:niters
        df = ForwardDiff.derivative(f, x)
        xnew = x - f(x)/df
        abs(x - xnew) < tol && return xnew
        x, xnew = xnew, x
    end
    xnew
end


text:
using Viznet
using Compose
using Random
Random.seed!(4)

ud = rand_unitdisk(20, 2.0)

function fancy_unitdisk(ud)
    unitc = compose(context(),
        circle(0.0, 0.0, ud.unit/2),
        stroke("gray"),
        fill("transparent"),
        strokedash([1mm,1mm]),
        )
    g1 = viz(ud)
    g2 = canvas() do
        unitc >> ud[1]
        unitc >> ud[7]
    end
    compose(context(), g2, g1)
end

fancy_unitdisk(ud) |> SVG("_unitdisk.svg")


text:
# Motion parameter utilities

export derivative1d_linop, derivative1d_motionpars_linop, interpolation1d_linop, interpolation1d_motionpars_linop

function derivative1d_linop(t::AbstractVector{T}, order::Integer) where {T<:Real}
    (order != 1) && (order != 2) && error("Order not supported (only 1 or 2)")
    nt = length(t)
    Δt = t[2:end]-t[1:end-1]
    if order == 1
        d_0 = [-T(1)./Δt; T(0)]
        d_p1 = T(1)./Δt;
        return spdiagm(nt, nt, 0 => d_0, 1 => d_p1)
    elseif order == 2
        Δt_m1 = Δt[1:end-1]
        Δt_p1 = Δt[2:end]
        Δt_mean = (Δt_m1+Δt_p1)/T(2)
        d_m1 = T(1)./(Δt_m1.*Δt_mean)
        d_p1 = [T(0); T(1)./(Δt_p1.*Δt_mean)]
        d_0  = [T(0); -d_m1-d_p1[2:end]]
        return spdiagm(nt, nt, -1 => d_m1, 0 => d_0, 1 => d_p1)
    end
end

function derivative1d_motionpars_linop(t::AbstractVector{T}, order::Integer; pars::NTuple{6,Bool}=(true, true, true, true, true, true)) where {T<:Real}
    D = derivative1d_linop(t, order)
    Id = sparse(I, length(t), length(t))
    A = []
    for i = 1:6
        pars[i] ? push!(A, D) : push!(A, Id)
    end
    return cat(A...; dims=(1,2))
end

function interpolation1d_linop(t::AbstractVector{T}, ti::AbstractVector{T}; interp::Symbol=:linear) where {T<:Real}
    (interp != :linear) && error("Only linear interpolation supported")
    nt = length(t)
    nti = length(ti)
    I = repeat(reshape(1:nti, :, 1); outer=(1,2))
    J = Array{Int64,2}(undef, nti, 2)
    V = Array{T,2}(undef, nti, 2)
    for i = 1:nti
        if     ti[i] < t[1]
            J[i, :] .= [1; 1]
            V[i, :] .= [T(0); T(0)]
        elseif ti[i] > t[end]
            J[i, :] .= [1; 1]
            V[i, :] .= [T(0); T(0)]
        elseif ti[i] == t[1]
            J[i, :] .= [1; 2]
            V[i, :] .= [T(1); T(0)]
        elseif ti[i] == t[end]
            J[i, :] .= [nt-1; nt]
            V[i, :] .= [T(0); T(1)]
        else
            idx_ = findall(t .< ti[i])
            if length(idx_) != 0
                idx = maximum(idx_)
                J[i,:] .= [idx; idx+1]
                Δt = t[idx+1]-t[idx]
                V[i,:] .= [(t[idx+1]-ti[i])/Δt; (ti[i]-t[idx])/Δt]
            else
                J[i, :] .= [1; 1]
                V[i, :] .= [T(0); T(0)]
            end
        end
    end
    return sparse(vec(I), vec(J), vec(V), nti, nt)
end

function interpolation1d_motionpars_linop(t::NTuple{6,AbstractVector{T}}, ti::NTuple{6,AbstractVector{T}}; interp::Symbol=:linear) where {T<:Real}
    Ip = []
    for i=1:6
        push!(Ip, interpolation1d_linop(t[i], ti[i]; interp=interp))
    end
    return cat(Ip...; dims=(1,2))
end

interpolation1d_motionpars_linop(t::AbstractVector{T}, ti::NTuple{6,AbstractVector{T}}; interp::Symbol=:linear) where {T<:Real} = interpolation1d_motionpars_linop((t,t,t,t,t,t), ti; interp=interp)

interpolation1d_motionpars_linop(t::NTuple{6,AbstractVector{T}}, ti::AbstractVector{T}; interp::Symbol=:linear) where {T<:Real} = interpolation1d_motionpars_linop(t, (ti,ti,ti,ti,ti,ti); interp=interp)

interpolation1d_motionpars_linop(t::AbstractVector{T}, ti::AbstractVector{T}; interp::Symbol=:linear) where {T<:Real} = interpolation1d_motionpars_linop((t,t,t,t,t,t), (ti,ti,ti,ti,ti,ti); interp=interp)

max_stars_repo_name: grizzuti/UtilitiesForMRI

max_issues_repo_name: grizzuti/UtilitiesForMRI

max_forks_repo_name: grizzuti/UtilitiesForMRI


text:
# TODO: Change verbosity on LightGBM's side after a booster is created.
"""
    predict(estimator, X; [predict_type = 0, num_iterations = -1, verbosity = 1,
    is_row_major = false])

Return an array with the labels that the `estimator` predicts for features data `X`.

# Arguments
* `estimator::LGBMEstimator`: the estimator to use in the prediction.
* `X::Matrix{T<:Real}`: the features data.
* `predict_type::Integer`: keyword argument that controls the prediction type. `0` for normal
    scores with transform (if needed), `1` for raw scores, `2` for leaf indices.
* `num_iterations::Integer`: keyword argument that sets the number of iterations of the model to
    use in the prediction. `< 0` for all iterations.
* `verbosity::Integer`: keyword argument that controls LightGBM's verbosity. `< 0` for fatal logs
    only, `0` includes warning logs, `1` includes info logs, and `> 1` includes debug logs.
* `is_row_major::Bool`: keyword argument that indicates whether or not `X` is row-major. `true`
    indicates that it is row-major, `false` indicates that it is column-major (Julia's default).
"""
function predict(estimator::LGBMEstimator, X::Matrix{TX}; predict_type::Integer = 0,
                           num_iterations::Integer = -1, verbosity::Integer = 1,
                           is_row_major::Bool = false) where TX<:Real
    @assert(estimator.booster.handle != C_NULL, "Estimator does not contain a fitted model.")
    log_debug(verbosity, "Started predicting\n")
    prediction = LGBM_BoosterPredictForMat(estimator.booster, X, predict_type, num_iterations,
                                           is_row_major)

    return prediction
end


text:
mutable struct EarthConstants
    mu::Float64
    JD_0::Float64
    igrf_Y_0::Float64
    igrf::Any
    EarthConstants(; mu = 3.986004418e14, JD_0 = 0., igrf_Y_0 = 0., igrf = nothing) = new(mu, JD_0, igrf_Y_0, igrf)
end

function LowEarthOrbit()

    # Tell the sim that Earths have no state (anarchy?) and that they produce
    # Gravitys.
    # function init(t, constants, state, draws)
    #     return state
    # end

    # Output all of the effects on the vehicle.
    function effects(t, constants, state, draws, effects, effects_bus)
        body, found = find_effect(effects, BodyStateEffect(zeros(3), zeros(3), zeros(3), zeros(3)))
        @assert found > 0 "No body effects were found."
        return (Gravity(point_mass_gravity(body.r_be_I, constants.mu)),)
    end

    # Create the IdealActuator as a DynamicalModel.
    DynamicalModel("low_earth_orbit",
                   #init = init,
                   effects = effects,
                   constants = EarthConstants())

end


text:
# RUN: julia -e "import Brutus; Brutus.lit(:emit_optimized)" --startup-file=no %s 2>&1 | FileCheck %s

struct A
    x
end
(a::A)(y) = a.x + y
a = A(10)
emit(a, Int64)

# CHECK: Core.MethodMatch(Tuple{Main.Main.A, Int64}, svec(), (a::Main.Main.A)(y) in Main.Main at /home/mccoy/Dev/brutus/test/Codegen/canonicalize/call_overload.jl:6, true)after translating to MLIR in JLIR dialect:module  {
# CHECK-NEXT:   func nested @"Tuple{Main.A, Int64}"(%arg0: !jlir.Main.A, %arg1: !jlir.Int64) -> !jlir.Any attributes {llvm.emit_c_interface} {
# CHECK-NEXT:     "jlir.goto"()[^bb1] : () -> ()
# CHECK-NEXT:   ^bb1:  // pred: ^bb0
# CHECK-NEXT:     %0 = "jlir.constant"() {value = #jlir.getfield} : () -> !jlir<"typeof(getfield)">
# CHECK-NEXT:     %1 = "jlir.constant"() {value = #jlir<":x">} : () -> !jlir.Symbol
# CHECK-NEXT:     %2 = "jlir.call"(%0, %arg0, %1) : (!jlir<"typeof(getfield)">, !jlir.Main.A, !jlir.Symbol) -> !jlir.Any
# CHECK-NEXT:     %3 = "jlir.constant"() {value = #jlir<"Base.:(+)">} : () -> !jlir<"typeof(Base.:(+))">
# CHECK-NEXT:     %4 = "jlir.call"(%3, %2, %arg1) : (!jlir<"typeof(Base.:(+))">, !jlir.Any, !jlir.Int64) -> !jlir.Any
# CHECK-NEXT:     "jlir.return"(%4) : (!jlir.Any) -> ()
# CHECK-NEXT:   }
# CHECK-NEXT: }

# CHECK:   func nested @"Tuple{Main.A, Int64}"(%arg0: !jlir.Main.A, %arg1: !jlir.Int64) -> !jlir.Any attributes {llvm.emit_c_interface} {
# CHECK-NEXT:     "jlir.goto"()[^bb1] : () -> ()
# CHECK-NEXT:   ^bb1:  // pred: ^bb0
# CHECK-NEXT:     %0 = "jlir.constant"() {value = #jlir<":x">} : () -> !jlir.Symbol
# CHECK-NEXT:     %1 = "jlir.getfield"(%arg0, %0) : (!jlir.Main.A, !jlir.Symbol) -> !jlir.Any
# CHECK-NEXT:     %2 = "jlir.constant"() {value = #jlir<"Base.:(+)">} : () -> !jlir<"typeof(Base.:(+))">
# CHECK-NEXT:     %3 = "jlir.call"(%2, %1, %arg1) : (!jlir<"typeof(Base.:(+))">, !jlir.Any, !jlir.Int64) -> !jlir.Any
# CHECK-NEXT:     "jlir.return"(%3) : (!jlir.Any) -> ()
# CHECK-NEXT:   }
# CHECK-NEXT: }

max_stars_repo_name: JuliaLabs/brutus

max_issues_repo_name: JuliaLabs/brutus

max_forks_repo_name: JuliaLabs/brutus


text:
module DDPG

using OpenAIGym
using DataStructures

# NOTE: to find other valid environment names, look at the universe code that registers them:
#   https://github.com/openai/universe/blob/master/universe/__init__.py
env = gym("flashgames.DuskDrive-v0")
# env = gym("Pong-v3")
# env = gym("wob.mini.CircleCenter-v0")
# env = gym("Breakout-v0")
@show actions(env, state(env))
# error()

# agent/policy
policy = RandomPolicy()

# -----------------------------------
# quick AR process for arbitrary vectors
# used for exploration policy of DDPG
using Distributions
mutable struct ARProcess{T}
    prev::T
    reversion
    noise
end
function Base.get(ar::ARProcess)
    ar.prev .= ar.reversion .* ar.prev .+ rand(noise)
end
# -----------------------------------

mutable struct DdpgPolicy
    ns; na; nϕ
    features
    actor
    actor_target
    critic
    critic_target
    experience
end

function build_actor_critic(env)
    ns = length(state(env))
    na = length(actions(env))

    # shared feature map: ϕ(s,a)
    nϕ = 10
    features = nnet(ns, nϕ, [10], :relu)

    # actor: μ(s | Θμ)
    actor = Chain(features, Affine(nϕ, na))
    actor_target = copy(actor)

    # critic: Q(s,a | ΘQ)
    critic = Chain(features, Concat(nϕ+na), Affine(nϕ+na, 1))
    critic_target = copy(critic)

    # experience replay buffer
    experience = CircularBuffer{Tuple}(100)

    DdpgPolicy(ns, na, nϕ, features, actor, actor_target, critic, critic_target, experience)
end

# main loop... run one episode, getting a tuple: (s, a, r, s′)
for sars in Episode(env, policy)
    # @show sars
    s,a,r,s′ = sars
    # @show a,r

    # push!(experience, sars)

    OpenAIGym.render(env)

    #= TODO
    - actor/critic share "feature transformation" ϕ(s,a)
    - update policy (agent) by sampling from the experience replay
    =#
end
# @show length(experience)

end # module


text:
# script to create all plots

for file in ("plot_intro_exact_figure_1a.jl",
             "plot_intro_approximate_figure_1b.jl",
             "plot_intro_recurrence_figure_1c.jl",
             "plot_homogenize_figure_2a_and_figure_8.jl",
             "plot_shrink_delta_figure_2b.jl",
             "plot_delta_figure_3.jl",
             "plot_X0_figure_4.jl",
             "plot_U_figure_5.jl",
             # figure 6 is generated from the experiment runtime data
             "plot_underapproximation_figure_7.jl")
    println("Running $file ...")
    include(file)
end

nothing


text:
@testset "EGM96 End to End Tests" begin
    @testset "EGM96 F477 Reference Implementation Tests" begin
        @test_approx_eq_eps undulation(38.6281550, 269.7791550) -31.629 1e-3
        @test_approx_eq_eps undulation(-14.6212170, 305.0211140) -2.966 1e-3
        @test_approx_eq_eps undulation(46.8743190, 102.4487290) -43.572 1e-3
        @test_approx_eq_eps undulation(-23.6174460, 133.8747120) 15.868 1e-3
        @test_approx_eq_eps undulation(38.6254730, 359.9995000) 50.065 1e-3
        @test_approx_eq_eps undulation(-.4667440, .0023000) 17.330 1e-3
    end
end


text:
#=
test.jl
=======
Uses .NonsmoothFwdAD to evaluate generalized derivative elements, replicating
the results of calculations from the following article:

[1]: KA Khan and PI Barton (2013), https://doi.org/10.1145/2491491.2491493
=#

include("../src/NonsmoothFwdAD.jl")

using .NonsmoothFwdAD

println("Replicating Example 6.1 from [1]:")

_, f1Grad = eval_gen_gradient(x -> abs(x[1]) - abs(x[1]), [0.0])
println("  Gen. gradient element ",
        "of (x -> abs(x) - abs(x)) at 0.0:")
println("    ", f1Grad)

_, g1Grad = eval_gen_gradient(x -> max(x[1], 0.0) - min(x[1], 0.0), [0.0])
println("  Gen. gradient element ",
        "of (x -> max(x, 0.0) - min(x, 0.0)) at 0.0:")
println("    ", g1Grad)

println("\nReplicating Example 6.2 from [1]:")

f2(x) = max(min(x[1], -x[2]), x[2] - x[1])

_, f2Grad = eval_gen_gradient(f2, [0.0, 0.0])
println("  Gen. gradient element ",
        "of (x -> max(min(x[1], -x[2]), x[2] - x[1]) ",
        "at [0.0, 0.0]:")
println("    ", f2Grad)

_, f2Compass = eval_compass_difference(f2, [0.0, 0.0])
println("  Compass difference ",
        "of this function ",
        "at [0.0, 0.0]:")
println("    ", f2Compass)

println("\nReplicating Example 6.3 from [1]:")

_, f3Grad = eval_gen_gradient([0.0, 0.0]) do x
    return (1.0 + abs(x[1] - x[2]))*(x[1] - x[2])
end
println("  Gen. gradient element ",
        "of (x -> (1.0 + abs(x[1] - x[2]))*(x[1] - x[2])) at [0.0, 0.0]:")
println("    ", f3Grad)

println("\nReplicating Example 6.5 from [1]:")

tInVal = [160.0, 170.0, 60.0, 116.0]
tOutVal = [93.0, 126.0, 160.0, 260.0]
t0Val = Float64[]
for (tI, tO) in zip(tInVal, tOutVal)
    append!(t0Val, [tI, tO])
end

qVal, qGrad = eval_gen_derivative(t0Val) do t
    # stream parameters
    fcP = [8.79, 10.55, 7.62, 6.08]
    deltaT = 10.0
    nStreams = length(fcP)

    # recover inlet/outlet temperatures from input
    tIn = [t[2*i - 1] for i in 1:nStreams]
    tOut = [t[2*i] for i in 1:nStreams]

    # subtract min. temperature approach from inlet/outlet temperatures of hot streams
    for j in 1:nStreams
        if tOut[j] < tIn[j]
            tIn[j] -= deltaT
            tOut[j] -= deltaT
        end
    end

    # compute minimum heating
    qH = 0.0
    for i in 1:nStreams
        zP = 0.0
        for j in 1:nStreams
            zP += fcP[j] * (max(0.0, tOut[j] - tIn[i]) - max(0.0, tIn[j] - tIn[i]))
        end
        qH = max(qH, zP)
    end

    # compute minimum cooling
    omega = 0.0
    for j in 1:nStreams
        omega += fcP[j] * (tIn[j] - tOut[j])
    end
    qC = qH + omega

    return [qH, qC]
end

println("  With Q^H and Q^C defined as in this example,")
println("  and with stream parameters and T defined in Table V in [1],")
println("  gen. derivative element of [Q^H, Q^C] at T:")
for row in eachrow(qGrad)
    println("    ", row)
end

max_stars_repo_name: kamilkhanlab/nonsmooth-forward-ad

max_issues_repo_name: kamilkhanlab/nonsmooth-forward-ad

max_forks_repo_name: kamilkhanlab/nonsmooth-forward-ad


text:
function animate_NN(xs, y, t, x_str, x_label=["" for i in length(xs)], filename=x_str, PATH=joinpath(pwd(), "Output"))
    anim = @animate for n in 1:size(xs[1], 2)
        x_max = maximum(maximum(x) for x in xs)
        x_min = minimum(minimum(x) for x in xs)
        @info "$x_str frame of $n/$(size(xs[1], 2))"
        fig = Plots.plot(xlim=(x_min, x_max), ylim=(minimum(y), maximum(y)), legend=:bottom)
        for i in 1:length(xs)
            Plots.plot!(fig, xs[i][:,n], y, label=x_label[i], title="t = $(round(t[n] / 86400, digits=2)) days")
        end
        Plots.xlabel!(fig, "$x_str")
        Plots.ylabel!(fig, "z")
    end
    # gif(anim, joinpath(PATH, "$(filename).gif"), fps=30)
    mp4(anim, joinpath(PATH, "$(filename).mp4"), fps=30)
end

function animate_profile(data, profile_type, FILE_PATH; dimensionless=true, fps=30, gif=false, mp4=true)
    truth_profile = data["truth_$profile_type"]
    test_profile = data["test_$profile_type"]

    profile_max = maximum([maximum(truth_profile), maximum(test_profile)])
    profile_min = minimum([minimum(truth_profile), minimum(test_profile)])

    t = data["t"]

    z_profile = data["depth_profile"]

    z_max = maximum(z_profile)
    z_min = minimum(z_profile)

    anim = @animate for i in 1:length(t)
        @info "Animating $profile_type frame $i/$(length(t))"
        fig = plot(truth_profile[:,i], z_profile, xlim=(profile_min, profile_max), ylim=(z_min, z_max), label="Truth", title="$(round(t[i]/86400, digits=2)) days")
        plot!(fig, test_profile[:,i], z_profile, label="NN")
        ylabel!(fig, "z /m")

        if dimensionless
            xlabel!(fig, profile_type)
        elseif profile_type == "T"
            xlabel!(fig, "T /K")
        else
            xlabel!(fig, "$profile_type /m s⁻¹")
        end

        if i == 1
            savefig(fig, "$FILE_PATH.pdf")
            savefig(fig, "$FILE_PATH.png")
        end
    end

    if gif
        Plots.gif(anim, "$FILE_PATH.gif", fps=fps)
    end

    if mp4
        Plots.mp4(anim, "$FILE_PATH.mp4", fps=fps)
    end
end

function animate_flux(data, flux_type, FILE_PATH; dimensionless=true, fps=30, gif=false, mp4=true)
    truth_flux = data["truth_$flux_type"]
    test_flux = data["test_$flux_type"]

    flux_max = maximum([maximum(truth_flux), maximum(test_flux)])
    flux_min = minimum([minimum(truth_flux), minimum(test_flux)])

    t = data["t"]

    z_flux = data["depth_flux"]

    z_max = maximum(z_flux)
    z_min = minimum(z_flux)

    anim = @animate for i in 1:length(t)
        @info "Animating $flux_type frame $i/$(length(t))"
        fig = plot(truth_flux[:,i], z_flux, xlim=(flux_min, flux_max), ylim=(z_min, z_max), label="Truth", title="$(round(t[i]/86400, digits=2)) days")
        plot!(fig, test_flux[:,i], z_flux, label = "NN")
        ylabel!(fig, "z /m")

        if dimensionless
            xlabel!(fig, flux_type)
        else
            xlabel!(fig, "$flux_type /m² s⁻²")
        end

    end

    if gif
        Plots.gif(anim, "$FILE_PATH.gif", fps=fps)
    end

    if mp4
        Plots.mp4(anim, "$FILE_PATH.mp4", fps=fps)
    end
end

function animate_profile_flux(data, profile_type, flux_type, FILE_PATH; dimensionless=true, fps=30, gif=false, mp4=true)
    truth_flux = data["truth_$flux_type"]
    test_flux = data["test_$flux_type"]

    truth_profile = data["truth_$profile_type"]
    test_profile = data["test_$profile_type"]

    flux_max = maximum([maximum(truth_flux), maximum(test_flux)])
    flux_min = minimum([minimum(truth_flux), minimum(test_flux)])

    profile_max = maximum([maximum(truth_profile), maximum(test_profile)])
    profile_min = minimum([minimum(truth_profile), minimum(test_profile)])

    t = data["t"]

    z_flux = data["depth_flux"]
    z_profile = data["depth_profile"]

    z_max = maximum([maximum(z_flux), maximum(z_profile)])
    z_min = minimum([minimum(z_flux), minimum(z_profile)])

    anim = @animate for i in 1:length(t)
        @info "Animating $flux_type/$profile_type frame $i/$(length(t))"
        l = @layout [a b]
        fig₁ = plot(truth_flux[:,i], z_flux, xlim=(flux_min, flux_max), ylim=(z_min, z_max), label="Truth")
        plot!(fig₁, test_flux[:,i], z_flux, label = "NN")
        ylabel!(fig₁, "z /m")
        if dimensionless
            xlabel!(fig₁, flux_type)
        elseif flux_type == "wT"
            xlabel!(fig₁, "$flux_type /m s⁻¹ °C")
        else
            xlabel!(fig₁, "$flux_type /m² s⁻²")
        end

        fig₂ = plot(truth_profile[:,i], z_profile, xlim=(profile_min, profile_max), ylim=(z_min, z_max), label="Truth", legend=:topleft)
        plot!(fig₂, test_profile[:,i], z_profile, label="NN")
        ylabel!(fig₂, "z /m")
        if dimensionless
            xlabel!(fig₂, profile_type)
        elseif profile_type == "T"
            xlabel!(fig₂, "T /°C")
        else
            xlabel!(fig₂, "$profile_type /m s⁻¹")
        end

        fig = plot(fig₁, fig₂, layout=l, title="$(round(t[i]/86400, digits=2)) days")
    end

    if gif
        Plots.gif(anim, "$FILE_PATH.gif", fps=fps)
    end

    if mp4
        Plots.mp4(anim, "$FILE_PATH.mp4", fps=fps)
    end
end

function animate_profiles(data, FILE_PATH; dimensionless=true, fps=30, gif=false, mp4=true)
    truth_u = data["truth_u"]
    truth_v = data["truth_v"]
    truth_T = data["truth_T"]

    test_u = data["test_u"]
    test_v = data["test_v"]
    test_T = data["test_T"]

    u_max = maximum([maximum(truth_u), maximum(test_u)])
    u_min = minimum([minimum(truth_u), minimum(test_u)])

    v_max = maximum([maximum(truth_v), maximum(test_v)])
    v_min = minimum([minimum(truth_v), minimum(test_v)])

    T_max = maximum([maximum(truth_T), maximum(test_T)])
    T_min = minimum([minimum(truth_T), minimum(test_T)])

    t = data["t"]

    z = data["depth_profile"]

    z_max = maximum(z)
    z_min = minimum(z)

    anim = @animate for i in 1:length(t)
        if i % 50 == 0
            @info "Animating frame $i/$(length(t))"
        end
        l = @layout [a b c]
        fig₁ = Plots.plot(truth_u[:,i], z, xlim=(u_min, u_max), ylim=(z_min, z_max), label="Truth", legend=:bottomright)
        Plots.plot!(fig₁, test_u[:,i], z, label = "NN")
        Plots.ylabel!(fig₁, "z /m")
        if dimensionless
            Plots.xlabel!(fig₁, "u")
        else
            Plots.xlabel!(fig₁, "u /m s⁻¹")
        end

        fig₂ = Plots.plot(truth_v[:,i], z, xlim=(v_min, v_max), ylim=(z_min, z_max), label="Truth", legend=:bottomleft)
        Plots.plot!(fig₂, test_v[:,i], z, label = "NN")
        Plots.ylabel!(fig₂, "z /m")
        if dimensionless
            Plots.xlabel!(fig₂, "v")
        else
            Plots.xlabel!(fig₂, "v /m s⁻¹")
        end

        fig₃ = Plots.plot(truth_T[:,i], z, xlim=(T_min, T_max), ylim=(z_min, z_max), label="Truth", legend=:bottomright)
        Plots.plot!(fig₃, test_T[:,i], z, label = "NN")
        Plots.ylabel!(fig₃, "z /m")
        if dimensionless
            Plots.xlabel!(fig₃, "T")
        else
            Plots.xlabel!(fig₃, "T /°C")
        end

        fig = Plots.plot(fig₁, fig₂, fig₃, layout=l, title="$(round(t[i]/86400, digits=2)) days")
    end

    if gif
        Plots.gif(anim, "$FILE_PATH.gif", fps=fps)
    end

    if mp4
        Plots.mp4(anim, "$FILE_PATH.mp4", fps=fps)
    end
end

function animate_profiles_fluxes(data, FILE_PATH; dimensionless=true, fps=30, gif=false, mp4=true, SIMULATION_NAME="")
    times = data["t"]

    frame = Node(1)

    truth_u = @lift data["truth_u"][:,$frame]
    truth_v = @lift data["truth_v"][:,$frame]
    truth_T = @lift data["truth_T"][:,$frame]

    test_u = @lift data["test_u"][:,$frame]
    test_v = @lift data["test_v"][:,$frame]
    test_T = @lift data["test_T"][:,$frame]

    truth_uw = @lift data["truth_uw"][:,$frame]
    truth_vw = @lift data["truth_vw"][:,$frame]
    truth_wT = @lift data["truth_wT"][:,$frame]

    test_uw = @lift data["test_uw"][:,$frame]
    test_vw = @lift data["test_vw"][:,$frame]
    test_wT = @lift data["test_wT"][:,$frame]

    truth_Ri = @lift clamp.(data["truth_Ri"][:,$frame], -1, 2)
    test_Ri = @lift clamp.(data["test_Ri"][:,$frame], -1, 2)

    u_max = maximum([maximum(data["truth_u"]), maximum(data["test_u"])])
    u_min = minimum([minimum(data["truth_u"]), minimum(data["test_u"])])

    v_max = maximum([maximum(data["truth_v"]), maximum(data["test_v"])])
    v_min = minimum([minimum(data["truth_v"]), minimum(data["test_v"])])

    T_max = maximum([maximum(data["truth_T"]), maximum(data["test_T"])])
    T_min = minimum([minimum(data["truth_T"]), minimum(data["test_T"])])

    uw_max = maximum([maximum(data["truth_uw"]), maximum(data["test_uw"])])
    uw_min = minimum([minimum(data["truth_uw"]), minimum(data["test_uw"])])

    vw_max = maximum([maximum(data["truth_vw"]), maximum(data["test_vw"])])
    vw_min = minimum([minimum(data["truth_vw"]), minimum(data["test_vw"])])

    wT_max = maximum([maximum(data["truth_wT"]), maximum(data["test_wT"])])
    wT_min = minimum([minimum(data["truth_wT"]), minimum(data["test_wT"])])

    plot_title = @lift "$SIMULATION_NAME: time = $(round(times[$frame]/86400, digits=2)) days, loss = $(round(data["loss"], sigdigits=3))"
    fig = Figure(resolution=(1920, 1080))
    colors=["navyblue", "hotpink2"]

    if dimensionless
        u_str = "u"
        v_str = "v"
        T_str = "T"
        uw_str = "uw"
        vw_str = "vw"
        wT_str = "wT"
    else
        u_str = "u / m s⁻¹"
        v_str = "v / m s⁻¹"
        T_str = "T / °C"
        uw_str = "uw / m² s⁻²"
        vw_str = "vw / m² s⁻²"
        wT_str = "wT / m s⁻¹ °C"
    end

    zc = data["depth_profile"]
    zf = data["depth_flux"]
    z_str = "z / m"

    ax_u = fig[1, 1] = Axis(fig, xlabel=u_str, ylabel=z_str)
    u_lines = [lines!(ax_u, truth_u, zc, linewidth=3, color=colors[1]), lines!(ax_u, test_u, zc, linewidth=3, color=colors[2])]
    CairoMakie.xlims!(ax_u, u_min, u_max)
    CairoMakie.ylims!(ax_u, minimum(zc), 0)

    ax_v = fig[1, 2] = Axis(fig, xlabel=v_str, ylabel=z_str)
    v_lines = [lines!(ax_v, truth_v, zc, linewidth=3, color=colors[1]), lines!(ax_v, test_v, zc, linewidth=3, color=colors[2])]
    CairoMakie.xlims!(ax_v, v_min, v_max)
    CairoMakie.ylims!(ax_v, minimum(zc), 0)

    ax_T = fig[1, 3] = Axis(fig, xlabel=T_str, ylabel=z_str)
    T_lines = [lines!(ax_T, truth_T, zc, linewidth=3, color=colors[1]), lines!(ax_T, test_T, zc, linewidth=3, color=colors[2])]
    CairoMakie.xlims!(ax_T, T_min, T_max)
    CairoMakie.ylims!(ax_T, minimum(zc), 0)

    ax_uw = fig[2, 1] = Axis(fig, xlabel=uw_str, ylabel=z_str)
    uw_lines = [lines!(ax_uw, truth_uw, zf, linewidth=3, color=colors[1]), lines!(ax_uw, test_uw, zf, linewidth=3, color=colors[2])]
    CairoMakie.xlims!(ax_uw, uw_min, uw_max)
    CairoMakie.ylims!(ax_uw, minimum(zf), 0)

    ax_vw = fig[2, 2] = Axis(fig, xlabel=vw_str, ylabel=z_str)
    vw_lines = [lines!(ax_vw, truth_vw, zf, linewidth=3, color=colors[1]), lines!(ax_vw, test_vw, zf, linewidth=3, color=colors[2])]
    CairoMakie.xlims!(ax_vw, vw_min, vw_max)
    CairoMakie.ylims!(ax_vw, minimum(zf), 0)

    ax_wT = fig[2, 3] = Axis(fig, xlabel=wT_str, ylabel=z_str)
    wT_lines = [lines!(ax_wT, truth_wT, zf, linewidth=3, color=colors[1]), lines!(ax_wT, test_wT, zf, linewidth=3, color=colors[2])]
    CairoMakie.xlims!(ax_wT, wT_min, wT_max)
    CairoMakie.ylims!(ax_wT, minimum(zf), 0)

    ax_Ri = fig[2, 4] = Axis(fig, xlabel="Ri", ylabel=z_str)
    Ri_lines = [lines!(ax_Ri, truth_Ri, zf, linewidth=3, color=colors[1]), lines!(ax_Ri, test_Ri, zf, linewidth=3, color=colors[2])]
    CairoMakie.xlims!(ax_Ri, -1, 2)
    CairoMakie.ylims!(ax_Ri, minimum(zf), 0)

    legend = fig[1, 4] = Legend(fig, u_lines, ["Oceananigans.jl LES", "NDE Prediction"])
    supertitle = fig[0, :] = Label(fig, plot_title, textsize=30)
    trim!(fig.layout)

    if gif
        CairoMakie.record(fig, "$FILE_PATH.gif", 1:length(times), framerate=fps) do n
            @info "Animating gif frame $n/$(length(times))..."
            frame[] = n
        end
    end

    if mp4
        CairoMakie.record(fig, "$FILE_PATH.mp4", 1:length(times), framerate=fps) do n
            @info "Animating mp4 frame $n/$(length(times))..."
            frame[] = n
        end
    end
end

function animate_profiles_fluxes_comparison(data_diffeq_explicit, data_diffeq_implicit, data_oceananigans, FILE_PATH; animation_type, n_trainings, training_types, fps=30, gif=false, mp4=true)
    times = data_diffeq_explicit["t"] ./ 86400

    frame = Node(1)

    time_point = @lift [times[$frame]]

    u_data = [
        data_diffeq_explicit["truth_u"],
        data_diffeq_explicit["test_u_modified_pacanowski_philander"],
        data_diffeq_explicit["test_u_kpp"],
        data_oceananigans["test_u"],
        data_diffeq_explicit["test_u"],
        data_diffeq_implicit["test_u"],
    ]

    v_data = [
        data_diffeq_explicit["truth_v"],
        data_diffeq_explicit["test_v_modified_pacanowski_philander"],
        data_diffeq_explicit["test_v_kpp"],
        data_oceananigans["test_v"],
        data_diffeq_explicit["test_v"],
        data_diffeq_implicit["test_v"],
    ]

    T_data = [
        data_diffeq_explicit["truth_T"],
        data_diffeq_explicit["test_T_modified_pacanowski_philander"],
        data_diffeq_explicit["test_T_kpp"],
        data_oceananigans["test_T"],
        data_diffeq_explicit["test_T"],
        data_diffeq_implicit["test_T"],
    ]

    uw_data = [
        data_diffeq_explicit["truth_uw"],
        data_diffeq_explicit["test_uw_modified_pacanowski_philander"],
        data_diffeq_explicit["test_uw_kpp"],
        data_oceananigans["test_uw"],
        data_diffeq_explicit["test_uw"],
        data_diffeq_implicit["test_uw"],
        data_diffeq_explicit["test_uw_NN_only"][2:end-1,:],
    ]

    vw_data = [
        data_diffeq_explicit["truth_vw"],
        data_diffeq_explicit["test_vw_modified_pacanowski_philander"],
        data_diffeq_explicit["test_vw_kpp"],
        data_oceananigans["test_vw"],
        data_diffeq_explicit["test_vw"],
        data_diffeq_implicit["test_vw"],
        data_diffeq_explicit["test_vw_NN_only"][2:end-1,:],
    ]

    wT_data = [
        data_diffeq_explicit["truth_wT"],
        data_diffeq_explicit["test_wT_modified_pacanowski_philander"],
        data_diffeq_explicit["test_wT_kpp"],
        data_oceananigans["test_wT"],
        data_diffeq_explicit["test_wT"],
        data_diffeq_implicit["test_wT"],
        data_diffeq_explicit["test_wT_NN_only"][2:end-1,:],
    ]

    Ri_data = [
        clamp.(data_diffeq_explicit["truth_Ri"], -1, 2),
        clamp.(data_diffeq_explicit["test_Ri_modified_pacanowski_philander"], -1, 2),
        clamp.(data_diffeq_explicit["test_Ri_kpp"], -1, 2),
        clamp.(data_oceananigans["test_Ri"], -1, 2),
        clamp.(data_diffeq_explicit["test_Ri"], -1, 2),
        clamp.(data_diffeq_implicit["test_Ri"], -1, 2),
    ]

    @inline function lowclamp(value, lo)
        if value >= lo
            return value
        else
            return lo
        end
    end

    losses_data = [
        lowclamp.(data_diffeq_explicit["losses_modified_pacanowski_philander"], 1f-5),
        lowclamp.(data_diffeq_explicit["losses_kpp"], 1f-5),
        lowclamp.(data_oceananigans["losses"], 1f-5),
        lowclamp.(data_diffeq_explicit["losses"], 1f-5),
        lowclamp.(data_diffeq_implicit["losses"], 1f-5),
        lowclamp.(data_diffeq_explicit["losses_modified_pacanowski_philander_gradient"], 1f-5),
        lowclamp.(data_diffeq_explicit["losses_kpp_gradient"], 1f-5),
        lowclamp.(data_oceananigans["losses_gradient"], 1f-5),
        lowclamp.(data_diffeq_explicit["losses_gradient"], 1f-5),
        lowclamp.(data_diffeq_implicit["losses_gradient"], 1f-5),
    ]

    @inline function add_ϵ!(losses)
        losses .= losses .+ (losses .== 0) .* eps(Float32)
    end

    add_ϵ!.(losses_data)

    u_frames = [@lift data[:,$frame] for data in u_data]
    v_frames = [@lift data[:,$frame] for data in v_data]
    T_frames = [@lift data[:,$frame] for data in T_data]

    uw_frames = [@lift data[:,$frame] for data in uw_data]
    vw_frames = [@lift data[:,$frame] for data in vw_data]
    wT_frames = [@lift data[:,$frame] for data in wT_data]

    Ri_frames = [@lift data[:,$frame] for data in Ri_data]

    losses_point_frames = [@lift [data[$frame]] for data in losses_data]

    @inline function find_lims(datasets)
        return maximum(maximum.(datasets)), minimum(minimum.(datasets))
    end

    u_max, u_min = find_lims(u_data)
    v_max, v_min = find_lims(v_data)
    T_max, T_min = find_lims(T_data)

    uw_max, uw_min = find_lims(uw_data)
    vw_max, vw_min = find_lims(vw_data)
    wT_max, wT_min = find_lims(wT_data)

    losses_max, losses_min = find_lims(losses_data)

    train_parameters = data_diffeq_explicit["train_parameters"]
    ν₀ = train_parameters.ν₀
    ν₋ = train_parameters.ν₋
    ΔRi = train_parameters.ΔRi
    Riᶜ = train_parameters.Riᶜ
    Pr = train_parameters.Pr
    loss_scalings = train_parameters.loss_scalings

    BC_str = @sprintf "Momentum Flux = %.1e m² s⁻², Buoyancy Flux = %.1e m² s⁻³" data_diffeq_explicit["truth_uw"][end, 1] data_diffeq_explicit["truth_wT"][end, 1]
    plot_title = @lift "$animation_type Data: $BC_str, Time = $(round(times[$frame], digits=2)) days"

    diffusivity_str = @sprintf "ν₀ = %.1e m² s⁻¹, ν₋ = %.1e m² s⁻¹, ΔRi = %.1e, Riᶜ = %.2f, Pr=%.1f" ν₀ ν₋ ΔRi Riᶜ Pr

    scaling_str = @sprintf "Loss Scalings: u = %.1e, v = %.1e, T = %.1e, ∂u∂z = %.1e, ∂v∂z = %.1e, ∂T∂z = %.1e" loss_scalings.u loss_scalings.v loss_scalings.T loss_scalings.∂u∂z loss_scalings.∂v∂z loss_scalings.∂T∂z

    plot_subtitle = "$n_trainings Training Simulations ($training_types): $diffusivity_str \n $scaling_str"

    fig = Figure(resolution=(1920, 1080))

    colors = distinguishable_colors(length(uw_data)+1, [RGB(1,1,1), RGB(0,0,0)], dropseed=true)

    u_str = "u / m s⁻¹"
    v_str = "v / m s⁻¹"
    T_str = "T / °C"
    uw_str = "uw / m² s⁻²"
    vw_str = "vw / m² s⁻²"
    wT_str = "wT / m s⁻¹ °C"

    zc = data_diffeq_explicit["depth_profile"]
    zf = data_diffeq_explicit["depth_flux"]
    zf_interior = zf[2:end-1]
    z_str = "z / m"

    alpha=0.5
    truth_linewidth = 7
    linewidth = 3

    ax_u = fig[1, 1] = Axis(fig, xlabel=u_str, ylabel=z_str)
    ax_v = fig[1, 2] = Axis(fig, xlabel=v_str, ylabel=z_str)
    ax_T = fig[1, 3:4] = Axis(fig, xlabel=T_str, ylabel=z_str)
    ax_Ri = fig[1, 5] = Axis(fig, xlabel="Ri", ylabel=z_str)
    ax_uw = fig[2, 1] = Axis(fig, xlabel=uw_str, ylabel=z_str)
    ax_vw = fig[2, 2] = Axis(fig, xlabel=vw_str, ylabel=z_str)
    ax_wT = fig[2, 3] = Axis(fig, xlabel=wT_str, ylabel=z_str)
    ax_losses = fig[2, 4] = Axis(fig, xlabel="Time / days", ylabel="Loss", yscale=CairoMakie.log10)

    legend_sublayout = GridLayout()
    fig[2, 5] = legend_sublayout

    CairoMakie.xlims!(ax_u, u_min, u_max)
    CairoMakie.xlims!(ax_v, v_min, v_max)
    CairoMakie.xlims!(ax_T, T_min, T_max)
    CairoMakie.xlims!(ax_uw, uw_min, uw_max)
    CairoMakie.xlims!(ax_vw, vw_min, vw_max)
    CairoMakie.xlims!(ax_wT, wT_min, wT_max)
    CairoMakie.xlims!(ax_Ri, -1, 2)
    CairoMakie.xlims!(ax_losses, times[1], times[end])

    CairoMakie.ylims!(ax_u, minimum(zc), 0)
    CairoMakie.ylims!(ax_v, minimum(zc), 0)
    CairoMakie.ylims!(ax_T, minimum(zc), 0)
    CairoMakie.ylims!(ax_uw, minimum(zf), 0)
    CairoMakie.ylims!(ax_vw, minimum(zf), 0)
    CairoMakie.ylims!(ax_wT, minimum(zf), 0)
    CairoMakie.ylims!(ax_Ri, minimum(zf), 0)
    CairoMakie.ylims!(ax_losses, losses_min, losses_max)

    u_lines = [
         lines!(ax_u, u_frames[1], zc, linewidth=truth_linewidth, color=(colors[1], alpha));
        [lines!(ax_u, u_frames[i], zc, linewidth=linewidth, color=colors[i]) for i in 2:length(u_data)]
    ]

    v_lines = [
         lines!(ax_v, v_frames[1], zc, linewidth=truth_linewidth, color=(colors[1], alpha));
        [lines!(ax_v, v_frames[i], zc, linewidth=linewidth, color=colors[i]) for i in 2:length(v_data)]
    ]

    T_lines = [
         lines!(ax_T, T_frames[1], zc, linewidth=truth_linewidth, color=(colors[1], alpha));
        [lines!(ax_T, T_frames[i], zc, linewidth=linewidth, color=colors[i]) for i in 2:length(T_data)]
    ]

    uw_lines = [
         lines!(ax_uw, uw_frames[1], zf, linewidth=truth_linewidth, color=(colors[1], alpha));
         lines!(ax_uw, uw_frames[end], zf_interior, linewidth=3, color=colors[end-1]);
        [lines!(ax_uw, uw_frames[i], zf, linewidth=linewidth, color=colors[i]) for i in 2:length(uw_data)-1]
    ]

   vw_lines = [
         lines!(ax_vw, vw_frames[1], zf, linewidth=truth_linewidth, color=(colors[1], alpha));
         lines!(ax_vw, vw_frames[end], zf_interior, linewidth=3, color=colors[end-1]);
        [lines!(ax_vw, vw_frames[i], zf, linewidth=linewidth, color=colors[i]) for i in 2:length(vw_data)-1]
    ]

    wT_lines = [
        lines!(ax_wT, wT_frames[1], zf, linewidth=truth_linewidth, color=(colors[1], alpha));
        lines!(ax_wT, wT_frames[end], zf_interior, linewidth=3, color=colors[end-1]);
       [lines!(ax_wT, wT_frames[i], zf, linewidth=linewidth, color=colors[i]) for i in 2:length(wT_data)-1]
   ]

    Ri_lines = [
         lines!(ax_Ri, Ri_frames[1], zf, linewidth=truth_linewidth, color=(colors[1], alpha));
        [lines!(ax_Ri, Ri_frames[i], zf, linewidth=linewidth, color=colors[i]) for i in 2:length(Ri_data)]
    ]

    losses_lines =  [
        [lines!(ax_losses, times, losses_data[i], linewidth=linewidth, color=colors[i+1]) for i in 1:Int(length(losses_data)/2)];
        [lines!(ax_losses, times, losses_data[Int(i+length(losses_data)/2)], linewidth=linewidth, color=colors[i+1], linestyle=:dot) for i in 1:Int(length(losses_data)/2)];
    ]

    losses_point = [CairoMakie.scatter!(ax_losses, time_point, point, color=colors[end]) for point in losses_point_frames]

    legend = Legend(fig, uw_lines, ["LES",
                                    "NN Only",
                                    "Modified Pac-Phil Only",
                                    "KPP",
                                    "Oceananigans, Implicit",
                                    "DiffEq, Explicit",
                                    "DiffEq, Implicit"
                                    ])

    legend_loss = Legend(fig, losses_lines, ["Profile Loss, Modified Pac-Phil Only",
                                             "Profile Loss, KPP",
                                             "Profile Loss, Oceananigans, Implicit",
                                             "Profile Loss, DiffEq, Explicit",
                                             "Profile Loss, DiffEq, Implicit",
                                             "Gradient Loss, Modified Pac-Phil Only",
                                             "Gradient Loss, KPP",
                                             "Gradient Loss, Oceananigans, Implicit",
                                             "Gradient Loss, DiffEq, Explicit",
                                             "Gradient Loss, DiffEq, Implicit"])
    legend_sublayout[:v] = [legend, legend_loss]

    supertitle = fig[0, :] = Label(fig, plot_title, textsize=25)
    subtitle = fig[end+1, :] = Label(fig, text=plot_subtitle, textsize=20)

    trim!(fig.layout)

    print_frame = maximum([1, Int(floor(length(times)/20))])

    function print_progress(n, n_total, print_frame, type)
        if n % print_frame == 0
            @info "Animating $(type) frame $n/$n_total"
        end
    end

    @info "Starting Animation"

    if gif
        CairoMakie.record(fig, "$FILE_PATH.gif", 1:length(times), framerate=fps) do n
            print_progress(n, length(times), print_frame, "gif")
            frame[] = n
        end
    end

    if mp4
        CairoMakie.record(fig, "$FILE_PATH.mp4", 1:length(times), framerate=fps) do n
            print_progress(n, length(times), print_frame, "mp4")
            frame[] = n
        end
    end
end

function animate_training_data_profiles_fluxes(train_files, FILE_PATH; fps=30, gif=false, mp4=true)
    all_data = [WindMixing.data(train_file, scale_type=ZeroMeanUnitVarianceScaling, enforce_surface_fluxes=false) for train_file in train_files]

    times = all_data[1].t ./ 86400

    frame = Node(1)

    time_point = @lift [times[$frame]]

    u_data = [data.u.coarse for data in all_data]
    v_data = [data.v.coarse for data in all_data]
    T_data = [data.T.coarse for data in all_data]

    us = [@lift u[:,$frame] for u in u_data]
    vs = [@lift v[:,$frame] for v in v_data]
    Ts = [@lift T[:,$frame] for T in T_data]

    uws = [@lift data.uw.coarse[:,$frame] for data in all_data]
    vws = [@lift data.vw.coarse[:,$frame] for data in all_data]
    wTs = [@lift data.wT.coarse[:,$frame] for data in all_data]

    u_max = maximum(maximum(data.u.coarse) for data in all_data)
    u_min = minimum(minimum(data.u.coarse) for data in all_data)

    v_max = maximum(maximum(data.v.coarse) for data in all_data)
    v_min = minimum(minimum(data.v.coarse) for data in all_data)

    T_max = maximum(maximum(data.T.coarse) for data in all_data)
    T_min = minimum(minimum(data.T.coarse) for data in all_data)

    uw_max = maximum(maximum(data.uw.coarse) for data in all_data)
    uw_min = minimum(minimum(data.uw.coarse) for data in all_data)

    vw_max = maximum(maximum(data.vw.coarse) for data in all_data)
    vw_min = minimum(minimum(data.vw.coarse) for data in all_data)

    wT_max = maximum(maximum(data.wT.coarse) for data in all_data)
    wT_min = minimum(minimum(data.wT.coarse) for data in all_data)

    Nz = all_data[1].grid_points - 1
    zc = all_data[1].u.z
    zf = all_data[1].uw.z

    D_face = Float32.(Dᶠ(Nz, zc[2] - zc[1]))

    @inline function ∂_∂z(profile)
        output = zeros(typeof(profile[1]), size(profile, 1) + 1, size(profile,2))
        for i in 1:size(profile,2)
            profile_col = @view profile[:,i]
            output_col = @view output[:,i]
            output_col .= D_face * profile_col
        end
        return output
    end

    ∂u∂zs = [∂_∂z(u) for u in u_data]
    ∂v∂zs = [∂_∂z(v) for v in v_data]
    ∂T∂zs = [∂_∂z(T) for T in T_data]

    @inline function local_richardson(∂u∂z, ∂v∂z, ∂T∂z, g, α)
        # ϵ = eps(typeof(∂u∂z))
        ϵ = 0
        ∂u∂z += ϵ
        ∂v∂z += ϵ
        ∂T∂z += ϵ
        Bz = g * α * ∂T∂z
        S² = ∂u∂z ^2 + ∂v∂z ^2
        return clamp.(Bz / S², -1, 2)
    end

    α = 1.67f-4
    g = 9.81f0

    Ris_data = [local_richardson.(∂u∂zs[i], ∂v∂zs[i], ∂T∂zs[i], g, α) for i in 1:length(∂u∂zs)]

    Ris = [@lift Ri[:,$frame] for Ri in Ris_data]

    plot_title = @lift "LES Simulations: Time = $(round(times[$frame], digits=2)) days"
    fig = Figure(resolution=(1920, 1080))
    color_palette = distinguishable_colors(length(all_data), [RGB(1,1,1), RGB(0,0,0)], dropseed=true)

    u_str = "u / m s⁻¹"
    v_str = "v / m s⁻¹"
    T_str = "T / °C"
    uw_str = "uw / m² s⁻²"
    vw_str = "vw / m² s⁻²"
    wT_str = "wT / m s⁻¹ °C"

    z_str = "z / m"

    uw_tops = [data.uw.coarse[end,1] for data in all_data]
    wT_tops = [data.wT.coarse[end,1] for data in all_data]

    BC_strs = [@sprintf "Momentum Flux = %.1e m² s⁻², Buoyancy Flux = %.1e m² s⁻³" uw_tops[i] wT_tops[i] for i in 1:length(uw_tops)]

    ax_u = fig[1, 1] = Axis(fig, xlabel=u_str, ylabel=z_str)
    u_lines = [lines!(ax_u, us[i], zc, linewidth=3, color=color_palette[i]) for i in 1:length(us)]
                # lines!(ax_u, test_u_NN_only, zc, linewidth=3, color=colors[4])]
    CairoMakie.xlims!(ax_u, u_min, u_max)
    CairoMakie.ylims!(ax_u, minimum(zc), 0)

    ax_v = fig[1, 2] = Axis(fig, xlabel=v_str, ylabel=z_str)
    v_lines = [lines!(ax_v, vs[i], zc, linewidth=3, color=color_palette[i]) for i in 1:length(vs)]
    CairoMakie.xlims!(ax_v, v_min, v_max)
    CairoMakie.ylims!(ax_v, minimum(zc), 0)

    ax_T = fig[1, 3] = Axis(fig, xlabel=T_str, ylabel=z_str)
    T_lines = [lines!(ax_T, Ts[i], zc, linewidth=3, color=color_palette[i]) for i in 1:length(Ts)]
    CairoMakie.xlims!(ax_T, T_min, T_max)
    CairoMakie.ylims!(ax_T, minimum(zc), 0)

    ax_uw = fig[2, 1] = Axis(fig, xlabel=uw_str, ylabel=z_str)
    uw_lines = [lines!(ax_uw, uws[i], zf, linewidth=3, color=color_palette[i]) for i in 1:length(uws)]
    CairoMakie.xlims!(ax_uw, uw_min, uw_max)
    CairoMakie.ylims!(ax_uw, minimum(zf), 0)

    ax_vw = fig[2, 2] = Axis(fig, xlabel=vw_str, ylabel=z_str)
    vw_lines = [lines!(ax_vw, vws[i], zf, linewidth=3, color=color_palette[i]) for i in 1:length(vws)]
    CairoMakie.xlims!(ax_vw, vw_min, vw_max)
    CairoMakie.ylims!(ax_vw, minimum(zf), 0)

    ax_wT = fig[2, 3] = Axis(fig, xlabel=wT_str, ylabel=z_str)
    wT_lines = [lines!(ax_wT, wTs[i], zf, linewidth=3, color=color_palette[i]) for i in 1:length(wTs)]
    CairoMakie.xlims!(ax_wT, wT_min, wT_max)
    CairoMakie.ylims!(ax_wT, minimum(zf), 0)

    ax_Ri = fig[2, 4] = Axis(fig, xlabel="Ri", ylabel=z_str)
    Ri_lines = [lines!(ax_Ri, Ris[i], zf, linewidth=3, color=color_palette[i]) for i in 1:length(Ris)]
                # lines!(ax_Ri, test_Ri_NN_only, zf, linewidth=3, color=colors[4])]
    CairoMakie.xlims!(ax_Ri, -1, 2)
    CairoMakie.ylims!(ax_Ri, minimum(zf), 0)

    legend = fig[1, 4] = Legend(fig, uw_lines, BC_strs)

    # legend = fig[1, 4] = Legend(fig, u_lines, ["Oceananigans.jl LES", "NN + Modified Pac-Phil", "Modified Pac-Phil Only"])
    supertitle = fig[0, :] = Label(fig, plot_title, textsize=30)
    trim!(fig.layout)

    if gif
        CairoMakie.record(fig, "$FILE_PATH.gif", 1:length(times), framerate=fps) do n
            @info "Animating gif frame $n/$(length(times))..."
            frame[] = n
        end
    end

    if mp4
        CairoMakie.record(fig, "$FILE_PATH.mp4", 1:length(times), framerate=fps) do n
            @info "Animating mp4 frame $n/$(length(times))..."
            frame[] = n
        end
    end
end

function generate_training_types_str(FILE_NAME)
    training_types = ""
    check_exists(str) = occursin(str, FILE_NAME)

    if check_exists("_wind_mixing_")
        training_types *= "Wind Mixing"
    end

    if check_exists("_cooling_")
        if training_types != ""
            training_types *= ", Cooling"
        else
            training_types *= "Cooling"
        end
    end

    if check_exists("_heating_")
        if training_types != ""
            training_types *= ", Heating"
        else
            training_types *= "Heating"
        end
    end

    if check_exists("_windcooling_")
        if training_types != ""
            training_types *= ", Wind + Cooling"
        else
            training_types *= "Wind + Cooling"
        end
    end

    if check_exists("_windheating_")
        if training_types != ""
            training_types *= ", Wind + Heating"
        else
            training_types *= "Wind + Heating"
        end
    end

    return training_types
end

function animate_training_results(test_files, FILE_NAME;
                                  EXTRACTED_DATA_DIR, OUTPUT_DIR,
                                  timestep=60, trange=1:1:1153,
                                  fps=30, gif=false, mp4=true,
                                  explicit_timestepper=ROCK4(), implicit_timestepper=RadauIIA5(autodiff=false),
                                  convective_adjustment=false)

    DATA_PATH = joinpath(EXTRACTED_DATA_DIR, "$(FILE_NAME)_extracted.jld2")
    OUTPUT_PATH = joinpath(OUTPUT_DIR, FILE_NAME)

    if !ispath(OUTPUT_PATH)
        mkdir(OUTPUT_PATH)
    end

    @info "Loading Data"
    file = jldopen(DATA_PATH, "r")
    if haskey(file, "losses/total")
        losses = (
            u = file["losses/u"],
            v = file["losses/v"],
            T = file["losses/T"],
            ∂u∂z = file["losses/∂u∂z"],
            ∂v∂z = file["losses/∂v∂z"],
            ∂T∂z = file["losses/∂T∂z"],
        )
        @info "Training Loss = $(minimum(file["losses/total"]))"
    else
        losses = file["losses"]
        @info "Training Loss = $(minimum(losses))"
    end

    train_files = file["training_info/train_files"]

    diurnal = occursin("diurnal", train_files[1])

    train_parameters = file["training_info/parameters"]
    uw_NN = file["neural_network/uw"]
    vw_NN = file["neural_network/vw"]
    wT_NN = file["neural_network/wT"]

    if haskey(file["training_info"], "loss_scalings")
        loss_scalings = file["training_info/loss_scalings"]
    elseif haskey(train_parameters, "gradient_scaling")
        gradient_scaling = train_parameters["gradient_scaling"]
        loss_scalings = (u=1f0, v=1f0, T=1f0, ∂u∂z=gradient_scaling, ∂v∂z=gradient_scaling, ∂T∂z=gradient_scaling)
    else
        loss_scalings = (u=1f0, v=1f0, T=1f0, ∂u∂z=1f0, ∂v∂z=1f0, ∂T∂z=1f0)
    end

    close(file)

    @info "Loading Training Data"
    𝒟train = WindMixing.data(train_files, scale_type=ZeroMeanUnitVarianceScaling, enforce_surface_fluxes=false)
    training_types = generate_training_types_str(FILE_NAME)

    ν₀ = train_parameters["ν₀"]
    ν₋ = train_parameters["ν₋"]
    ΔRi = train_parameters["ΔRi"]
    Riᶜ = train_parameters["Riᶜ"]
    Pr = train_parameters["Pr"]

    @info "Plotting Loss..."
    plot_loss(losses, joinpath(OUTPUT_PATH, "loss.pdf"))

    @info "Solving NDE: Oceananigans"
    solve_oceananigans_modified_pacanowski_philander_nn(test_files, DATA_PATH, OUTPUT_PATH; timestep=timestep, convective_adjustment=convective_adjustment)

    for test_file in test_files
        @info "Generating Data: $test_file"
        𝒟test = WindMixing.data(test_file, scale_type=ZeroMeanUnitVarianceScaling, enforce_surface_fluxes=false)

        @info "Processing $test_file solution"
        if test_file in train_files
            SOL_DIR = joinpath(OUTPUT_PATH, "train_$test_file")
            animation_type = "Training"
        else
            SOL_DIR = joinpath(OUTPUT_PATH, "test_$test_file")
            animation_type = "Testing"
        end

        @info "Solving NDE: $test_file, Explicit timestepper"
        plot_data_explicit = NDE_profile(uw_NN, vw_NN, wT_NN, test_file, 𝒟test, 𝒟train, trange,
                                modified_pacanowski_philander=train_parameters["modified_pacanowski_philander"],
                                ν₀=ν₀, ν₋=ν₋, ΔRi=ΔRi, Riᶜ=Riᶜ, Pr=Pr,
                                convective_adjustment=false,
                                smooth_NN=train_parameters["smooth_NN"], smooth_Ri=train_parameters["smooth_Ri"],
                                zero_weights=train_parameters["zero_weights"],
                                loss_scalings=loss_scalings,
                                timestepper=explicit_timestepper,
                                OUTPUT_PATH=joinpath(SOL_DIR, "solution_diffeq_explicit.jld2"))

        @info "Solving NDE: $test_file, Implicit timestepper"
        plot_data_implicit = NDE_profile(uw_NN, vw_NN, wT_NN, test_file, 𝒟test, 𝒟train, trange,
                                modified_pacanowski_philander=train_parameters["modified_pacanowski_philander"],
                                ν₀=ν₀, ν₋=ν₋, ΔRi=ΔRi, Riᶜ=Riᶜ, Pr=Pr,
                                convective_adjustment=false,
                                smooth_NN=train_parameters["smooth_NN"], smooth_Ri=train_parameters["smooth_Ri"],
                                zero_weights=train_parameters["zero_weights"],
                                loss_scalings=loss_scalings,
                                timestepper=implicit_timestepper,
                                OUTPUT_PATH=joinpath(SOL_DIR, "solution_diffeq_implicit.jld2"))

        @info "Solving NDE: $test_file, Oceananigans"

        plot_data_oceananigans = NDE_profile_oceananigans(SOL_DIR, train_files, [test_file],
                                            ν₀=ν₀, ν₋=ν₋, ΔRi=ΔRi, Riᶜ=Riᶜ, Pr=Pr,
                                            loss_scalings=loss_scalings,
                                            OUTPUT_PATH=joinpath(SOL_DIR, "solution_oceananigans.jld2"))

        if test_file in train_files
            animation_type = "Training"
        else
            animation_type = "Testing"
        end
        n_trainings = length(train_files)

        if animation_type == "Training"
            VIDEO_NAME = "train_$test_file"
        else
            VIDEO_NAME = "test_$test_file"
        end

        VIDEO_PATH = joinpath(SOL_DIR, "$VIDEO_NAME")

        @info "Animating $test_file Video"

        animate_profiles_fluxes_comparison(plot_data_explicit, plot_data_implicit, plot_data_oceananigans, VIDEO_PATH, fps=fps, gif=gif, mp4=mp4,
                                                animation_type=animation_type, n_trainings=n_trainings, training_types=training_types)
        @info "$test_file Animation Completed"
    end
end

function animate_profiles_fluxes_final(data, axis_images, FILE_PATH; animation_type, n_trainings, training_types, fps=30, gif=false, mp4=true)
    times = data["t"] ./ 86400

    frame = Node(1)

    time_point = @lift [times[$frame]]

    u_data = [
        data["truth_u"],
        data["test_u_modified_pacanowski_philander"],
        data["test_u_kpp"],
        data["test_u"],
    ]

    v_data = [
        data["truth_v"],
        data["test_v_modified_pacanowski_philander"],
        data["test_v_kpp"],
        data["test_v"],
    ]

    T_data = [
        data["truth_T"],
        data["test_T_modified_pacanowski_philander"],
        data["test_T_kpp"],
        data["test_T"],
    ]

    uw_data = [
        data["truth_uw"],
        data["test_uw_modified_pacanowski_philander"],
        data["test_uw_kpp"],
        data["test_uw"],
    ]

    vw_data = [
        data["truth_vw"],
        data["test_vw_modified_pacanowski_philander"],
        data["test_vw_kpp"],
        data["test_vw"],
    ]

    wT_data = [
        data["truth_wT"],
        data["test_wT_modified_pacanowski_philander"],
        data["test_wT_kpp"],
        data["test_wT"],
    ]

    uw_data .*= 1f4
    vw_data .*= 1f4
    wT_data .*= 1f5

    Ri_data = [
        clamp.(data["truth_Ri"], -1, 2),
        clamp.(data["test_Ri_modified_pacanowski_philander"], -1, 2),
        clamp.(data["test_Ri_kpp"], -1, 2),
        clamp.(data["test_Ri"], -1, 2),
    ]

    # @inline function lowclamp(value, lo)
    #     if value >= lo
    #         return value
    #     else
    #         return lo
    #     end
    # end

    # losses_data = [
    #     lowclamp.(data["losses_modified_pacanowski_philander"] .+ data["losses_modified_pacanowski_philander_gradient"], 1f-5),
    #     lowclamp.(data["losses_kpp"] .+ data["losses_kpp_gradient"], 1f-5),
    #     lowclamp.(data["losses"] .+ data["losses_gradient"], 1f-5),
    # ]

    u_frames = [@lift data[:,$frame] for data in u_data]
    v_frames = [@lift data[:,$frame] for data in v_data]
    T_frames = [@lift data[:,$frame] for data in T_data]

    uw_frames = [@lift data[:,$frame] for data in uw_data]
    vw_frames = [@lift data[:,$frame] for data in vw_data]
    wT_frames = [@lift data[:,$frame] for data in wT_data]

    Ri_frames = [@lift data[:,$frame] for data in Ri_data]

    # losses_point_frames = [@lift [data[$frame]] for data in losses_data]

    @inline function find_lims(datasets)
        return maximum(maximum.(datasets)), minimum(minimum.(datasets))
    end

    u_max, u_min = find_lims(u_data)
    v_max, v_min = find_lims(v_data)
    T_max, T_min = find_lims(T_data)

    uw_max, uw_min = find_lims(uw_data)
    vw_max, vw_min = find_lims(vw_data)
    wT_max, wT_min = find_lims(wT_data)

    # losses_max, losses_min = find_lims(losses_data)

    train_parameters = data["train_parameters"]
    ν₀ = train_parameters.ν₀
    ν₋ = train_parameters.ν₋
    ΔRi = train_parameters.ΔRi
    Riᶜ = train_parameters.Riᶜ
    Pr = train_parameters.Pr
    loss_scalings = train_parameters.loss_scalings

    BC_str = @sprintf "Momentum Flux = %.1e m² s⁻², Temperature Flux = %.1e m s⁻¹ °C" data["truth_uw"][end, 1] maximum(data["truth_wT"][end, :])
    plot_title = @lift "$animation_type Data: $BC_str, Time = $(round(times[$frame], digits=2)) days"

    diffusivity_str = @sprintf "ν₀ = %.2e m² s⁻¹, ν₋ = %.2e m² s⁻¹, ΔRi = %.2e, Riᶜ = %.3f, Pr=%.2f" ν₀ ν₋ ΔRi Riᶜ Pr

    # scaling_str = @sprintf "Loss Scalings: u = %.1e, v = %.1e, T = %.1e, ∂u∂z = %.1e, ∂v∂z = %.1e, ∂T∂z = %.1e" loss_scalings.u loss_scalings.v loss_scalings.T loss_scalings.∂u∂z loss_scalings.∂v∂z loss_scalings.∂T∂z
    # plot_subtitle = "$n_trainings Training Simulations ($training_types): $diffusivity_str \n $scaling_str"

    plot_subtitle = "$n_trainings Training Simulations ($training_types): $diffusivity_str"

    # fig = Figure(resolution=(1920, 1080))
    fig = Figure(resolution=(1920, 960))

    # colors = distinguishable_colors(length(uw_data)+1, [RGB(1,1,1), RGB(0,0,0)], dropseed=true)
    colors = distinguishable_colors(length(uw_data), [RGB(1,1,1), RGB(0,0,0)], dropseed=true)

    temp_color = colors[2]
    colors[2] = colors[4]
    colors[4] = temp_color
    # colors[4] = RGB(0,0,0)

    u_img = axis_images.u
    v_img = axis_images.v
    T_img = axis_images.T
    uw_img = axis_images.uw
    vw_img = axis_images.uw
    wT_img = axis_images.uw
    Ri_img = axis_images.uw
    z_img = axis_images.z

    zc = data["depth_profile"]
    zf = data["depth_flux"]
    zf_interior = zf[2:end-1]

    rel_size = 30
    # aspect = 1920 / 1080
    aspect = 2

    ax_u = fig[1, 2] = Axis(fig)
    ax_v = fig[1, 4] = Axis(fig)

    T_layout = fig[1:4, 5] = GridLayout()
    colsize!(fig.layout, 5, CairoMakie.Relative(0.4))

    ax_T = T_layout[1, 2] = Axis(fig)
    y_ax_T = T_layout[1,1] = CairoMakie.Axis(fig, aspect=DataAspect())
    x_ax_T = T_layout[2,2] = CairoMakie.Axis(fig, aspect=DataAspect())

    ax_Ri = fig[1, 7] = Axis(fig)
    ax_uw = fig[3, 2] = Axis(fig)
    ax_vw = fig[3, 4] = Axis(fig)
    ax_wT = fig[3, 7] = Axis(fig)

    y_ax_u = CairoMakie.Axis(fig[1,1], aspect=DataAspect())
    y_ax_v = CairoMakie.Axis(fig[1,3], aspect=DataAspect())
    y_ax_Ri = CairoMakie.Axis(fig[1,6], aspect=DataAspect())
    y_ax_uw = CairoMakie.Axis(fig[3,1], aspect=DataAspect())
    y_ax_vw = CairoMakie.Axis(fig[3,3], aspect=DataAspect())
    y_ax_wT = CairoMakie.Axis(fig[3,6], aspect=DataAspect())

    x_ax_u = CairoMakie.Axis(fig[2,2], aspect=DataAspect())
    x_ax_v = CairoMakie.Axis(fig[2,4], aspect=DataAspect())
    x_ax_Ri = CairoMakie.Axis(fig[2,7], aspect=DataAspect())
    x_ax_uw = CairoMakie.Axis(fig[4,2], aspect=DataAspect())
    x_ax_vw = CairoMakie.Axis(fig[4,4], aspect=DataAspect())
    x_ax_wT = CairoMakie.Axis(fig[4,7], aspect=DataAspect())

    hidedecorations!(y_ax_u)
    hidedecorations!(y_ax_v)
    hidedecorations!(y_ax_Ri)
    hidedecorations!(y_ax_uw)
    hidedecorations!(y_ax_vw)
    hidedecorations!(y_ax_wT)

    hidedecorations!(x_ax_u)
    hidedecorations!(x_ax_v)
    hidedecorations!(x_ax_Ri)
    hidedecorations!(x_ax_uw)
    hidedecorations!(x_ax_vw)
    hidedecorations!(x_ax_wT)

    hidespines!(y_ax_u)
    hidespines!(y_ax_v)
    hidespines!(y_ax_T)
    hidespines!(y_ax_Ri)
    hidespines!(y_ax_uw)
    hidespines!(y_ax_vw)
    hidespines!(y_ax_wT)

    hidespines!(x_ax_u)
    hidespines!(x_ax_v)
    hidespines!(x_ax_T)
    hidespines!(x_ax_Ri)
    hidespines!(x_ax_uw)
    hidespines!(x_ax_vw)
    hidespines!(x_ax_wT)

    image!(x_ax_u, axis_images.u)
    image!(x_ax_v, axis_images.v)
    image!(x_ax_Ri, axis_images.Ri)
    image!(x_ax_uw, axis_images.uw)
    image!(x_ax_vw, axis_images.vw)
    image!(x_ax_wT, axis_images.wT)

    image!(y_ax_u, axis_images.z)
    image!(y_ax_v, axis_images.z)
    image!(y_ax_Ri, axis_images.z)
    image!(y_ax_uw, axis_images.z)
    image!(y_ax_vw, axis_images.z)
    image!(y_ax_wT, axis_images.z)

    hidedecorations!(y_ax_T)
    hidedecorations!(x_ax_T)
    image!(x_ax_T, axis_images.T)
    image!(y_ax_T, axis_images.z)

    colsize!(T_layout, 1, CairoMakie.Relative(1 / rel_size * 1.5))
    rowsize!(T_layout, 2, CairoMakie.Relative(1 / rel_size / aspect * 1.75))
    colgap!(T_layout, Relative(1 / rel_size / aspect / 2))
    rowgap!(T_layout, Relative(1 / rel_size / aspect))

    rowsize!(fig.layout, 2, CairoMakie.Relative(1 / rel_size))
    rowsize!(fig.layout, 4, CairoMakie.Relative(1 / rel_size))
    colsize!(fig.layout, 1, CairoMakie.Relative(1 / rel_size / aspect))
    colsize!(fig.layout, 3, CairoMakie.Relative(1 / rel_size / aspect))
    colsize!(fig.layout, 6, CairoMakie.Relative(1 / rel_size / aspect))

    colgap!(fig.layout, Relative(1 / rel_size / aspect / 2))
    rowgap!(fig.layout, Relative(1 / rel_size / aspect))

    alpha=0.5
    truth_linewidth =8
    linewidth = 3

    CairoMakie.xlims!(ax_u, u_min, u_max)
    CairoMakie.xlims!(ax_v, v_min, v_max)
    CairoMakie.xlims!(ax_T, T_min, T_max)
    CairoMakie.xlims!(ax_uw, uw_min, uw_max)
    CairoMakie.xlims!(ax_vw, vw_min, vw_max)
    CairoMakie.xlims!(ax_wT, wT_min, wT_max)
    CairoMakie.xlims!(ax_Ri, -1, 2)
    # CairoMakie.xlims!(ax_losses, times[1], times[end])

    CairoMakie.ylims!(ax_u, minimum(zc), 0)
    CairoMakie.ylims!(ax_v, minimum(zc), 0)
    CairoMakie.ylims!(ax_T, minimum(zc), 0)
    CairoMakie.ylims!(ax_uw, minimum(zf), 0)
    CairoMakie.ylims!(ax_vw, minimum(zf), 0)
    CairoMakie.ylims!(ax_wT, minimum(zf), 0)
    CairoMakie.ylims!(ax_Ri, minimum(zf), 0)
    # CairoMakie.ylims!(ax_losses, losses_min, losses_max)

    u_lines = [
         lines!(ax_u, u_frames[1], zc, linewidth=truth_linewidth, color=(colors[1], alpha));
        [lines!(ax_u, u_frames[i], zc, linewidth=linewidth, color=colors[i]) for i in 2:length(u_data)]
    ]

    v_lines = [
         lines!(ax_v, v_frames[1], zc, linewidth=truth_linewidth, color=(colors[1], alpha));
        [lines!(ax_v, v_frames[i], zc, linewidth=linewidth, color=colors[i]) for i in 2:length(v_data)]
    ]

    T_lines = [
         lines!(ax_T, T_frames[1], zc, linewidth=truth_linewidth, color=(colors[1], alpha));
        [lines!(ax_T, T_frames[i], zc, linewidth=linewidth, color=colors[i]) for i in 2:length(T_data)]
    ]

    uw_lines = [
         lines!(ax_uw, uw_frames[1], zf, linewidth=truth_linewidth, color=(colors[1], alpha));
        [lines!(ax_uw, uw_frames[i], zf, linewidth=linewidth, color=colors[i]) for i in 2:length(uw_data)]
    ]

   vw_lines = [
         lines!(ax_vw, vw_frames[1], zf, linewidth=truth_linewidth, color=(colors[1], alpha));
        [lines!(ax_vw, vw_frames[i], zf, linewidth=linewidth, color=colors[i]) for i in 2:length(vw_data)]
    ]

    wT_lines = [
        lines!(ax_wT, wT_frames[1], zf, linewidth=truth_linewidth, color=(colors[1], alpha));
       [lines!(ax_wT, wT_frames[i], zf, linewidth=linewidth, color=colors[i]) for i in 2:length(wT_data)]
   ]

    Ri_lines = [
         lines!(ax_Ri, Ri_frames[1], zf, linewidth=truth_linewidth, color=(colors[1], alpha));
        [lines!(ax_Ri, Ri_frames[i], zf, linewidth=linewidth, color=colors[i]) for i in 2:length(Ri_data)]
    ]

    axislegend(ax_T, T_lines, ["Oceananigans.jl Large Eddy Simulation", "Ri-based Diffusivity Only", "K-Profile Parameterisation", "NN Embedded in Oceananigans.jl"], "Data Type", position = :rb)

    # legend = fig[3, 8] = Legend(fig, u_lines, ["LES",
    #                                     "Ri-based Diffusivity Only",
    #                                     "KPP",
    #                                     "Oceananigans.jl",
    #                                     ])

    supertitle = fig[0, :] = Label(fig, plot_title, textsize=25)
    subtitle = fig[end+1, :] = Label(fig, text=plot_subtitle, textsize=20)

    trim!(fig.layout)

    print_frame = maximum([1, Int(floor(length(times)/20))])

    function print_progress(n, n_total, print_frame, type)
        if n % print_frame == 0
            @info "Animating $(type) frame $n/$n_total"
        end
    end

    @info "Starting Animation"

    if gif
        CairoMakie.record(fig, "$FILE_PATH.gif", 1:length(times), framerate=fps) do n
            print_progress(n, length(times), print_frame, "gif")
            frame[] = n
        end
    end

    if mp4
        CairoMakie.record(fig, "$FILE_PATH.mp4", 1:length(times), framerate=fps, compression=1) do n
            print_progress(n, length(times), print_frame, "mp4")
            frame[] = n
        end
    end

end

function animate_LES_3D(FILE_DIR, OUTPUT_PATH, axis_images; num_frames, simulation_str, fps=30, gif=false, mp4=true, colorscheme=:ice, rev=false)
    xy_file = jldopen(joinpath(FILE_DIR, "xy_slice.jld2"))
    xz_file = jldopen(joinpath(FILE_DIR, "xz_slice.jld2"))
    yz_file = jldopen(joinpath(FILE_DIR, "yz_slice.jld2"))
    instantaneous_statistics = jldopen(joinpath(FILE_DIR, "instantaneous_statistics.jld2"))

    iterations = keys(instantaneous_statistics["timeseries/t"])
    times = [instantaneous_statistics["timeseries/t/$iter"] for iter in iterations]
    xC = xz_file["grid/xC"][4:end-3]
    yC = xz_file["grid/yC"][4:end-3]
    zC = xz_file["grid/zC"][4:end-3]

    x_xz = fill(xC[1], 128)
    y_xz = yC

    z_xz = zeros(length(x_xz), length(y_xz))
    for i in 1:size(z_xz, 1)
        z_xz[i,:] .= zC[i]
    end

    T_xzs = [transpose(hcat([xz_file["timeseries/T/$iter"][:, :, i] for i in 1:length(zC)]... )) for iter in iterations]

    x_yz = xC
    y_yz = fill(yC[1], 128)

    z_yz = zeros(length(x_yz), length(y_yz))
    for i in 1:size(z_yz, 2)
        z_yz[:,i] .= zC[i]
    end

    @inline function obtain_T_yz(iter)
        T_yz = similar(z_yz)
        for i in 1:size(T_yz, 2)
            T_yz[:,i] = yz_file["timeseries/T/$iter"][:, :, i]
        end
        return T_yz
    end
    T_yzs = [obtain_T_yz(iter) for iter in iterations]

    x_xy = xC
    y_xy = yC
    z_xy = zeros(length(y_xy), length(x_xy))

    T_xys = [xy_file["timeseries/T/$iter"][:,:,1] for iter in iterations]

    us = [instantaneous_statistics["timeseries/u/$iter"][:] for iter in iterations]
    vs = [instantaneous_statistics["timeseries/v/$iter"][:] for iter in iterations]
    Ts = [instantaneous_statistics["timeseries/T/$iter"][:] for iter in iterations]

    @inline function find_lims(profiles)
        return minimum(minimum.(profiles)), maximum(maximum.(profiles))
    end

    u_min, u_max = find_lims(us)
    v_min, v_max = find_lims(vs)
    T_min, T_max = find_lims(Ts)

    color_range_max = maximum([maximum(maximum.(T_xzs)), maximum(maximum.(T_yzs)), maximum(maximum.(T_xys))])
    color_range_min = minimum([minimum(minimum.(T_xzs)), minimum(minimum.(T_yzs)), minimum(minimum.(T_xys))])

    color_range = (color_range_min, color_range_max)
    colormap = cgrad(colorscheme, scale=:log10, rev=rev)

    close(xy_file)
    close(xz_file)
    close(yz_file)
    close(instantaneous_statistics)

    frame = Node(1)
    iteration = @lift iterations[$frame]

    T_xz = @lift T_xzs[$frame]
    T_yz = @lift T_yzs[$frame]
    T_xy = @lift T_xys[$frame]

    u = @lift us[$frame]
    v = @lift vs[$frame]
    T = @lift Ts[$frame]

    fig = Figure(resolution=(1920, 1080))
    ax = fig[1,1] = CairoMakie.Axis3(fig, aspect=(1, 1, 0.5), xlabel="x /m", ylabel="y /m", zlabel="z /m")

    xz_surface = CairoMakie.surface!(ax, x_xz, y_xz, z_xz, color=T_xz, colormap=colormap, colorrange=color_range)

    yz_surface = CairoMakie.surface!(ax, x_yz, y_yz, z_yz, color=T_yz, colormap=colormap, colorrange=color_range)

    xy_surface = CairoMakie.surface!(ax, x_xy, y_xy, z_xy, color=T_xy, colormap=colormap, colorrange=color_range)

    ax_T_3D = fig[2,1] = CairoMakie.Axis(fig, aspect=DataAspect())

    rel_size = 40
    aspect = 1 / 4

    hidedecorations!(ax_T_3D)
    hidespines!(ax_T_3D)
    image!(ax_T_3D, axis_images.T_3D)
    rowsize!(fig.layout, 2, CairoMakie.Relative(1 / rel_size))

    colorbar = CairoMakie.Colorbar(fig[3, 1], xz_surface, vertical=false)

    plots_sublayout = fig[:,2] = GridLayout()

    colsize!(fig.layout, 2, CairoMakie.Relative(aspect))
    rowgap!(fig.layout, Relative(1 / 50 / 3))

    y_ax_u = plots_sublayout[1,1] = CairoMakie.Axis(fig, aspect=DataAspect())
    y_ax_v = plots_sublayout[3,1] = CairoMakie.Axis(fig, aspect=DataAspect())
    y_ax_T = plots_sublayout[5,1] = CairoMakie.Axis(fig, aspect=DataAspect())

    x_ax_u = plots_sublayout[2,2] = CairoMakie.Axis(fig, aspect=DataAspect())
    x_ax_v = plots_sublayout[4,2] = CairoMakie.Axis(fig, aspect=DataAspect())
    x_ax_T = plots_sublayout[6,2] = CairoMakie.Axis(fig, aspect=DataAspect())

    ax_u = plots_sublayout[1,2] = CairoMakie.Axis(fig)
    ax_v = plots_sublayout[3,2] = CairoMakie.Axis(fig)
    ax_T = plots_sublayout[5,2] = CairoMakie.Axis(fig)

    hidedecorations!(y_ax_u)
    hidedecorations!(y_ax_v)
    hidedecorations!(y_ax_T)

    hidedecorations!(x_ax_u)
    hidedecorations!(x_ax_v)
    hidedecorations!(x_ax_T)

    hidespines!(y_ax_u)
    hidespines!(y_ax_v)
    hidespines!(y_ax_T)

    hidespines!(x_ax_u)
    hidespines!(x_ax_v)
    hidespines!(x_ax_T)

    image!(x_ax_u, axis_images.u)
    image!(x_ax_v, axis_images.v)
    image!(x_ax_T, axis_images.T)

    image!(y_ax_u, axis_images.z)
    image!(y_ax_v, axis_images.z)
    image!(y_ax_T, axis_images.z)

    rowsize!(plots_sublayout, 2, CairoMakie.Relative(1 / rel_size))
    rowsize!(plots_sublayout, 4, CairoMakie.Relative(1 / rel_size))
    rowsize!(plots_sublayout, 6, CairoMakie.Relative(1 / rel_size))
    colsize!(plots_sublayout, 1, CairoMakie.Aspect(2, 1))

    colgap!(plots_sublayout, 1 / 50)
    rowgap!(plots_sublayout, Relative(1 / 50 / 3))

    u_line = CairoMakie.lines!(ax_u, u, zC)
    v_line = CairoMakie.lines!(ax_v, v, zC)
    T_line = CairoMakie.lines!(ax_T, T, zC)

    CairoMakie.xlims!(ax_u, u_min, u_max)
    CairoMakie.xlims!(ax_v, v_min, v_max)
    CairoMakie.xlims!(ax_T, T_min, T_max)

    CairoMakie.ylims!(ax_u, minimum(zC), 0)
    CairoMakie.ylims!(ax_v, minimum(zC), 0)
    CairoMakie.ylims!(ax_T, minimum(zC), 0)

    plot_title = @lift "$(simulation_str), Time = $(round(times[$frame]/86400, digits=2)) days"

    supertitle = fig[0, :] = Label(fig, plot_title, textsize=25)

    trim!(fig.layout)

    function print_progress(n, n_total, type)
        @info "Animating $(type) frame $n/$n_total"
    end

    @info "Starting Animation"

    if num_frames == 0
        if gif
            CairoMakie.record(fig, "$OUTPUT_PATH.gif", 1:length(times), framerate=fps, compression=1) do n
                print_progress(n, length(times), "gif")
                frame[] = n
            end
        end

        if mp4
            CairoMakie.record(fig, "$OUTPUT_PATH.mp4", 1:length(times), framerate=fps, compression=1) do n
                print_progress(n, length(times), "mp4")
                frame[] = n
            end
        end
    else
        if gif
            CairoMakie.record(fig, "$OUTPUT_PATH.gif", 1:num_frames, framerate=fps, compression=1) do n
                print_progress(n, num_frames, "gif")
                frame[] = n
            end
        end

        if mp4
            CairoMakie.record(fig, "$OUTPUT_PATH.mp4", 1:num_frames, framerate=fps, compression=1) do n
                print_progress(n, num_frames, "mp4")
                frame[] = n
            end
        end
    end
end


text:
export Codecov

"""
Coverage.Codecov Module

This module provides functionality to push coverage information to the CodeCov.io
web service. It exports the `submit` and `submit_token` methods.
"""

module Codecov
    using Requests
    using Coverage
    using JSON
    using Compat

    export submit, submit_token

    #=
    JSON structure for Codecov.io
    https://codecov.io/api#post-report

    {
      "coverage": {
        "path/to/file.py": [null, 1, 0, null, true, 0, 0, 1, 1],
        "path/to/other.py": [null, 0, 1, 1, "1/3", null]
      },
      "messages": {
        "path/to/other.py": {
          "1": "custom message for line 1"
        }
      }
    }
    =#

    # Turn vector of filename : coverage pairs into a dictionary
    function to_json(fcs::Vector{FileCoverage})
        cov = Dict()
        for fc in fcs
            cov[fc.filename] = vcat(nothing, fc.coverage)
        end
        return @compat Dict("coverage" => cov)
    end

    """
        submit(fcs::Vector{FileCoverage})

    Take a vector of file coverage results (produced by `process_folder`),
    and submits them to Codecov.io. Assumes that this code is being run
    on TravisCI. If running locally, use `submit_token`.
    """
    function submit(fcs::Vector{FileCoverage})
        branch  = ENV["TRAVIS_BRANCH"]
        pr      = ENV["TRAVIS_PULL_REQUEST"]
        job     = ENV["TRAVIS_JOB_ID"]
        slug    = ENV["TRAVIS_REPO_SLUG"]
        build   = ENV["TRAVIS_JOB_NUMBER"]
        commit  = ENV["TRAVIS_COMMIT"]
        uri_str = "https://codecov.io/upload/v2?service=travis-org" *
                    "&branch=$(branch)&commit=$(commit)" *
                    "&build=$(build)&pull_request=$(pr)&" *
                    "job=$(job)&slug=$(slug)"
        println("Codecov.io API URL:")
        println(uri_str)

        heads   = @compat Dict("Content-Type" => "application/json")
        data    = to_json(fcs)
        req     = Requests.post(URI(uri_str); json = data, headers = heads)
        println("Result of submission:")
        dump(req.data)
    end

    import Base.Git

    """
        submit_token(fcs::Vector{FileCoverage},
                     commit=Git.readchomp(`rev-parse HEAD`, dir=""),
                     branch=Git.branch(dir=""))

    Take a `Vector` of file coverage results (produced by `process_folder`),
    and submits them to Codecov.io. For submissions not from TravisCI.
    """
    function submit_token(fcs::Vector{FileCoverage},
                            commit=Git.readchomp(`rev-parse HEAD`, dir=""),
                            branch=Git.branch(dir=""))
        repo_token = ENV["REPO_TOKEN"]
        uri_str = "https://codecov.io/upload/v2?&token=$(repo_token)&commit=$(commit)&branch=$(branch)"
        heads   = @compat Dict("Content-Type" => "application/json")
        data    = to_json(fcs)
        req     = Requests.post(URI(uri_str); json = data, headers = heads)
        println("Result of submission:")
        dump(req.data)
    end
end  # module Codecov


text:
#=
tests for PolyUtils module
=#

using Test
import Random
using LinearAlgebra
import Hypatia.PolyUtils
import Hypatia.PolyUtils.interpolate

function test_interp_domain(T::Type{<:Real})
    Random.seed!(1)
    @testset "domain interpolation n = $n, halfdeg = $halfdeg" for (n, halfdeg) in
        ((1, 5), (2, 2), (3, 1))
        free_dom = PolyUtils.FreeDomain{T}(n)
        box_dom = PolyUtils.BoxDomain{T}(-ones(T, n), ones(T, n))
        ball_dom = PolyUtils.BallDomain{T}(ones(T, n), one(T))
        Q = rand(T, n, n)
        Q = Symmetric(Q * Q')
        ellip_dom = PolyUtils.EllipsoidDomain{T}(ones(T, n), Q)

        for sample in (true, false)
            interp_options = (; sample = sample, sample_factor = 10)

            free = interpolate(free_dom, halfdeg; interp_options...)
            box = interpolate(box_dom, halfdeg; interp_options...)
            @test length(free.Ps) == 1
            @test length(box.Ps) == 1 + n
            @test free.U == box.U
            @test size(box.pts) == size(free.pts)
            @test size(box.Ps[1]) == size(free.Ps[1])
            @test norm(box.Ps[1]) ≈ norm(free.Ps[1]) atol=1e-1 rtol=1e-1

            if sample
                ball = interpolate(ball_dom, halfdeg; interp_options...)
                ellip = interpolate(ellip_dom, halfdeg; interp_options...)
                @test length(ball.Ps) == 2
                @test length(ellip.Ps) == 2
                @test free.U == ball.U == ellip.U
                @test size(ball.pts) == size(ellip.pts)
            end
        end
    end
end

function test_complex_interp(T::Type{<:Real})
    @testset "complex interp n = $n, halfdeg = $halfdeg, use_qr = $use_qr" for
        (n, halfdeg) in ((1, 3), (2, 2), (3, 1)), use_qr in (false, true)
        gs = [z -> 1 - sum(abs2, z)]
        (pts, Ps) = interpolate(Complex{T}, halfdeg, n, gs, [1], use_qr = use_qr)
        @test length(Ps) == 2
        L = binomial(n + halfdeg, n)
        U = L^2
        @test length(pts) == U
        @test size(Ps[1]) == (U, L)
    end
end

function test_quadrature(T::Type{<:Real})
    # univariate
    @testset "quadrature n = 1, halfdeg = $halfdeg" for halfdeg in 1:4
        free_dom = PolyUtils.FreeDomain{T}(1)
        free = interpolate(free_dom, halfdeg, sample = false, get_quadr = true)
        pts_polys = [sum(free.pts[i, 1] ^ d for d in
            0:(2 * halfdeg)) for i in 1:free.U]
        @test dot(free.w, pts_polys) ≈ sum(2 / (2i + 1) for i in 0:halfdeg)
    end

    # multivariate box
    @testset "quadrature n = $n, halfdeg = $halfdeg" for (n, halfdeg) in
        ((1, 5), (2, 2), (3, 1))
        box_dom = PolyUtils.BoxDomain{T}(fill(-1, n), fill(3, n))
        box = interpolate(box_dom, halfdeg, sample = false, get_quadr = true)
        @test sum(box.w) ≈ 2^n
    end
end


text:
# Julia wrapper for header: /usr/local/opt/curl/include/curl/curl.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

function curl_strequal(s1, s2)
    ccall((:curl_strequal, libcurl), Cint, (Ptr{UInt8}, Ptr{UInt8}), s1, s2)
end

function curl_strnequal(s1, s2, n)
    ccall((:curl_strnequal, libcurl), Cint, (Ptr{UInt8}, Ptr{UInt8}, Csize_t), s1, s2, n)
end

function curl_mime_init(easy)
    ccall((:curl_mime_init, libcurl), Ptr{curl_mime}, (Ptr{CURL},), easy)
end

function curl_mime_free(mime)
    ccall((:curl_mime_free, libcurl), Cvoid, (Ptr{curl_mime},), mime)
end

function curl_mime_addpart(mime)
    ccall((:curl_mime_addpart, libcurl), Ptr{curl_mimepart}, (Ptr{curl_mime},), mime)
end

function curl_mime_name(part, name)
    ccall((:curl_mime_name, libcurl), CURLcode, (Ptr{curl_mimepart}, Ptr{UInt8}), part, name)
end

function curl_mime_filename(part, filename)
    ccall((:curl_mime_filename, libcurl), CURLcode, (Ptr{curl_mimepart}, Ptr{UInt8}), part, filename)
end

function curl_mime_type(part, mimetype)
    ccall((:curl_mime_type, libcurl), CURLcode, (Ptr{curl_mimepart}, Ptr{UInt8}), part, mimetype)
end

function curl_mime_encoder(part, encoding)
    ccall((:curl_mime_encoder, libcurl), CURLcode, (Ptr{curl_mimepart}, Ptr{UInt8}), part, encoding)
end

function curl_mime_data(part, data, datasize)
    ccall((:curl_mime_data, libcurl), CURLcode, (Ptr{curl_mimepart}, Ptr{UInt8}, Csize_t), part, data, datasize)
end

function curl_mime_filedata(part, filename)
    ccall((:curl_mime_filedata, libcurl), CURLcode, (Ptr{curl_mimepart}, Ptr{UInt8}), part, filename)
end

function curl_mime_data_cb(part, datasize, readfunc, seekfunc, freefunc, arg)
    ccall((:curl_mime_data_cb, libcurl), CURLcode, (Ptr{curl_mimepart}, curl_off_t, curl_read_callback, curl_seek_callback, curl_free_callback, Ptr{Cvoid}), part, datasize, readfunc, seekfunc, freefunc, arg)
end

function curl_mime_subparts(part, subparts)
    ccall((:curl_mime_subparts, libcurl), CURLcode, (Ptr{curl_mimepart}, Ptr{curl_mime}), part, subparts)
end

function curl_mime_headers(part, headers, take_ownership)
    ccall((:curl_mime_headers, libcurl), CURLcode, (Ptr{curl_mimepart}, Ptr{Cvoid}, Cint), part, headers, take_ownership)
end

function curl_formget(form, arg, append)
    ccall((:curl_formget, libcurl), Cint, (Ptr{Cvoid}, Ptr{Cvoid}, curl_formget_callback), form, arg, append)
end

function curl_formfree(form)
    ccall((:curl_formfree, libcurl), Cvoid, (Ptr{Cvoid},), form)
end

function curl_getenv(variable)
    ccall((:curl_getenv, libcurl), Ptr{UInt8}, (Ptr{UInt8},), variable)
end

function curl_version()
    ccall((:curl_version, libcurl), Ptr{UInt8}, ())
end

function curl_easy_escape(handle, string, length)
    ccall((:curl_easy_escape, libcurl), Ptr{UInt8}, (Ptr{CURL}, Ptr{UInt8}, Cint), handle, string, length)
end

function curl_escape(string, length)
    ccall((:curl_escape, libcurl), Ptr{UInt8}, (Ptr{UInt8}, Cint), string, length)
end

function curl_easy_unescape(handle, string, length, outlength)
    ccall((:curl_easy_unescape, libcurl), Ptr{UInt8}, (Ptr{CURL}, Ptr{UInt8}, Cint, Ptr{Cint}), handle, string, length, outlength)
end

function curl_unescape(string, length)
    ccall((:curl_unescape, libcurl), Ptr{UInt8}, (Ptr{UInt8}, Cint), string, length)
end

function curl_free(p)
    ccall((:curl_free, libcurl), Cvoid, (Ptr{Cvoid},), p)
end

function curl_global_init(flags)
    ccall((:curl_global_init, libcurl), CURLcode, (Clong,), flags)
end

function curl_global_init_mem(flags, m, f, r, s, c)
    ccall((:curl_global_init_mem, libcurl), CURLcode, (Clong, curl_malloc_callback, curl_free_callback, curl_realloc_callback, curl_strdup_callback, curl_calloc_callback), flags, m, f, r, s, c)
end

function curl_global_cleanup()
    ccall((:curl_global_cleanup, libcurl), Cvoid, ())
end

function curl_global_sslset(id, name, avail)
    ccall((:curl_global_sslset, libcurl), CURLsslset, (curl_sslbackend, Ptr{UInt8}, Ptr{Ptr{Ptr{curl_ssl_backend}}}), id, name, avail)
end

function curl_slist_append(arg1, arg2)
    ccall((:curl_slist_append, libcurl), Ptr{Cvoid}, (Ptr{Cvoid}, Ptr{UInt8}), arg1, arg2)
end

function curl_slist_free_all(arg1)
    ccall((:curl_slist_free_all, libcurl), Cvoid, (Ptr{Cvoid},), arg1)
end

function curl_getdate(p, unused)
    ccall((:curl_getdate, libcurl), time_t, (Ptr{UInt8}, Ptr{time_t}), p, unused)
end

function curl_share_init()
    ccall((:curl_share_init, libcurl), Ptr{CURLSH}, ())
end

function curl_share_setopt(handle, opt, param)
    ccall((:curl_share_setopt, libcurl), CURLSHcode, (Ptr{CURLSH}, CURLSHoption, Any...), handle, opt, param)
end

function curl_share_cleanup(arg1)
    ccall((:curl_share_cleanup, libcurl), CURLSHcode, (Ptr{CURLSH},), arg1)
end

function curl_version_info(arg1)
    ccall((:curl_version_info, libcurl), Ptr{curl_version_info_data}, (CURLversion,), arg1)
end

function curl_easy_strerror(arg1)
    ccall((:curl_easy_strerror, libcurl), Ptr{UInt8}, (CURLcode,), arg1)
end

function curl_share_strerror(arg1)
    ccall((:curl_share_strerror, libcurl), Ptr{UInt8}, (CURLSHcode,), arg1)
end

function curl_easy_pause(handle, bitmask)
    ccall((:curl_easy_pause, libcurl), CURLcode, (Ptr{CURL}, Cint), handle, bitmask)
end

function curl_easy_init()
    ccall((:curl_easy_init, libcurl), Ptr{CURL}, ())
end

function curl_easy_setopt(handle, opt, param)
    ccall((:curl_easy_setopt, libcurl), CURLcode, (Ptr{CURL}, CURLoption, Any...), handle, opt, param)
end

function curl_easy_perform(curl)
    ccall((:curl_easy_perform, libcurl), CURLcode, (Ptr{CURL},), curl)
end

function curl_easy_cleanup(curl)
    ccall((:curl_easy_cleanup, libcurl), Cvoid, (Ptr{CURL},), curl)
end

function curl_easy_getinfo(handle, info, arg)
    ccall((:curl_easy_getinfo, libcurl), CURLcode, (Ptr{CURL}, CURLINFO, Any...), handle, info, arg)
end

function curl_easy_duphandle(curl)
    ccall((:curl_easy_duphandle, libcurl), Ptr{CURL}, (Ptr{CURL},), curl)
end

function curl_easy_reset(curl)
    ccall((:curl_easy_reset, libcurl), Cvoid, (Ptr{CURL},), curl)
end

function curl_easy_recv(curl, buffer, buflen, n)
    ccall((:curl_easy_recv, libcurl), CURLcode, (Ptr{CURL}, Ptr{Cvoid}, Csize_t, Ptr{Csize_t}), curl, buffer, buflen, n)
end

function curl_easy_send(curl, buffer, buflen, n)
    ccall((:curl_easy_send, libcurl), CURLcode, (Ptr{CURL}, Ptr{Cvoid}, Csize_t, Ptr{Csize_t}), curl, buffer, buflen, n)
end

function curl_multi_init()
    ccall((:curl_multi_init, libcurl), Ptr{CURLM}, ())
end

function curl_multi_add_handle(multi_handle, curl_handle)
    ccall((:curl_multi_add_handle, libcurl), CURLMcode, (Ptr{CURLM}, Ptr{CURL}), multi_handle, curl_handle)
end

function curl_multi_remove_handle(multi_handle, curl_handle)
    ccall((:curl_multi_remove_handle, libcurl), CURLMcode, (Ptr{CURLM}, Ptr{CURL}), multi_handle, curl_handle)
end

function curl_multi_fdset(multi_handle, read_fd_set, write_fd_set, exc_fd_set, max_fd)
    ccall((:curl_multi_fdset, libcurl), CURLMcode, (Ptr{CURLM}, Ptr{fd_set}, Ptr{fd_set}, Ptr{fd_set}, Ptr{Cint}), multi_handle, read_fd_set, write_fd_set, exc_fd_set, max_fd)
end

function curl_multi_wait(multi_handle, extra_fds, extra_nfds, timeout_ms, ret)
    ccall((:curl_multi_wait, libcurl), CURLMcode, (Ptr{CURLM}, Ptr{Cvoid}, UInt32, Cint, Ptr{Cint}), multi_handle, extra_fds, extra_nfds, timeout_ms, ret)
end

function curl_multi_perform(multi_handle, running_handles)
    ccall((:curl_multi_perform, libcurl), CURLMcode, (Ptr{CURLM}, Ptr{Cint}), multi_handle, running_handles)
end

function curl_multi_cleanup(multi_handle)
    ccall((:curl_multi_cleanup, libcurl), CURLMcode, (Ptr{CURLM},), multi_handle)
end

function curl_multi_info_read(multi_handle, msgs_in_queue)
    ccall((:curl_multi_info_read, libcurl), Ptr{CURLMsg}, (Ptr{CURLM}, Ptr{Cint}), multi_handle, msgs_in_queue)
end

function curl_multi_strerror(arg1)
    ccall((:curl_multi_strerror, libcurl), Ptr{UInt8}, (CURLMcode,), arg1)
end

function curl_multi_socket(multi_handle, s, running_handles)
    ccall((:curl_multi_socket, libcurl), CURLMcode, (Ptr{CURLM}, curl_socket_t, Ptr{Cint}), multi_handle, s, running_handles)
end

function curl_multi_socket_action(multi_handle, s, ev_bitmask, running_handles)
    ccall((:curl_multi_socket_action, libcurl), CURLMcode, (Ptr{CURLM}, curl_socket_t, Cint, Ptr{Cint}), multi_handle, s, ev_bitmask, running_handles)
end

function curl_multi_socket_all(multi_handle, running_handles)
    ccall((:curl_multi_socket_all, libcurl), CURLMcode, (Ptr{CURLM}, Ptr{Cint}), multi_handle, running_handles)
end

function curl_multi_timeout(multi_handle, milliseconds)
    ccall((:curl_multi_timeout, libcurl), CURLMcode, (Ptr{CURLM}, Ptr{Clong}), multi_handle, milliseconds)
end

function curl_multi_setopt(multi_handle, opt, param)
    ccall((:curl_multi_setopt, libcurl), CURLMcode, (Ptr{CURLM}, CURLMoption, Any...), multi_handle, opt, param)
end

function curl_multi_assign(multi_handle, sockfd, sockp)
    ccall((:curl_multi_assign, libcurl), CURLMcode, (Ptr{CURLM}, curl_socket_t, Ptr{Cvoid}), multi_handle, sockfd, sockp)
end

function curl_pushheader_bynum(h, num)
    ccall((:curl_pushheader_bynum, libcurl), Ptr{UInt8}, (Ptr{Cvoid}, Csize_t), h, num)
end

function curl_pushheader_byname(h, name)
    ccall((:curl_pushheader_byname, libcurl), Ptr{UInt8}, (Ptr{Cvoid}, Ptr{UInt8}), h, name)
end