text:
import MitosisStochasticDiffEq as MSDE
using Mitosis
using Test, Random
using LinearAlgebra

# set true/estimated model parameters
# p = [-0.1,0.2,0.9]

"""
    backwardfilter() -> ps, p0, c
Backward filtering using the Euler method, starting with `N(ν, P)` prior
and integration scale c2 between observations

"""
function backwardfilter((c, ν, P)::NamedTuple{(:logscale, :μ, :Σ)}, p, s)
    ps = [[ν, P, c]]
    B, β, σ̃ = p
    for i in eachindex(s)[end-1:-1:1]
        dt = s[i+1] - s[i]
        H = inv(P)
        F = H*ν

        P = P - dt*(B*P + P*B' - σ̃*σ̃'*I)
        ν = ν - dt*(B*ν .+ β)
        c = c - dt*tr(B)

        push!(ps, [ν, P, c])
    end
    ps, [ν, P, c]
end

@testset "backward filtering tests" begin

  # define SDE function
  f(u,p,t) = p[1]*u .+ p[2]
  g(u,p,t) = p[3]

  # time span
  tstart = 0.0
  tend = 1.0
  dt = 0.02
  trange = tstart:dt:tend

  # set of linear parameters Eq.~(2.2)
  B, β, σ̃ = -0.1, 0.2, 1.3
  plin = [B, β, σ̃]
  kernel = MSDE.SDEKernel(f,g,trange,plin)

  # initial values for ODE
  mynames = (:logscale, :μ, :Σ);
  myvalues = [0.0, 0.0, 10.0];
  NT = NamedTuple{mynames}(myvalues)

  message, solend = MSDE.backwardfilter(kernel, NT)

  kernel2 = MSDE.SDEKernel(Mitosis.AffineMap(B, β), Mitosis.ConstantMap(σ̃), trange, plin)
  message2, solend2 = MSDE.backwardfilter(kernel2, NT)

  @test isapprox(solend, solend2, rtol=1e-15)
  @test isapprox(Array(message.sol), Array(message2.sol), rtol=1e-15)

  message2, solend2 = backwardfilter(NT, plin, message.ts)

  @test isapprox(solend, solend2, rtol=1e-15)
  @test isapprox(Array(message.sol), reduce(hcat, message2), rtol=1e-15)

  # multivariate tests
  dim = 5
  Random.seed!(123)
  logscale = randn()
  μ = randn(dim)
  Σ = randn(dim,dim)
  myvalues = [logscale, μ, Σ];
  NT = NamedTuple{mynames}(myvalues)

  kernel = MSDE.SDEKernel(f,g,trange,plin)
  message, solend = MSDE.backwardfilter(kernel, NT)
  message2, solend2 = backwardfilter(NT, plin, message.ts)

  @test isapprox(solend.x[1], solend2[1], rtol=1e-15)
  @test isapprox(solend.x[2], solend2[2], rtol=1e-15)
  @test isapprox(solend.x[3][1], solend2[3], rtol=1e-15)

  kernel2 = MSDE.SDEKernel(Mitosis.AffineMap(B, β), Mitosis.ConstantMap(σ̃), trange, plin)
  message2, solend2 = MSDE.backwardfilter(kernel2, NT)
  @test isapprox(solend, solend2, rtol=1e-15)
  @test isapprox(Array(message.sol), Array(message2.sol), rtol=1e-15)

  # test inplace version
  message2, solend2 = MSDE.backwardfilter(kernel, NT, inplace=true)
  @test isapprox(solend, solend2, rtol=1e-15)
  @test isapprox(Array(message.sol), Array(message2.sol), rtol=1e-15)

  message2, solend2 = MSDE.backwardfilter(kernel2, NT, inplace=true)
  @test isapprox(solend, solend2, rtol=1e-15)
  @test isapprox(Array(message.sol), Array(message2.sol), rtol=1e-15)

  m = 3 # some number of Brownian processes
  plin = [randn(dim,dim), randn(dim), randn(dim,m)] # B, β, σtil

  kernel = MSDE.SDEKernel(f,g,trange,plin)
  message, solend = MSDE.backwardfilter(kernel, NT)
  message2, solend2 = backwardfilter(NT, plin, message.ts)

  @test isapprox(solend.x[1], solend2[1], rtol=1e-14)
  @test isapprox(solend.x[2], solend2[2], rtol=1e-14)
  @test isapprox(solend.x[3][1], solend2[3], rtol=1e-14)

  kernel2 = MSDE.SDEKernel(Mitosis.AffineMap(B, β), Mitosis.ConstantMap(σ̃), trange, plin)
  message2, solend2 = MSDE.backwardfilter(kernel2, NT)
  @test isapprox(solend, solend2, rtol=1e-15)
  @test isapprox(Array(message.sol), Array(message2.sol), rtol=1e-15)

  # test inplace version
  message2, solend2 = MSDE.backwardfilter(kernel, NT, inplace=true)
  @test isapprox(solend, solend2, rtol=1e-15)
  @test isapprox(Array(message.sol), Array(message2.sol), rtol=1e-15)

  # test symmetric matrix
  plin = [Symmetric(randn(dim,dim)), randn(dim), randn(dim,m)] # B, β, σtil
  kernel = MSDE.SDEKernel(f,g,trange,plin)
  kernel2 = MSDE.SDEKernel(Mitosis.AffineMap(B, β), Mitosis.ConstantMap(σ̃), trange, plin)
  message, solend  = MSDE.backwardfilter(kernel, NT)
  message2, solend2 = MSDE.backwardfilter(kernel2, NT)
  @test isapprox(solend, solend2, rtol=1e-15)
  @test isapprox(Array(message.sol), Array(message2.sol), rtol=1e-15)

  plin = [Array(plin[1]), plin[2], plin[3]] # B, β, σtil
  kernel = MSDE.SDEKernel(f,g,trange,plin)
  message2, solend2  = MSDE.backwardfilter(kernel, NT)

  @test isapprox(solend, solend2, rtol=1e-15)
  @test isapprox(Array(message.sol), Array(message2.sol), rtol=1e-15)
end

@testset "backward filtering timechange tests" begin
  # define SDE function
  f(u,p,t) = p[1]*u .+ p[2]
  g(u,p,t) = p[3]

  # time span
  tstart = 0.0
  tend = 1.0
  dt = 0.02
  trange = tstart:dt:tend

  # set of linear parameters Eq.~(2.2)
  B, β, σ̃ = -0.1, 0.2, 1.3
  plin = [B, β, σ̃]

  # initial values for ODE
  mynames = (:logscale, :μ, :Σ);
  myvalues = [0.0, 0.0, 10.0];
  NT = NamedTuple{mynames}(myvalues)

  kernel = MSDE.SDEKernel(Mitosis.AffineMap(B, β), Mitosis.ConstantMap(σ̃), trange, plin)
  message, solend = MSDE.backwardfilter(kernel, NT, apply_timechange=true)

  @test length(message.ts) == length(trange)
  @test message.ts != trange
  @test message.ts == MSDE.timechange(trange)

  message2, solend2 = backwardfilter(NT, plin, message.ts)

  @test isapprox(solend, solend2, rtol=1e-15)
  @test isapprox(Array(message.sol), reduce(hcat, message2), rtol=1e-15)

end

@testset "backward filtering adaptive tests" begin
  # define SDE function
  using OrdinaryDiffEq

  f(u,p,t) = p[1]*u .+ p[2]
  g(u,p,t) = p[3]

  # time span
  tstart = 0.0
  tend = 1.0
  dt = 0.05 # more coarse grained dt
  trange = tstart:dt:tend

  # set of linear parameters Eq.~(2.2)
  B, β, σ̃ = -0.1, 0.2, 1.3
  plin = [B, β, σ̃]

  # initial values for ODE
  mynames = (:logscale, :μ, :Σ);
  myvalues = [0.0, 0.0, 10.0];
  NT = NamedTuple{mynames}(myvalues)

  tildekernel = MSDE.SDEKernel(Mitosis.AffineMap(B, β), Mitosis.ConstantMap(σ̃), trange, plin)
  message, solend = MSDE.backwardfilter(tildekernel, NT, alg=Tsit5(), apply_timechange=true)

  # adaptive solution with large tolerances follows tstops
  @test length(message.ts) == length(trange)
  @test message.ts != trange
  @test message.ts == MSDE.timechange(trange)

  message2, solend2 = MSDE.backwardfilter(tildekernel, NT, alg=Tsit5(),
    abstol=1e-12, reltol=1e-12, apply_timechange=true)

  # adaptive solution with small tolerances steps to tstops but also takes additional substeps
  @test length(message.ts) != length(message2.ts)
  @test minimum(message.ts .∈ [message2.ts])
end


text:
#function bodies
@createOpenGLFun glDebugMessageEnableAMD(category::GLenum, severity::GLenum, count::GLsizei, ids::Ptr{GLuint}, enabled::GLboolean)::Void
@createOpenGLFun glGetPerfMonitorGroupStringAMD(group::GLuint, bufSize::GLsizei, length::Ptr{GLsizei}, groupString::Ptr{GLchar})::Void
@createOpenGLFun glGetPerfMonitorCounterInfoAMD(group::GLuint, counter::GLuint, pname::GLenum, data::Ptr{Void})::Void
@createOpenGLFun glBlendEquationSeparateIndexedAMD(buf::GLuint, modeRGB::GLenum, modeAlpha::GLenum)::Void
@createOpenGLFun glIsNameAMD(identifier::GLenum, name::GLuint)::Bool
@createOpenGLFun glGenNamesAMD(identifier::GLenum, num::GLuint, names::Ptr{GLuint})::Void
@createOpenGLFun glGetPerfMonitorCounterDataAMD(monitor::GLuint, pname::GLenum, dataSize::GLsizei, data::Ptr{GLuint}, bytesWritten::Ptr{GLint})::Void
@createOpenGLFun glBlendFuncSeparateIndexedAMD(buf::GLuint, srcRGB::GLenum, dstRGB::GLenum, srcAlpha::GLenum, dstAlpha::GLenum)::Void
@createOpenGLFun glGenPerfMonitorsAMD(n::GLsizei, monitors::Ptr{GLuint})::Void
@createOpenGLFun glSelectPerfMonitorCountersAMD(monitor::GLuint, enable::GLboolean, group::GLuint, numCounters::GLint, counterList::Ptr{GLuint})::Void
@createOpenGLFun glEndPerfMonitorAMD(monitor::GLuint)::Void
@createOpenGLFun glSetMultisamplefvAMD(pname::GLenum, index::GLuint, val::Ptr{GLfloat})::Void
@createOpenGLFun glBeginPerfMonitorAMD(monitor::GLuint)::Void
@createOpenGLFun glGetPerfMonitorCountersAMD(group::GLuint, numCounters::Ptr{GLint}, maxActiveCounters::Ptr{GLint}, counterSize::GLsizei, counters::Ptr{GLuint})::Void
@createOpenGLFun glGetPerfMonitorCounterStringAMD(group::GLuint, counter::GLuint, bufSize::GLsizei, length::Ptr{GLsizei}, counterString::Ptr{GLchar})::Void
@createOpenGLFun glTexStorageSparseAMD(target::GLenum, internalFormat::GLenum, width::GLsizei, height::GLsizei, depth::GLsizei, layers::GLsizei, flags::GLbitfield)::Void
@createOpenGLFun glMultiDrawArraysIndirectAMD(mode::GLenum, indirect::Ptr{Void}, primcount::GLsizei, stride::GLsizei)::Void
@createOpenGLFun glTessellationFactorAMD(factor::GLfloat)::Void
@createOpenGLFun glBlendEquationIndexedAMD(buf::GLuint, mode::GLenum)::Void
@createOpenGLFun glGetDebugMessageLogAMD(count::GLuint, bufsize::GLsizei, categories::Ptr{GLenum}, severities::Ptr{GLuint}, ids::Ptr{GLuint}, lengths::Ptr{GLsizei}, message::Ptr{GLchar})::Cuint
@createOpenGLFun glTextureStorageSparseAMD(texture::GLuint, target::GLenum, internalFormat::GLenum, width::GLsizei, height::GLsizei, depth::GLsizei, layers::GLsizei, flags::GLbitfield)::Void
@createOpenGLFun glGetPerfMonitorGroupsAMD(numGroups::Ptr{GLint}, groupsSize::GLsizei, groups::Ptr{GLuint})::Void
@createOpenGLFun glDeleteNamesAMD(identifier::GLenum, num::GLuint, names::Ptr{GLuint})::Void
@createOpenGLFun glDeletePerfMonitorsAMD(n::GLsizei, monitors::Ptr{GLuint})::Void
@createOpenGLFun glDebugMessageInsertAMD(category::GLenum, severity::GLenum, id::GLuint, length::GLsizei, buf::Ptr{GLchar})::Void
@createOpenGLFun glDebugMessageCallbackAMD(callback::Ptr{Void}, userParam::Ptr{Void})::Void
@createOpenGLFun glStencilOpValueAMD(face::GLenum, value::GLuint)::Void
@createOpenGLFun glTessellationModeAMD(mode::GLenum)::Void
@createOpenGLFun glBlendFuncIndexedAMD(buf::GLuint, src::GLenum, dst::GLenum)::Void
@createOpenGLFun glMultiDrawElementsIndirectAMD(mode::GLenum, type_::GLenum, indirect::Ptr{Void}, primcount::GLsizei, stride::GLsizei)::Void
#constants
const GL_VIRTUAL_PAGE_SIZE_X_AMD                                      = 0x9195
const GL_PERFMON_RESULT_SIZE_AMD                                      = 0x8BC5
const GL_TEXTURE_STORAGE_SPARSE_BIT_AMD                               = 0x00000001
const GL_MAX_DEBUG_LOGGED_MESSAGES_AMD                                = 0x9144
const GL_PERFORMANCE_MONITOR_AMD                                      = 0x9152
const GL_DEBUG_CATEGORY_OTHER_AMD                                     = 0x9150
const GL_PERFMON_RESULT_AMD                                           = 0x8BC6
const GL_TESSELLATION_FACTOR_AMD                                      = 0x9005
const GL_MIN_LOD_WARNING_AMD                                          = 0x919C
const GL_SUBSAMPLE_DISTANCE_AMD                                       = 0x883F
const GL_CONTINUOUS_AMD                                               = 0x9007
const GL_QUERY_BUFFER_BINDING_AMD                                     = 0x9193
const GL_COUNTER_TYPE_AMD                                             = 0x8BC0
const GL_MAX_DEBUG_MESSAGE_LENGTH_AMD                                 = 0x9143
const GL_VIRTUAL_PAGE_SIZE_Y_AMD                                      = 0x9196
const GL_COUNTER_RANGE_AMD                                            = 0x8BC1
const GL_DEPTH_CLAMP_FAR_AMD                                          = 0x901F
const GL_FACTOR_MIN_AMD                                               = 0x901C
const GL_DEBUG_CATEGORY_UNDEFINED_BEHAVIOR_AMD                        = 0x914C
const GL_DEBUG_CATEGORY_WINDOW_SYSTEM_AMD                             = 0x914A
const GL_DEPTH_CLAMP_NEAR_AMD                                         = 0x901E
const GL_DEBUG_SEVERITY_MEDIUM_AMD                                    = 0x9147
const GL_SET_AMD                                                      = 0x874A
const GL_TEXTURE_CUBE_MAP_SEAMLESS                                    = 0x884F
const GL_DEBUG_LOGGED_MESSAGES_AMD                                    = 0x9145
const GL_DEBUG_SEVERITY_LOW_AMD                                       = 0x9148
const GL_DEBUG_CATEGORY_API_ERROR_AMD                                 = 0x9149
const GL_TESSELLATION_MODE_AMD                                        = 0x9004
const GL_DEBUG_CATEGORY_APPLICATION_AMD                               = 0x914F
const GL_MAX_SPARSE_ARRAY_TEXTURE_LAYERS                              = 0x919A
const GL_DEBUG_SEVERITY_HIGH_AMD                                      = 0x9146
const GL_MAX_SPARSE_TEXTURE_SIZE_AMD                                  = 0x9198
const GL_SAMPLER_BUFFER_AMD                                           = 0x9001
const GL_PERFMON_RESULT_AVAILABLE_AMD                                 = 0x8BC4
const GL_MIN_SPARSE_LEVEL_AMD                                         = 0x919B
const GL_INT_SAMPLER_BUFFER_AMD                                       = 0x9002
const GL_UNSIGNED_INT_SAMPLER_BUFFER_AMD                              = 0x9003
const GL_UNSIGNED_INT64_AMD                                           = 0x8BC2
const GL_VERTEX_ARRAY_OBJECT_AMD                                      = 0x9154
const GL_FACTOR_MAX_AMD                                               = 0x901D
const GL_DEBUG_CATEGORY_DEPRECATION_AMD                               = 0x914B
const GL_STENCIL_OP_VALUE_AMD                                         = 0x874C
const GL_MAX_SPARSE_3D_TEXTURE_SIZE_AMD                               = 0x9199
const GL_DEBUG_CATEGORY_SHADER_COMPILER_AMD                           = 0x914E
const GL_DATA_BUFFER_AMD                                              = 0x9151
const GL_VIRTUAL_PAGE_SIZE_Z_AMD                                      = 0x9197
const GL_PERCENTAGE_AMD                                               = 0x8BC3
const GL_SAMPLER_OBJECT_AMD                                           = 0x9155
const GL_STENCIL_BACK_OP_VALUE_AMD                                    = 0x874D
const GL_DISCRETE_AMD                                                 = 0x9006
const GL_QUERY_OBJECT_AMD                                             = 0x9153
const GL_QUERY_BUFFER_AMD                                             = 0x9192
const GL_QUERY_RESULT_NO_WAIT_AMD                                     = 0x9194
const GL_DEBUG_CATEGORY_PERFORMANCE_AMD                               = 0x914D
const GL_REPLACE_VALUE_AMD                                            = 0x874B
const GL_EXTERNAL_VIRTUAL_MEMORY_BUFFER_AMD                           = 0x9160
# Export everything!
export glDebugMessageEnableAMD, glGetPerfMonitorGroupStringAMD, glGetPerfMonitorCounterInfoAMD, glBlendEquationSeparateIndexedAMD, glIsNameAMD, glGenNamesAMD, glGetPerfMonitorCounterDataAMD, glBlendFuncSeparateIndexedAMD, glGenPerfMonitorsAMD, glSelectPerfMonitorCountersAMD, glEndPerfMonitorAMD, glSetMultisamplefvAMD, glBeginPerfMonitorAMD, glGetPerfMonitorCountersAMD, glGetPerfMonitorCounterStringAMD, glTexStorageSparseAMD, glMultiDrawArraysIndirectAMD, glTessellationFactorAMD, glBlendEquationIndexedAMD, glGetDebugMessageLogAMD, glTextureStorageSparseAMD, glGetPerfMonitorGroupsAMD, glDeleteNamesAMD, glDeletePerfMonitorsAMD, glDebugMessageInsertAMD, glDebugMessageCallbackAMD, glStencilOpValueAMD, glTessellationModeAMD, glBlendFuncIndexedAMD, glMultiDrawElementsIndirectAMD, GL_VIRTUAL_PAGE_SIZE_X_AMD, GL_PERFMON_RESULT_SIZE_AMD, GL_TEXTURE_STORAGE_SPARSE_BIT_AMD, GL_MAX_DEBUG_LOGGED_MESSAGES_AMD, GL_PERFORMANCE_MONITOR_AMD, GL_DEBUG_CATEGORY_OTHER_AMD, GL_PERFMON_RESULT_AMD, GL_TESSELLATION_FACTOR_AMD, GL_MIN_LOD_WARNING_AMD, GL_SUBSAMPLE_DISTANCE_AMD, GL_CONTINUOUS_AMD, GL_QUERY_BUFFER_BINDING_AMD, GL_COUNTER_TYPE_AMD, GL_MAX_DEBUG_MESSAGE_LENGTH_AMD, GL_VIRTUAL_PAGE_SIZE_Y_AMD, GL_COUNTER_RANGE_AMD, GL_DEPTH_CLAMP_FAR_AMD, GL_FACTOR_MIN_AMD, GL_DEBUG_CATEGORY_UNDEFINED_BEHAVIOR_AMD, GL_DEBUG_CATEGORY_WINDOW_SYSTEM_AMD, GL_DEPTH_CLAMP_NEAR_AMD, GL_DEBUG_SEVERITY_MEDIUM_AMD, GL_SET_AMD, GL_TEXTURE_CUBE_MAP_SEAMLESS, GL_DEBUG_LOGGED_MESSAGES_AMD, GL_DEBUG_SEVERITY_LOW_AMD, GL_DEBUG_CATEGORY_API_ERROR_AMD, GL_TESSELLATION_MODE_AMD, GL_DEBUG_CATEGORY_APPLICATION_AMD, GL_MAX_SPARSE_ARRAY_TEXTURE_LAYERS, GL_DEBUG_SEVERITY_HIGH_AMD, GL_MAX_SPARSE_TEXTURE_SIZE_AMD, GL_SAMPLER_BUFFER_AMD, GL_PERFMON_RESULT_AVAILABLE_AMD, GL_MIN_SPARSE_LEVEL_AMD, GL_INT_SAMPLER_BUFFER_AMD, GL_UNSIGNED_INT_SAMPLER_BUFFER_AMD, GL_UNSIGNED_INT64_AMD, GL_VERTEX_ARRAY_OBJECT_AMD, GL_FACTOR_MAX_AMD, GL_DEBUG_CATEGORY_DEPRECATION_AMD, GL_STENCIL_OP_VALUE_AMD, GL_MAX_SPARSE_3D_TEXTURE_SIZE_AMD, GL_DEBUG_CATEGORY_SHADER_COMPILER_AMD, GL_DATA_BUFFER_AMD, GL_VIRTUAL_PAGE_SIZE_Z_AMD, GL_PERCENTAGE_AMD, GL_SAMPLER_OBJECT_AMD, GL_STENCIL_BACK_OP_VALUE_AMD, GL_DISCRETE_AMD, GL_QUERY_OBJECT_AMD, GL_QUERY_BUFFER_AMD, GL_QUERY_RESULT_NO_WAIT_AMD, GL_DEBUG_CATEGORY_PERFORMANCE_AMD, GL_REPLACE_VALUE_AMD, GL_EXTERNAL_VIRTUAL_MEMORY_BUFFER_AMD

ext: jl


text:
module TaskGraphs

using Parameters
using LightGraphs
using GraphUtils
using Random
using LinearAlgebra
using DataStructures
using JuMP, MathOptInterface
using GLPK
using TOML
using CRCBS
using SparseArrays
using Reexport

export
    DEBUG_PATH,
    EXPERIMENT_DIR,
    ENVIRONMENT_DIR,
    PROBLEM_DIR,
    RESULTS_DIR,
    VIDEO_DIR

global DEBUG_PATH       = joinpath(dirname(pathof(TaskGraphs)),"..","debug")
global EXPERIMENT_DIR   = joinpath(dirname(pathof(TaskGraphs)),"..","experiments")
global ENVIRONMENT_DIR  = joinpath(EXPERIMENT_DIR,"environments")
global PROBLEM_DIR      = joinpath(EXPERIMENT_DIR,"problem_instances")
global RESULTS_DIR      = joinpath(EXPERIMENT_DIR,"results")
global VIDEO_DIR        = joinpath(EXPERIMENT_DIR,"videos")

include("JuMP_interface.jl")
include("planning_predicates.jl")
include("task_graphs_core.jl")
include("utils.jl")
include("task_assignment_solvers.jl")
include("path_planning.jl")
include("pc_tapf_solvers.jl")
include("milp_formulation.jl")
include("disturbances.jl")
include("path_planners/pibt_planner.jl")
include("replanning.jl")
include("helpers/problem_instances.jl")
include("helpers/profiling.jl")
include("experiments/replanning_experiments.jl")

@reexport using CRCBS

# set GLPK as the default optimizer
set_default_milp_optimizer!(GLPK.Optimizer)
set_default_optimizer_attributes!(MOI.TimeLimitSec()=>100, MOI.Silent()=>true)

end


text:
using Documenter, FrankWolfe
using SparseArrays
using LinearAlgebra

using Literate, Test

const _EXAMPLE_DIR = joinpath(@__DIR__, "src", "examples")

"""
    _include_sandbox(filename)
Include the `filename` in a temporary module that acts as a sandbox. (Ensuring
no constants or functions leak into other files.)
"""
function _include_sandbox(filename)
    mod = @eval module $(gensym()) end
    return Base.include(mod, filename)
end

function _file_list(full_dir, relative_dir, extension)
    return map(
        file -> joinpath(relative_dir, file),
        filter(file -> endswith(file, extension), sort(readdir(full_dir))),
    )
end

function literate_directory(dir)
    rm.(_file_list(dir, dir, ".md"))
    for filename in _file_list(dir, dir, ".jl")
        # `include` the file to test it before `#src` lines are removed. It is
        # in a testset to isolate local variables between files.
        @testset "$(filename)" begin
            _include_sandbox(filename)
        end
        Literate.markdown(
            filename,
            dir;
            documenter = true,
        )
    end
    return
end

literate_directory(_EXAMPLE_DIR)

ENV["GKSwstype"] = "100"

generated_path = joinpath(@__DIR__, "src")
base_url = "https://github.com/ZIB-IOL/FrankWolfe.jl/blob/master/"
isdir(generated_path) || mkdir(generated_path)

open(joinpath(generated_path, "contributing.md"), "w") do io
    # Point to source license file
    println(
        io,
        """
        ```@meta
        EditURL = "$(base_url)CONTRIBUTING.md"
        ```
        """,
    )
    # Write the contents out below the meta block
    for line in eachline(joinpath(dirname(@__DIR__), "CONTRIBUTING.md"))
        println(io, line)
    end
end

makedocs(
    modules=[FrankWolfe],
    sitename="FrankWolfe.jl",
    format=Documenter.HTML(prettyurls=get(ENV, "CI", nothing) == "true"),
    pages=[
        "Home" => "index.md",
        "Examples" => [
        joinpath("examples", f) for f in readdir(_EXAMPLE_DIR) if endswith(f, ".md")
        ],
        "Contributing" => "contributing.md",
        "References" => "reference.md",
        "Index" => "indexlist.md",
    ],
)

deploydocs(repo="github.com/ZIB-IOL/FrankWolfe.jl.git", push_preview=true)


text:
struct GradientCache{CacheType1, CacheType2, CacheType3, fdtype, returntype, inplace}
    fx :: CacheType1
    c1 :: CacheType2
    c2 :: CacheType3
end

function GradientCache(
    df         :: Union{<:Number, AbstractArray{<:Number}},
    x          :: Union{<:Number, AbstractArray{<:Number}},
    fdtype     :: Type{T1} = Val{:central},
    returntype :: Type{T2} = eltype(df),
    inplace    :: Type{Val{T3}} = Val{true}) where {T1,T2,T3}

    if typeof(x)<:AbstractArray # the vector->scalar case
        if fdtype!=Val{:complex} # complex-mode FD only needs one cache, for x+eps*im
            if typeof(x)<:StridedVector
                if eltype(df)<:Complex && !(eltype(x)<:Complex)
                    _c1 = fill(zero(Complex{eltype(x)}), size(x))
                    _c2 = nothing
                else
                    _c1 = nothing
                    _c2 = nothing
                end
            else
                _c1 = similar(x)
                _c2 = fill(zero(real(eltype(x))), size(x))
            end
        else
            if !(returntype<:Real)
                fdtype_error(returntype)
            else
                _c1 = x .+ 0*im
                _c2 = nothing
            end
        end
    else # the scalar->vector case
        # need cache arrays for fx1 and fx2, except in complex mode, which needs one complex array
        if fdtype != Val{:complex}
            _c1 = similar(df)
            _c2 = similar(df)
        else
            _c1 = fill(zero(Complex{eltype(x)}), size(df))
            _c2 = nothing
        end
    end

    GradientCache{Nothing,typeof(_c1),typeof(_c2),fdtype,
                  returntype,inplace}(nothing,_c1,_c2)

end

function GradientCache(
    c1         :: Union{Nothing,AbstractArray{<:Number}},
    c2         :: Union{Nothing,AbstractArray{<:Number}},
    fx         :: Union{Nothing,<:Number,AbstractArray{<:Number}} = nothing,
    fdtype     :: Type{T1} = Val{:central},
    returntype :: Type{T2} = eltype(c1),
    inplace    :: Type{Val{T3}} = Val{true}) where {T1,T2,T3}

    if fdtype!=Val{:forward} && typeof(fx)!=Nothing
        @warn("Pre-computed function values are only useful for fdtype == Val{:forward}.")
        _fx = nothing
    else
        # more runtime sanity checks?
        _fx = fx
    end

    if typeof(x)<:AbstractArray # the vector->scalar case
        # need cache arrays for x1 (c1) and epsilon (c2) (both only if non-StridedArray)
        if fdtype!=Val{:complex} # complex-mode FD only needs one cache, for x+eps*im
            if typeof(x)<:StridedVector
                if eltype(df)<:Complex && !(eltype(x)<:Complex)
                    _c1 = fill(zero(Complex{eltype(x)}), size(x))
                    _c2 = nothing
                else
                    _c1 = nothing
                    _c2 = nothing
                    if typeof(c1)!=Nothing || typeof(c2)!=Nothing
                        @warn("For StridedVectors, neither c1 nor c2 are necessary.")
                    end
                end
            else
                _c1 = c1
                _c2 = c2
            end
        else
            if !(returntype<:Real)
                fdtype_error(returntype)
            else
                _c1 = x + 0*im
                _c2 = nothing
            end
        end

    else # the scalar->vector case
        # need cache arrays for fx1 and fx2, except in complex mode, which needs one complex array
        if fdtype != Val{:complex}
            _c1 = c1
            _c2 = c2
        else
            _c1 = c1
            _c2 = nothing
        end
    end
    GradientCache{typeof(_fx),typeof(_c1),typeof(_c2),fdtype,returntype,inplace}(_fx,_c1,_c2)
end

function finite_difference_gradient(f, x, fdtype::Type{T1}=Val{:central},
    returntype::Type{T2}=eltype(x), inplace::Type{Val{T3}}=Val{true},
    fx::Union{Nothing,AbstractArray{<:Number}}=nothing,
    c1::Union{Nothing,AbstractArray{<:Number}}=nothing,
    c2::Union{Nothing,AbstractArray{<:Number}}=nothing) where {T1,T2,T3}

    if typeof(x) <: AbstractArray
        df = fill(zero(returntype), size(x))
    else
        if inplace == Val{true}
            if typeof(fx)==Nothing && typeof(c1)==Nothing && typeof(c2)==Nothing
                error("In the scalar->vector in-place map case, at least one of fx, c1 or c2 must be provided, otherwise we cannot infer the return size.")
            else
                if     c1 != nothing    df = similar(c1)
                elseif fx != nothing    df = similar(fx)
                elseif c2 != nothing    df = similar(c2)
                end
            end
        else
            df = similar(f(x))
        end
    end
    cache = GradientCache(df,x,fdtype,returntype,inplace)
    finite_difference_gradient!(df,f,x,cache)
end

function finite_difference_gradient!(df, f, x, fdtype::Type{T1}=Val{:central},
    returntype::Type{T2}=eltype(df), inplace::Type{Val{T3}}=Val{true},
    fx::Union{Nothing,AbstractArray{<:Number}}=nothing,
    c1::Union{Nothing,AbstractArray{<:Number}}=nothing,
    c2::Union{Nothing,AbstractArray{<:Number}}=nothing,
    ) where {T1,T2,T3}

    cache = GradientCache(df,x,fdtype,returntype,inplace)
    finite_difference_gradient!(df,f,x,cache)
end

function finite_difference_gradient(f,x,
    cache::GradientCache{T1,T2,T3,fdtype,returntype,inplace}) where {T1,T2,T3,fdtype,returntype,inplace}

    if typeof(x) <: AbstractArray
        df = fill(zero(returntype), size(x))
    else
        df = zero(cache.c1)
    end
    finite_difference_gradient!(df,f,x,cache)
    df
end

# vector of derivatives of a vector->scalar map by each component of a vector x
# this ignores the value of "inplace", because it doesn't make much sense
function finite_difference_gradient!(df::AbstractArray{<:Number}, f, x::AbstractArray{<:Number},
    cache::GradientCache{T1,T2,T3,fdtype,returntype,inplace}) where {T1,T2,T3,fdtype,returntype,inplace}

    # NOTE: in this case epsilon is a vector, we need two arrays for epsilon and x1
    # c1 denotes x1, c2 is epsilon
    fx, c1, c2 = cache.fx, cache.c1, cache.c2
    if fdtype != Val{:complex}
        epsilon_factor = compute_epsilon_factor(fdtype, eltype(x))
        @. c2 = compute_epsilon(fdtype, x, epsilon_factor)
        copyto!(c1,x)
    end
    if fdtype == Val{:forward}
        @inbounds for i ∈ eachindex(x)
            epsilon = c2[i]
            c1_old = c1[i]
            c1[i] += epsilon
            if typeof(fx) != Nothing
                dfi = (f(c1) - fx) / epsilon
            else
                fx0 = f(x)
                dfi = (f(c1) - fx0) / epsilon
            end
            df[i] = real(dfi)
            c1[i] = c1_old
            if eltype(df)<:Complex
                c1[i] += im * epsilon
                if typeof(fx) != Nothing
                    dfi = (f(c1) - fx) / (im*epsilon)
                else
                    dfi = (f(c1) - fx0) / (im*epsilon)
                end
                c1[i] = c1_old
                df[i] -= im * imag(dfi)
            end
        end
    elseif fdtype == Val{:central}
        @inbounds for i ∈ eachindex(x)
            epsilon = c2[i]
            c1_old = c1[i]
            c1[i] += epsilon
            x_old  = x[i]
            x[i]  -= epsilon
            df[i]  = real((f(c1) - f(x)) / (2*epsilon))
            c1[i]  = c1_old
            x[i]   = x_old
            if eltype(df)<:Complex
                c1[i] += im*epsilon
                x[i]  -= im*epsilon
                df[i] -= im*imag( (f(c1) - f(x)) / (2*im*epsilon) )
                c1[i] = c1_old
                x[i] = x_old
            end
        end
    elseif fdtype == Val{:complex} && returntype <: Real
        copyto!(c1,x)
        epsilon_complex = eps(real(eltype(x)))
        # we use c1 here to avoid typing issues with x
        @inbounds for i ∈ eachindex(x)
            c1_old = c1[i]
            c1[i] += im*epsilon_complex
            df[i]  = imag(f(c1)) / epsilon_complex
            c1[i]  = c1_old
        end
    else
        fdtype_error(returntype)
    end
    df
end

function finite_difference_gradient!(df::StridedVector{<:Number}, f, x::StridedVector{<:Number},
    cache::GradientCache{T1,T2,T3,fdtype,returntype,inplace}) where {T1,T2,T3,fdtype,returntype,inplace}

    # c1 is x1 if we need a complex copy of x, otherwise Nothing
    # c2 is Nothing
    fx, c1, c2 = cache.fx, cache.c1, cache.c2
    if fdtype != Val{:complex}
        epsilon_factor = compute_epsilon_factor(fdtype, eltype(x))
        if eltype(df)<:Complex && !(eltype(x)<:Complex)
            copyto!(c1,x)
        end
    end
    if fdtype == Val{:forward}
        for i ∈ eachindex(x)
            epsilon = compute_epsilon(fdtype, x[i], epsilon_factor)
            x_old = x[i]
            if typeof(fx) != Nothing
                x[i] += epsilon
                dfi = (f(x) - fx) / epsilon
                x[i] = x_old
            else
                fx0 = f(x)
                x[i] += epsilon
                dfi = (f(x) - fx0) / epsilon
                x[i] = x_old
            end

            df[i] = real(dfi)
            if eltype(df)<:Complex
                if eltype(x)<:Complex
                    x[i] += im * epsilon
                    if typeof(fx) != Nothing
                        dfi = (f(x) - fx) / (im*epsilon)
                    else
                        dfi = (f(x) - fx0) / (im*epsilon)
                    end
                    x[i] = x_old
                else
                    c1[i] += im * epsilon
                    if typeof(fx) != Nothing
                        dfi = (f(c1) - fx) / (im*epsilon)
                    else
                        dfi = (f(c1) - fx0) / (im*epsilon)
                    end
                    c1[i] = x_old
                end
                df[i] -= im * imag(dfi)
            end
        end
    elseif fdtype == Val{:central}
        @inbounds for i ∈ eachindex(x)
            epsilon = compute_epsilon(fdtype, x[i], epsilon_factor)
            x_old = x[i]
            x[i] += epsilon
            dfi = f(x)
            x[i] = x_old - epsilon
            dfi -= f(x)
            x[i] = x_old
            df[i] = real(dfi / (2*epsilon))
            if eltype(df)<:Complex
                if eltype(x)<:Complex
                    x[i] += im*epsilon
                    dfi = f(x)
                    x[i] = x_old - im*epsilon
                    dfi -= f(x)
                    x[i] = x_old
                else
                    c1[i] += im*epsilon
                    dfi = f(c1)
                    c1[i] = x_old - im*epsilon
                    dfi -= f(c1)
                    c1[i] = x_old
                end
                df[i] -= im*imag(dfi / (2*im*epsilon))
            end
        end
    elseif fdtype==Val{:complex} && returntype<:Real && eltype(df)<:Real && eltype(x)<:Real
        copyto!(c1,x)
        epsilon_complex = eps(real(eltype(x)))
        # we use c1 here to avoid typing issues with x
        @inbounds for i ∈ eachindex(x)
            c1_old = c1[i]
            c1[i] += im*epsilon_complex
            df[i]  = imag(f(c1)) / epsilon_complex
            c1[i]  = c1_old
        end
    else
        fdtype_error(returntype)
    end
    df
end

# vector of derivatives of a scalar->vector map
# this is effectively a vector of partial derivatives, but we still call it a gradient
function finite_difference_gradient!(df::AbstractArray{<:Number}, f, x::Number,
    cache::GradientCache{T1,T2,T3,fdtype,returntype,inplace}) where {T1,T2,T3,fdtype,returntype,inplace}

    # NOTE: in this case epsilon is a scalar, we need two arrays for fx1 and fx2
    # c1 denotes fx1, c2 is fx2, sizes guaranteed by the cache constructor
    fx, c1, c2 = cache.fx, cache.c1, cache.c2

    if fdtype == Val{:forward}
        epsilon_factor = compute_epsilon_factor(fdtype, eltype(x))
        epsilon = compute_epsilon(Val{:forward}, x, epsilon_factor)
        if inplace == Val{true}
            f(c1, x+epsilon)
        else
            c1 .= f(x+epsilon)
        end
        if typeof(fx) != Nothing
            @. df = (c1 - fx) / epsilon
        else
            if inplace == Val{true}
                f(c2, x)
            else
                c2 .= f(x)
            end
            @. df = (c1 - c2) / epsilon
        end
    elseif fdtype == Val{:central}
        epsilon_factor = compute_epsilon_factor(fdtype, eltype(x))
        epsilon = compute_epsilon(Val{:central}, x, epsilon_factor)
        if inplace == Val{true}
            f(c1, x+epsilon)
            f(c2, x-epsilon)
        else
            c1 .= f(x+epsilon)
            c2 .= f(x-epsilon)
        end
        @. df = (c1 - c2) / (2*epsilon)
    elseif fdtype == Val{:complex} && returntype <: Real
        epsilon_complex = eps(real(eltype(x)))
        if inplace == Val{true}
            f(c1, x+im*epsilon_complex)
        else
            c1 .= f(x+im*epsilon_complex)
        end
        @. df = imag(c1) / epsilon_complex
    else
        fdtype_error(returntype)
    end
    df
end


text:
include("../preprocess/boostweights.jl")
include("tree_base.jl")
using Random
using CUDA
using DataFrames
struct TREECLASS end
struct Result{T, S}
    sc::Vector{Node{T}}
    d::Vector{S}
end
function fit(::TREECLASS, X, y, rng = Random.GLOBAL_RNG, max_depth = 6, min_node_records = 1,
        n_features_per_node = Int(floor(sqrt(size(X, 2)))), n_trees = 100)
    if n_features_per_node > size(X, 2)
        n_features_per_node = size(X, 2)
    end
    out_classes = unique(y)
    n_classes = length(out_classes)
    in_classes = collect(1:n_classes)
    d = Dict(zip(out_classes, in_classes))
    target = map(z -> d[z], y)

    T = eltype(X)
    sc = Node{T}[]
    nrow = size(X, 1)
    for i in 1:n_trees
        root = Node{T}()
        push!(sc, root)
        ids = sample(1:nrow, nrow)
        X1 = X[ids, :]
        target1 = target[ids]
        dtc = DecisionTreeContainer{T}(root, n_features_per_node, n_classes, max_depth, min_node_records)
        process_node(dtc, root, X1, target1, rng)
    end
    return Result(sc, out_classes)
end

function rf_predict(scr::Result, X)
    d = scr.d

    n = length(d)
    vals = Vector{Int}(undef, n)
    res = Vector{eltype(d)}(undef, size(X, 1))
    for i in axes(X, 1)
        vals .= 0
        for tree in scr.sc
            vals[predict(tree, X, i)] += 1
        end
        res[i] = d[argmax(vals)]
    end
    return res
end
"""
    ## Random Forest Classifier
    ### Description
      The Random Forest Classifier uses a multitude of decision trees to
      solve classification problems.\n
      --------------------\n
    ### Input
      RandomForestClassifier(x, y, rng ;
      max_depth = 6, min_node_records = 1, n_trees = 100)\n
      --------------------\n
      #### Positional Arguments
      Array{Any} - X:: Array of x's for which the model will use to
      predict y.\n
      Array{Any} - Y:: Array of y's for which the x's are used to predict.\n
      RandomNumberGenerator - rng:: Determines the seed for the given model.\n
      #### Key-word Arguments
      Int64 - max_depth:: Determines the max depth which the tree should use
      as a stop parameter.\n
      Int64 - min_node_records:: Determines the minimum number of nodes a
      constructed node is allowed to have.\n
      Int64 - n_trees:: Determines how many decision trees should be trained.
      --------------------\n
     ### Output
     model:: A Lathe Model.\n
     ---------------------\n
     ### Functions
     Model.predict(xt) :: Predicts a new y based on the data provided as xt and
      the weights obtained from X.\n
     ---------------------\n
     ### Data
     storedata :: A tree node type that contains the weights and their
     corresponding values.
       """
mutable struct RandomForestClassifier{P} <: Classifier
    predict::P
    storedata::Result
function RandomForestClassifier(X::Array, Y::Array, rng = Random.GLOBAL_RNG; max_depth = 6,
     min_node_records = 1,
    n_features_per_node = Int(floor(sqrt(size(X, 2)))), n_trees = 100, cuda = false)
    if cuda == true
        X = CuArray(X)
        Y = CuArray(Y)
    end
    storedata = fit(TREECLASS(), X, Y, rng, max_depth, min_node_records,
        Int(floor(sqrt(size(X, 2)))), n_trees)
    predict(xt::Array) = rf_predict(storedata, xt)
    new{typeof(predict)}(predict, storedata)
end
function RandomForestClassifier(X::DataFrame, Y::Array, rng = Random.GLOBAL_RNG;
    max_depth = 6,
     min_node_records = 1,
     weights = NoWeights(Dict()), cuda = false,
    n_features_per_node = Int(floor(sqrt(size(X, 2)))),
     n_trees = 100)
    classifiers = []
    treec = 0
    n_features = size(X)[1]
    divamount = n_trees / n_features
    if cuda == true
        Y = CuArray(Y)
    end
    for data in eachcol(X)
        if cuda == true
            data = CuArray(data)
        end
        mdl = RandomForestClassifier(data, Y, n_trees = divamount)
        push!(classifiers, mdl)
    end
    predict(xt) = _compare_predCat(classifiers, xt)
    new{typeof(predict)}(predict, result)
end
end
"""
    ## Decision Tree Classifier
    ### Description
      The decision tree classifier is a model ideal for solving most
          classification problems.\n
      --------------------\n
    ### Input
      DecisionTreeClassifier(x, y, rng ;
      max_depth = 6, min_node_records = 1)\n
      --------------------\n
      #### Positional Arguments
      Array{Any} - X:: Array of x's for which the model will use to
      predict y.\n
      Array{Any} - Y:: Array of y's for which the x's are used to predict.\n
      RandomNumberGenerator - rng:: Determines the seed for the given model.\n
      #### Key-word Arguments
      Int64 - max_depth:: Determines the max depth which the tree should use
      as a stop parameter.\n
      Int64 - min_node_records:: Determines the minimum number of nodes a
      constructed node is allowed to have.\n
      --------------------\n
     ### Output
     model:: A Lathe Model.\n
     ---------------------\n
     ### Functions
     Model.predict(xt) :: Predicts a new y based on the data provided as xt and
      the weightsz obtained from X.\n
     ---------------------\n
     ### Data
     storedata :: A tree node type that contains the weights and their
     corresponding values.
       """
function DecisionTreeClassifier(X, Y, rng = Random.GLOBAL_RNG; max_depth = 6,
     min_node_records = 1,
    n_features_per_node = Int(floor(sqrt(size(X, 2)))))
    storedata = fit(TREECLASS(), X, Y, rng, max_depth, min_node_records,
        Int(floor(sqrt(size(X, 2)))), 1)
    predict(xt) = rf_predict(storedata, xt)
    (var)->(predict;storedata)
end


text:
module BangBang

# Use README as the docstring of the module:
@doc read(joinpath(dirname(@__DIR__), "README.md"), String) BangBang

export @!,
       @set!!,
       Empty,
       append!!,
       delete!!,
       deleteat!!,
       empty!!,
       lmul!!,
       materialize!!,
       mul!!,
       pop!!,
       popfirst!!,
       prefermutation,
       push!!,
       pushfirst!!,
       rmul!!,
       setindex!!,
       setproperties!!,
       setproperty!!,
       singletonof,
       splice!!

using Base.Broadcast: materialize!
using Compat: hasproperty
using InitialValues
using LinearAlgebra
using Requires
using Tables: Tables
using ZygoteRules: @adjoint

# Used in NoBang:
function ismutable end
function push!! end

include("NoBang/NoBang.jl")
using .NoBang: Empty, ImmutableContainer, singletonof

include("core.jl")
include("base.jl")
include("linearalgebra.jl")
include("initials.jl")
include("macro.jl")
include("dataframes_impl.jl")
include("zygote.jl")

include("setfield.jl")
using .SetfieldImpl: @set!!, prefermutation

function __init__()
    @require StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" begin
        include("staticarrays.jl")
    end
    @require StructArrays = "09ab397b-f2b6-538f-b94a-2f83cf4a842a" begin
        include("structarrays.jl")
    end
    @require DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" begin
        include("dataframes.jl")
    end
end

end # module


text: const ComputableType = Union{AbstractFloat,Complex}


text:
using ExamplePkg

# To avoid compiling the command-line parsing, here we moved that into the module, as well.
x, mu, sigma = ExamplePkg.parse_args_for_demo(ARGS)

result = pdf_normal(x,mu,sigma)
println("prob(x | mu, sigma) = ", result)


text:
"""
"""
struct RestrictedTriangulation{Dc,Dp,G} <: Triangulation{Dc,Dp}
  oldtrian::G
  cell_to_oldcell::Vector{Int}
  oldcell_to_cell::Vector{Int}
  void_to_oldcell::Vector{Int}

  function RestrictedTriangulation(
    oldtrian::Triangulation{Dc,Dp},
    cell_to_oldcell::Vector{Int},
    oldcell_to_cell::Vector{Int},
    void_to_oldcell::Vector{Int}) where {Dc,Dp}

    new{Dc,Dp,typeof(oldtrian)}(
      oldtrian,cell_to_oldcell,oldcell_to_cell,void_to_oldcell)
  end

end

function RestrictedTriangulation(
  oldtrian::Triangulation{Dc,Dp},cell_to_oldcell::Vector{Int}) where {Dc,Dp}
  n_oldcells = num_cells(oldtrian)
  oldcell_to_cell = fill(Int(UNSET),n_oldcells)
  oldcell_to_cell[cell_to_oldcell] .= 1:length(cell_to_oldcell)
  void_to_oldcell = findall(oldcell_to_cell .== UNSET)
  oldcell_to_cell[void_to_oldcell] .= -(1:length(void_to_oldcell))
  RestrictedTriangulation(oldtrian,cell_to_oldcell,oldcell_to_cell,void_to_oldcell)
end

function RestrictedTriangulation(
  oldtrian::Triangulation{Dc,Dp},oldcell_to_mask::Vector{Bool}) where {Dc,Dp}
  cell_to_oldcell = findall(oldcell_to_mask)
  RestrictedTriangulation(oldtrian,cell_to_oldcell)
end

function get_reffes(trian::RestrictedTriangulation)
  get_reffes(trian.oldtrian)
end

function get_cell_type(trian::RestrictedTriangulation)
  reindex(get_cell_type(trian.oldtrian),trian.cell_to_oldcell)
end

function get_cell_coordinates(trian::RestrictedTriangulation)
  reindex(get_cell_coordinates(trian.oldtrian),trian.cell_to_oldcell)
end

function restrict(f::AbstractArray,trian::RestrictedTriangulation)
  reindex(f,trian)
end

function get_cell_id(trian::RestrictedTriangulation)
  trian.cell_to_oldcell
end

function get_cell_map(trian::RestrictedTriangulation)
  cell_map = get_cell_map(trian.oldtrian)
  reindex(cell_map,trian.cell_to_oldcell)
end


text:
"""
__author__ = "Gil Miranda"
__copyright__ = "MIT"
__credits__ = ["Gil Miranda", "Ricardo R. Rosa", "Beatriz Farah"]
__license__ = "MIT"
__version__ = "0.1"
__email__ = "gil@matematica.ufrj.br"
__status__ = "Production
"""

include("neuralactivations.jl")
include("costfunctions.jl")

mutable struct NeuralNetwork
    """
    Mutable Struct for the Artificial Neural Network

    Fields:
        layer_dimensions: The number of neurons on each layer.
        θ: The parameters of the NN, it contains a dictionary storing all weights and biases.
        ϕ: A dictionary storing all activation functions for each layer.
    """

    layer_dimensions::Array
    θ::Dict
    ϕ::Dict
end

  function NeuralNetwork(layer_dimensions, all_activations = "relu")
    """
    Constructor for the NeuralNetwork struct.

    Args:
        layer_dimensions: The number of neurons on each layer.
        all_activations: Set all layers to this activation function.

    Returns:
        A mutable struct NeuralNetwork

    Raises:
        KeyError: Raises an exception.
    """

    @assert(all_activations ∈ ACTV_FUNCS, "Please, use a valid activation function from $ACTV_FUNCS")

    parameters, activations = init_param(layer_dimensions, all_activations)
    NeuralNetwork([l[1] for l in layer_dimensions], parameters, activations)
end

function init_param(layer_dimensions::Vector{Int64}, all_activations = "relu")
    """
    Function that generates the dictionary of parameters for the NeuralNetwork struct.
    This function will be called if all activation functions for each layer must be the same.

    Args:
        layer_dimensions: Array of `Int` with number of neurons on each layer
        all_activations: Set all layers to this activation function.

    Returns:
       param: Dict with weights and biases
       actv: Dict with activation functions
    """

    param = Dict()
    actv = Dict()
    for l = 2:length(layer_dimensions)
        param["W_$(l-1)"] = rand(layer_dimensions[l], layer_dimensions[l-1])*.1
        param["b_$(l-1)"] = zeros(layer_dimensions[l],1)
        actv["ϕ_$(l-1)"] = all_activations
    end
    return param, actv
end

function init_param(layer_dimensions::Vector{Any}, all_activations = "relu")
    """
    Function that generates the dictionary of parameters for the NeuralNetwork struct.
    This function will be called if some activation functions is set differently from others.

    Args:
        layer_dimensions: Array of `Int` or `Tuple` of (Int, String) with number of neurons on each layer
        all_activations: Set all default layers to this activation function.

    Returns:
       param: Dict with weights and biases
       actv: Dict with activation functions
    """

    param = Dict()
    actv = Dict()
    @assert(all_activations ∈ ACTV_FUNCS, "Please, use a valid activation function from $ACTV_FUNCS")
    for l = 2:length(layer_dimensions)
        if length(layer_dimensions[l-1]) > 1
            @assert(layer_dimensions[l-1][2] ∈ ACTV_FUNCS, "Please, use a valid activation function from $ACTV_FUNCS")

            param["W_$(l-1)"] = rand(layer_dimensions[l][1], layer_dimensions[l-1][1])*.1
            param["b_$(l-1)"] = zeros(layer_dimensions[l][1],1)
            actv["ϕ_$(l-1)"] = layer_dimensions[l-1][2]
        else
            param["W_$(l-1)"] = rand(layer_dimensions[l][1], layer_dimensions[l-1][1])*.1
            param["b_$(l-1)"] = zeros(layer_dimensions[l][1],1)
            actv["ϕ_$(l-1)"] = all_activations
        end
    end
    return param, actv
end

function predict(nn::NeuralNetwork, X)
    """
    Model prediction function.
    Makes one forward pass and returns the output of the NeuralNetwork struct, given an input

    Args:
        nn: Artificial Neural Network
        X: Input data to feed the Neural Network

    Returns:
       activation: Output layer for the Neural Network
    """
    activation = X
    for l in 1:div(length(nn.θ),2)
        activation, _  = forward_activation(activation, nn.θ["W_$l"], nn.θ["b_$l"], nn.ϕ["ϕ_$l"])
    end
    return activation
end

function forward_linear(A,w,b)
    """
    Affine transformation on the neuron input data.
    Given an input A, calculates Z = wA+b

    Args:
        A: Last layer output to feed current layer input
        w: Array of weights for the current
        b: Array of biases for the current layer

    Returns:
       Z: Affine transformation on the input data
    """

    Z = w*A .+ b
    return Z
end

function forward_activation(A_pre, W, b, func = "sigmoid")
    """
    Makes one step of a forward pass.
    Given an input A, call forward_linear for Z = wA+b then apply ϕ(Z)

    Args:
        A_pre: Last layer output to feed current layer input
        W: Array of weights for the current
        b: Array of biases for the current layer
        func = Activation function for the current layer

    Returns:
       A: Activation function ϕ applied on Z
       Z: Affine transformation on the input data
    """

    @assert(func ∈ ACTV_FUNCS, "Please, use a valid activation function from $ACTV_FUNCS")
    if func == "sigmoid" || func == "σ"
        activation_func = σ
    elseif func == "relu"
        activation_func = relu
    end
    Z = forward_linear(A_pre,W,b)
    A = activation_func(Z)

    return A, Z
end

# function loss(Ŷ, Y)
#     return 1/length(Y) * sum(abs2, Y .- Ŷ)
# end

function d_activation(z, actv)
    """
    Call function to calculate ϕ'(z) based on current layer activations.

    Args:
       z: Affine input data from last layer.
    actv: Current layer activation function.

    Returns:
       deriv: Derivative of activation function ϕ applied on z with respect to z
    """

    @assert(actv ∈ ACTV_FUNCS, "Please, use a valid activation function from $ACTV_FUNCS")

    if actv == "σ" || actv == "sigmoid"
        deriv = dσ(z)
    elseif actv == "relu"
        deriv = drelu(z)
    elseif actv == "tanh"
        deriv = dtanh(z)
    else
        deriv = 0
    end
    return deriv
end

function GradDescent!(nn::NeuralNetwork, grads, η, batch_size = 1)
    """
    Function that makes a Stochastic Gradient Descent step.

    Args:
            nn: NeuralNetwork
         grads: Dict of gradients calculates with backpropagation
             η: SGD learning rate
    batch_size: Size of a batch of learning examples

    Returns:
       This function does not return anything, it simply change the values on θ field for the NN struct
    """
    for l=1:div(length(nn.θ),2)
        nn.θ["W_$l"] -= (η/batch_size) * grads["W_$l"]
        nn.θ["b_$l"] -= (η/batch_size) * grads["b_$l"]
    end
end

function backprop(nn::NeuralNetwork, X, Y, cost_function = loss().mse)
    """
    Backpropagation Algorithm.
    Reverse mode of Automatic Differentiation for Neural Networks.

    Args:
    nn: NeuralNetwork
     X: Input data for the Neural Network
     Y: Output labels for the Neural Network (supervised learning)

    Returns:
       grads: Dict with the gradient for Cost function with respect to the Neural Network parameters
    """

    grad = Dict()

    activation = X
    activation_cache = Any[X]
    z_cache = Any[]
    L = div(length(nn.θ),2)

    for l in 1:L
        activation, z  = forward_activation(activation, nn.θ["W_$l"], nn.θ["b_$l"], nn.ϕ["ϕ_$l"])
        push!(z_cache, z)
        push!(activation_cache, activation)
    end
        delta = cost_function(activation_cache[end], Y, true) .* d_activation(z_cache[end], nn.ϕ["ϕ_$L"])

    grad["b_$L"] = sum(delta,dims=2)
    grad["W_$L"] = delta*activation_cache[end-1]'
    for l in L:-1:2
        z = z_cache[l-1]
        delta = nn.θ["W_$(l)"]'*delta .* d_activation(z, nn.ϕ["ϕ_$l"])
        grad["b_$(l-1)"] = sum(delta,dims=2)
        grad["W_$(l-1)"] = (delta*activation_cache[l-1]')
    end
    return grad
end

function train!(nn::NeuralNetwork, X, Y; η = 0.1, epochs = 100, verbose = true, batch_size = 2, cost_function = loss().mse)
    """
    Function that executes the training for the Neural Network
    """
    N = size(Y)[2]

    costs = []
    epch = []
    accuracy = []
    for k in 1:epochs

        randomized_index = shuffle(1:N)
        xrandom = X[:, randomized_index]
        yrandom = Y[:, randomized_index]
        batch_index = []
        for k in 1:batch_size:div(N,batch_size)*batch_size
            push!(batch_index, (k,k+batch_size-1))
        end
        if N%batch_size != 0
            push!(batch_index, (div(N,batch_size)*batch_size, div(N,batch_size)*batch_size + N%batch_size))
        end
        for i in batch_index
            grads = []
            grads = backprop(nn, xrandom[:,i[1]:i[2]], yrandom[:,i[1]:i[2]], cost_function)
            ## Fazer a media
            GradDescent!(nn, grads, η, batch_size)
        end

        cost = cost_function(predict(nn,xrandom), yrandom, false)
        acc = check_accuracy(nn, xrandom, yrandom)

        push!(costs, cost)
        push!(epch, k)
        push!(accuracy, acc)

        if verbose == true
            println("Loss at epoch $k = $cost")
        end
    end

    return (costs=costs, epochs=epch, accuracy=accuracy)
end


text:
import Base: convert, size, eltype, ~
export Variable

struct Variable{T, N, A <: AbstractArray{T,N}} <: AbstractExpression
	x::A
end

# constructors
"""
	Variable([T::Type,] dims...)

Returns a `Variable` of dimension `dims` initialized with an array of all zeros.

`Variable(x::AbstractArray)`

Returns a `Variable` of dimension `size(x)` initialized with `x`

"""
function Variable(T::Type, args::Vararg{I,N}) where {I <: Integer,N}
	Variable{T,N,Array{T,N}}(zeros(T, args...))
end

function Variable(args::Vararg{I}) where {I <: Integer}
  Variable(zeros(args...))
end

# Utils

function Base.show(io::IO, x::Variable)
  print(io, "Variable($(eltype(x.x)), $(size(x.x)))")
end

"""
	~(x::Variable)

Returns the `Array` of the variable `x`
"""
~(x::Variable) = x.x
~(x::Tuple{Variable}) = (~)(x[1])
~(x::NTuple{N,Variable}) where {N} = ArrayPartition((~).(x))

"""
size(x::Variable, [dim...])

Like `size(A::AbstractArray, [dims...])` returns the tuple containing the dimensions of the variable `x`.
"""
size(x::Variable) = size(x.x)
size(x::Variable, dim::I) where { I <: Integer} = size(x.x, dim)

"""
eltype(x::Variable)

Like `eltype(x::AbstractArray)` returns the type of the elements of `x`.
"""
eltype(x::Variable) = eltype(x.x)


text:
"""
`vstorent!` (non-temporal store) requires data to be aligned.
`alignstores!` will align `y` in preparation for the non-temporal maps.
"""
function setup_vmap!(
    f::F, y::AbstractArray{T}, ::Val{true},
    args::Vararg{AbstractArray,A}
) where {F, T <: Base.HWReal, A}
    N = length(y)
    ptry = VectorizationBase.zstridedpointer(y)
    ptrargs = VectorizationBase.zstridedpointer.(args)
    V = pick_vector_width(T)
    W = unwrap(V)
    zero_index = MM{W}(Static(0))
    uintptry = reinterpret(UInt, pointer(ptry))
    @assert iszero(uintptry & (sizeof(T) - 1)) "The destination vector (`dest`) must be aligned to `sizeof(eltype(dest)) == $(sizeof(T))` bytes."
    alignment = uintptry & (register_size() - 1)
    if alignment > 0
        i = reinterpret(Int, W - (alignment >>> VectorizationBase.intlog2(sizeof(T))))
        m = mask(T, i)
        if N < i
            m &= mask(T, N & (W - 1))
        end
        vstore!(ptry, f(vload.(ptrargs, ((zero_index,),), m)...), (zero_index,), m, False(), True(), False(), register_size())
        gesp(ptry, (i,)), gesp.(ptrargs, ((i,),)), N - i
    else
        ptry, ptrargs, N
    end
end

@inline function setup_vmap!(f, y, ::Val{false}, args::Vararg{AbstractArray,A}) where {A}
    N = length(y)
    ptry = VectorizationBase.zstridedpointer(y)
    ptrargs = VectorizationBase.zstridedpointer.(args)
    ptry, ptrargs, N
end

function vmap_singlethread!(
    f::F, y::AbstractArray{T},
    ::Val{NonTemporal},
    args::Vararg{AbstractArray,A}
) where {F,T <: Base.HWReal, A, NonTemporal}
    ptry, ptrargs, N = setup_vmap!(f, y, Val{NonTemporal}(), args...)
    vmap_singlethread!(f, ptry, Zero(), N, Val{NonTemporal}(), ptrargs)
    nothing
end
function vmap_singlethread!(
    f::F, ptry::AbstractStridedPointer{T},
    start, N, ::Val{NonTemporal},
    ptrargs::Tuple{Vararg{AbstractStridedPointer,A}}
) where {F, T, NonTemporal, A}
    i = convert(Int, start)
    V = VectorizationBase.pick_vector_width(T)
    W = unwrap(V)
    st = VectorizationBase.static_sizeof(T)
    UNROLL = 4
    LOG2UNROLL = 2
    while i < N - ((W << LOG2UNROLL) - 1)
        index = VectorizationBase.Unroll{1,W,UNROLL,1,W,0x0000000000000000}((i,))
        v = f(vload.(ptrargs, index)...)
        if NonTemporal
            vstore!(ptry, v, index, True(), True(), True(), register_size())
        else
            vstore!(ptry, v, index, False(), True(), False(), register_size())
        end
        i = vadd_fast(i, StaticInt{UNROLL}() * W)
    end
    # if Base.libllvm_version ≥ v"11" # this seems to be slower
    #     Nm1 = vsub_fast(N, 1)
    #     while i < N # stops at 16 when
    #         m = mask(V, i, Nm1)
    #         vnoaliasstore!(ptry, f(vload.(ptrargs, ((MM{W}(i),),), m)...), (MM{W}(i,),), m)
    #         i = vadd_fast(i, W)
    #     end
    # else
    while i < N - (W - 1) # stops at 16 when
        vᵣ = f(vload.(ptrargs, ((MM{W}(i),),))...)
        if NonTemporal
            vstore!(ptry, vᵣ, (MM{W}(i),), True(), True(), True(), register_size())
        else
            vstore!(ptry, vᵣ, (MM{W}(i),), False(), True(), False(), register_size())
        end
        i = vadd_fast(i, W)
    end
    if i < N
        m = mask(T, N & (W - 1))
        if NonTemporal
            vstore!(ptry, f(vload.(ptrargs, ((MM{W}(i),),), m)...), (MM{W}(i,),), m, True(), True(), False(), register_size())
        else
            vstore!(ptry, f(vload.(ptrargs, ((MM{W}(i),),), m)...), (MM{W}(i,),), m, False(), True(), False(), register_size())
        end
    end
    # end
    nothing
end

abstract type AbstractVmapClosure{NonTemporal,F,D,N,A<:Tuple{Vararg{StridedPointer,N}}} <: Function end
struct VmapClosure{NonTemporal,F,D,N,A} <: AbstractVmapClosure{NonTemporal,F,D,N,A}
    f::F
    function VmapClosure{NonTemporal}(f::F, ::D, ::A) where {NonTemporal,F,D,N,A<:Tuple{Vararg{StridedPointer,N}}}
        new{NonTemporal,F,D,N,A}(f)
    end
end
# struct VmapKnownClosure{NonTemporal,F,D,N,A} <: AbstractVmapClosure{NonTemporal,F,D,N,A} end

@inline function _vmap_thread_call!(
    f::F, p::Ptr{UInt}, ::Type{D}, ::Type{A}, ::Val{NonTemporal}
) where {F,D,A,NonTemporal}
    (offset, dest) = ThreadingUtilities.load(p, D, 1)
    (offset, args) = ThreadingUtilities.load(p, A, offset)

    (offset, start) = ThreadingUtilities.load(p, Int, offset)
    (offset, stop ) = ThreadingUtilities.load(p, Int, offset)

    vmap_singlethread!(f, dest, start, stop, Val{NonTemporal}(), args)
    nothing
end
# @generated function (::VmapKnownClosure{NonTemporal,F,D,N,A})(p::Ptr{UInt})  where {NonTemporal,F,D,N,A}
#     :(_vmap_thread_call!($(F.instance), p, $D, $A, Val{$NonTemporal}()))
# end
function (m::VmapClosure{NonTemporal,F,D,N,A})(p::Ptr{UInt}) where {NonTemporal,F,D,N,A}
    _vmap_thread_call!(m.f, p, D, A, Val{NonTemporal}())
end

@inline function _get_fptr(cfunc::Base.CFunction)
    Base.unsafe_convert(Ptr{Cvoid}, cfunc)
end
# @generated function _get_fptr(cfunc::F) where {F<:VmapKnownClosure}
#     precompile(F(), (Ptr{UInt},))
#     quote
#         $(Expr(:meta,:inline))
#         @cfunction($(F()), Cvoid, (Ptr{UInt},))
#     end
# end

@inline function setup_thread_vmap!(
    p, cfunc, ptry, ptrargs, start, stop
)
    fptr = _get_fptr(cfunc)
    offset = ThreadingUtilities.store!(p, fptr, 0)
    offset = ThreadingUtilities.store!(p, ptry, offset)
    offset = ThreadingUtilities.store!(p, ptrargs, offset)
    offset = ThreadingUtilities.store!(p, start, offset)
    offset = ThreadingUtilities.store!(p, stop, offset)
    nothing
end
@inline function launch_thread_vmap!(tid, cfunc, ptry, ptrargs, start, stop)
    p = ThreadingUtilities.taskpointer(tid)
    while true
        if ThreadingUtilities._atomic_cas_cmp!(p, ThreadingUtilities.SPIN, ThreadingUtilities.STUP)
            setup_thread_vmap!(p, cfunc, ptry, ptrargs, start, stop)
            @assert ThreadingUtilities._atomic_cas_cmp!(p, ThreadingUtilities.STUP, ThreadingUtilities.TASK)
            return
        elseif ThreadingUtilities._atomic_cas_cmp!(p, ThreadingUtilities.WAIT, ThreadingUtilities.STUP)
            setup_thread_vmap!(p, cfunc, ptry, ptrargs, start, stop)
            @assert ThreadingUtilities._atomic_cas_cmp!(p, ThreadingUtilities.STUP, ThreadingUtilities.LOCK)
            ThreadingUtilities.wake_thread!(tid % UInt)
            return
        end
        ThreadingUtilities.pause()
    end
end

@inline function vmap_closure(f::F, ptry::D, ptrargs::A, ::Val{NonTemporal}) where {F,D<:StridedPointer,N,A<:Tuple{Vararg{StridedPointer,N}},NonTemporal}
    vmc = VmapClosure{NonTemporal}(f, ptry, ptrargs)
    @cfunction($vmc, Cvoid, (Ptr{UInt},))
end
# @inline function _cfunc_closure(f, ptry, ptrargs, ::Val{NonTemporal}) where {NonTemporal}
#     vmc = VmapClosure{NonTemporal}(f, ptry, ptrargs)
#     @cfunction($vmc, Cvoid, (Ptr{UInt},))
# end
# @generated function vmap_closure(f::F, ptry::D, ptrargs::A, ::Val{NonTemporal}) where {F,D<:StridedPointer,N,A<:Tuple{Vararg{StridedPointer,N}},NonTemporal}
#     # fsym = get(FUNCTIONSYMBOLS, F, Symbol("##NOTFOUND##"))
#      # fsym === Symbol("##NOTFOUND##")
#     if false# iszero(sizeof(F))
#         quote
#             $(Expr(:meta,:inline))
#             VmapKnownClosure{$NonTemporal,$F,$D,$N,$A}()
#         end
#     else
#         quote
#             $(Expr(:meta,:inline))
#             _cfunc_closure(f, ptry, ptrargs, Val{$NonTemporal}())
#         end
#     end
# end

function vmap_multithread!(
    f::F,
    y::AbstractArray{T},
    ::Val{NonTemporal},
    args::Vararg{AbstractArray,A}
) where {F,T,A,NonTemporal}
    W, Wshift = VectorizationBase.pick_vector_width_shift(T)
    ptry, ptrargs, N = setup_vmap!(f, y, Val{NonTemporal}(), args...)
    # nt = min(Threads.nthreads(), VectorizationBase.SYS_CPU_THREADS, N >> (Wshift + 3))
    nt = min(Threads.nthreads(), VectorizationBase.num_cores(), N >> (Wshift + 5))

    # if !((nt > 1) && iszero(ccall(:jl_in_threaded_region, Cint, ())))
    if nt < 2
        vmap_singlethread!(f, ptry, Zero(), N, Val{NonTemporal}(), ptrargs)
        return
    end

    cfunc = vmap_closure(f, ptry, ptrargs, Val{NonTemporal}())
    vmc = VmapClosure{NonTemporal}(f, ptry, ptrargs)
    Nveciter = (N + (W-1)) >> Wshift
    Nd, Nr = divrem(Nveciter, nt)
    NdW = Nd << Wshift
    NdWr = NdW + W
    GC.@preserve cfunc begin
        start = 0
        for tid ∈ 1:nt-1
            stop = start + ifelse(tid ≤ Nr, NdWr, NdW)
            launch_thread_vmap!(tid, cfunc, ptry, ptrargs, start, stop)
            start = stop
        end
        vmap_singlethread!(f, ptry, start, N, Val{NonTemporal}(), ptrargs)
        for tid ∈ 1:nt-1
            ThreadingUtilities.__wait(tid)
        end
    end
    nothing
end
@generated function gc_preserve_vmap!(f::F,
    y::AbstractArray,
    ::Val{NonTemporal},
    ::Val{Threaded},
    args::Vararg{AbstractArray,A}
) where {F,A,NonTemporal,Threaded}
    m = Threaded ? :vmap_multithread! : :vmap_singlethread!
    call = Expr(:call, m, :f, :y, Expr(:call, Expr(:curly, :Val, NonTemporal)))
    q = Expr(:block, Expr(:meta, :inline))
    gcpres = Expr(:gc_preserve, call)
    for a ∈ 1:A
        arg = Symbol(:arg_,a)
        parg = Symbol(:parg_,a)
        push!(q.args, Expr(:(=), arg, :(@inbounds args[$a])))#Expr(:ref, :args, a)))
        push!(q.args, Expr(:(=), parg, Expr(:call, :preserve_buffer, arg)))
        push!(call.args, arg)
        push!(gcpres.args, parg)
    end
    push!(q.args, gcpres, :y)
    q
end

@inline _all_dense(t::Tuple{ArrayInterface.True}) = true
@inline _all_dense(t::Tuple{ArrayInterface.True,ArrayInterface.True,Vararg}) = _all_dense(Base.tail(t))
@inline _all_dense(t::Tuple{ArrayInterface.True,ArrayInterface.False,Vararg}) = false
@inline _all_dense(t::Tuple{ArrayInterface.False,Vararg}) = false
@inline all_dense() = true
@inline all_dense(A::Array) = true
@inline all_dense(A::AbstractArray) = _all_dense(ArrayInterface.dense_dims(A))
@inline all_dense(A::AbstractArray, B::AbstractArray, C::Vararg{AbstractArray,K}) where {K} = all_dense(A) && all_dense(B, C...)

"""
    vmap!(f, destination, a::AbstractArray)
    vmap!(f, destination, a::AbstractArray, b::AbstractArray, ...)
Vectorized-`map!`, applying `f` to each element of `a` (or paired elements of `a`, `b`, ...)
and storing the result in `destination`.
"""
function vmap!(
    f::F, y::AbstractArray, arg1::AbstractArray, arg2::AbstractArray, args::Vararg{AbstractArray,A}
) where {F,A}
    if check_args(y, arg1, arg2, args...) && all_dense(y, arg1, arg2, args...)
        gc_preserve_vmap!(f, y, Val{false}(), Val{false}(), arg1, arg2, args...)
    else
        map!(f, y, arg1, arg2, args...)
    end
end
function vmap!(f::F, y::AbstractArray, arg::AbstractArray) where {F}
    if check_args(y, arg) && all_dense(y, arg)
        gc_preserve_vmap!(f, y, Val{false}(), Val{false}(), arg)
    else
        map!(f, y, arg)
    end
end

"""
    vmapt!(::Function, dest, args...)
Like `vmap!` (see `vmap!`), but uses `Threads.@threads` for parallel execution.
"""
function vmapt!(
    f::F, y::AbstractArray, args::Vararg{AbstractArray,A}
) where {F,A}
    if check_args(y, args...) && all_dense(y, args...)
        gc_preserve_vmap!(f, y, Val{false}(), Val{true}(), args...)
    else
        map!(f, y, args...)
    end
end

"""
    vmapnt!(::Function, dest, args...)
This is a vectorized map implementation using nontemporal store operations. This means that the write operations to the destination will not go to the CPU's cache.
If you will not immediately be reading from these values, this can improve performance because the writes won't pollute your cache. This can especially be the case if your arguments are very long.
```julia
julia> using LoopVectorization, BenchmarkTools
julia> x = rand(10^8); y = rand(10^8); z = similar(x);
julia> f(x,y) = exp(-0.5abs2(x - y))
f (generic function with 1 method)
julia> @benchmark map!(f, \$z, \$x, \$y)
BenchmarkTools.Trial:
  memory estimate:  0 bytes
  allocs estimate:  0
  --------------
  minimum time:     439.613 ms (0.00% GC)
  median time:      440.729 ms (0.00% GC)
  mean time:        440.695 ms (0.00% GC)
  maximum time:     441.665 ms (0.00% GC)
  --------------
  samples:          12
  evals/sample:     1
julia> @benchmark vmap!(f, \$z, \$x, \$y)
BenchmarkTools.Trial:
  memory estimate:  0 bytes
  allocs estimate:  0
  --------------
  minimum time:     178.147 ms (0.00% GC)
  median time:      178.381 ms (0.00% GC)
  mean time:        178.430 ms (0.00% GC)
  maximum time:     179.054 ms (0.00% GC)
  --------------
  samples:          29
  evals/sample:     1
julia> @benchmark vmapnt!(f, \$z, \$x, \$y)
BenchmarkTools.Trial:
  memory estimate:  0 bytes
  allocs estimate:  0
  --------------
  minimum time:     144.183 ms (0.00% GC)
  median time:      144.338 ms (0.00% GC)
  mean time:        144.349 ms (0.00% GC)
  maximum time:     144.641 ms (0.00% GC)
  --------------
  samples:          35
  evals/sample:     1
```
"""
function vmapnt!(
    f::F, y::AbstractArray, args::Vararg{AbstractArray,A}
) where {F,A}
    if check_args(y, args...) && all_dense(y, args...)
        gc_preserve_vmap!(f, y, Val{true}(), Val{false}(), args...)
    else
        map!(f, y, args...)
    end
end

"""
    vmapntt!(::Function, dest, args...)
Like `vmapnt!` (see `vmapnt!`), but uses `Threads.@threads` for parallel execution.
"""
function vmapntt!(
    f::F, y::AbstractArray, args::Vararg{AbstractArray,A}
) where {F,A}
    if check_args(y, args...) && all_dense(y, args...)
        gc_preserve_vmap!(f, y, Val{true}(), Val{true}(), args...)
    else
        map!(f, y, args...)
    end
end

# generic fallbacks
@inline vmap!(f, args...) = map!(f, args...)
@inline vmapt!(f, args...) = map!(f, args...)
@inline vmapnt!(f, args...) = map!(f, args...)
@inline vmapntt!(f, args...) = map!(f, args...)

function vmap_call(f::F, vm!::V, args::Vararg{Any,N}) where {V,F,N}
    T = Base._return_type(f, Base.Broadcast.eltypes(args))
    dest = similar(first(args), T)
    vm!(f, dest, args...)
end

"""
    vmap(f, a::AbstractArray)
    vmap(f, a::AbstractArray, b::AbstractArray, ...)
SIMD-vectorized `map`, applying `f` to each element of `a` (or paired elements of `a`, `b`, ...)
and returning a new array.
"""
vmap(f::F, args::Vararg{Any,N}) where {F,N} = vmap_call(f, vmap!, args...)

"""
    vmapt(f, a::AbstractArray)
    vmapt(f, a::AbstractArray, b::AbstractArray, ...)
A threaded variant of [`vmap`](@ref).
"""
vmapt(f::F, args::Vararg{Any,N}) where {F,N} = vmap_call(f, vmapt!, args...)

"""
    vmapnt(f, a::AbstractArray)
    vmapnt(f, a::AbstractArray, b::AbstractArray, ...)
A "non-temporal" variant of [`vmap`](@ref). This can improve performance in cases where
`destination` will not be needed soon.
"""
vmapnt(f::F, args::Vararg{Any,N}) where {F,N} = vmap_call(f, vmapnt!, args...)

"""
    vmapntt(f, a::AbstractArray)
    vmapntt(f, a::AbstractArray, b::AbstractArray, ...)
A threaded variant of [`vmapnt`](@ref).
"""
vmapntt(f::F, args::Vararg{Any,N}) where {F,N} = vmap_call(f, vmapntt!, args...)


text:
# Testing reference values and precisions
# Each test block of varr and parr should be followed by an append to refVals, refPrecs arrays.
# e.g.
#   refVals=[]
#   refPrecs=[]
#
#   varr = ..........
#   par  = ..........
#
#   append!(refVals ,[ varr ] )
#   append!(refPrecs,[ parr ] )
#
#   varr = ..........
#   par  = ..........
#
#   append!(refVals ,[ varr ] )
#   append!(refPrecs,[ parr ] )
#
#   varr = ..........
#   par  = ..........
#
#   append!(refVals ,[ varr ] )
#   append!(refPrecs,[ parr ] )
#
#   etc.....
#
#   Now for real!
#
#! format: off

    refVals=[]
    refPrecs=[]
# SC ========== Test number 1 reference values and precision match template. =======
# SC ========== /Users/chrishill/projects/clima/cm/test/Ocean/HydrostaticBoussinesq/test_ocean_gyre.jl test reference values ======================================
# BEGIN SCPRINT
# varr - reference values (from reference run)
# parr - digits match precision (hand edit as needed)
#
# [
#  [ MPIStateArray Name, Field Name, Maximum, Minimum, Mean, Standard Deviation ],
#  [         :                :          :        :      :          :           ],
# ]
varr = [
 [       "Q",   "u[1]", -2.18427865357219835873e-02,  4.54905817704320050709e-02,  2.91799084468619303323e-03,  9.86739717788472255056e-03 ],
 [       "Q",   "u[2]", -6.47928098321387119229e-02,  7.44690631237251432495e-02, -1.82439830964300966215e-03,  1.02482752029024841434e-02 ],
 [       "Q",       :η, -6.35241759974374819997e-01,  6.25677877693153861038e-01, -8.60068718305024041901e-04,  2.24583755068219675932e-01 ],
 [       "Q",       :θ,  9.03726335555428482264e-05,  9.03968461560725344839e+00,  2.49953258705220227043e+00,  2.19711762918947606238e+00 ],
 [   "s_aux",       :y,  0.00000000000000000000e+00,  4.00000000000000046566e+06,  2.00000000000000000000e+06,  1.18025270281469495967e+06 ],
 [   "s_aux",       :w, -4.50036075461459223737e-05,  4.09915220804811021062e-05,  3.06722602595587297903e-07,  8.22999489798410955438e-06 ],
 [   "s_aux",    :pkin, -9.00231747934270076783e-01,  0.00000000000000000000e+00, -3.32080874523538072118e-01,  2.56218992025050162908e-01 ],
 [   "s_aux",     :wz0, -2.22688579102610652283e-05,  1.52895210358102938127e-05, -2.78589672590168561413e-08,  7.72520757253795920232e-06 ],
 [ "s_gflux", "ν∇u[1]", -1.84915426791050691892e-03,  1.81823541139391206892e-03,  4.41670652876216560588e-08,  1.43003334026535461722e-04 ],
 [ "s_gflux", "ν∇u[2]", -2.18923622988041743920e-03,  1.83543406028729679316e-03,  3.98744828991449433208e-06,  3.02040505665129290977e-04 ],
 [ "s_gflux", "ν∇u[3]", -1.31369743560772895408e-05,  6.60039477945146449384e-06, -1.12119891567145960325e-08,  9.06617548358093396365e-07 ],
 [ "s_gflux", "ν∇u[4]", -3.29092487463896604671e-03,  3.71582885889340866492e-03, -2.53728229769815408283e-06,  2.97060576409856780689e-04 ],
 [ "s_gflux", "ν∇u[5]", -3.81669587050501284836e-03,  7.74606451320633592264e-04, -5.95763172118155189140e-06,  2.09071039491200258605e-04 ],
 [ "s_gflux", "ν∇u[6]", -7.41189897881484784109e-06,  1.26149616846000803353e-05,  1.41878743800085061125e-08,  1.12070628757683659611e-06 ],
 [ "s_gflux", "κ∇θ[1]", -5.86968584224644812338e-05,  4.34652887121820283387e-05, -2.82719382567508372835e-06,  8.15301391557321884090e-06 ],
 [ "s_gflux", "κ∇θ[2]", -3.70308432284375849564e-05,  3.11721545144766831389e-03,  9.81447035540739714035e-04,  8.31872571479874375368e-04 ],
 [ "s_gflux", "κ∇θ[3]", -2.15846132870457054810e-05,  2.00382267962942220408e-05, -5.17436031104969784121e-07,  1.79150072112032684306e-06 ],
]
parr = [
 [       "Q",   "u[1]",    12,    12,    12,    12 ],
 [       "Q",   "u[2]",    12,    12,    12,    12 ],
 [       "Q",       :η,    12,    12,    12,    12 ],
 [       "Q",       :θ,    12,    12,    12,    12 ],
 [   "s_aux",       :y,    12,    12,    12,    12 ],
 [   "s_aux",       :w,    12,    12,    12,    12 ],
 [   "s_aux",    :pkin,    12,    12,    12,    12 ],
 [   "s_aux",     :wz0,    12,    12,    12,    12 ],
 [ "s_gflux", "ν∇u[1]",    12,    12,    12,    12 ],
 [ "s_gflux", "ν∇u[2]",    12,    12,    12,    12 ],
 [ "s_gflux", "ν∇u[3]",    12,    12,    12,    12 ],
 [ "s_gflux", "ν∇u[4]",    12,    12,    12,    12 ],
 [ "s_gflux", "ν∇u[5]",    12,    12,    12,    12 ],
 [ "s_gflux", "ν∇u[6]",    12,    12,    12,    12 ],
 [ "s_gflux", "κ∇θ[1]",    12,    12,    12,    12 ],
 [ "s_gflux", "κ∇θ[2]",    12,    12,    12,    12 ],
 [ "s_gflux", "κ∇θ[3]",    12,    12,     8,    12 ],
]
# END SCPRINT

    append!(refVals ,[ varr ] )
    append!(refPrecs,[ parr ] )

# SC ========== Test number 2 reference values and precision match template. =======
# SC ========== /Users/chrishill/projects/clima/cm/test/Ocean/HydrostaticBoussinesq/test_ocean_gyre.jl test reference values ======================================
# BEGIN SCPRINT
# varr - reference values (from reference run)
# parr - digits match precision (hand edit as needed)
#
# [
#  [ MPIStateArray Name, Field Name, Maximum, Minimum, Mean, Standard Deviation ],
#  [         :                :          :        :      :          :           ],
# ]
varr = [
 [       "Q",   "u[1]", -2.14403161060213384714e-02,  4.51092912722281316751e-02,  2.91883577365915263327e-03,  9.85543151135863208789e-03 ],
 [       "Q",   "u[2]", -6.47849908268208068973e-02,  7.44377079062340796245e-02, -1.82528856707940708402e-03,  1.02465626974943338490e-02 ],
 [       "Q",       :η, -6.35241717401733962944e-01,  6.25675643032199912952e-01, -8.59997171611655590248e-04,  2.24578663698960928619e-01 ],
 [       "Q",       :θ,  1.21624650983828366162e-04,  9.05348870526675320036e+00,  2.49953463844110279624e+00,  2.19723475068906504148e+00 ],
 [   "s_aux",       :y,  0.00000000000000000000e+00,  4.00000000000000046566e+06,  2.00000000000000000000e+06,  1.18025270281469495967e+06 ],
 [   "s_aux",       :w, -4.48750663401728297594e-05,  4.07425796079212936126e-05,  3.05580091823959372610e-07,  8.20754136359875149074e-06 ],
 [   "s_aux",    :pkin, -9.00231018468760080253e-01,  0.00000000000000000000e+00, -3.32080814803801249724e-01,  2.56218661345592568779e-01 ],
 [   "s_aux",     :wz0, -2.24168211854045644460e-05,  1.54464254869175437296e-05, -2.71999565492576589218e-08,  7.76164707621753630128e-06 ],
 [ "s_gflux", "ν∇u[1]", -1.85093041460835305498e-03,  1.81880838666876047845e-03,  4.56526424711389789871e-08,  1.42643910766357428534e-04 ],
 [ "s_gflux", "ν∇u[2]", -2.20819366719671563229e-03,  1.84128108384750657951e-03,  3.98074575776293245019e-06,  3.03797417171225480129e-04 ],
 [ "s_gflux", "ν∇u[3]", -1.32036367846675870225e-05,  6.56099008285145167548e-06, -1.24087351933463812712e-08,  9.05500396408119435374e-07 ],
 [ "s_gflux", "ν∇u[4]", -3.28419921150994737619e-03,  3.70704618050693088249e-03, -2.54827541969113706709e-06,  2.96228406753828749521e-04 ],
 [ "s_gflux", "ν∇u[5]", -3.80879251696829247351e-03,  7.74834565383260352560e-04, -5.93513695165015881065e-06,  2.08479602570907525654e-04 ],
 [ "s_gflux", "ν∇u[6]", -7.39839704966590930855e-06,  1.26017207509591005658e-05,  1.46129449482955510494e-08,  1.11733219690338497597e-06 ],
 [ "s_gflux", "κ∇θ[1]", -5.81601257251870891614e-05,  4.70677319813339798652e-05, -2.82659431063543911260e-06,  8.16581551142350108366e-06 ],
 [ "s_gflux", "κ∇θ[2]", -3.76891063975906055088e-05,  3.10874693302757654570e-03,  9.82036601802479310747e-04,  8.31189786431030404844e-04 ],
 [ "s_gflux", "κ∇θ[3]", -2.42394160392341078493e-06,  8.68196014640797840416e-08, -5.04918136536487334668e-07,  3.01622429419160335704e-07 ],
]
parr = [
 [       "Q",   "u[1]",    12,    12,    12,    12 ],
 [       "Q",   "u[2]",    12,    12,    12,    12 ],
 [       "Q",       :η,    12,    12,    12,    12 ],
 [       "Q",       :θ,    12,    12,    12,    12 ],
 [   "s_aux",       :y,    12,    12,    12,    12 ],
 [   "s_aux",       :w,    12,    12,    12,    12 ],
 [   "s_aux",    :pkin,    12,    12,    12,    12 ],
 [   "s_aux",     :wz0,    12,    12,    12,    12 ],
 [ "s_gflux", "ν∇u[1]",    12,    12,    12,    12 ],
 [ "s_gflux", "ν∇u[2]",    12,    12,    12,    12 ],
 [ "s_gflux", "ν∇u[3]",    12,    12,    12,    12 ],
 [ "s_gflux", "ν∇u[4]",    12,    12,    12,    12 ],
 [ "s_gflux", "ν∇u[5]",    12,    12,    12,    12 ],
 [ "s_gflux", "ν∇u[6]",    12,    12,    12,    12 ],
 [ "s_gflux", "κ∇θ[1]",    12,    12,    12,    12 ],
 [ "s_gflux", "κ∇θ[2]",    12,    12,    12,    12 ],
 [ "s_gflux", "κ∇θ[3]",    12,     8,    12,    12 ],
]

    append!(refVals ,[ varr ] )
    append!(refPrecs,[ parr ] )

# SC ========== Test number 3 reference values and precision match template. =======
# SC ========== /Users/chrishill/projects/clima/cm/test/Ocean/HydrostaticBoussinesq/test_ocean_gyre.jl test reference values ======================================
# BEGIN SCPRINT
# varr - reference values (from reference run)
# parr - digits match precision (hand edit as needed)
#
# [
#  [ MPIStateArray Name, Field Name, Maximum, Minimum, Mean, Standard Deviation ],
#  [         :                :          :        :      :          :           ],
# ]
varr = [
 [       "Q",   "u[1]", -1.32150247347356725180e-01,  1.64031203423611798842e-01,  4.86341857180362744512e-03,  1.48290239846339700486e-02 ],
 [       "Q",   "u[2]", -6.36612598877540902809e-02,  7.44560153894480020176e-02, -2.47624053318063700596e-03,  1.31303104153666864301e-02 ],
 [       "Q",       :η, -6.17793649659350796455e-01,  6.06966232332265209770e-01, -1.29090689600396902117e-03,  2.25023233970694097605e-01 ],
 [       "Q",       :θ,  8.59363830767245827184e-05,  9.03903412174035025828e+00,  2.49947426680983575409e+00,  2.19708020497952460914e+00 ],
 [   "s_aux",       :y,  0.00000000000000000000e+00,  4.00000000000000046566e+06,  2.00000000000000000000e+06,  1.18025270281469495967e+06 ],
 [   "s_aux",       :w, -4.88109548294750840891e-05,  5.07779930498253729947e-05,  3.71427081618204492464e-07,  9.36953994676888197318e-06 ],
 [   "s_aux",    :pkin, -9.00072683420330310611e-01,  0.00000000000000000000e+00, -3.32075420474812310623e-01,  2.56210240101286457204e-01 ],
 [   "s_aux",     :wz0, -2.37163738698016553002e-05,  1.61967572887668093411e-05, -4.64606200633942812947e-08,  7.86254275820167091446e-06 ],
 [ "s_gflux", "ν∇u[1]", -6.67376825614485365279e-03,  6.97979211528434052153e-03,  1.01948640467980356077e-07,  2.29183131022463946721e-04 ],
 [ "s_gflux", "ν∇u[2]", -7.69310192750550202179e-03,  1.37445482375157909845e-03, -3.35376028298812290740e-05,  5.28895496282001859913e-04 ],
 [ "s_gflux", "ν∇u[3]", -1.70129300483239544585e-05,  1.30704917153533933284e-05, -4.30738980896587680255e-08,  1.59157632234111246367e-06 ],
 [ "s_gflux", "ν∇u[4]", -4.05738110420349878948e-03,  4.20907825808012042018e-03,  2.59008171395293793600e-06,  4.12110908407640325318e-04 ],
 [ "s_gflux", "ν∇u[5]", -4.15559862073769641783e-03,  1.11805092953295753984e-03, -7.36416840356856789053e-06,  2.39482266912490057102e-04 ],
 [ "s_gflux", "ν∇u[6]", -7.43240009691175527648e-06,  2.23076045923041977150e-05,  2.75166886765055948798e-08,  1.43297925696224871009e-06 ],
 [ "s_gflux", "κ∇θ[1]", -7.28126206584962977356e-05,  6.51042269918416155244e-05, -3.06828578003712164183e-06,  9.88453091533718813987e-06 ],
 [ "s_gflux", "κ∇θ[2]", -5.64043404057886982382e-05,  3.11699982926543850761e-03,  9.81313644046314260266e-04,  8.32013174688638057486e-04 ],
 [ "s_gflux", "κ∇θ[3]", -2.16059510376949300253e-05,  2.05358490375227276093e-05, -5.17225691641381943229e-07,  1.79077714232301003982e-06 ],
]
parr = [
 [       "Q",   "u[1]",    12,    12,    12,    12 ],
 [       "Q",   "u[2]",    12,    12,    12,    12 ],
 [       "Q",       :η,    12,    12,    12,    12 ],
 [       "Q",       :θ,    12,    12,    12,    12 ],
 [   "s_aux",       :y,    12,    12,    12,    12 ],
 [   "s_aux",       :w,    12,    12,    12,    12 ],
 [   "s_aux",    :pkin,    12,    12,    12,    12 ],
 [   "s_aux",     :wz0,    12,    12,    12,    12 ],
 [ "s_gflux", "ν∇u[1]",    12,    12,    12,    12 ],
 [ "s_gflux", "ν∇u[2]",    12,    12,    12,    12 ],
 [ "s_gflux", "ν∇u[3]",    12,    12,    12,    12 ],
 [ "s_gflux", "ν∇u[4]",    12,    12,    12,    12 ],
 [ "s_gflux", "ν∇u[5]",    12,    12,    12,    12 ],
 [ "s_gflux", "ν∇u[6]",    12,    12,    12,    12 ],
 [ "s_gflux", "κ∇θ[1]",    12,    12,    12,    12 ],
 [ "s_gflux", "κ∇θ[2]",    12,    12,    12,    12 ],
 [ "s_gflux", "κ∇θ[3]",    12,    12,    12,    12 ],
]
# END SCPRINT
# SC ====================================================================================

    append!(refVals ,[ varr ] )
    append!(refPrecs,[ parr ] )

# SC ========== Test number 4 reference values and precision match template. =======
# SC ========== /Users/chrishill/projects/clima/cm/test/Ocean/HydrostaticBoussinesq/test_ocean_gyre.jl test reference values ======================================
# BEGIN SCPRINT
# varr - reference values (from reference run)
# parr - digits match precision (hand edit as needed)
#
# [
#  [ MPIStateArray Name, Field Name, Maximum, Minimum, Mean, Standard Deviation ],
#  [         :                :          :        :      :          :           ],
# ]
varr = [
 [       "Q",   "u[1]", -1.32522061560302495398e-01,  1.64403805704383537689e-01,  4.86426609935955867742e-03,  1.48297097518645653452e-02 ],
 [       "Q",   "u[2]", -6.36533081952946194759e-02,  7.44240286788056537581e-02, -2.47705914433662445495e-03,  1.31287947739917917195e-02 ],
 [       "Q",       :η, -6.17792146223299143415e-01,  6.06964210100974344009e-01, -1.29086897946112909297e-03,  2.25017541338622450997e-01 ],
 [       "Q",       :θ,  1.18930697909414831350e-04,  9.05193063653280205472e+00,  2.49947692790997599843e+00,  2.19719805926361999582e+00 ],
 [   "s_aux",       :y,  0.00000000000000000000e+00,  4.00000000000000046566e+06,  2.00000000000000000000e+06,  1.18025270281469495967e+06 ],
 [   "s_aux",       :w, -4.86355979078713029592e-05,  5.05698562144389179547e-05,  3.70131833873677610995e-07,  9.34929694591741612729e-06 ],
 [   "s_aux",    :pkin, -9.00073156073595503912e-01,  0.00000000000000000000e+00, -3.32075359710396145196e-01,  2.56209888633661342361e-01 ],
 [   "s_aux",     :wz0, -2.38160453952588820092e-05,  1.63295706147758099968e-05, -4.63837554347600732956e-08,  7.90783852649017554050e-06 ],
 [ "s_gflux", "ν∇u[1]", -6.68560102683574044441e-03,  6.98986179365855032908e-03,  1.00393045025810862916e-07,  2.29127846391889087837e-04 ],
 [ "s_gflux", "ν∇u[2]", -7.69995573676401749708e-03,  1.37839162179005680388e-03, -3.35981134949950310548e-05,  5.29965375899606985950e-04 ],
 [ "s_gflux", "ν∇u[3]", -1.70723021047607414786e-05,  1.30922704334956021486e-05, -4.44618327893026796269e-08,  1.59336218473792970699e-06 ],
 [ "s_gflux", "ν∇u[4]", -4.05946426934079233201e-03,  4.21052254828468888043e-03,  2.59778123025053437326e-06,  4.12726714325575502491e-04 ],
 [ "s_gflux", "ν∇u[5]", -4.14814071152546075955e-03,  1.11726677539041160600e-03, -7.34353368996733004402e-06,  2.38907876999768631733e-04 ],
 [ "s_gflux", "ν∇u[6]", -7.41889880459672764628e-06,  2.22738260107801488963e-05,  2.78343092988202756278e-08,  1.42972197860415487431e-06 ],
 [ "s_gflux", "κ∇θ[1]", -7.25681212396385411870e-05,  7.16011471902109551376e-05, -3.07008530876072845326e-06,  9.89983434920985542473e-06 ],
 [ "s_gflux", "κ∇θ[2]", -5.79706743991101401507e-05,  3.10853728745289804858e-03,  9.81905426072233198168e-04,  8.31329782146290126484e-04 ],
 [ "s_gflux", "κ∇θ[3]", -3.24307611623586362264e-06,  1.48205642064209931587e-06, -5.05122832861326802479e-07,  3.05070684012206338078e-07 ],
]
parr = [
 [       "Q",   "u[1]",    12,    12,    12,    12 ],
 [       "Q",   "u[2]",    12,    12,    12,    12 ],
 [       "Q",       :η,    12,    12,    12,    12 ],
 [       "Q",       :θ,    12,    12,    12,    12 ],
 [   "s_aux",       :y,    12,    12,    12,    12 ],
 [   "s_aux",       :w,    12,    12,    12,    12 ],
 [   "s_aux",    :pkin,    12,    12,    12,    12 ],
 [   "s_aux",     :wz0,    12,    12,    12,    12 ],
 [ "s_gflux", "ν∇u[1]",    12,    12,    12,    12 ],
 [ "s_gflux", "ν∇u[2]",    12,    12,    12,    12 ],
 [ "s_gflux", "ν∇u[3]",    12,    12,    12,    12 ],
 [ "s_gflux", "ν∇u[4]",    12,    12,    12,    12 ],
 [ "s_gflux", "ν∇u[5]",    12,    12,    12,    12 ],
 [ "s_gflux", "ν∇u[6]",    12,    12,    12,    12 ],
 [ "s_gflux", "κ∇θ[1]",    12,    12,    12,    12 ],
 [ "s_gflux", "κ∇θ[2]",    12,    12,    12,    12 ],
 [ "s_gflux", "κ∇θ[3]",    12,    12,    12,    12 ],
]
# END SCPRINT
# SC ====================================================================================

    append!(refVals ,[ varr ] )
    append!(refPrecs,[ parr ] )

#! format: on

max_stars_repo_name: leios/CLIMA

max_issues_repo_name: leios/CLIMA

max_forks_repo_name: leios/CLIMA


text:
function SiteOp(o::Op)
  return SiteOp(Ops.which_op(o), Ops.sites(o), Ops.params(o))
end

function MPOTerm(o::α{∏{Op}})
  return MPOTerm(coefficient(o), [SiteOp(oₙ) for oₙ in o])
end

function OpSum(o::∑{<:α{∏{Op}}})
  return OpSum([MPOTerm(oₙ) for oₙ in o])
end

function OpSum(o::Union{Op,Applied})
  return OpSum(∑{<:α{∏{Op}}}(o))
end

# Conversions from other formats
function MPO(o::Union{Op,Applied}, s::Vector{<:Index}; kwargs...)
  return MPO(OpSum(o), s; kwargs...)
end


text:
mpSim = ManpowerSimulation( "sim" )

# Setting Attributes.
attribute1 = Attribute( "level" )
setPossibleAttributeValues!( attribute1, ["Junior", "Senior", "Master"] )
setInitialAttributeValues!( attribute1, Dict( "Junior" => 1.0 ) )
attribute2 = Attribute( "branch" )
setPossibleAttributeValues!( attribute2, ["A", "B", "reserve", "none"] )
setInitialAttributeValues!( attribute2, Dict( "none" => 1.0 ) )
attribute3 = Attribute( "isCareer" )
setPossibleAttributeValues!( attribute3, ["yes", "no"] )
setInitialAttributeValues!( attribute3, Dict( "no" => 1.0 ) )
setSimulationAttributes!( mpSim, [attribute1, attribute2, attribute3] )

# Setting Base Nodes.
node1 = BaseNode( "A junior" )
setNodeRequirements!( node1, ("level", "Junior"), ("branch", "A"),
    ("isCareer", "no") )
node2 = BaseNode( "B junior" )
setNodeRequirements!( node2, ("level", "Junior"), ("branch", "B"),
    ("isCareer", "no") )
node3 = BaseNode( "Reserve junior" )
setNodeRequirements!( node3, ("level", "Junior"), ("branch", "reserve"),
    ("isCareer", "no") )
node4 = BaseNode( "A senior" )
setNodeRequirements!( node4, ("level", "Senior"), ("branch", "A"),
    ("isCareer", "yes") )
node5 = BaseNode( "B senior" )
setNodeRequirements!( node5, ("level", "Senior"), ("branch", "B"),
    ("isCareer", "yes") )
node6 = BaseNode( "Master" )
setNodeRequirements!( node6, ("level", "Master"), ("branch", "none"),
    ("isCareer", "yes") )
setNodeTarget!.( [node1, node2, node3, node4, node5, node6],
    [-1, -1, -1, 30, 30, 25] )
setSimulationBaseNodes!( mpSim, [node1, node2, node3, node4, node5, node6] )
setSimulationBaseNodeOrder!( mpSim, Dict( "Master" => 1, "A senior" => 2,
    "B senior" => 2, "A junior" => 3, "B junior" => 3 ) )

# Setting Compound Nodes.
cnode1 = CompoundNode( "Branch A" )
setCompoundNodeComponents!( cnode1, ["A junior", "A senior"] )
cnode2 = CompoundNode( "Branch B" )
setCompoundNodeComponents!( cnode2, ["B junior", "B senior"] )
cnode3 = CompoundNode( "Junior" )
setCompoundNodeComponents!( cnode3,
    ["A junior", "B junior", "Reserve junior"] )
cnode4 = CompoundNode( "Career" )
setCompoundNodeComponents!( cnode4, ["A senior", "B senior", "Master"] )
cnode5 = CompoundNode( "Empty" )
setSimulationCompoundNodes!( mpSim,
    [cnode1, cnode2, cnode3, cnode4 , cnode5] )

# Setting recruitment.
rec1 = Recruitment( "EW" )
setRecruitmentTarget!( rec1, "A junior" )
rec2 = Recruitment( "EW" )
setRecruitmentTarget!( rec2, "B junior" )
rec3 = Recruitment( "EW" )
setRecruitmentTarget!( rec3, "Reserve junior" )

setRecruitmentSchedule!.( [rec1, rec2, rec3], 12 )
setRecruitmentFixed!.( [rec1, rec2, rec3], 10 )
setRecruitmentAgeFixed!.( [rec1, rec2, rec3], 240 )
addSimulationRecruitment!( mpSim, rec1, rec2, rec3 )

# Setting THROUGH Transitions.
ttrans1 = Transition( "Promotion", "A junior", "A senior" )
addTransitionCondition!( ttrans1, MPcondition( "time in node", >=, 24 ) )
ttrans2 = Transition( "Promotion", "B junior", "B senior" )
addTransitionCondition!( ttrans2, MPcondition( "time in node", >=, 24 ) )
ttrans3 = Transition( "Promotion", "A senior", "Master" )
addTransitionCondition!( ttrans3, MPcondition( "time in node", >=, 36 ) )
ttrans4 = Transition( "Promotion", "B senior", "Master" )
addTransitionCondition!( ttrans4, MPcondition( "time in node", >=, 36 ) )
ttrans5 = Transition( "Reserve", "Reserve junior", "A senior" )
addTransitionCondition!( ttrans5, MPcondition( "time in node", >=, 24 ) )
ttrans6 = Transition( "Reserve", "Reserve junior", "B senior" )
addTransitionCondition!( ttrans6, MPcondition( "time in node", >=, 24 ) )

setTransitionSchedule!.( [ttrans1, ttrans2, ttrans3, ttrans4, ttrans5,
    ttrans6], 12 )
setTransitionMaxAttempts!.( [ttrans1, ttrans2, ttrans3, ttrans4, ttrans5,
    ttrans6], 1 )
setTransitionFluxLimits!.( [ttrans1, ttrans2, ttrans3, ttrans4, ttrans5,
    ttrans6], [8, 8, 0, 0, 0, 0], [8, 8, 5, 5, 4, 4] )
addSimulationTransition!( mpSim, ttrans1, ttrans2, ttrans3, ttrans4, ttrans5,
    ttrans6 )

# Setting OUT Transitions.
otrans1 = Transition( "B-", "A junior" )
addTransitionCondition!( otrans1, MPcondition( "time in node", >=, 24 ) )
otrans2 = Transition( "B-", "B junior" )
addTransitionCondition!( otrans2, MPcondition( "time in node", >=, 24 ) )
otrans3 = Transition( "B-", "Reserve junior" )
addTransitionCondition!( otrans3, MPcondition( "time in node", >=, 24 ) )
otrans4 = Transition( "PE", "A senior" )
addTransitionCondition!( otrans4, MPcondition( "tenure", >=, 120 ) )
otrans5 = Transition( "PE", "B senior" )
addTransitionCondition!( otrans5, MPcondition( "tenure", >=, 120 ) )
otrans6 = Transition( "PE", "Master" )
addTransitionCondition!( otrans6, MPcondition( "tenure", >=, 120 ) )

setTransitionMaxAttempts!.( [otrans1, otrans2, otrans3, otrans4, otrans5,
    otrans6], 1 )
setTransitionFluxLimits!.( [otrans1, otrans2, otrans3, otrans4, otrans5,
    otrans6], 0, -1 )
addSimulationTransition!( mpSim, otrans1, otrans2, otrans3, otrans4, otrans5,
    otrans6 )
setSimulationTransitionTypeOrder!( mpSim, Dict( "PE" => 1, "Promotion" => 2,
    "Reserve" => 3, "B-" => 4 ) )

# Miscellaneous configs.
setSimulationLength!( mpSim, 300 )

@test verifySimulation!( mpSim )

@testset "function saveSimulationConfiguration" begin
    setSimulationDatabase!( mpSim, joinpath( "reports", "simDBcopy" ) )
    saveSimulationConfiguration( mpSim )
    setSimulationDatabase!( mpSim, joinpath( "reports", "simDB" ) )
end  # @test "function saveSimulationConfiguration"