text:
using Base.Threads

using ArgParse
using CUDA
using Dates
using LinearAlgebra
using Logging
using MPI
using Printf
using Random

using CLIMAParameters

using ..Atmos
using ..Callbacks
using ..Checkpoint
using ..SystemSolvers
using ..ConfigTypes
using ..Diagnostics
using ..DiagnosticsMachine
using ..DGMethods
using ..BalanceLaws
using ..DGMethods: remainder_DGModel, SpaceDiscretization
using ..DGMethods.NumericalFluxes
using ..DGMethods.FVReconstructions

using ..Mesh.Grids
using ..Mesh.Topologies
using ..Mesh.Filters
using Thermodynamics
using ..MPIStateArrays
using ..ODESolvers
using ..TicToc
using ..VariableTemplates
using ..VTK

function _init_array(::Type{CuArray})
    comm = MPI.COMM_WORLD
    # allocate GPUs among MPI ranks
    local_comm =
        MPI.Comm_split_type(comm, MPI.MPI_COMM_TYPE_SHARED, MPI.Comm_rank(comm))
    # we intentionally oversubscribe GPUs for testing: may want to disable this for production
    CUDA.device!(MPI.Comm_rank(local_comm) % length(devices()))
    CUDA.allowscalar(false)
    return nothing
end

_init_array(::Type{Array}) = nothing

function gpu_allowscalar(val::Bool)
    CUDA.allowscalar(val)
    return
end

# Note that the initial values specified here are overwritten by the
# command line argument defaults in `parse_commandline()`.
Base.@kwdef mutable struct ClimateMachine_Settings
    disable_gpu::Bool = false
    show_updates::String = "60secs"
    start_datetime::DateTime = DateTime(2000, 1, 1, 12)
    diagnostics::String = "never"
    no_overwrite::Bool = false
    vtk::String = "never"
    vtk_number_sample_points::Int = 0
    monitor_timestep_duration::String = "never"
    monitor_courant_numbers::String = "never"
    adapt_timestep::String = "never"
    checkpoint::String = "never"
    checkpoint_keep_one::Bool = true
    checkpoint_at_end::Bool = false
    checkpoint_on_crash::Bool = false
    checkpoint_dir::String = "checkpoint"
    restart_from_num::Int = -1
    fix_rng_seed::Bool = false
    log_level::String = "INFO"
    disable_custom_logger::Bool = false
    output_dir::String = "output"
    debug_init::Bool = false
    integration_testing::Bool = false
    array_type::Type = Array
    sim_time::Float64 = NaN
    fixed_number_of_steps::Int = -1
    degree::NTuple{2, Int} = (-1, -1)
    cutoff_degree::NTuple{2, Int} = (-1, -1)
    nelems::NTuple{3, Int} = (-1, -1, -1)
    domain_height::Float64 = -1
    resolution::NTuple{3, Float64} = (-1, -1, -1)
    domain_min::NTuple{3, Float64} = (-1, -1, -1)
    domain_max::NTuple{3, Float64} = (-1, -1, -1)
end

const Settings = ClimateMachine_Settings()

"""
    ClimateMachine.array_type()

Return the array type used by ClimateMachine. This defaults to (CPU-based)
`Array` and is only correctly set (based on choice from the command
line, from an environment variable, or from experiment code) after
`ClimateMachine.init()` is called.
"""
array_type() = Settings.array_type

"""
    ClimateMachine.datetime(solver::AbstractODESolver)

Return the current simulation date and time.
"""
function datetime(solver::AbstractODESolver)
    simtime = gettime(solver)
    return Settings.start_datetime + Second(round(Int, simtime))
end

"""
    get_setting(setting_name::Symbol, settings, defaults)

Define fallback behavior for driver settings, first accessing overloaded
`settings` if defined, followed by constructed global ENV variable
`CLIMATEMACHINE_SETTINGS_<SETTING_NAME>`, then (global) defaults.

Returns setting value.
"""
function get_setting(setting_name::Symbol, settings, defaults)
    if !haskey(defaults, setting_name)
        error("setting $setting_name is not defined in `defaults`")
    end
    setting_type = typeof(defaults[setting_name])
    setting_env = "CLIMATEMACHINE_SETTINGS_" * uppercase(String(setting_name))
    if haskey(settings, setting_name)
        return convert(setting_type, settings[setting_name])
    elseif haskey(ENV, setting_env)
        env_val = ENV[setting_env]
        if setting_type == String
            v = env_val
        elseif setting_type == NTuple{2, Int} ||
               setting_type == NTuple{3, Float64} ||
               setting_type == DateTime
            v = ArgParse.parse_item(setting_type, env_val)
        else
            v = tryparse(setting_type, env_val)
        end
        if isnothing(v)
            error("cannot parse ENV $setting_env value $env_val, to setting type $setting_type")
        end
        return v
    elseif haskey(defaults, setting_name)
        return defaults[setting_name]
    else
        error("setting $setting_name is not contained in either settings or defaults")
    end
end

function get_gpu_setting(setting_name::Symbol, settings, defaults)
    # do not override disable_gpu keyword argument setting if it exists
    if !haskey(settings, setting_name)
        # if old GPU ENV keyword exists, overwrite the settings variable if not defined
        if haskey(ENV, "CLIMATEMACHINE_GPU")
            settings[setting_name] = ENV["CLIMATEMACHINE_GPU"] == "false"
        end
    end
    # fallback behavior
    return get_setting(setting_name, settings, defaults)
end

"""
    parse_commandline(
        defaults::Dict{Symbol, Any},
        global_defaults::Dict{Symbol, Any},
        custom_clargs::Union{Nothing, ArgParseSettings} = nothing,
    )

Parse process command line ARGS values. If `defaults` is specified, it
overrides default values for command line argument defaults. If
`custom_clargs` is specified, it is added to the parsing step.

Returns a `Dict` containing parsed process ARGS values.
"""
function parse_commandline(
    defaults::Dict{Symbol, Any},
    global_defaults::Dict{Symbol, Any},
    custom_clargs::Union{Nothing, ArgParseSettings} = nothing,
)
    exc_handler = ArgParse.default_handler
    if Base.isinteractive()
        exc_handler = ArgParse.debug_handler
    end
    s = ArgParseSettings(
        prog = PROGRAM_FILE,
        description = "Climate Machine: an Earth System Model that automatically learns from data\n",
        preformatted_description = true,
        epilog = """
            Any <interval> unless otherwise stated may be specified as:
                - 2hours or 10mins or 30secs => wall-clock time
                - 9.5smonths or 3.3sdays or 1.5shours => simulation time
                - 1000steps => simulation steps
                - never => disable
                - default => use experiment specified interval (only for diagnostics at present)
            """,
        preformatted_epilog = true,
        version = string(CLIMATEMACHINE_VERSION),
        exc_handler = exc_handler,
        autofix_names = true,     # switches --flag-name to 'flag_name'
        error_on_conflict = true, # don't allow custom_clargs' settings to override these
    )
    add_arg_group!(s, "ClimateMachine")

    @add_arg_table! s begin
        "--disable-gpu"
        help = "do not use the GPU"
        action = :store_const
        constant = true
        default = get_gpu_setting(:disable_gpu, defaults, global_defaults)
        "--show-updates"
        help = "interval at which to show simulation updates"
        metavar = "<interval>"
        arg_type = String
        default = get_setting(:show_updates, defaults, global_defaults)
        "--start-datetime"
        help = "specify simulation start date/time as 'yyyy-mm-ddTHH:MM:SS'"
        metavar = "<date/time>"
        arg_type = DateTime
        default = get_setting(:start_datetime, defaults, global_defaults)
        "--diagnostics"
        help = "interval at which to collect diagnostics"
        metavar = "<interval>"
        arg_type = String
        default = get_setting(:diagnostics, defaults, global_defaults)
        "--no-overwrite"
        help = "throw an error if an output file would be overwritten"
        action = :store_const
        constant = true
        default = get_setting(:no_overwrite, defaults, global_defaults)
        "--vtk"
        help = "interval at which to output VTK"
        metavar = "<interval>"
        arg_type = String
        default = get_setting(:vtk, defaults, global_defaults)
        "--vtk-number-sample-points"
        help = "number of sampling points in each element for VTK output"
        metavar = "<number>"
        arg_type = Int
        default =
            get_setting(:vtk_number_sample_points, defaults, global_defaults)
        "--monitor-timestep-duration"
        help = "interval in time-steps at which to output wall-clock time per time-step"
        metavar = "<interval>"
        arg_type = String
        default =
            get_setting(:monitor_timestep_duration, defaults, global_defaults)
        "--monitor-courant-numbers"
        help = "interval at which to output acoustic, advective, and diffusive Courant numbers"
        metavar = "<interval>"
        arg_type = String
        default =
            get_setting(:monitor_courant_numbers, defaults, global_defaults)
        "--adapt-timestep"
        help = "interval at which to update the timestep"
        metavar = "<interval>"
        arg_type = String
        default = get_setting(:adapt_timestep, defaults, global_defaults)
        "--checkpoint"
        help = "interval at which to create a checkpoint"
        metavar = "<interval>"
        arg_type = String
        default = get_setting(:checkpoint, defaults, global_defaults)
        "--checkpoint-keep-all"
        help = "keep all checkpoints (instead of just the most recent)"
        action = :store_const
        constant = true
        default = !get_setting(:checkpoint_keep_one, defaults, global_defaults)
        "--checkpoint-at-end"
        help = "create a checkpoint at the end of the simulation"
        action = :store_const
        constant = true
        default = get_setting(:checkpoint_at_end, defaults, global_defaults)
        "--checkpoint-on-crash"
        help = "create a checkpoint on a kernel crash (hurts performance!)"
        action = :store_const
        constant = true
        default = get_setting(:checkpoint_on_crash, defaults, global_defaults)
        "--checkpoint-dir"
        help = "the directory in which to store checkpoints"
        metavar = "<path>"
        arg_type = String
        default = get_setting(:checkpoint_dir, defaults, global_defaults)
        "--restart-from-num"
        help = "checkpoint number from which to restart (in <checkpoint-dir>)"
        metavar = "<number>"
        arg_type = Int
        default = get_setting(:restart_from_num, defaults, global_defaults)
        "--fix-rng-seed"
        help = "set RNG seed to a fixed value for reproducibility"
        action = :store_const
        constant = true
        default = get_setting(:fix_rng_seed, defaults, global_defaults)
        "--disable-custom-logger"
        help = "do not use a custom logger"
        action = :store_const
        constant = true
        default = get_setting(:disable_custom_logger, defaults, global_defaults)
        "--log-level"
        help = "set the log level to one of debug/info/warn/error"
        metavar = "<level>"
        arg_type = String
        default = uppercase(get_setting(:log_level, defaults, global_defaults))
        "--output-dir"
        help = "directory for output data"
        metavar = "<path>"
        arg_type = String
        default = get(defaults, :output_dir) do
            get(ENV, "CLIMATEMACHINE_OUTPUT_DIR") do
                get(ENV, "CLIMATEMACHINE_SETTINGS_OUTPUT_DIR", Settings.output_dir)
            end
        end
        "--debug-init"
        help = "fill state arrays with NaNs and dump them post-initialization"
        action = :store_const
        constant = true
        default = get_setting(:debug_init, defaults, global_defaults)
        "--integration-testing"
        help = "enable integration testing"
        action = :store_const
        constant = true
        default = get_setting(:integration_testing, defaults, global_defaults)
        "--sim-time"
        help = "run for the specified time (in simulation seconds)"
        metavar = "<number>"
        arg_type = Float64
        default = get_setting(:sim_time, defaults, global_defaults)
        "--fixed-number-of-steps"
        help = "if `≥0` perform specified number of steps"
        metavar = "<number>"
        arg_type = Int
        default = get_setting(:fixed_number_of_steps, defaults, global_defaults)
        "--degree"
        help = "tuple of horizontal and vertical polynomial degrees for spatial discretization order (no space before/after comma)"
        metavar = "<horizontal>,<vertical>"
        arg_type = NTuple{2, Int}
        default = get_setting(:degree, defaults, global_defaults)
        "--cutoff-degree"
        help = "tuple of horizontal and vertical polynomial degrees to keep when applying a cutoff filter (no space before/after comma)"
        metavar = "<horizontal>,<vertical>"
        arg_type = NTuple{2, Int}
        default = get_setting(:cutoff_degree, defaults, global_defaults)
        "--nelems"
        help = "number of elements in each direction: 3 for Ocean GCM, 2 for Atmos GCM or 1 for Atmos single-stack (no space before/after comma)"
        metavar = "<nelem_1>[,<nelem_2>[,<nelem_3>]]"
        arg_type = NTuple{3, Int}
        default = get_setting(:nelems, defaults, global_defaults)
        "--domain-height"
        help = "domain height (in meters) for GCM or single-stack configurations"
        metavar = "<number>"
        arg_type = Float64
        default = get_setting(:domain_height, defaults, global_defaults)
        "--resolution"
        help = "tuple of three element resolutions (in meters) for LES and MultiColumnLandModel configurations"
        metavar = "<Δx>,<Δy>,<Δz>"
        arg_type = NTuple{3, Float64}
        default = get_setting(:resolution, defaults, global_defaults)
        "--domain-min"
        help = "tuple of three minima for the domain size (in meters) for LES and MultiColumnLandModel configurations"
        metavar = "<xmin>,<ymin>,<zmin>"
        arg_type = NTuple{3, Float64}
        default = get_setting(:domain_min, defaults, global_defaults)
        "--domain-max"
        help = "tuple of three maxima for the domain size (in meters) for LES and MultiColumnLandModel configurations"
        metavar = "<xmax>,<ymax>,<zmax>"
        arg_type = NTuple{3, Float64}
        default = get_setting(:domain_max, defaults, global_defaults)
    end
    # add custom cli argparse settings if provided
    if !isnothing(custom_clargs)
        import_settings!(s, custom_clargs)
    end
    return parse_args(s)
end

"""
    ClimateMachine.init(;
        parse_clargs::Bool = false,
        custom_clargs::Union{Nothing, ArgParseSettings} = nothing,
        init_driver::Bool = true,
        keyword_args...,
    )

Initialize the ClimateMachine. If `parse_clargs` is set, parse command line
arguments (additional driver-specific arguments can be added by specifying
`custom_clargs`).

Setting `init_driver = false` will set up the `ClimateMachine.Settings`
singleton values without initializing the ClimateMachine runtime. Otherwise,
the runtime will be initialized (see `init_runtime()`).

Finally, key-value pairs can be supplied to `ClimateMachine.init()` to set
system default settings -- the final settings are decided as follows (in
order of precedence):
1. Command line arguments (if `parse_clargs = true`).
2. Environment variables.
3. Keyword arguments to `init()`.
4. Defaults (in `ClimateMachine_Settings`).

Recognized keyword arguments are:
- `disable_gpu::Bool = false`:
        do not use the GPU
- `show_updates::String = "60secs"`:
        interval at which to show simulation updates
- `start_datetime::DateTime = DateTime(2000, 1, 1, 12)`:
        date/time at which the simulation starts
- `diagnostics::String = "never"`:
        interval at which to collect diagnostics"
- `no_overwrite::Bool = false`:
        throw an error if an output file would be overwritten
- `vtk::String = "never"`:
        interval at which to write simulation vtk output
- `vtk-number-sample-points::Int` = 0:
        the number of sampling points in each element for VTK output
- `monitor_timestep_duration::String = "never"`:
        interval in time-steps at which to output wall-clock time per time-step
- `monitor_courant_numbers::String = "never"`:
        interval at which to output acoustic, advective, and diffusive Courant numbers
- `adapt-timestep::String = "never"`:
        interval at which to update the timestep
- `checkpoint::String = "never"`:
        interval at which to write a checkpoint
- `checkpoint_keep_one::Bool = true`: (interval)
        keep all checkpoints (instead of just the most recent)
- `checkpoint_at_end::Bool = false`:
        create a checkpoint at the end of the simulation
- `checkpoint_on_crash::Bool = false`:
    create a checkpoint on a kernel crash (hurts performance!)
- `checkpoint_dir::String = "checkpoint"`:
        absolute or relative path to checkpoint directory
- `restart_from_num::Int = -1`:
        checkpoint number from which to restart (in `checkpoint_dir`)
- `fix_rng_seed::Bool = false`:
        set RNG seed to a fixed value for reproducibility
- `log_level::String = "INFO"`:
        log level for ClimateMachine global default runtime logger
- `disable_custom_logger::String = false`:
        disable using a global custom logger for ClimateMachine
- `output_dir::String = "output"`: (path)
        absolute or relative path to output data directory
- `debug_init::Bool = false`:
        fill state arrays with NaNs and dump them post-initialization
- `integration_testing::Bool = false`:
        enable integration_testing
- `sim_time::Float64 = NaN`:
        run for the specified time (in simulation seconds)
- `fixed_number_of_steps::Int = -1`:
        if `≥0` perform specified number of steps
- `degree::NTuple{2, Int} = (-1, -1)`:
        tuple of horizontal and vertical polynomial degrees for spatial discretization order
- `cutoff_degree::NTuple{2, Int} = (-1, -1)`:
        tuple of horizontal and vertical polynomial degrees for cutoff filter
- `nelems::NTuple{3, Int} = (-1, -1, -1)`:
        tuple of number of elements in each direction: 3 for Ocean, 2 for GCM or 1 for single-stack
- `domain_height::Float64 = -1`:
        domain height (in meters) for GCM or single-stack configurations
- `resolution::NTuple{3, Float64} = (-1, -1, -1)`:
        tuple of three element resolutions (in meters) for LES and MultiColumnLandModel configurations
- `domain_min::NTuple{3, Float64} = (-1, -1, -1)`:
        tuple of three minima for the domain size (in meters) for LES and MultiColumnLandModel configurations
- `domain_max::NTuple{3, Float64} = (-1, -1, -1)`:
        tuple of three maxima for the domain size (in meters) for LES and MultiColumnLandModel configurations

Returns `nothing`, or if `parse_clargs = true`, returns parsed command line
arguments.
"""
function init(;
    parse_clargs::Bool = false,
    custom_clargs::Union{Nothing, ArgParseSettings} = nothing,
    init_driver::Bool = true,
    keyword_args...,
)
    # `Settings` contains global defaults
    global_defaults = Dict{Symbol, Any}(
        (n, getproperty(Settings, n)) for n in propertynames(Settings)
    )

    # keyword arguments must be applicable to `Settings`
    all_args = Dict{Symbol, Any}(keyword_args)
    for kwarg in keys(all_args)
        if get(global_defaults, kwarg, nothing) === nothing
            throw(ArgumentError(string(kwarg)))
        end
    end

    # if command line arguments should be processed, do so and override
    # keyword arguments
    cl_args = nothing
    if parse_clargs
        cl_args = parse_commandline(all_args, global_defaults, custom_clargs)

        # We need to munge the parsed arg dict a bit as parsed arg keys
        # and initialization keywords are not 1:1
        cl_args["checkpoint_keep_one"] = !cl_args["checkpoint_keep_all"]

        # parsed command line arguments override keyword arguments
        all_args = merge(
            all_args,
            Dict{Symbol, Any}((Symbol(k), v) for (k, v) in cl_args),
        )
    end

    # TODO: also add validation for initialization values

    # Here, `all_args` contains command line arguments and keyword arguments.
    # They must be applied to `Settings`.
    #
    # special cases for backward compatibility
    if haskey(all_args, :disable_gpu)
        Settings.disable_gpu = all_args[:disable_gpu]
    elseif haskey(ENV, "CLIMATEMACHINE_GPU")
        @warn(
            "CLIMATEMACHINE_GPU will be deprecated; " *
            "use CLIMATEMACHINE_SETTINGS_DISABLE_GPU"
        )
        Settings.disable_gpu = ENV["CLIMATEMACHINE_GPU"] == "false"
    elseif haskey(ENV, "CLIMATEMACHINE_SETTINGS_DISABLE_GPU")
        Settings.disable_gpu =
            parse(Bool, ENV["CLIMATEMACHINE_SETTINGS_DISABLE_GPU"])
    end

    if haskey(all_args, :output_dir)
        Settings.output_dir = all_args[:output_dir]
    elseif haskey(ENV, "CLIMATEMACHINE_OUTPUT_DIR")
        @warn(
            "CLIMATEMACHINE_OUTPUT_DIR will be deprecated; " *
            "use CLIMATEMACHINE_SETTINGS_OUTPUT_DIR"
        )
        Settings.output_dir = ENV["CLIMATEMACHINE_OUTPUT_DIR"]
    elseif haskey(ENV, "CLIMATEMACHINE_SETTINGS_OUTPUT_DIR")
        Settings.output_dir = ENV["CLIMATEMACHINE_SETTINGS_OUTPUT_DIR"]
    end

    # all other settings
    for n in propertynames(Settings)
        # skip over the special backwards compat cases defined above
        if n == :disable_gpu || n == :output_dir
            continue
        end
        setproperty!(Settings, n, get_setting(n, all_args, global_defaults))
    end

    # set up the array type appropriately depending on whether we're using GPUs
    if !Settings.disable_gpu && CUDA.has_cuda_gpu()
        Settings.array_type = CUDA.CuArray
    else
        Settings.array_type = Array
    end

    # initialize the runtime
    if init_driver
        init_runtime(Settings)
    end
    return cl_args
end

"""
    init_runtime(settings::ClimateMachine_Settings)

Initialize the ClimateMachine runtime: initialize MPI, set the default
array type, set RNG seeds for all threads, create output directory and
set up logging.
"""
function init_runtime(settings::ClimateMachine_Settings)
    # set up timing mechanism
    tictoc()

    # initialize MPI
    if !MPI.Initialized()
        MPI.Init()
    end

    # initialize the array GPU backend if appropriate
    _init_array(settings.array_type)

    # fix the RNG seeds if requested
    if settings.fix_rng_seed
        rank = MPI.Comm_rank(MPI.COMM_WORLD)
        for tid in 1:nthreads()
            s = 1000 * rank + tid
            Random.seed!(Random.default_rng(tid), s)
        end
    end

    # TODO: write a better MPI logging back-end and also integrate Dlog
    # for large scale

    # set up logging
    log_level_str = uppercase(settings.log_level)
    loglevel =
        log_level_str == "DEBUG" ? Logging.Debug :
        log_level_str == "WARN" ? Logging.Warn :
        log_level_str == "ERROR" ? Logging.Error : Logging.Info
    if !settings.disable_custom_logger
        # cannot use `NullLogger` here because MPI collectives may be
        # used in logging calls!
        logger_stream = MPI.Comm_rank(MPI.COMM_WORLD) == 0 ? stderr : devnull
        prev_logger = global_logger(ConsoleLogger(logger_stream, loglevel))
        atexit() do
            global_logger(prev_logger)
        end
    end
    return
end

include("driver_configs.jl")
include("solver_configs.jl")
include("diagnostics_configs.jl")

"""
    ClimateMachine.ConservationCheck

Pass a tuple of these to `ClimateMachine.invoke!` to perform a
conservation check of each `varname` at the specified `interval`. This
computes `Σv = weightedsum(Q.varname)` and `δv = (Σv - Σv₀) / Σv`.
`invoke!` throws an error if `abs(δv)` exceeds `error_threshold. If
`show`, `δv` is displayed.
"""
struct ConservationCheck{FT}
    varname::String
    interval::String
    error_threshold::FT
    show::Bool
end
ConservationCheck(varname::String, interval::String) =
    ConservationCheck(varname, interval, Inf, true)
ConservationCheck(
    varname::String,
    interval::String,
    error_threshold::FT,
) where {FT} = ConservationCheck(varname, interval, error_threshold, true)

"""
    ClimateMachine.invoke!(
        solver_config::SolverConfiguration;
        adjustfinalstep = false,
        diagnostics_config = nothing,
        user_callbacks = (),
        user_info_callback = () -> nothing,
        check_cons = (),
        check_euclidean_distance = false,
    )

Run the simulation defined by `solver_config`.

Keyword Arguments:

The value of 'adjustfinalstep` is passed to the ODE solver; see
[`solve!`](@ref ODESolvers.solve!).

The `user_callbacks` are passed to the ODE solver as callback functions;
see [`solve!`](@ref ODESolvers.solve!).

The function `user_info_callback` is called after the default info
callback (which is called every `Settings.show_updates` interval). The
single input argument `init` is `true` when the callback is called for
initialization (before time stepping begins) and `false` when called
during the actual ODE solve; see [`GenericCallbacks`](@ref) and
[`solve!`](@ref ODESolvers.solve!).

If conservation checks are to be performed, `check_cons` must be a
tuple of [`ConservationCheck`](@ref).

If `check_euclidean_distance` is `true, then the Euclidean distance
between the final solution and initial condition function evaluated with
`solver_config.timeend` is reported.
"""
function invoke!(
    solver_config::SolverConfiguration;
    adjustfinalstep = false,
    diagnostics_config = nothing,
    user_callbacks = (),
    user_info_callback = () -> nothing,
    check_cons = (),
    check_euclidean_distance = false,
)
    mpicomm = solver_config.mpicomm
    dg = solver_config.dg
    bl = dg.balance_law
    Q = solver_config.Q
    FT = eltype(Q)
    timeend = solver_config.timeend
    init_on_cpu = solver_config.init_on_cpu
    init_args = solver_config.init_args
    solver = solver_config.solver

    # create the output directories if needed on delegated rank
    if MPI.Comm_rank(MPI.COMM_WORLD) == 0
        if Settings.diagnostics != "never" || Settings.vtk != "never"
            mkpath(Settings.output_dir)
        end
        if Settings.checkpoint != "never" ||
           Settings.checkpoint_at_end ||
           Settings.checkpoint_on_crash
            mkpath(Settings.checkpoint_dir)
        end
    end
    MPI.Barrier(MPI.COMM_WORLD)

    # set up callbacks
    callbacks = ()

    # info callback
    cb_updates = Callbacks.show_updates(
        Settings.show_updates,
        solver_config,
        user_info_callback,
    )
    if !isnothing(cb_updates)
        callbacks = (callbacks..., cb_updates)
    end

    # diagnostics callback(s)
    if Settings.diagnostics != "never" && !isnothing(diagnostics_config)
        dgn_starttime = replace(string(now()), ":" => ".")
        Diagnostics.init(
            mpicomm,
            solver_config.param_set,
            dg,
            Q,
            dgn_starttime,
            Settings.output_dir,
            Settings.no_overwrite,
        )
        DiagnosticsMachine.init(
            mpicomm,
            solver_config.param_set,
            dg,
            Q,
            dgn_starttime,
            Settings.output_dir,
            Settings.no_overwrite,
        )

        dgncbs = Callbacks.diagnostics(
            Settings.diagnostics,
            solver_config,
            dgn_starttime,
            diagnostics_config,
        )
        callbacks = (callbacks..., dgncbs...)
    end

    # vtk callback
    cb_vtk = Callbacks.vtk(
        Settings.vtk,
        solver_config,
        Settings.output_dir,
        Settings.vtk_number_sample_points,
    )
    if !isnothing(cb_vtk)
        callbacks = (callbacks..., cb_vtk)
    end

    # timestep duration monitor
    cb_mtd = Callbacks.monitor_timestep_duration(
        Settings.monitor_timestep_duration,
        Settings.array_type,
        mpicomm,
    )
    if !isnothing(cb_mtd)
        callbacks = (callbacks..., cb_mtd)
    end

    # Courant number monitor
    cb_mcn = Callbacks.monitor_courant_numbers(
        Settings.monitor_courant_numbers,
        solver_config,
    )
    if !isnothing(cb_mcn)
        callbacks = (callbacks..., cb_mcn)
    end

    # Timestep adapter
    cb_adp = Callbacks.adapt_timestep(Settings.adapt_timestep, solver_config)
    if !isnothing(cb_adp)
        callbacks = (callbacks..., cb_adp)
    end

    # checkpointing callback
    cb_checkpoint = Callbacks.checkpoint(
        Settings.checkpoint,
        Settings.checkpoint_keep_one,
        solver_config,
        Settings.checkpoint_dir,
    )
    if !isnothing(cb_checkpoint)
        callbacks = (callbacks..., cb_checkpoint)
    end

    # conservation callbacks
    cccbs = Callbacks.check_cons(check_cons, solver_config)
    callbacks = (callbacks..., cccbs...)

    # user callbacks
    callbacks = (user_callbacks..., callbacks...)

    # initial condition norm
    eng0 = norm(Q)
    @info @sprintf(
        """
        %s %s
            dt              = %.5e
            timeend         = %8.2f
            number of steps = %d
            norm(Q)         = %.16e""",
        Settings.restart_from_num > 0 ? "Restarting" : "Starting",
        solver_config.name,
        solver_config.dt,
        solver_config.timeend,
        solver_config.numberofsteps,
        eng0
    )

    # run the simulation
    try
        @tic solve!
        solve!(
            Q,
            solver;
            timeend = timeend,
            callbacks = callbacks,
            adjustfinalstep = adjustfinalstep,
        )
        @toc solve!
    catch exc
        if Settings.checkpoint_on_crash
            Callbacks.write_checkpoint(
                solver_config,
                Settings.checkpoint_dir,
                solver_config.name,
                mpicomm,
                solver_config.numberofsteps,
            )
        end
        rethrow()
    end

    # write end checkpoint if requested
    if Settings.checkpoint_at_end
        Callbacks.write_checkpoint(
            solver_config,
            Settings.checkpoint_dir,
            solver_config.name,
            mpicomm,
            solver_config.numberofsteps,
        )
    end

    engf = norm(Q)
    @info @sprintf(
        """
        Finished
            norm(Q)            = %.16e
            norm(Q) / norm(Q₀) = %.16e
            norm(Q) - norm(Q₀) = %.16e""",
        engf,
        engf / eng0,
        engf - eng0
    )

    if check_euclidean_distance
        Qe =
            init_ode_state(dg, timeend, init_args...; init_on_cpu = init_on_cpu)
        engfe = norm(Qe)
        errf = euclidean_distance(Q, Qe)
        @info @sprintf(
            """
            Euclidean distance
                norm(Q - Qe)            = %.16e
                norm(Q - Qe) / norm(Qe) = %.16e""",
            errf,
            errf / engfe
        )
    end

    return engf / eng0
end


text:
import Dance

include("./utils/main.jl")

Dance.start_project("demo")
project_settings()
routes("static")

# Create file for static file route
mkdir("files")
cp("../sample/static/shipping_containers.jpg", "files/image.jpg")

# Ensure route paths cannot go higher than project route in directory structure
@test_logs (:error, "Route paths cannot go higher than project route in directory structure") include(joinpath(abspath(pwd()), "routes.jl"))

lines_array = []
open("routes.jl", "r") do io
    global lines_array = readlines(io)
    for (idx, line) in enumerate(lines_array)
        if line=="route(\"/../files/image.jpg\", output_file_as_string; method=GET, endpoint=STATIC)"
            global lines_array
            deleteat!(lines_array, idx)
        end
    end
end
lines_string = ""
for line in lines_array
    global lines_string
    lines_string = lines_string * line * "\n"
end
open("routes.jl", "w") do io
    write(io, lines_string)
end
Dance.Router.delete_routes!()

## Test all routes with pending slash, to ensure is removed ##
@testset "HTTP.listen" begin
    @async include(joinpath(abspath(pwd()), "dance.jl"))
    sleep(1)

    # Favicon
    make_and_test_request_get("/favicon.ico/", 200, Dict("Content-Type" => "image/x-icon"), 15406, false, read("../../files/html/favicon.ico"))

    # Single static file
    make_and_test_request_get("/files/image.jpg/", 200, Dict("Content-Type" => "image/jpeg"), 84668, false, read("files/image.jpg"))

    # Static dir
    for item in [
        ["shipping_containers.jpg", 84668],
        ["office/coffee_and_laptop.jpg", 2126200],
        ["outdoors/scenary/hills_in_mist.jpg", 3181705],
        ["outdoors/wildlife/goat_and_sheep.jpg", 3175701]
    ]
        path::String = item[1]
        content_length::Int64 = item[2]
        make_and_test_request_get("/static/$path/", 200, Dict("Content-Type" => "image/jpeg"), content_length, false, read(joinpath("../sample/static", path)))
    end
end

delete_project()

ext: jl


text:
round_random(x)::Int = round(x) + (rand() < mod(x, 1))
# d6(n) = floor.(Int, rand(n) .* 6) .+ 1

function fire(;evade, HPPercent, ShooterHardAtk, ShooterSoftAtk, VictimHardness,
               ShooterPierce, ShooterArmor, VictimPierce, VictimArmor)
    hardness = VictimHardness
    firepower = ShooterHardAtk * VictimHardness + ShooterSoftAtk * (1-VictimHardness)
    if ShooterPierce < VictimArmor
        firepower /= 2
    end

    fire_hit = round_random(firepower / 10)
    evade_hit = round_random(evade / 10)

    fire_heavy = max(fire_hit - evade_hit, 0)
    fire_light = min(fire_hit, evade_hit)

    hit = sum(rand(fire_heavy) .< 0.4) + sum(rand(fire_light) .< 0.1)

    HP_dice_size = 2
    if ShooterArmor > VictimPierce
        Org_dice_size = 6
    else
        Org_dice_size = 4
    end

    HP_loss =  sum(floor.(rand(hit) .*  HP_dice_size) .+ 1) * 0.05
    Org_loss = sum(floor.(rand(hit) .* Org_dice_size) .+ 1) * 0.05

    if ShooterPierce < VictimArmor
        HP_loss /= 2
        Org_loss /= 2
    end

    p = floor(HPPercent, digits=1)
    HP_loss *= p
    Org_loss *= p

    return HP_loss, Org_loss
end

"""
    fire(shooter, victin, victim_pos)

shooter == attacker, victim == defender in wiki jargon, but I prefer to not use words atk/def
since they make confusion with "region" attacker and defender, which is indicated by victim_pos.
victim ∈ :attacker, :defender

Note: I misunderstood the last argument as shooter_pos many times! -_-
"""
function fire(shooter::Division, victim::Division, victim_pos::Symbol)
    if victim_pos == :attacker
        evade = victim.T.Breakthr
    elseif victim_pos == :defender
        evade = victim.T.Defense
    else
        error("Unknown victim_pos")
    end

    HPPercent = shooter.HP / shooter.T.HP

    fire(evade=evade, HPPercent=HPPercent, ShooterHardAtk=shooter.T.HardAtk,
         ShooterSoftAtk=shooter.T.SoftAtk, VictimHardness=victim.T.Hardness,
         ShooterPierce=shooter.T.Pierce, ShooterArmor=shooter.T.Armor,
         VictimPierce=victim.T.Pierce, VictimArmor=victim.T.Armor)
end

function fire!(shooter::Division, victim::Division, victim_pos::Symbol)
    HP_loss, Org_loss = fire(shooter, victim, victim_pos)
    victim.HP = max(victim.HP - HP_loss, 0)
    victim.Org = max(victim.Org - Org_loss, 0)
end

"Used to test only, real game resolve damage in meantime."
function duel(div_atk::Division, div_def::Division)
    HP_atk = [div_atk.HP]
    HP_def = [div_def.HP]
    Org_atk = [div_atk.Org]
    Org_def = [div_def.Org]
    while !defeated(div_atk) && !defeated(div_def)
        fire!(div_def, div_atk, :defender)
        fire!(div_atk, div_def, :attacker)
        push!(HP_atk, div_atk.HP)
        push!(HP_def, div_def.HP)
        push!(Org_atk, div_atk.Org)
        push!(Org_def, div_def.Org)
    end
    return (HP_atk=HP_atk, HP_def=HP_def, Org_atk=Org_atk, Org_def=Org_def)
end


text:
# ---
# title: 477. Total Hamming Distance
# id: problem477
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Bit Manipulation
# link: <https://leetcode.com/problems/total-hamming-distance/description/>
# hidden: true
# ---
#
# The [Hamming distance](https://en.wikipedia.org/wiki/Hamming_distance) between
# two integers is the number of positions at which the corresponding bits are
# different.
#
# Now your job is to find the total Hamming distance between all pairs of the
# given numbers.
#
# **Example:**
#
#
#
#     Input: 4, 14, 2
#
#     Output: 6
#
#     Explanation: In binary representation, the 4 is 0100, 14 is 1110, and 2 is 0010 (just
#     showing the four bits relevant in this case). So the answer will be:
#     HammingDistance(4, 14) + HammingDistance(4, 2) + HammingDistance(14, 2) = 2 + 2 + 2 = 6.
#
#
# **Note:**
#
#   1. Elements of the given array are in the range of `0 ` to `10^9`
#   2. Length of the array will not exceed `10^4`.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
const READSTATE_INIT = 0
const READSTATE_STARTED = 1
const READSTATE_CANCELLED = 2

startreply!(state::Threads.Atomic{Int}) =
    Threads.atomic_cas!(state, READSTATE_INIT, READSTATE_STARTED) == READSTATE_INIT

cancelrequest!(state::Threads.Atomic{Int}) =
    Threads.atomic_cas!(state, READSTATE_INIT, READSTATE_CANCELLED) == READSTATE_INIT

struct ReadBytesReq
    out::Vector{UInt8}
    nbytes::Int
    state::Threads.Atomic{Int}
    waitwrite::typeof(Promise{Nothing}())
end

struct ReadUntilVectorReq
    target::Vector{UInt8}
    keep::Bool
    out::Vector{UInt8}
    state::Threads.Atomic{Int}
    waitwrite::typeof(Promise{Nothing}())
end

struct UnsafeWriteReq
    out::Vector{UInt8}
    nbytes::UInt
end

struct IOManager{Raw<:IO}
    raw::Raw
    buffer::Vector{UInt8}
    request_send::typeof(Reagents.channel(Any, Nothing)[1])
    request_receive::typeof(Reagents.channel(Any, Nothing)[2])
    shutdown::typeof(Promise{Nothing}())
    sent::typeof(Ref(0))
end

struct IOWrapper{Raw<:IO} <: IO
    manager::IOManager{Raw}
    task::Task
    close_raw::Bool
end

function Julio.open(raw::IO; close = false)
    request_send, request_receive = Reagents.channel(Any, Nothing)
    shutdown = Promise{Nothing}()
    manager = IOManager(
        raw,
        empty!(Vector{UInt8}(undef, 1024)),
        request_send,
        request_receive,
        shutdown,
        Ref(0),
    )
    task = Threads.@spawn manage!(manager)
    return IOWrapper(manager, task, close)
end

function Base.close(io::IOWrapper)
    Reagents.trysync!(tryputting(io.manager.shutdown), nothing)
    if io.close_raw
        close(io.manager.raw)
        wait(io.task)
    elseif istaskdone(io.task)
        wait(io.task)  # let it throw
    else
        errormonitor(io.task)  # TODO: Find a less unstructured handling.
    end
end

function manage!(manager::IOManager)
    (; request_receive, shutdown) = manager
    shutdown_receive = _fetching(shutdown)
    # shutdown_receive = shutdown_receive ⨟ PostCommit(_ -> @info "shutting down")
    while true
        try_receive = request_receive | shutdown_receive
        (; request, reply, abort) = @something(try_receive(), break)
        ans = try
            Some(exec!(request, manager))
        catch err
            # @error "`exec!`" exception = (err, catch_backtrace())
            Error(err)
        end
        try_reply = (
            (Return(ans) ⨟ reply ⨟ Return(true)) |
            (abort ⨟ Return(false)) |
            shutdown_receive
        )
        @trace(label = :iomanager_start_try_reply)
        if @something(try_reply(), break)
            commit!(manager)
        else
            @trace(label = :iomanager_aborted)
        end
    end
    @trace(label = :iomanager_done)
    return
end

function commit!(manager::IOManager)
    (; buffer, sent) = manager
    # TODO: lazily move elements
    deleteat!(buffer, 1:sent[])
    @trace(label = :iomanager_commit, sent = sent[], buffer_length = length(buffer))
    sent[] = 0
end

function send!(manager::IOManager, request)
    (; request_send) = manager
    return Reagents.WithNack() do abort
        reply, receive = Reagents.channel(Any, Nothing)
        request_send((; request, reply, abort))
        return receive
    end
end

function readbytes_append!(io, buffer, nbytes)
    readbytes!(io, resizabletail(buffer), nbytes - length(buffer))
    return buffer
end

event(
    ::typeof(readbytes!),
    io::IOWrapper,
    out::Vector{UInt8},
    nbytes::Integer,
    state::Threads.Atomic{Int},
    waitwrite::Promise{Nothing},
) = send!(io.manager, ReadBytesReq(out, nbytes, state, waitwrite))

function exec!(request::ReadBytesReq, manager::IOManager)
    (; raw, buffer, sent) = manager
    (; out, nbytes, state, waitwrite) = request
    if nbytes <= length(buffer)
        startreply!(state) || return 0
        try
            sent[] = nbytes
            resize!(out, nbytes)
            copyto!(out, @view buffer[1:nbytes])
            @trace(label = :readbytes_use_buffer, sent = sent[])
        finally
            waitwrite[] = nothing
        end
    else
        # TODO: chunk large `nbytes` and check for cancellation
        readbytes_append!(raw, buffer, nbytes)
        startreply!(state) || return 0
        try
            sent[] = length(buffer)
            resize!(out, length(buffer))
            copyto!(out, buffer)
            @trace(label = :readbytes_read_raw, sent = sent[])
        finally
            waitwrite[] = nothing
        end
    end
    return sent[]
end

event(
    ::typeof(Base.readuntil_vector!),
    io::IOWrapper,
    target::Vector{UInt8},
    keep::Bool,
    out::Vector{UInt8},
    state::Threads.Atomic{Int},
    waitwrite::Promise{Nothing},
) = send!(io.manager, ReadUntilVectorReq(target, keep, out, state, waitwrite))

function exec!(request::ReadUntilVectorReq, manager::IOManager)
    (; raw, buffer, sent) = manager
    (; target, keep, out, state, waitwrite) = request
    iobuf = IOBuffer(buffer)
    if has_subseq(buffer, target)
        startreply!(state) || return
        try
            Base.readuntil_vector!(iobuf, target, keep, out)
        finally
            waitwrite[] = nothing
        end
    else
        Base.readuntil_vector!(raw, target, true, resizabletail(buffer))
        if keep
            result = @view buffer[begin:end]
        else
            result = @view buffer[begin:end-length(target)]
        end
        startreply!(state) || return
        try
            resize!(out, length(result))
            copyto!(out, result)
        finally
            waitwrite[] = nothing
        end
    end
    sent[] = length(buffer)
    return
end

function event(::typeof(unsafe_write), io::IOWrapper, ptr::Ptr{UInt8}, nbytes::UInt)
    # For cancel-safety, we need to copy the input now:
    out = Vector{UInt8}(undef, nbytes)
    copyto!(out, unsafe_wrap(Array, ptr, nbytes))
    send!(io.manager, UnsafeWriteReq(out, nbytes))
end

function exec!(request::UnsafeWriteReq, manager::IOManager)
    (; raw) = manager
    (; out, nbytes) = request
    GC.@preserve out begin
        ans = unsafe_write(raw, pointer(out), nbytes)
    end
    return ans
end

function event(::typeof(read), io::IOWrapper, nbytes::Integer = typemax(Int))
    state = Threads.Atomic{Int}()
    waitwrite = Promise{Nothing}()
    out = UInt8[]
    return event(readbytes!, io, out, nbytes, state, waitwrite) ⨟ Return(Some(out))
end

function event(::typeof(readline), io::IOWrapper; keep::Bool = false)
    if Sys.iswindows()
        target = collect(transcode(UInt8, "\r\n"))
    else
        target = collect(transcode(UInt8, "\n"))
    end
    return Computed() do _
        state = Threads.Atomic{Int}()
        waitwrite = Promise{Nothing}()
        out = UInt8[]
        return event(Base.readuntil_vector!, io, target, keep, out, state, waitwrite) ⨟
               Map() do _
            Some(String(out))
        end
    end
end

# ## Blocking interface

function Base.readbytes!(io::IOWrapper, out::Vector{UInt8}, nbytes::Integer = typemax(Int))
    state = Threads.Atomic{Int}()
    waitwrite = Promise{Nothing}()
    try
        return Julio.sync(event(readbytes!, io, out, nbytes, state, waitwrite))
    finally
        cancelrequest!(state) || fetch(waitwrite)
        # TODO: automate this (`out` ownership/locking handling)?
    end
end

function Base.readuntil_vector!(
    io::IOWrapper,
    target::Vector{UInt8},
    keep::Bool,
    out::Vector{UInt8},
)
    state = Threads.Atomic{Int}()
    waitwrite = Promise{Nothing}()
    try
        return Julio.sync(
            event(Base.readuntil_vector!, io, target, keep, out, state, waitwrite),
        )
    finally
        cancelrequest!(state) || fetch(waitwrite)
    end
end

function Base.readuntil(io::IOWrapper, delim::UInt8; keep::Bool = false)
    out = UInt8[]
    Base.readuntil_vector!(io, [delim], keep, out)
    return out
end

Base.unsafe_write(io::IOWrapper, ptr::Ptr{UInt8}, nbytes::UInt) =
    Julio.sync(event(unsafe_write, io, ptr, nbytes))

Julio.open(path::AbstractString; kwargs...) =
    Julio.open(open(path; kwargs...); close = true)

Julio.open(cmd::Base.AbstractCmd, args...; kwargs...) =
    Julio.open(open(cmd, args...; kwargs...); close = true)

function withopen(f, cmd::Base.AbstractCmd, args...; kwargs...)
    open(cmd, args...; kwargs...) do io
        withopen(f, io)
    end
end


text:
using Revise
using NNFEM
using Statistics

nodes, elements = meshread("crack.msh")

visualize_mesh(nodes, elements)

ids = Int64[]
for i = 1:size(nodes, 1)
    x, y = nodes[i,:]
    if  abs(y-0.5*(x-0.5))<0.01 && (0.5 < x <= 1.55)
        push!(ids, i)
        @info x, y
        plot([x], [y], ".g")
    end
end

new_nodes = []
adic = Dict{Int64, Int64}()
for i = 1:size(elements, 1)
    e = elements[i,:]
    if !any([j in ids for j in e])
        continue
    end

    @info i
    x, y = mean(nodes[e,:], dims=1)[:]
    if y-0.5*(x-0.5)<0
        plot([x], [y], ".r")
    else
        continue
    end

    for j = 1:4
        if e[j] in ids
            if haskey(adic,e[j])
                elements[i,j] = adic[e[j]]
            else
                push!(new_nodes, nodes[e[j],:] + 0.02*rand(2))
                elements[i,j] = size(new_nodes,1) + size(nodes,1)
                adic[e[j]] = size(new_nodes,1) + size(nodes,1)
            end
        end
    end

end

new_nodes = hcat(new_nodes...)
nodes = [nodes;new_nodes']

close("all")
visualize_mesh(nodes, elements)


text:
#!/usr/bin/env julia

# Because we've been having so many errors with launching, run a check script.
# Load config, print the parameters that we will be sampling and the mean starting position of the walkers.

using ArgParse

s = ArgParseSettings()
@add_arg_table s begin
    "config"
    help = "a YAML configuration file"
    default = "config.yaml"
end

parsed_args = parse_args(ARGS, s)

import YAML
config = YAML.load(open(parsed_args["config"]))

using Statistics
using DiskJockey.model
using DiskJockey.visibilities
using DiskJockey.constants
using NPZ

species = config["species"]
transition = config["transition"]
lam0 = lam0s[species * transition]
model = config["model"]
params = config["parameters"]
fix_params = config["fix_params"]
reg_params = registered_params[model]

println("Running model $model with the $species $transition transition.")
println()

println("We are fixing the following parameters to these values:")
for param in fix_params
  println("$param : ", params[param])
end

println()

fit_params = filter(x->∉(x, fix_params), reg_params)
n_fit = length(fit_params)

pos0 = npzread(config["pos0"])
ndim, nwalkers = size(pos0)
means = mean(pos0, dims = 2) # Find the mean walker position for each parameter

@assert ndim == n_fit "Number of specified parameters ($n_fit) does not match the number of dimensions in the walker file ($ndim)."

println("We are fitting with the following $n_fit parameters:")
for (par, mu) in zip(fit_params, means)
  println("$par : $mu")
end

println("Using $nwalkers walkers in the MCMC sampling.")

# load data and figure out how many channels
dvarr = DataVis(config["data_file"])
nchan = length(dvarr)

if haskey(config, "exclude")
    exclude = config["exclude"]
    # which channels of the dset to fit
    # keylist = filter(x->(!in(x, exclude)), Int[i for i=1:nchan])

    lam0 = lam0s[config["species"] * config["transition"]]
    # calculate the velocities corresponding to dvarr
    lams = Float64[dv.lam for dv in dvarr]
    vels = c_kms * (lams .- lam0) / lam0
    # get the mask
    vel_mask = generate_vel_mask(exclude, vels)
else
    # keylist = Int[i for i=1:nchan]
    vel_mask = trues(nchan)
end

println("Using the following velocity channels: ", vels[vel_mask])


text:
@testset "Linear Algebra" begin
    @testset "lu" begin
        m = AxisArray([1.0 2; 3 4], (2:3, 3:4))
        x = @inferred(lu(m))

        @test m[x.p, :] == parent(m)[parent(x).p,:]
        @test x.factors == lu(parent(m)).factors
    end

    # FIXME inferrence for eigen
    @testset "eigen" begin
        m = AxisArray([1.0 2; 3 4], (2:3, 3:4))
        x = eigen(m)
        @test keys.(axes(x.vectors)) == (2:3, 3:4)
        @test eigvals(m) == eigvals(parent(m))
    end

    @testset "svd" begin
        m = AxisArray([1.0 2; 3 4], (2:3, 3:4))
        x = svd(m)
        u, s, v = x
        @test size(x) == (2, 2)
        @test size(x) == (2, 2)
        @test size(x, 1) == 2
        @test keys.(axes(x.S)) == (1:2,)

        F = svd(AxisArray([1.0 2; 3 4], (Axis(2:3 => Base.OneTo(2)), Axis(3:4 => Base.OneTo(2)))));
        io = IOBuffer()
        show(io, F)
        str = String(take!(io))
        @test str[1:7] == "AxisSVD"
    end

    @testset "LQ" begin
        m = AxisArray([1.0 2; 3 4], (2:3, 3:4))
        x = @inferred(lq(m))
        @test x.factors == lq(parent(m)).factors
    end

    @testset "qr" begin
        for pivot in (true, false)
            for data in ([1.0 2; 3 4], [big"1.0" 2; 3 4])
                m = AxisArray(data, (2:3, 3:4))
                x = @inferred(qr(m, Val(pivot)));

                @test keys.(axes(x.Q)) == (2:3, 1:2)
                @test keys.(axes(x.R)) == (1:2, 3:4)
                @test keys.(axes(x.Q * x.R)) == (2:3, 3:4)
                @test x.factors == qr(parent(m), Val(pivot)).factors

                pivot && @testset "pivoted" begin
                    @test x isa QRPivoted
                    @test x.p == [3, 2]
                    @test keys.(axes(x.P)) == (2:3, 2:3)
                    @test keys.(axes(x.P * m)) == (2:3, 3:4)
                    @test m[x.p,:] == parent(m)[parent(x).p,:]
                end
            end
        end
    end

    @testset "diagonal" begin
        m = AxisArray([1 2 3; 4 5 6; 7 8 9], ["a", "b", "c"], nothing);
        @test @inferred(diag(m)) == diag(parent(m))
        @test @inferred(diag(m, 1)) == diag(parent(m), 1)
        @test @inferred(diag(m, 1; dim=Val(2))) == diag(parent(m), 1)
    end
end

#= TODO figure allocations
@testset "allocations: matmul names" begin
    @test 0 == @allocated (() -> matrix_prod_names((:foo, :bar), (:bar,)))()
    @test 0 == @allocated (() -> symmetric_names((:foo, :bar), 1))()
end
=#

@testset "Mutmul with special types" begin
    a = AxisArray(ones(5,5), (2:6, 3:7))
    @testset "$T" for T in (Diagonal, Symmetric, Tridiagonal, SymTridiagonal, BitArray,)
        x = T(ones(5,5))
        @test keys.(axes(x * a)) == (1:5, 3:7)
        @test keys.(axes(a * x)) == (2:6, 1:5)
    end
end

@testset "inv" begin
    a = AxisArray([1.0 2; 3 4], (2:3, 4:5));
    @test keys.(axes(inv(a))) == (4:5, 2:3)
    @test a * inv(a) ≈ [1.0 0; 0 1]
    @test keys.(axes(a * inv(a))) == (2:3, 2:3)

    @test inv(a) * a ≈ [1.0 0; 0 1]
    @test keys.(axes(inv(a) * a)) == (4:5, 4:5)
end

@testset "cov/cor" begin
    @testset "$f" for f in (cov, cor)
        @testset "matrix input, matrix result" begin
            A = rand(3, 5)
            a = AxisArray(A, (2:4, 2:6))
            @test f(a; dims=1) == f(A, dims=1)
            @test keys.(axes(f(a; dims=1))) == (2:6, 2:6)
            @test keys.(axes(f(a, dims=2))) == (2:4, 2:4)
            @test keys.(axes(f(a, dims=3))) == (2:4, 2:4)
            @test keys.(axes(f(a[:,2:4], dims=3))) == (2:4, 2:4)
        end
        @testset "vector input, scalar result" begin
            v = rand(4)
            av = AxisArray(v, (2:5,))
            @test f(av) isa Number
            @test f(av) == f(v)
        end
    end
    @testset "cov corrected=$bool" for bool in (true, false)
        # test that kwargs get passed on correctly
        A = rand(2, 4)
        a = AxisArray(A)
        @test cov(a; corrected=bool) == cov(A; corrected=bool)
        @test cov(a; corrected=bool, dims=2)  == cov(A; corrected=bool, dims=2)
    end
end


text:
## Overridden behaviours
"""
    WrapperLogType

Abstract wrapper type to override the default backend to which a type would be
logged in TensorBoard.

Every concrete implementation should implement the method `content` for the
concrete type, returning the content.
"""
abstract type WrapperLogType end

"""
    content(x::WrapperLogType)

Returns the wrapped data inside of `x`.
"""
content(x::WrapperLogType) = @error "You should define `content($(typeof(x)))`"

# When we hit a wrapped type, we dispatch based on it's type.
preprocess(name, val::WrapperLogType, data) = push!(data, name=>val)

# When logging to text or console, unwrap the LogType to fix formatting.
Base.show(io::IO, mime::AbstractString, x::WrapperLogType) =
    Base.show(io, mime, content(x))

########## For things going to LogImage ########################
"""
    TBImage(data, format)

Forces `data` to be serialized as an Image to TensorBoard.
"""
struct TBImages <:WrapperLogType
    data::AbstractArray
    format::ImageFormat
end
content(x::TBImages) = x.data
function preprocess(name, val::TBImages, data)
    for (n, img) in enumerate(val.data)
        preprocess(name*"/$n", TBImage(img, val.format), data)
    end
end

struct TBImage <: WrapperLogType
    data::AbstractArray
    format::ImageFormat
end
content(x::TBImage) = x.data
function preprocess(name, val::TBImage, data)
    imgArray = val.data
    format = val.format
    imgArray = channelview(imgArray)
    dims = ndims(imgArray)
    @assert dims == expected_ndims(format)
    obsdim = obs_dim(format)
    if iszero(obsdim)
        push!(data, name=>TBImage(imgArray, format))
    else
        format = strip_obs[format]
        index = collect("[:"* ",:"^(dims-1) *"]")
        index[2*obsdim] = 'g'
        index = join(index)
        global gimgArray = imgArray
        nth_img = Meta.parse("gimgArray$(index)")
        for n in 1:size(imgArray, obsdim)
            global g = n
            push!(data, name*"/$n"=>TBImage(eval(nth_img), format))
        end
    end
    data
end
summary_impl(name, val::TBImage) = image_summary(name, val.data, val.format)

########## For things going to LogText ##############################
"""
    TBText(data)

Forces `data` to be serialized as text to TensorBoard.
"""
struct TBText <: WrapperLogType data; end
content(x::TBText) = x.data
summary_impl(name, val::TBText) = text_summary(name, val.data)

"""
    RealArrayWrapperLogType{T,N} <: WrapperLogType

A wrapper type that accepts only arrays and converts complex-valued arrays
into two real value ones
"""
abstract type RealArrayWrapperLogType{T,N} <: WrapperLogType end

########## For things going to LogHistograms ########################
"""
    TBHistogram(data)

Forces `data` to be serialized as an histogram to TensorBoard.
"""
struct TBHistogram{T,N} <: RealArrayWrapperLogType{T,N} data::AbstractArray{T,N}; end
content(x::TBHistogram) = x.data
preprocess(name, val::TBHistogram{T,N}, data) where {T<:Complex,N} =
    push!(data, name*"/re"=>TBHistogram(real.(content(val))), name*"/im"=>TBHistogram(imag.(content(val))))
function summary_impl(name, val::TBHistogram)
    hvals = fit(Histogram, val.data)
    return histogram_summary(name, collect(hvals.edges[1]), hvals.weights)
end

"""
    TBVector(data)

Forces `data` to be serialized as a vector in the histogram backend of
TensorBoard.
"""
struct TBVector{T,N} <: RealArrayWrapperLogType{T,N} data::AbstractArray{T,N}; end
content(x::TBVector) = x.data
preprocess(name, val::TBVector{T,N}, data) where {T<:Complex,N} =
    push!(data, name*"/re"=>TBVector(real.(content(val))), name*"/im"=>TBVector(imag.(content(val))))
summary_impl(name, val::TBVector) = histogram_summary(name, collect(0:length(val.data)),
                                                            val.data)


text:
using Test
using Quaternions

include("helpers.jl")
include("Quaternion.jl")
include("Octonion.jl")
include("DualQuaternion.jl")


text:
export basemeasure

import MeasureTheory

function MeasureTheory.basemeasure(c::ConditionalModel{A,B,M}, x=NamedTuple()) where {A,B,M}
    _basemeasure(M, Model(c), argvals(c), x)  | observations(c)
end

export sourceBasemeasure

sourceBasemeasure(m::AbstractModel) = sourceBasemeasure()(Model(m))

function sourceBasemeasure()
    function(_m::Model)
        proc(_m, st :: Assign)     = :($(st.x) = $(st.rhs))
        proc(_m, st :: Return)     = nothing
        proc(_m, st :: LineNumber) = nothing
        function proc(_m, st :: Sample)
            x = st.x
            xname = QuoteNode(x)
            rhs = st.rhs
            @q begin
                _bm = merge(_bm, NamedTuple{($xname,)}((basemeasure($rhs),)))
                $x = Soss.testvalue($rhs)
            end
        end

        wrap(kernel) = @q begin
            _bm = (;)
            $kernel
            return Soss.productmeasure(_bm)
        end

        buildSource(_m, proc, wrap) |> MacroTools.flatten
    end
end

@gg function _basemeasure(M::Type{<:TypeLevel}, _m::Model, _args, _pars)
    body = type2model(_m) |> sourceBasemeasure() |> loadvals(_args, _pars)
    @under_global from_type(_unwrap_type(M)) @q let M
        $body
    end
end


text:
#=
inputs matrices A,B
outputs a matrix C
Steps
  check compatibility
  create nxn matrix C
  partition the three matrices into n/2 x n/2 submatrices (O(1) work/span)
  create 10 matrices S1..S10 of size n/2 x n/2 and sum or difference of two partitions
    via doubly nested parallel for loops
  recursively spawn the computation of seven n/2 x n/2 matrix products
  compute submatrices of C by adding subtracting various of the above
    via doubly nested parallel for loops (Theta(n^2) work, Theta(logn) span)

max_stars_repo_name: K-P-S/Julia-Scientific-Computing

max_issues_repo_name: K-P-S/Julia-Scientific-Computing

max_forks_repo_name: K-P-S/Julia-Scientific-Computing


text:
export PolyhedralStartSolutionsIterator, update_cell!, mixed_volume

import MixedSubdivisions: MixedCell, MixedCellIterator, RegenerationTraverser, mixed_volume
import ElasticArrays: ElasticArray

#####################
## START SOLUTIONS ##
#####################

"""
    PolyhedralStartSolutionsIterator(f)

Creates an iterator over all mixed cells and their solutions corresponding to the polynomial
system `f`.
The iterator returns a tuple `(cell, x)` where `cell` is a `MixedCell` and
`x` is a vector matrix containing  a solution of the binomial system corresponding to
the mixed cell.

    PolyhedralStartSolutionsIterator(support)

Construct a polyhedral start solutions iterator by picking a random lifting
and random start coefficients.
"""
mutable struct PolyhedralStartSolutionsIterator{CT<:CoreTracker}
    support::Vector{Matrix{Int32}}
    start_coefficients::Vector{Vector{ComplexF64}}
    lifting::Union{Nothing,Vector{Vector{Int32}}}
    mixed_cells::Union{Nothing,Vector{MixedCell}} # all mixed_cells
    X::ElasticArray{ComplexF64,2,1} # solutions
    # Cache
    D::Matrix{Int} # binomial system lhs
    b::Vector{ComplexF64} # binomial system rhs
    H::Matrix{Int} # hermite normal form
    U::Matrix{Int} # trafo matrix for hnf
    γ::Vector{Float64} # angle of the elements in b
    μ::Vector{Float64} # γ^U
    S::ElasticArray{Int,2,1}
    αs::Vector{Float64} # partial angles in the triangular solve
    DT::Matrix{Float64} # transposed of D to solve coords
    # binomial homotopy fallback solution
    γU::Vector{Int}
    x::Vector{ComplexF64}
    αs_rational::Vector{Rational{Int}}
    binomial_tracker::CT
end

function PolyhedralStartSolutionsIterator(f::MPPolys)
    PolyhedralStartSolutionsIterator(MixedSubdivisions.support(f))
end
function PolyhedralStartSolutionsIterator(support::Vector{<:Matrix{<:Integer}})
    # Intuitively we want to sample the coefficients from the normal distribution
    # because we like this one very much.
    # But this yields to a bad surpise for the start solutions in the case that
    # the ratios of the coefficients in the binomial system is very different in between
    # the equations
    start_coefficients = map(A -> rand_approx_unit(size(A, 2)), support)
    PolyhedralStartSolutionsIterator(
        convert(Vector{Matrix{Int32}}, support),
        start_coefficients,
    )
end

"""
    rand_approx_unit(n::Integer)::Vector{ComplexF64}

This samples uniformly from the rectangle ``[0.9,1.1] × [0,2π]`` and transforms the sampled
values with the map ``(r, φ) ↦ r * e^{i φ}``.
"""
rand_approx_unit(n) = map(_ -> (0.9 + 0.2 * rand()) * cis(2π * rand()), 1:n)

function PolyhedralStartSolutionsIterator(
    support::Vector{Matrix{Int32}},
    start_coefficients::Vector{Vector{ComplexF64}},
)
    n = size(first(support), 1)
    D = zeros(Int, n, n)
    b = ones(ComplexF64, n)
    H = zeros(Int, n, n)
    U = zeros(Int, n, n)
    γ = ones(Float64, n)
    μ = zeros(Float64, n)
    X = ElasticArray{ComplexF64}(undef, n, 0)
    S = ElasticArray{Int}(undef, n, 0)
    αs = zeros(Float64, n)
    DT = zeros(Float64, n, n)

    γU = zeros(Int, n)
    x = randn(ComplexF64, n)
    αs_rational = Vector{Rational{Int64}}(undef, n)
    binomial_tracker = coretracker(BinomialHomotopy(D, b, γ), x; affine_tracking = true)

    PolyhedralStartSolutionsIterator(
        support,
        start_coefficients,
        nothing,
        nothing,
        X,
        D,
        b,
        H,
        U,
        γ,
        μ,
        S,
        αs,
        DT,
        γU,
        x,
        αs_rational,
        binomial_tracker,
    )
end

Base.IteratorSize(::Type{<:PolyhedralStartSolutionsIterator}) = Base.HasLength()
Base.IteratorEltype(::Type{<:PolyhedralStartSolutionsIterator}) = Base.HasEltype()
Base.eltype(iter::PolyhedralStartSolutionsIterator) = Tuple{MixedCell,Vector{ComplexF64}}
function Base.length(iter::PolyhedralStartSolutionsIterator)
    compute_mixed_cells!(iter)
    sum(MixedSubdivisions.volume, iter.mixed_cells)
end

lifting_sampler(nterms) = rand(Int32(-2^12):Int32(2^12), nterms)

function compute_mixed_cells!(iter::PolyhedralStartSolutionsIterator)
    if isnothing(iter.mixed_cells) || isnothing(iter.lifting)
        res = MixedSubdivisions.fine_mixed_cells(iter.support, lifting_sampler)
        if isnothing(res) || isempty(res[1])
            throw(OverflowError("Cannot compute a start system due."))
        end
        mixed_cells, lifting = res
        iter.mixed_cells = mixed_cells
        iter.lifting = lifting
    end
    iter
end

function Base.iterate(iter::PolyhedralStartSolutionsIterator)
    compute_mixed_cells!(iter)

    isempty(iter.mixed_cells) && return nothing

    cell = iter.mixed_cells[1]
    cell_solutions!(iter, cell)

    el = (cell, iter.X[:, 1])
    next_state = size(iter.X, 2) == 1 ? (2, 1) : (1, 2)
    el, next_state
end

function Base.iterate(iter::PolyhedralStartSolutionsIterator, (i, j)::Tuple{Int,Int})
    i > length(iter.mixed_cells) && return nothing

    cell = iter.mixed_cells[i]
    j == 1 && cell_solutions!(iter, cell)
    el = (cell, iter.X[:, j])
    next_state = size(iter.X, 2) == j ? (i + 1, 1) : (i, j + 1)
    return el, next_state
end

start_solution_sample(iter::PolyhedralStartSolutionsIterator) =
    randn(ComplexF64, size(iter.X, 1))

###########################
## SOLVE BINOMIAL SYSTEM ##
###########################
function cell_solutions!(iter::PolyhedralStartSolutionsIterator, cell::MixedCell)
    @unpack support, start_coefficients, D, b, H, U = iter
    n = length(iter.support)

    # setup system and compute hermite normal form
    construct_binomial_system!(D, b, support, start_coefficients, cell.indices)

    # This beauty tries first to compute the hnf using Int64,
    # then falls back to Int128 and finally gives up and just uses
    # BigInt
    try
        hnf!(H, U, D)
        cell_solutions!(iter, cell, H, U)
    catch e
        isa(e, OverflowError) || rethrow(e)
        try
            H128, U128 = hnf(D, Int128)
            cell_solutions!(iter, cell, H128, U128)
        catch e128
            isa(e128, OverflowError) || rethrow(e128)

            H_big, U_big = hnf(D, BigInt)
            cell_solutions!(iter, cell, H_big, U_big)
        end
    end

    nothing
end

function cell_solutions!(iter::PolyhedralStartSolutionsIterator, cell::MixedCell, H, U)
    @unpack D, b, γ, μ, X, S, DT = iter
    n = size(H, 1)
    # resize memory
    d̂ = cell.volume
    resize!(S, n, d̂)
    resize!(X, n, d̂)

    # setup a table of all combinations of the types of roots of unity we need
    fill_unit_roots_combinations!(S, H, d̂)

    cell_solutions_direct!(iter, cell, H, U)
    if !verify_solutions(iter)
        cell_solutions_homotopy!(iter, cell, H, U)
    end
end

function verify_solutions(iter)
    @unpack D, b, X = iter
    n = size(D, 1)
    for k = 1:size(X, 2), i = 1:n
        yᵢ = one(X[1, k])
        for j = 1:n
            if D[j, i] != 0
                yᵢ *= X[j, k]^D[j, i]
            end
        end
        abs(yᵢ - b[i]) < 1e-10 || return false
    end
    true
end

function cell_solutions_direct!(
    iter::PolyhedralStartSolutionsIterator,
    cell::MixedCell,
    H,
    U,
)
    @unpack D, b, μ, X, DT = iter
    n = size(D, 1)
    # We split the computation of the solution in 2 stages
    # 1) Compute solutions for the angular part -> d̂ many solutions
    compute_angular_part!(iter, cell, H, U)
    # 2) Solve for the absolute value
    μ .= log.(abs.(b))
    LinearAlgebra.transpose!(DT, D)
    LinearAlgebra.ldiv!(LinearAlgebra.lu!(DT), μ)
    μ .= exp.(μ)
    for j = 1:cell.volume, i = 1:n
        X[i, j] *= μ[i]
    end

    nothing
end

function fill_unit_roots_combinations!(S, H, d̂)
    n = size(H, 1)
    d, e = d̂, 1
    @inbounds for I = 1:n
        dᵢ = convert(Int64, H[I, I])
        d = d ÷ dᵢ
        k = 1
        for _ = 1:e, j = 0:(dᵢ-1), _ = 1:d
            S[I, k] = j
            k += 1
        end
        e *= dᵢ
    end
end

function compute_angular_part!(
    iter::PolyhedralStartSolutionsIterator,
    cell::MixedCell,
    H::Matrix{I},
    U::Matrix{I},
) where {I<:Union{Int64,Int128}}
    @unpack γ, b, μ, X, S, αs = iter
    n = length(iter.support)
    γ .= angle.(b) ./ (2π)
    apply_trafo_angles!(μ, γ, U)
    @inbounds for i = 1:cell.volume, j = n:-1:1
        α = μ[j] + S[j, i]
        for k = n:-1:(j + 1)
            α -= αs[k] * H[k, j]
        end
        α /= H[j, j]
        X[j, i] = cis(2π * α)
        αs[j] = α
    end
end

function compute_angular_part!(
    iter::PolyhedralStartSolutionsIterator,
    cell::MixedCell,
    H::Matrix{BigInt},
    U::Matrix{BigInt},
)
    @unpack b, γ, X, S = iter
    n = length(iter.support)
    μ = Vector{BigFloat}(undef, n)
    αs = Vector{BigFloat}(undef, n)
    γ .= angle.(b) ./ (2π)
    apply_trafo_angles!(μ, γ, U)
    @inbounds for i = 1:cell.volume, j = n:-1:1
        # This is always small in our applications
        α = μ[j] + S[j, i]
        for k = n:-1:(j + 1)
            α -= αs[k] * H[k, j]
        end
        α /= H[j, j]
        X[j, i] = cis(2π * convert(Float64, mod(α, 1.0)))
        αs[j] = α
    end
end

function apply_trafo_angles!(μ, γ, U)
    n = length(γ)
    μ .= zero(eltype(μ))
    @inbounds for j = 1:n, i = 1:n
        μ[j] += U[i, j] * γ[i]
    end
    μ
end

function construct_binomial_system!(D, b, support, coeffs, indices)
    n = length(indices)
    @inbounds for i = 1:n
        aᵢ, bᵢ = indices[i]
        for j = 1:n
            D[j, i] = support[i][j, aᵢ] - support[i][j, bᵢ]
        end
        b[i] = -coeffs[i][bᵢ] / coeffs[i][aᵢ]
    end
    nothing
end

function cell_solutions_homotopy!(
    iter::PolyhedralStartSolutionsIterator,
    cell::MixedCell,
    H,
    U,
)
    @unpack D, b, γ, μ, X, γU, x, S, αs_rational, binomial_tracker = iter

    γ .= sign.(real.(b))
    apply_trafo!(γU, γ, U)

    n = size(D, 1)
    for i = 1:cell.volume
        αs_denom = 1
        for j = n:-1:1
            α = S[j, i]
            denom = lcm(αs_denom, convert(Int, H[j, j]))
            if γU[j] == -1
                α += 1 // 2
                denom *= 2
            end
            for k = n:-1:(j + 1)
                α -= αs_rational[k] * convert(Int, H[k, j] % denom)
            end
            α //= H[j, j]
            αs_denom = lcm(αs_denom, denominator(α))
            αs_rational[j] = (numerator(α) % denominator(α)) // denominator(α)
            x[j] = cis(2π * αs_rational[j])
        end

        ret = track!(binomial_tracker, x, 1.0, 0.0)
        x̄ = current_x(binomial_tracker)
        for j = 1:n
            X[j, i] = x̄[j]
        end
    end

    nothing
end

function apply_trafo!(γU, γ, U)
    for i = 1:length(γ)
        γUᵢ = 1
        for j = 1:size(U, 2)
            if γ[j] == -1 && isodd(U[j, i])
                γUᵢ = -γUᵢ
            end
        end
        γU[i] = γUᵢ
    end
    γU
end

########################
## POLYHEDRAL TRACKER ##
########################

"""
    PolyhedralTracker{CT<:CoreTracker, PT<:PathTracker, H<:ToricHomotopy}

The polyhedral tracker combines the tracking from toric infinity toward the target system
by a two-stage approach.
"""
struct PolyhedralTracker{
    CT<:CoreTracker,
    PT<:PathTracker,
    H<:ToricHomotopy,
} <: AbstractPathTracker
    toric_homotopy::H
    toric_tracker::CT
    generic_tracker::PT
    # state
    s₀::Base.RefValue{Float64}
    curr_mixed_cell::Base.RefValue{Union{Nothing,MixedCell}}
end

function construct_tracker(
    prob::PolyhedralProblem,
    startsolutions::PolyhedralStartSolutionsIterator;
    kwargs...,
)
    x = randn(ComplexF64, size(prob.toric_homotopy)[2])
    toric_tracker = coretracker(
        prob.toric_homotopy,
        [x];
        seed = prob.seed,
        affine_tracking = true,
        logarithmic_time_scale = true,
        from_infinity = true,
        track_cond = false,
    )
    generic_prob = Problem(
        prob.tracking_type,
        prob.generic_homotopy,
        prob.vargroups,
        prob.seed,
        prob.startsolutions_need_reordering,
        prob.regauging_factors,
    )
    generic_tracker = construct_tracker(generic_prob, startsolutions; kwargs...)

    PolyhedralTracker(
        prob.toric_homotopy,
        toric_tracker,
        generic_tracker,
        Ref(NaN),
        Base.RefValue{Union{Nothing,MixedCell}}(nothing),
    )
end

seed(PT::PolyhedralTracker) = seed(PT.generic_tracker)
function PathResult(
    PT::PolyhedralTracker,
    (cell, x)::Tuple{MixedCell,Vector{ComplexF64}},
    path_number = nothing;
    kwargs...,
)
    PathResult(PT.generic_tracker, nothing, path_number; kwargs...)
end

result_type(PT::PolyhedralTracker) = result_type(PT.generic_tracker)
solution(tracker::PolyhedralTracker) = solution(tracker.generic_tracker)
LA.norm(tracker::PolyhedralTracker) = LA.norm(tracker.generic_tracker)

function prepare!(PT::PolyhedralTracker, S::PolyhedralStartSolutionsIterator)
    if S.lifting === nothing
        compute_mixed_cells!(S)
    end
    update_lifting!(PT.toric_homotopy, S.lifting)
    PT
end
function track!(
    PT::PolyhedralTracker,
    (cell, x∞)::Tuple{MixedCell,Vector{ComplexF64}};
    accuracy::Union{Nothing,Float64} = nothing,
    max_corrector_iters::Union{Nothing,Int} = nothing,
)
    update_cell!(PT, cell)
    retcode = track!(PT.toric_tracker, x∞, PT.s₀[], 0.0)
    if is_invalid_startvalue(retcode)
        track!(PT.toric_tracker, x∞, 2 * PT.s₀[], 0.0)
    end
    track!(
        PT.generic_tracker,
        current_x(PT.toric_tracker);
        accuracy = accuracy,
        max_corrector_iters = max_corrector_iters,
    )
end

function track(
    tracker::PolyhedralTracker,
    cell_x::Tuple{MixedCell,Vector{ComplexF64}};
    details::Symbol = :default,
    path_number::Union{Int,Nothing} = nothing,
    accuracy::Union{Nothing,Float64} = nothing,
    max_corrector_iters::Union{Nothing,Int} = nothing,
)
    track!(tracker, cell_x; accuracy = accuracy, max_corrector_iters = max_corrector_iters)
    PathResult(tracker, cell_x, path_number; details = details)
end

"""
    update_cell!(PT::PolyhedralTracker, cell::MixedCell)

Update the polyhedral tracker `PT` to track paths coming from the mixed cell `cell`.
"""
function update_cell!(PT::PolyhedralTracker, cell::MixedCell)
    curr_cell = PT.curr_mixed_cell[]
    if curr_cell !== nothing
        curr_cell === cell && return PT
    end

    PT.curr_mixed_cell[] = cell
    min_weight = update_cell!(PT.toric_homotopy, cell)
    # We construct s₀ such that
    #   min_{1 ≤ i ≤ n} min_{aᵢ ∈ Aᵢ} exp(s₀ * w(aᵢ)+aᵢ⋅γ-βᵢ) = 10^-8
    # From this follows:
    PT.s₀[] = -8.0 * log(10) / min_weight
    PT
end


text:
######################################################################
# kalman.jl
# helper functions for ekf, ukf, etc
# Each initializer should have:
#  length
#  mu
#  Sigma
######################################################################

abstract type GaussianFilter <: AbstractFilter end
abstract type Initializer end

type NaiveInitializer <: Initializer
	length::Float64
	mu::Vector{Float64}
	Sigma::Matrix{Float64}

	function NaiveInitializer(m::SearchDomain)
		mu = m.length/2 * ones(2)
		S = 1e7 * eye(2)
		return new(m.length, mu, S)
	end
end

initialize(ni::NaiveInitializer, x::Vehicle, o::Float64) = true

type LSInitializer <: Initializer
	length::Float64
	mu::Vector{Float64}
	Sigma::Matrix{Float64}
	obs::Vector{Float64}
	pos::Vector{NTuple{2,Float64}}
	min_obs_num::Int

	function LSInitializer(m::SearchDomain)
		return LSInitializer(m.length)
	end
	function LSInitializer(length::Float64)
		S = 3e3*eye(2)
		obs = Array(Float64,0)
		pos = Array(NTuple{2,Float64},0)
		min_obs_num = 5
		return new(length, zeros(2), S, obs, pos, min_obs_num)
	end
end
function initialize(lsi::LSInitializer, x::Vehicle, o::Float64)
	push!(lsi.pos, (x.x,x.y))
	push!(lsi.obs, o)
	if length(lsi.obs) < lsi.min_obs_num
		return false
	end
	# ok, do our triangulation/least squares stuff
	best_sum = Inf
	for xj in linspace(0.0, lsi.length)
		for yj in linspace(0.0, lsi.length)
			square_sum = 0.0
			for i = 1:length(lsi.pos)
				ot = true_bearing(lsi.pos[i], (xj,yj))
				square_sum += fit_180(ot - lsi.obs[i])^2
			end
			if square_sum < best_sum
				best_sum = square_sum
				lsi.mu[1] = xj
				lsi.mu[2] = yj
			end
		end
	end
	return true
end

function check_initialization(gf::GaussianFilter, x::Vehicle, o::Float64)
	gf.initialized = initialize(gf.initializer, x, o)
	if gf.initialized
		gf.mu = gf.initializer.mu
		gf.Sigma = gf.initializer.Sigma
	end
	return gf.initialized
end


text:
plt = Observable{Any}()
ui = dom"div"(plt)
# Ploting recon
map!(t-> begin
        @manipulate for slice = 1:size(image,3)
            aux = abs.(image) * prod(size(image)[1:2])
            plot_image(aux[:,:,slice],zmin=minimum(aux[:]),zmax=maximum(aux[:]);darkmode,title="Reconstruction ($slice/$(size(image,3)))")
        end
    end
    , plt, img_obs)
content!(w, "div#content", ui)


text:
export pplot, pplot!, yplot, yplot!

PLOTSIZE = 400
ALPHA = 0.4

pointsize(x, sz=1) = 44 .* sqrt.(sz .* weights(x)) .+ 1

using RecipesBase

@recipe function f(x::Weighted{<:Matrix, <:Vector}; sz=1) ## exlude both 1D and Tracked
    if size(x,1)==1
        out = Weighted(vec(x.array), x.weights, x.opt) ## just to use 1D type signature
    elseif size(x,1)>=4
        out = wPCA(x,2)(x)
    else

        legend --> :false
        lab    --> " "
        if x.opt.clamp
            extr = [x.opt.lo - 0.05, x.opt.hi + 0.05]
            xlim --> extr
            ylim --> extr
            # zlim --> [x.opt.lo, x.opt.hi] ## leaves too much space, in GR
        end

        if endswith(x.opt.aname,"PC")
            size   --> (1.5*PLOTSIZE, PLOTSIZE)
            ylim   --> pcaylim(x, 1.5)
        elseif size(x,1)==2
            size   --> (1.1*PLOTSIZE, PLOTSIZE)
        elseif size(x,1)==3
            size   --> (1.3*PLOTSIZE, 1.2*PLOTSIZE)
        end

        seriestype := :scatter
        markeralpha --> ALPHA ## nickname alpha=0.2 doesn't overwrite this
        #markerstrokewidth --> 0
        markersize := pointsize(x, sz)

        if size(x,1)==2
            xaxis --> grsafe(x.opt.aname)*"_1"
            yaxis --> grsafe(x.opt.aname)*"_2"

            out = array(x)[1,:], array(x)[2,:]

        elseif size(x,1)==3
            xaxis --> xaxis_grbug(x.opt.aname)
            yaxis --> yaxis_grbug(x.opt.aname)
            zaxis --> grsafe(x.opt.aname)*"_3"

            zcolor --> array(x)[3,:]

            out = array(x)[1,:], array(x)[2,:], array(x)[3,:]
        end
    end
    out
end

@recipe f(x::Weighted{<:Base.ReshapedArray}; sz=1) = copy(x) ## these were otherwise not caught
@recipe f(x::Weighted{<:Base.ReinterpretArray}; sz=1) = copy(x)

# using Requires
# @require Plots begin

    ## Hack to work around bug:  https://github.com/JuliaPlots/Plots.jl/issues/743
    function xaxis_grbug(s::Union{Symbol,String})
        # Plots.backend_name() != :gr && return string(s)*"_1"
        return "        "*grsafe(s)*"\\_1                                       "*grsafe(s)*"\\_2"
    end
    function yaxis_grbug(s::Union{Symbol,String})
        # Plots.backend_name() != :gr && return string(s)*"_2"
        return grsafe(s)*"\\_3"
    end

    ## GR doesn't like unicode, but likes latex style greek
    function grsafe(s::Union{Symbol,String})
        # Plots.backend_name() != :gr && return string(s)
        return reduce(*, grsafe(c) for c in string(s))
    end

# end

function grsafe(c::Char)
    c=='θ' && return "\\theta"
    c=='ϕ' && return "\\phi"
    c=='σ' && return "\\sigma"
    c=='π' && return "\\pi"
    c=='_' && return "\\_"
    c=='(' && return "\\("
    c==')' && return "\\)"
    string(c)
end

function pcaylim(x::Weighted, ratio=1.5)
    mat = array(x)

    ex1 = [extrema(mat[1,:])...]
    diff1 = ex1[2] - ex1[1]

    ex2 = [extrema(mat[2,:])...]
    diff2 = ex2[2] - ex2[1]

    if diff1 > ratio * diff2
        ex2 .*= ( diff1 / (ratio * diff2) ) ## expand the yilm
    end

    return ex2
end

## 1D with shadow
@recipe function f(x::Weighted{<:Vector}, shadow="yes"; sz=0.8) ## sz doesn't get passed to here :(
    size   --> (1.2*PLOTSIZE, 0.8*PLOTSIZE)
    legend --> :false
    if x.opt.clamp && x.opt.hi < Inf
        xlim --> [x.opt.lo - 0.05, x.opt.hi + 0.05]
    end
    xaxis --> grsafe(x.opt.aname) #"\\theta"
    yaxis --> grsafe(x.opt.wname)
    ylim  --> [0, 1.4*maximum(weights(x)) ]
    yticks --> [0, round(maximum(weights(x)),digits=2)]

    if shadow != "only"
        @series begin ## plot the points
            seriestype := :scatter
            markersize := pointsize(x, sz)
            markeralpha --> ALPHA
            # markerstrokewidth --> 0

            array(x), weights(x)
        end
    end
    if shadow != "no"
        @series begin ## plot the shadow
            seriestype := :line
            lab   := ""
            fill  := 0
            color := :black
            alpha := 0.5*ALPHA

            shadowxy(array(x), weights(x))
        end
    end
end

function shadowxy(x::Vector, y::Vector, smooth=(maximum(x)-minimum(x))/100, numpoints=321; fixmax=true)
    @assert length(x)==length(y)
    xlist = range(minimum(x)-10*smooth, stop=maximum(x)+10*smooth, length=numpoints)
    ylist = zeros(numpoints)
    invtst = 1/(2*smooth^2)
    for i=1:numpoints
        for j=1:length(x)
            ylist[i] += y[j] * exp( -(xlist[i]-x[j])^2 * invtst )
        end
    end
    scale = fixmax ? maximum(y) / maximum(ylist) : 1
    return xlist, scale .* ylist
end

@recipe function f(x::Weighted, fun::Function)
    unique!(fun(x))
end
@recipe function f(fun::Function, x::Weighted)
    unique!(fun(x))
end

@recipe function f(x::Weighted, λ::Number)
    rmul!(copy(x), λ)
end

using Requires

function __init__()
@require Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" begin
    # Am not sure whether these can easily be made into recipies.
    # import WeightedArrays: pplot, pplot!, yplot, yplot!, rPCA, sPCA, rmul!

    """
        pplot(x::Weighted)
        pplot!(x)
    The first calls `sPCA`, the second calls `rPCA` to plot on the same axes.

        pplot(x, f)
        pplot(x, λ)
    Ditto, but first applies function `f`, or scales by `λ`.
    """
    pplot(x::Weighted, f::Function=identity; kw...) = Plots.plot(sPCA(unique!(f(x))); kw...)
    pplot(x::Weighted, λ::Number; kw...) = Plots.plot(rmul!(sPCA(x),λ); kw...)

    pplot!(x::Weighted, f::Function=identity; kw...) = Plots.plot!(rPCA(unique!(f(x))); kw...)
    pplot!(x::Weighted, λ::Number; kw...) = Plots.plot!(rmul!(rPCA(x),λ); kw...)

    pplot(f::Function, x::Weighted; kw...) = pplot(x, f; kw...) # standard order
    pplot!(f::Function, x::Weighted; kw...) = pplot!(x, f; kw...)

    """
        yplot(y, x::Weighted, t)

    Time-series plotting. If `y(Π, 1:3)` is what we observe (up to noise)
    then `yplot(y, Π, 1:3)` will plot a line for each column of `Π`, for times `range(0, tmax, length=200)` points; `pts=200` and `tmax=3.3` are keywords.
    Their opacity by default represents weight, unless you override with e.g. `alpha=1/20`.
    """
    function yplot(args...; kw...)
        Plots.plot(xaxis = "time")#, yaxis="y")
        yplot!(args...; kw...)
    end

    function yplot!(yfun::Function, prior::Weighted, times::AbstractVector;
            pts=200, c=:blue, m=:uptriangle, alpha=:auto, tmax=1.1*maximum(times), lab="", kw...)
        if alpha == :auto
            alpha = sqrt.(prior.weights ./ maximum(prior.weights)) |> transpose
        end
        ptimes = range(0, tmax, length=pts)
        res = yfun(prior, ptimes).array

        # case of y_react etc, when y(θ) should be a matrix
        if :syms in fieldnames(typeof(yfun))
            dout, k = size(res)
            tensor = reshape(res, length(yfun.syms), :,k)
            # alpha = alpha ./ 2
            for (i,s) in enumerate(yfun.syms)
                Plots.plot!(0:0, 0:0, lab=lab * string(s), c=i, alpha=1, kw...)
                Plots.plot!(ptimes, tensor[i, :, :], lab="", c=i, alpha=alpha, kw...)
            end

        # simple yexp etc.
        else
            Plots.plot!(0:0, 0:0, lab=lab, c=c, alpha=1, kw...)
            Plots.plot!(ptimes, res, lab="", c=c, alpha=alpha, kw...)
        end

        Plots.scatter!(times, zero(times), c=:grey, m=m, lab="")
    end

end # @require
end #__init__


text:
using TmuxDisplays
using Test

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


text:
using RCall,PyCall,DataFrames
R"library(lhs)"
len = length

rows = 10

constrain= readtable("constrained_ics.csv",separator = ',')
changing = readtable("changing_ics.csv",separator = ',')
columns = len(changing[1])
[constrain[Symbol("x$i")] = constrain[:x1] for i in 1:rows]
R" lhsgrid = optimumLHS(n=$(columns), k=$(rows), maxSweeps=20, eps=.1, verbose=FALSE)"
@rget lhsgrid
#concat df vcat or hcat  # .operation = appply to all elements
lhsgrid =  ((changing[:max]-changing[:min]) .* lhsgrid) .+ changing[:min]
lhsgrid =  DataFrame(Array(hcat( changing[:species],lhsgrid)))
indexarray = [Symbol("x$(i)") for i in 0:rows]
indexarray[1] = "species"
names!(lhsgrid,indexarray)
ics = vcat(constrain,lhsgrid)

writetable("latin.csv", ics, separator = ',', header = false,quotemark=' ')

#changing[:species]
print("ahhh",ics)


text:
"""
	anyoverlapping(a::GenomicFeatures.Interval{T},
	b::GenomicFeatures.IntervalCollection{S})

Extend GenomicFeatures.isoverlapping to linearly check if interval a.
overlaps collection. Return true or false

```jldoctest
using GenomicFeatures
a = GenomicFeatures.Interval("chr1", 0, 10)
b = GenomicFeatures.IntervalCollection([
	GenomicFeatures.Interval("chr1", 5, 15),
	GenomicFeatures.Interval("chr1", 20, 25)
	])
anyoverlapping(a, b)

# output

true
```

See Also [`GenomePermutations.alloverlapping`](@ref), [`GenomePermutations.countoverlapping`](@ref).
"""
function anyoverlapping(a::GenomicFeatures.Interval{T}, b::GenomicFeatures.IntervalCollection{S}) where {T, S}
	for interval in b
		GenomicFeatures.isoverlapping(a, interval) && return(true::Bool)
	end
	return(false::Bool)
end

"""
	alloverlapping(a::GenomicFeatures.IntervalCollection{T}, b::GenomicFeatures.IntervalCollection{S})

Linearly check if all intervals in collection a overlap collection b.
Return true or false

```jldoctest
using GenomicFeatures
a = GenomicFeatures.IntervalCollection([
	GenomicFeatures.Interval("chr1", 14, 17),
	GenomicFeatures.Interval("chr1", 20, 25)
	])
b = GenomicFeatures.IntervalCollection([
	GenomicFeatures.Interval("chr1", 5, 10),
	GenomicFeatures.Interval("chr1", 20, 25)
	])
alloverlapping(a, b)

# output

false
```

See Also [`GenomePermutations.anyoverlapping`](@ref), [`GenomePermutations.countoverlapping`](@ref).
"""
function alloverlapping(a::GenomicFeatures.IntervalCollection{T}, b::GenomicFeatures.IntervalCollection{S}) where {T,S}
	for interval in a
		anyoverlapping(interval, b) || return(false::Bool)
	end
	return(true::Bool)
end

"""
	countoverlapping(a::GenomicFeatures.IntervalCollection{T} , b::GenomicFeatures.IntervalCollection{T})

Linearly count how many intervals in collection a overlap with collection b.

```jldoctest

using GenomicFeatures
a = GenomicFeatures.IntervalCollection([
	GenomicFeatures.Interval("chr1", 14, 17),
	GenomicFeatures.Interval("chr1", 20, 50)
	])
b = GenomicFeatures.IntervalCollection([
	GenomicFeatures.Interval("chr1", 5, 10),
	GenomicFeatures.Interval("chr1", 25, 35 ),
	GenomicFeatures.Interval("chr1", 40, 65)
	])

[countoverlapping(a, b), countoverlapping(b, a)]

# output

2-element Vector{Int64}:
 1
 2
```

#Note
As seen above, countoverlapping(a,b) != countoverlapping(b,a) as one interval in a
can overlap with >1 intervals in b.

See Also [`GenomePermutations.anyoverlapping`](@ref), [`GenomePermutations.alloverlapping`](@ref).
"""
function countoverlapping(
	a::GenomicFeatures.IntervalCollection{T},
	b::GenomicFeatures.IntervalCollection{S}
) where {T,S}

	overlaps::Int = 0
	for interval in a
		anyoverlapping(interval, b) && (overlaps += 1)
	end
	return(overlaps)
end

"""
	isin(a::GenomicFeatures.Interval{S}, GenomicFeatures.Interval{T})

Check if interval a is fully contained in interval b.

```jldoctest
using GenomicFeatures
a = GenomicFeatures.Interval("chr1", 5, 10)
b = GenomicFeatures.Interval("chr1", 1, 15)

isin(a, b)

# output

true
```

```jldoctest
using GenomicFeatures
a = GenomicFeatures.Interval("chr1", 20, 30)
b = GenomicFeatures.Interval("chr1", 1, 15)

isin(a, b)

# output

false
```

```jldoctest
using GenomicFeatures
a = GenomicFeatures.Interval("chr1", 10, 20)
b = GenomicFeatures.Interval("chr1", 1, 15)

isin(a, b)

# output

false
```
"""
function isin(a::GenomicFeatures.Interval{S},b::GenomicFeatures.Interval{T}) where {S,T}
	return (a.seqname == b.seqname && a.first >= b.first && a.last <= b.last)
end

"""
	isin(a::GenomicFeatures.Interval{T}, b::GenomicFeatures.IntervalCollection{S})

Linearly check if interval a is fully contained in collection b.
"""
function isin(a::GenomicFeatures.Interval{S},
				 b::GenomicFeatures.IntervalCollection{T}) where {S,T}

	for interval in b
		isin(a, interval) && return(true::Bool)
	end
	return(false::Bool)
end

"""
	isin(a:::GenomicFeatures.IntervalCollection{T}, b:::GenomicFeatures.IntervalCollection{S})

Linearly check if all intervals in collection a are contained in collection b.

# Note

At this point there are 2 layers of linear search so the time complexity is n^2.
"""
function isin(a::GenomicFeatures.IntervalCollection{T},
				 b::GenomicFeatures.IntervalCollection{S}) where {T,S}

	for interval in a
		isin(interval, b) || return(false::Bool)
	end
	return(true::Bool)
end


text:
CuMesh(::MetaGraph)::CuMesh
CuMesh(::MetaDiGraph)::CuMesh

MetaGraph(::CuMesh)::MetaGraph
MetaDiGraph(::CuMesh)::MetaDiGraph


text:
function global_align(x, y;
                      nsplits=2, initblock=nothing,
                      spacefun=UncertaintyRegion, objfun=alignment_objective, tformfun=AffineMap,
                      atol=0.1, rtol=0, maxblocks=5e8, maxsplits=Inf, maxevals=Inf, maxstagnant=Inf, threads=false)
    if isodd(nsplits)
        throw(ArgumentError("`nsplits` must be even"))
    end
    if dims(gmmx) != dims(gmmy)
        throw(ArgumentError("Dimensionality of the point sets must be equal"))
    end

    # initialization
    if isnothing(initblock)

        initblock = blockfun()
    end
    ndims = dims(initblock)
    ub, bestloc = initblock.upperbound, initblock.center    # best-so-far objective value and transformation
    t = promote_type(numbertype(gmmx), numbertype(gmmy))
    lb = typemin(t)
    pq = PriorityQueue{Block{ndims, t}, Tuple{t,t}}()
    enqueue!(pq, initblock, (initblock.lowerbound, initblock.upperbound))

    # split cubes until convergence
    ndivisions = 0
    sinceimprove = 0
    evalsperdiv = length(gmmx)*length(gmmy)*nsplits^ndims
    while !isempty(pq)
        if (length(pq) > maxblocks) || (ndivisions*evalsperdiv > maxevals) || (sinceimprove > maxstagnant) || (ndivisions > maxsplits)
            break
        end
        ndivisions += 1
        sinceimprove += 1

        # take the block with the lowest lower bound
        bl, (lb, blub) = dequeue_pair!(pq)

        # if the best solution so far is close enough to the best possible solution, end
        if abs((ub - lb)/lb) < rtol || abs(ub-lb) < atol
            return GMAlignmentResult(gmmx, gmmy, ub, lb, tformfun(bestloc...), bestloc, ndivisions*evalsperdiv, ndivisions, length(pq), sinceimprove)
        end

        # split up the block into `nsplits` smaller blocks across each dimension
        subrngs = subranges(bl.ranges, nsplits)
        sblks = fill(Block{ndims,t}(), nsplits^ndims)
        for i=1:length(subrngs)
            sblks[i] = blockfun(gmmx, gmmy, subrngs[i], pσ, pϕ, rot, trl)
        end

        # reset the upper bound if appropriate
        minub, ubidx = findmin([sblk.upperbound for sblk in sblks])
        if minub < ub
            ub, bestloc = local_align(gmmx, gmmy, sblks[ubidx], pσ, pϕ; objfun=objfun, rot=rot, trl=trl)
            sinceimprove = 0
        end

        # only add sub-blocks to the queue if they present possibility for improvement
        for sblk in sblks
            if sblk.lowerbound < ub
                enqueue!(pq, sblk, (sblk.lowerbound, sblk.upperbound))
            end
        end
    end
    if isempty(pq)
        return GMAlignmentResult(gmmx, gmmy, ub, lb, tformfun(bestloc...), bestloc, ndivisions*evalsperdiv, ndivisions, length(pq), sinceimprove)
    else
        return GMAlignmentResult(gmmx, gmmy, ub, dequeue_pair!(pq)[2][1], tformfun(bestloc...), bestloc, ndivisions*evalsperdiv, ndivisions, length(pq), sinceimprove)
    end
end