text:
const SIMULATION_SERIALIZATION_FILENAME = "simulation.bin"
const SIMULATION_LOG_FILENAME = "simulation.log"
const REQUIRED_RECORDERS = [:simulation_status]

mutable struct SimulationInternal
    logs_dir::String
    models_dir::String
    raw_dir::String
    recorder_dir::String
    results_dir::String
    stages_count::Int
    run_count::Dict{Int, Dict{Int, Int}}
    date_ref::Dict{Int, Dates.DateTime}
    date_range::NTuple{2, Dates.DateTime} #Inital Time of the first forecast and Inital Time of the last forecast
    current_time::Dates.DateTime
    reset::Bool
    compiled_status::BUILD_STATUS
    simulation_cache::Dict{<:CacheKey, AbstractCache}
    recorders::Vector{Symbol}
    console_level::Base.CoreLogging.LogLevel
    file_level::Base.CoreLogging.LogLevel
end

function SimulationInternal(
    steps::Int,
    stages_keys::Base.KeySet,
    sim_dir,
    name;
    output_dir = nothing,
    recorders = [],
    console_level = Logging.Error,
    file_level = Logging.Info,
)
    count_dict = Dict{Int, Dict{Int, Int}}()

    for s in 1:steps
        count_dict[s] = Dict{Int, Int}()
        for st in stages_keys
            count_dict[s][st] = 0
        end
    end

    base_dir = joinpath(sim_dir, name)
    mkpath(base_dir)

    output_dir = _get_output_dir_name(base_dir, output_dir)
    simulation_dir = joinpath(base_dir, output_dir)
    if isdir(simulation_dir)
        error("$simulation_dir already exists. Delete it or pass a different output_dir.")
    end

    logs_dir = joinpath(simulation_dir, "logs")
    models_dir = joinpath(simulation_dir, "models_json")
    raw_dir = joinpath(simulation_dir, "raw_output")
    recorder_dir = joinpath(simulation_dir, "recorder")
    results_dir = joinpath(simulation_dir, "results")

    for path in (simulation_dir, logs_dir, models_dir, raw_dir, recorder_dir, results_dir)
        mkpath(path)
    end

    unique_recorders = Set(REQUIRED_RECORDERS)
    foreach(x -> push!(unique_recorders, x), recorders)

    init_time = Dates.now()
    return SimulationInternal(
        logs_dir,
        models_dir,
        raw_dir,
        recorder_dir,
        results_dir,
        length(stages_keys),
        count_dict,
        Dict{Int, Dates.DateTime}(),
        (init_time, init_time),
        init_time,
        true,
        false,
        Dict{CacheKey, AbstractCache}(),
        collect(unique_recorders),
        console_level,
        file_level,
    )
end

function configure_logging(internal::SimulationInternal, file_mode)
    return IS.configure_logging(
        console = true,
        console_stream = stderr,
        console_level = internal.console_level,
        file = true,
        filename = joinpath(internal.logs_dir, SIMULATION_LOG_FILENAME),
        file_level = internal.file_level,
        file_mode = file_mode,
        tracker = nothing,
        set_global = false,
    )
end

function register_recorders!(internal::SimulationInternal, file_mode)
    for name in internal.recorders
        IS.register_recorder!(name; mode = file_mode, directory = internal.recorder_dir)
    end
end

function unregister_recorders!(internal::SimulationInternal)
    for name in internal.recorders
        IS.unregister_recorder!(name)
    end
end

function _add_initial_condition_caches(
    sim::SimulationInternal,
    stage::Stage,
    caches::Union{Nothing, Vector{<:AbstractCache}},
)
    initial_conditions = stage.internal.psi_container.initial_conditions
    for (ic_key, init_conds) in initial_conditions.data
        _create_cache(ic_key, caches)
    end
    return
end

function _create_cache(ic_key::ICKey, caches::Union{Nothing, Vector{<:AbstractCache}})
    return
end

function _create_cache(
    ic_key::ICKey{TimeDurationON, T},
    caches::Union{Nothing, Vector{<:AbstractCache}},
) where {T <: PSY.Device}

    cache_keys = CacheKey.(caches)
    if isempty(cache_keys) || !in(CacheKey(TimeStatusChange, T), cache_keys)
        cache = TimeStatusChange(T, ON)
        push!(caches, cache)
    end
    return
end

function _create_cache(
    ic_key::ICKey{TimeDurationOFF, T},
    caches::Vector{<:AbstractCache},
) where {T <: PSY.Device}

    cache_keys = CacheKey.(caches)
    if isempty(cache_keys) || !in(CacheKey(TimeStatusChange, T), cache_keys)
        cache = TimeStatusChange(T, ON)
        push!(caches, cache)
    end
    return
end

function _create_cache(
    ic_key::ICKey{EnergyLevel, T},
    caches::Vector{<:AbstractCache},
) where {T <: PSY.Device}

    cache_keys = CacheKey.(caches)
    if isempty(cache_keys) || !in(CacheKey(StoredEnergy, T), cache_keys)
        cache = StoredEnergy(T, ENERGY)
        push!(caches, cache)
    end
    return
end

function _set_internal_caches(
    internal::SimulationInternal,
    stage::Stage,
    caches::Vector{<:AbstractCache},
)
    for c in caches
        cache_key = CacheKey(c)
        push!(stage.internal.caches, cache_key)
        if !haskey(internal.simulation_cache, cache_key)
            @debug "Cache $(cache_key) added to the simulation"
            internal.simulation_cache[cache_key] = c
        end
        internal.simulation_cache[cache_key].value = get_initial_cache(c, stage)
    end
    return
end

function _get_output_dir_name(path, output_dir)
    if !isnothing(output_dir)
        # The user wants a custom name.
        return output_dir
    end

    # Return the next highest integer.
    output_dir = 1
    for name in readdir(path)
        if occursin(r"^\d+$", name)
            num = parse(Int, name)
            if num >= output_dir
                output_dir = num + 1
            end
        end
    end

    return string(output_dir)
end

# TODO: Add DocString
@doc raw"""
    Simulation(steps::Int
                stages::Dict{String, Stage{<:AbstractOperationsProblem}}
                sequence::Union{Nothing, SimulationSequence}
                simulation_folder::String
                name::String
                internal::Union{Nothing, SimulationInternal}
                )

"""
mutable struct Simulation
    steps::Int
    stages::Dict{String, Stage{<:AbstractOperationsProblem}}
    initial_time::Union{Nothing, Dates.DateTime}
    sequence::SimulationSequence
    simulation_folder::String
    name::String
    internal::Union{Nothing, SimulationInternal}

    function Simulation(;
        stages_sequence::SimulationSequence,
        name::String,
        steps::Int,
        stages = Dict{String, Stage{AbstractOperationsProblem}}(),
        simulation_folder::String,
        initial_time = nothing,
    )
        new(steps, stages, initial_time, stages_sequence, simulation_folder, name, nothing)
    end
end

"""
    Simulation(directory::AbstractString)

Constructs Simulation from a serialized directory. Callers should pass any kwargs here that
they passed to the original Simulation.

# Arguments
- `directory::AbstractString`: the directory returned from the call to serialize
- `stage_info::Dict`: Two-level dictionary containing stage parameters that cannot be
  serialized. The outer dict should be keyed by the stage name. The inner dict must contain
  'optimizer' and may contain 'jump_model'. These should be the same values used for the
  original simulation.
"""
function Simulation(directory::AbstractString, stage_info::Dict)
    obj = deserialize(Simulation, directory, stage_info)
end

################# accessor functions ####################
get_initial_time(s::Simulation) = s.initial_time
get_sequence(s::Simulation) = s.sequence
get_steps(s::Simulation) = s.steps
get_date_range(s::Simulation) = s.internal.date_range

function get_base_powers(s::Simulation)
    base_powers = Dict()
    for (k, v) in s.stages
        base_powers[k] = PSY.get_base_power(v.sys)
    end
    return base_powers
end

function get_stage(s::Simulation, name::String)
    stage = get(s.stages, name, nothing)
    isnothing(stage) && throw(ArgumentError("Stage $(name) not present in the simulation"))
    return stage
end

get_stage_interval(s::Simulation, name::String) = get_stage_interval(s.sequence, name)

function get_stage(s::Simulation, number::Int)
    name = get(s.sequence.order, number, nothing)
    isnothing(name) && throw(ArgumentError("Stage with $(number) not defined"))
    return get_stage(s, name)
end

get_stages_quantity(s::Simulation) = s.internal.stages_count

function get_simulation_time(s::Simulation, stage_number::Int)
    return s.internal.date_ref[stage_number]
end

get_ini_cond_chronology(s::Simulation) = s.sequence.ini_cond_chronology
get_stage_name(s::Simulation, stage::Stage) = get_stage_name(s.sequence, stage)
IS.get_name(s::Simulation) = s.name
get_simulation_folder(s::Simulation) = s.simulation_folder
get_execution_order(s::Simulation) = s.sequence.execution_order
get_current_execution_index(s::Simulation) = s.sequence.current_execution_index
get_logs_folder(s::Simulation) = s.internal.logs_dir
get_recorder_folder(s::Simulation) = s.internal.recorder_dir

function get_stage_cache_definition(s::Simulation, stage::String)
    caches = s.sequence.cache
    cache_ref = Array{AbstractCache, 1}()
    for stage_names in keys(caches)
        if stage in stage_names
            push!(cache_ref, caches[stage_names])
        end
    end
    return cache_ref
end

function get_cache(simulation_cache::Dict{<:CacheKey, AbstractCache}, key::CacheKey)
    c = get(simulation_cache, key, nothing)
    isnothing(c) && @debug("Cache with key $(key) not present in the simulation")
    return c
end

function get_cache(
    simulation_cache::Dict{<:CacheKey, AbstractCache},
    ::Type{T},
    ::Type{D},
) where {T <: AbstractCache, D <: PSY.Device}
    return get_cache(simulation_cache, CacheKey(T, D))
end

function get_cache(
    s::Simulation,
    ::Type{T},
    ::Type{D},
) where {T <: AbstractCache, D <: PSY.Device}
    return get_cache(s.internal.simulation_cache, CacheKey(T, D))
end

function _check_forecasts_sequence(sim::Simulation)
    for (stage_number, stage_name) in sim.sequence.order
        stage = get_stage(sim, stage_name)
        resolution = get_resolution(stage)
        horizon = get_stage_horizon(get_sequence(sim), stage_name)
        interval = get_stage_interval(get_sequence(sim), stage_name)
        horizon_time = resolution * horizon
        if horizon_time < interval
            throw(IS.ConflictingInputsError("horizon ($horizon_time) is
                                shorter than interval ($interval) for $stage_name"))
        end
    end
end

function _check_feedforward_chronologies(sim::Simulation)
    if isempty(sim.sequence.feedforward_chronologies)
        @info("No Feedforward Chronologies defined")
    end
    for (key, chron) in sim.sequence.feedforward_chronologies
        check_chronology!(sim, key, chron)
    end
    return
end

function _assign_feedforward_chronologies(sim::Simulation)
    for (key, chron) in sim.sequence.feedforward_chronologies
        destination_stage = get_stage(sim, key.second)
        destination_stage_interval =
            IS.time_period_conversion(get_stage_interval(sim, key.second))
        source_stage_number = find_key_with_value(sim.sequence.order, key.first)
        if isempty(source_stage_number)
            throw(ArgumentError("Stage $(key.first) not specified in the order dictionary"))
        end
        for stage_number in source_stage_number
            destination_stage.internal.chronolgy_dict[stage_number] = chron
            source_stage = get_stage(sim, stage_number)
            source_stage_resolution =
                IS.time_period_conversion(PSY.get_forecasts_resolution(source_stage.sys))
            # This line keeps track of the executions of a stage relative to other stages.
            # This might be needed in the future to run multiple stages. For now it is disabled
            #destination_stage.internal.synchronized_executions[stage_number] =
            #Int(source_stage_resolution / destination_stage_interval)
        end
    end
    return
end

function _get_simulation_initial_times!(sim::Simulation)
    k = keys(sim.sequence.order)
    k_size = length(k)
    @assert k_size == maximum(k)

    stage_initial_times = Dict{Int, Vector{Dates.DateTime}}()
    time_range = Vector{Dates.DateTime}(undef, 2)
    sim_ini_time = get_initial_time(sim)
    for (stage_number, stage_name) in sim.sequence.order
        stage_system = sim.stages[stage_name].sys
        PSY.check_forecast_consistency(stage_system)
        interval = PSY.get_forecasts_interval(stage_system)
        horizon = get_stage_horizon(get_sequence(sim), stage_name)
        seq_interval = get_stage_interval(get_sequence(sim), stage_name)
        if PSY.are_forecasts_contiguous(stage_system)
            stage_initial_times[stage_number] =
                PSY.generate_initial_times(stage_system, seq_interval, horizon)
            if isempty(stage_initial_times[stage_number])
                throw(IS.ConflictingInputsError("Simulation interval ($seq_interval) and
                        forecast interval ($interval) definitions are not compatible"))
            end
        else
            stage_initial_times[stage_number] = PSY.get_forecast_initial_times(stage_system)
            interval = PSY.get_forecasts_interval(stage_system)
            if interval != seq_interval
                throw(IS.ConflictingInputsError("Simulation interval ($seq_interval) and
                        forecast interval ($interval) definitions are not compatible"))
            end
            for (ix, element) in enumerate(stage_initial_times[stage_number][1:(end - 1)])
                if !(element + interval == stage_initial_times[stage_number][ix + 1])
                    throw(IS.ConflictingInputsError("The sequence of forecasts is invalid"))
                end
            end
        end
        if !isnothing(sim_ini_time) &&
           !mapreduce(x -> x == sim_ini_time, |, stage_initial_times[stage_number])
            throw(IS.ConflictingInputsError("The specified simulation initial_time $sim_ini_time isn't contained in stage $stage_number.
            Manually provided initial times have to be compatible with the specified interval and horizon in the stages."))
        end
        stage_number == 1 && (time_range[1] = stage_initial_times[stage_number][1])
        (
            stage_number == k_size &&
            (time_range[end] = stage_initial_times[stage_number][end])
        )
    end
    sim.internal.date_range = Tuple(time_range)

    if isnothing(get_initial_time(sim))
        sim.initial_time = stage_initial_times[1][1]
        @debug("Initial Simulation Time will be infered from the data.
               Initial Simulation Time set to $(sim.initial_time)")
    end

    sim.internal.current_time = sim.initial_time
    return stage_initial_times
end

function _attach_feedforward!(sim::Simulation, stage_name::String)
    stage = get(sim.stages, stage_name, nothing)
    feedforward = filter(p -> (p.first[1] == stage_name), sim.sequence.feedforward)
    for (key, ff) in feedforward
        # Note: key[1] = Stage name, key[2] = template field name, key[3] = device model key
        field_dict = getfield(stage.template, key[2])
        device_model = get(field_dict, key[3], nothing)
        isnothing(device_model) &&
            throw(IS.ConflictingInputsError("Device model $(key[3]) not found in stage $(stage_name)"))
        device_model.feedforward = ff
    end
    return
end

function _check_steps(
    sim::Simulation,
    stage_initial_times::Dict{Int, Vector{Dates.DateTime}},
)
    for (stage_number, stage_name) in sim.sequence.order
        stage = sim.stages[stage_name]
        execution_counts = get_executions(stage)
        @assert length(findall(x -> x == stage_number, sim.sequence.execution_order)) ==
                execution_counts
        forecast_count = length(stage_initial_times[stage_number])
        if get_steps(sim) * execution_counts > forecast_count
            throw(IS.ConflictingInputsError("The number of available time series ($(forecast_count)) is not enough to perform the
            desired amount of simulation steps ($(sim.steps*stage.internal.execution_count))."))
        end
    end
    return
end

function _check_required_ini_cond_caches(sim::Simulation)
    for (stage_number, stage_name) in sim.sequence.order
        stage = get_stage(sim, stage_name)
        for (k, v) in iterate_initial_conditions(stage.internal.psi_container)
            # No cache needed for the initial condition -> continue
            isnothing(v[1].cache_type) && continue
            c = get_cache(sim, v[1].cache_type, k.device_type)
            if isnothing(c)
                throw(ArgumentError("Cache $(v[1].cache_type) not defined for initial condition $(k) in stage $stage_name"))
            end
            @debug "found cache $(v[1].cache_type) for initial condition $(k) in stage $(stage_name)"
        end
    end
    return
end

function _populate_caches!(sim::Simulation, stage_name::String)
    caches = get_stage_cache_definition(sim, stage_name)
    stage = get_stage(sim, stage_name)
    _add_initial_condition_caches(sim.internal, stage, caches)
    _set_internal_caches(sim.internal, stage, caches)
    return
end

function _build_stages!(sim::Simulation)
    for (stage_number, stage_name) in sim.sequence.order
        TimerOutputs.@timeit BUILD_SIMULATION_TIMER "Build Stage $(stage_name)" begin
            @info("Building Stage $(stage_number)-$(stage_name)")
            horizon = get_stage_horizon(get_sequence(sim), stage_name)
            stage = get_stage(sim, stage_name)
            stage_interval = get_stage_interval(get_sequence(sim), stage_name)
            initial_time = get_initial_time(sim)
            stage.internal.stage_path =
                joinpath(sim.internal.models_dir, "stage_$(stage_name)_model")
            mkpath(stage.internal.stage_path)
            build!(stage, initial_time, horizon, stage_interval)
            _populate_caches!(sim, stage_name)
            sim.internal.date_ref[stage_number] = initial_time
        end
    end
    _check_required_ini_cond_caches(sim)
    return
end

function _check_folder(sim::Simulation)
    folder = get_simulation_folder(sim)
    !isdir(folder) && throw(IS.ConflictingInputsError("Specified folder is not valid"))
    try
        mkdir(joinpath(folder, "fake"))
        rm(joinpath(folder, "fake"))
    catch e
        throw(IS.ConflictingInputsError("Specified folder does not have write access [$e]"))
    end
end

"""
    build!(sim::Simulation)

Build the Simulation and all stages.

# Arguments
- `sim::Simulation`: simulation object
- `output_dir = nothing`: If nothing then generate a unique name.
- `recorders::Vector{Symbol} = []`: recorder names to register
- `console_level = Logging.Error`:
- `file_level = Logging.Info`:

Throws an exception if label is passed and the directory already exists.
"""
function build!(
    sim::Simulation;
    output_dir = nothing,
    recorders = [],
    console_level = Logging.Error,
    file_level = Logging.Info,
)
    TimerOutputs.reset_timer!(BUILD_SIMULATION_TIMER)
    TimerOutputs.@timeit BUILD_SIMULATION_TIMER "Build Simulation" begin
        _check_forecasts_sequence(sim)
        _check_feedforward_chronologies(sim)
        _check_folder(sim)
        sim.internal = SimulationInternal(
            sim.steps,
            keys(sim.sequence.order),
            get_simulation_folder(sim),
            get_name(sim);
            output_dir = output_dir,
            recorders = recorders,
            console_level = console_level,
            file_level = file_level,
        )

        file_mode = "w"
        logger = configure_logging(sim.internal, file_mode)
        register_recorders!(sim.internal, file_mode)
        try
            Logging.with_logger(logger) do
                _build!(sim)
                sim.internal.compiled_status = BUILT
                @info "\n$(BUILD_SIMULATION_TIMER)\n"
            end
        finally
            unregister_recorders!(sim.internal)
            close(logger)
        end
    end
    return
end

function _build!(sim::Simulation)
    sim.internal.compiled_status = IN_PROGRESS
    stage_initial_times = _get_simulation_initial_times!(sim)
    for (stage_number, stage_name) in sim.sequence.order
        stage = get_stage(sim, stage_name)
        if isnothing(stage)
            throw(IS.ConflictingInputsError("Stage $(stage_name) not found in the stages definitions"))
        end
        stage_interval = get_stage_interval(sim, stage_name)
        stage.internal.executions = Int(get_step_resolution(sim.sequence) / stage_interval)
        stage.internal.number = stage_number
        _attach_feedforward!(sim, stage_name)
    end
    _assign_feedforward_chronologies(sim)
    _check_steps(sim, stage_initial_times)
    _build_stages!(sim)
end

#Defined here because it requires Stage and Simulation to defined

#############################Interfacing Functions##########################################
# These are the functions that the user will have to implement to update a custom IC Chron #
# or custom InitialConditionType #

""" Updates the initial conditions of the stage"""
function initial_condition_update!(
    stage::Stage,
    ini_cond_key::ICKey,
    initial_conditions::Vector{InitialCondition},
    ::IntraStageChronology,
    sim::Simulation,
)
    # TODO: Replace this convoluted way to get information with access to data store
    simulation_cache = sim.internal.simulation_cache
    for ic in initial_conditions
        name = device_name(ic)
        interval_chronology =
            get_stage_interval_chronology(sim.sequence, get_stage_name(sim, stage))
        var_value =
            get_stage_variable(interval_chronology, (stage => stage), name, ic.update_ref)
        # We pass the simulation cache instead of the whole simulation to avoid definition dependencies. All the inputs to calculate_ic_quantity are defined before the simulation object
        quantity = calculate_ic_quantity(ini_cond_key, ic, var_value, simulation_cache)
        previous_value = get_condition(ic)
        PJ.fix(ic.value, quantity)
        IS.@record :simulation InitialConditionUpdateEvent(
            sim.internal.current_time,
            ini_cond_key,
            ic,
            quantity,
            previous_value,
            get_number(stage),
        )
    end
end

""" Updates the initial conditions of the stage"""
function initial_condition_update!(
    stage::Stage,
    ini_cond_key::ICKey,
    initial_conditions::Vector{InitialCondition},
    ::InterStageChronology,
    sim::Simulation,
)
    # TODO: Replace this convoluted way to get information with access to data store
    simulation_cache = sim.internal.simulation_cache
    execution_index = get_execution_order(sim)
    execution_count = get_execution_count(stage)
    for ic in initial_conditions
        name = device_name(ic)
        current_ix = get_current_execution_index(sim)
        source_stage_ix = current_ix == 1 ? last(execution_index) : current_ix - 1
        source_stage = get_stage(sim, execution_index[source_stage_ix])
        source_stage_name = get_stage_name(sim, source_stage)

        # If the stage that ran before is lower in the order of execution the chronology needs to grab the first result as the initial condition
        if get_number(source_stage) >= get_number(stage)
            interval_chronology =
                get_stage_interval_chronology(sim.sequence, source_stage_name)
        elseif get_number(source_stage) < get_number(stage)
            interval_chronology = RecedingHorizon()
        end

        var_value = get_stage_variable(
            interval_chronology,
            (source_stage => stage),
            name,
            ic.update_ref,
        )
        quantity = calculate_ic_quantity(ini_cond_key, ic, var_value, simulation_cache)
        previous_value = get_condition(ic)
        PJ.fix(ic.value, quantity)
        IS.@record :simulation InitialConditionUpdateEvent(
            sim.internal.current_time,
            ini_cond_key,
            ic,
            quantity,
            previous_value,
            get_number(stage),
        )
    end
    return
end

function _update_caches!(sim::Simulation, stage::Stage)
    for cache in stage.internal.caches
        update_cache!(sim, cache, stage)
    end
    return
end

################################Cache Update################################################
# TODO: Need to be careful here if 2 stages modify the same cache. This function might need
# dispatch on the Statge{OpModel} to assign different actions. e.g. HAUC and DAUC
"""
TimeStatusChange Cache is updated by calculating the final status of the stage before running again. It only stores the final status and time at status once the remaining stages have been completed.
"""
function update_cache!(
    sim::Simulation,
    ::CacheKey{TimeStatusChange, D},
    stage::Stage,
) where {D <: PSY.Device}
    # TODO: Remove debug statements and use recorder here
    c = get_cache(sim, TimeStatusChange, D)
    increment = get_increment(sim, stage, c)
    variable = get_variable(stage.internal.psi_container, c.ref)
    t_range = 1:get_end_of_interval_step(stage)
    for name in variable.axes[1]
        #Store the initial condition
        c.value[name][:series] = Vector{Float64}(undef, length(t_range) + 1)
        c.value[name][:series][1] = c.value[name][:status]
        c.value[name][:current] = 2
        for t in t_range
            # Implemented this way because JuMPDenseAxisArrays doesn't support pasing a ranges Array[name, 1:n]
            device_status = JuMP.value.(variable[name, t])
            c.value[name][:series][t + 1] = device_status
            if c.value[name][:status] == device_status
                c.value[name][:count] += increment
                @debug("Cache value TimeStatus for device $name set to $device_status and count to $(c.value[name][:count])")
            else
                c.value[name][:count] = increment
                c.value[name][:status] = device_status
                @debug("Cache value TimeStatus for device $name set to $device_status and count to 1.0")
            end
        end
    end

    return
end

function get_increment(sim::Simulation, stage::Stage, cache::TimeStatusChange)
    units = cache.units
    stage_name = get_stage_name(sim, stage)
    stage_interval = IS.time_period_conversion(get_stage_interval(sim, stage_name))
    horizon = get_stage_horizon(sim.sequence, stage_name)
    stage_resolution = stage_interval / horizon
    return float(stage_resolution / units)
end

function update_cache!(
    sim::Simulation,
    ::CacheKey{StoredEnergy, D},
    stage::Stage,
) where {D <: PSY.Device}
    c = get_cache(sim, StoredEnergy, D)
    variable = get_variable(stage.internal.psi_container, c.ref)
    t = get_end_of_interval_step(stage)
    for name in variable.axes[1]
        device_energy = JuMP.value(variable[name, t])
        @debug name, device_energy
        c.value[name] = device_energy
        @debug("Cache value StoredEnergy for device $name set to $(c.value[name])")
    end

    return
end

#########################TimeSeries Data Updating###########################################
function update_parameter!(
    param_reference::UpdateRef{T},
    container::ParameterContainer,
    stage::Stage,
    sim::Simulation,
) where {T <: PSY.Component}
    components = get_available_components(T, stage.sys)
    initial_forecast_time = get_simulation_time(sim, get_number(stage))
    horizon = length(model_time_steps(stage.internal.psi_container))
    for d in components
        # RECORDER TODO: Parameter Update from forecast
        forecast = PSY.get_forecast(
            PSY.Deterministic,
            d,
            initial_forecast_time,
            get_accessor_func(param_reference),
            horizon,
        )
        ts_vector = TS.values(PSY.get_data(forecast))
        component_name = PSY.get_name(d)
        for (ix, val) in enumerate(get_parameter_array(container)[component_name, :])
            value = ts_vector[ix]
            JuMP.fix(val, value)
        end
    end

    return
end

function update_parameter!(
    param_reference::UpdateRef{T},
    container::ParameterContainer,
    stage::Stage,
    sim::Simulation,
) where {T <: PSY.Service}
    # RECORDER TODO: Parameter Update from forecast
    components = get_available_components(T, stage.sys)
    initial_forecast_time = get_simulation_time(sim, get_number(stage))
    horizon = length(model_time_steps(stage.internal.psi_container))
    param_array = get_parameter_array(container)
    for ix in axes(param_array)[1]
        service = PSY.get_component(T, stage.sys, ix)
        forecast = PSY.get_forecast(
            PSY.Deterministic,
            service,
            initial_forecast_time,
            get_accessor_func(param_reference),
            horizon,
        )
        ts_vector = TS.values(PSY.get_data(forecast))
        for (jx, value) in enumerate(ts_vector)
            JuMP.fix(get_parameter_array(container)[ix, jx], value)
        end
    end

    return
end

"""Updates the forecast parameter value"""
function update_parameter!(
    param_reference::UpdateRef{JuMP.VariableRef},
    container::ParameterContainer,
    stage::Stage,
    sim::Simulation,
)
    param_array = get_parameter_array(container)
    for (k, chronology) in stage.internal.chronolgy_dict
        source_stage = get_stage(sim, k)
        feedforward_update!(
            stage,
            source_stage,
            chronology,
            param_reference,
            param_array,
            sim.internal.current_time,
        )
    end

    return
end

function _update_initial_conditions!(stage::Stage, sim::Simulation)
    ini_cond_chronology = sim.sequence.ini_cond_chronology
    for (k, v) in iterate_initial_conditions(stage.internal.psi_container)
        initial_condition_update!(stage, k, v, ini_cond_chronology, sim)
    end
    return
end

function _update_parameters(stage::Stage, sim::Simulation)
    for container in iterate_parameter_containers(stage.internal.psi_container)
        update_parameter!(container.update_ref, container, stage, sim)
    end
    return
end

function _apply_warm_start!(stage::Stage)
    for variable in values(stage.internal.psi_container.variables)
        for e in variable
            current_solution = JuMP.value(e)
            JuMP.set_start_value(e, current_solution)
        end
    end
    return
end

""" Required update stage function call"""
function _update_stage!(stage::Stage, sim::Simulation)
    _update_parameters(stage, sim)
    _update_initial_conditions!(stage, sim)
    return
end

#############################Interfacing Functions##########################################
## These are the functions that the user will have to implement to update a custom stage ###
""" Generic Stage update function for most problems with no customization"""
function update_stage!(
    stage::Stage{M},
    sim::Simulation,
) where {M <: PowerSimulationsOperationsProblem}
    _update_stage!(stage, sim)
    return
end

"""
    execute!(sim::Simulation; kwargs...)

Solves the simulation model for sequential Simulations
and populates a nested folder structure created in Simulation()
with a dated folder of featherfiles that contain the results for
each stage and step.

# Arguments
- `sim::Simulation=sim`: simulation object created by Simulation()

# Example
```julia
sim = Simulation("Test", 7, stages, "/Users/folder")
execute!(sim::Simulation; kwargs...)
```
"""

function execute!(sim::Simulation; kwargs...)
    file_mode = "a"
    logger = configure_logging(sim.internal, file_mode)
    register_recorders!(sim.internal, file_mode)
    try
        return Logging.with_logger(logger) do
            _execute!(sim; kwargs...)
        end
    finally
        unregister_recorders!(sim.internal)
        close(logger)
    end
end

function _execute!(sim::Simulation; kwargs...)
    if sim.internal.reset
        sim.internal.reset = false
    elseif sim.internal.reset == false
        error("Re-build the simulation")
    end
    system_to_file = get(kwargs, :system_to_file, true)
    isnothing(sim.internal) &&
        error("Simulation not built, build the simulation to execute")
    TimerOutputs.reset_timer!(RUN_SIMULATION_TIMER)
    TimerOutputs.@timeit RUN_SIMULATION_TIMER "Execute Simulation" begin
        execution_order = get_execution_order(sim)
        for step in 1:get_steps(sim)
            TimerOutputs.@timeit RUN_SIMULATION_TIMER "Execution Step $(step)" begin
                println("Executing Step $(step)")
                IS.@record :simulation_status SimulationStepEvent(
                    sim.internal.current_time,
                    step,
                    "start",
                )
                for (ix, stage_number) in enumerate(execution_order)
                    IS.@record :simulation_status SimulationStageEvent(
                        sim.internal.current_time,
                        step,
                        stage_number,
                        "start",
                    )
                    TimerOutputs.@timeit RUN_SIMULATION_TIMER "Execution Stage $(stage_number)" begin
                        # TODO: implement some efficient way of indexing with stage name.
                        stage = get_stage(sim, stage_number)
                        stage_name = get_stage_name(sim, stage)
                        !stage_built(stage) &&
                            error("Stage $(stage_name) status is not BUILT")
                        stage_interval = get_stage_interval(sim, stage_name)
                        run_name = "step-$(step)-stage-$(stage_name)"
                        sim.internal.current_time = sim.internal.date_ref[stage_number]
                        sim.sequence.current_execution_index = ix
                        @info "Starting run $run_name $(sim.internal.current_time)"
                        raw_results_path = joinpath(
                            sim.internal.raw_dir,
                            run_name,
                            replace_chars("$(sim.internal.current_time)", ":", "-"),
                        )
                        mkpath(raw_results_path)
                        # Is first run of first stage? Yes -> don't update stage
                        TimerOutputs.@timeit RUN_SIMULATION_TIMER "Update Stage $(stage_number)" begin
                            !(step == 1 && ix == 1) && update_stage!(stage, sim)
                        end
                        TimerOutputs.@timeit RUN_SIMULATION_TIMER "Run Stage $(stage_number)" begin
                            run_stage(stage, sim.internal.current_time, raw_results_path)
                            sim.internal.run_count[step][stage_number] += 1
                        end
                        TimerOutputs.@timeit RUN_SIMULATION_TIMER "Update Cache $(stage_number)" begin
                            _update_caches!(sim, stage)
                        end
                        if warm_start_enabled(stage)
                            TimerOutputs.@timeit RUN_SIMULATION_TIMER "Warm Start $(stage_number)" begin
                                _apply_warm_start!(stage)
                            end
                        end
                    end
                    IS.@record :simulation_status SimulationStageEvent(
                        sim.internal.current_time,
                        step,
                        stage_number,
                        "done",
                    )
                    sim.internal.date_ref[stage_number] += stage_interval
                end
                IS.@record :simulation_status SimulationStepEvent(
                    sim.internal.current_time,
                    step,
                    "done",
                )
            end
        end
        sim_results = SimulationResultsReference(sim)
    end

    @info ("\n$(RUN_SIMULATION_TIMER)\n")
    serialize_sim_output(sim_results)
    return sim_results
end

struct SimulationSerializationWrapper
    steps::Int
    stages::Dict{String, StageSerializationWrapper}
    initial_time::Union{Nothing, Dates.DateTime}
    sequence::Union{Nothing, SimulationSequence}
    simulation_folder::String
    name::String
end

"""
    serialize(simulation::Simulation, path = ".")

Serialize the simulation to a directory in path.

Return the serialized simulation directory name that is created.

# Arguments
- `simulation::Simulation`: simulation to serialize
- `path = "."`: path in which to create the serialzed directory
- `force = false`: If true, delete the directory if it already exists. Otherwise, it will
   throw an exception.
"""
function serialize(simulation::Simulation; path = ".", force = false)
    directory = joinpath(path, "simulation-$(get_name(simulation))")
    stages = Dict{String, StageSerializationWrapper}()

    orig = pwd()
    if isdir(directory) || ispath(directory) && !force
        throw(ArgumentError("$directory already exists. Please delete it or pass force = true."))
    end
    rm(directory, recursive = true, force = true)
    mkdir(directory)
    cd(directory)

    try
        for (key, stage) in simulation.stages
            if isnothing(stage.internal)
                throw(ArgumentError("stage $(stage.internal.number) has not been built"))
            end
            sys_filename = "system-$(IS.get_uuid(stage.sys)).json"
            # Skip serialization if multiple stages have the same system.
            if !ispath(sys_filename)
                PSY.to_json(stage.sys, sys_filename)
            end
            stages[key] = StageSerializationWrapper(
                stage.template,
                sys_filename,
                stage.internal.psi_container.settings_copy,
                typeof(stage),
            )
        end
    finally
        cd(orig)
    end

    filename = joinpath(directory, SIMULATION_SERIALIZATION_FILENAME)
    obj = SimulationSerializationWrapper(
        simulation.steps,
        stages,
        simulation.initial_time,
        simulation.sequence,
        simulation.simulation_folder,
        get_name(simulation),
    )
    Serialization.serialize(filename, obj)
    @info "Serialized simulation" get_name(simulation) directory
    return directory
end

function deserialize(::Type{Simulation}, directory::AbstractString, stage_info::Dict)
    orig = pwd()
    cd(directory)

    try
        filename = SIMULATION_SERIALIZATION_FILENAME
        if !ispath(filename)
            throw(ArgumentError("$filename does not exist"))
        end

        obj = Serialization.deserialize(filename)
        if !(obj isa SimulationSerializationWrapper)
            throw(IS.DataFormatError("deserialized object has incorrect type $(typeof(obj))"))
        end

        stages = Dict{String, Stage{<:AbstractOperationsProblem}}()
        for (key, wrapper) in obj.stages
            sys_filename = wrapper.sys
            if !ispath(sys_filename)
                throw(ArgumentError("stage PowerSystems serialized file $sys_filename does not exist"))
            end
            sys = PSY.System(sys_filename)
            if !haskey(stage_info[key], "optimizer")
                throw(ArgumentError("stage_info must define 'optimizer'"))
            end
            stages[key] = wrapper.stage_type(
                wrapper.template,
                sys,
                restore_from_copy(
                    wrapper.settings;
                    optimizer = stage_info[key]["optimizer"],
                ),
                get(stage_info[key], "jump_model", nothing),
            )
        end

        sim = Simulation(;
            name = obj.name,
            steps = obj.steps,
            stages = stages,
            stages_sequence = obj.sequence,
            simulation_folder = obj.simulation_folder,
        )
        return sim
    finally
        cd(orig)
    end
end


text:
# syntax: proto2
using Compat
using ProtoBuf
import ProtoBuf.meta

# service methods for ContainerManagementProtocolService
const _ContainerManagementProtocolService_methods = MethodDescriptor[
        MethodDescriptor("startContainers", 1, hadoop.yarn.StartContainersRequestProto, hadoop.yarn.StartContainersResponseProto),
        MethodDescriptor("stopContainers", 2, hadoop.yarn.StopContainersRequestProto, hadoop.yarn.StopContainersResponseProto),
        MethodDescriptor("getContainerStatuses", 3, hadoop.yarn.GetContainerStatusesRequestProto, hadoop.yarn.GetContainerStatusesResponseProto)
    ] # const _ContainerManagementProtocolService_methods
const _ContainerManagementProtocolService_desc = ServiceDescriptor("hadoop.yarn.ContainerManagementProtocolService", 1, _ContainerManagementProtocolService_methods)

ContainerManagementProtocolService(impl::Module) = ProtoService(_ContainerManagementProtocolService_desc, impl)

mutable struct ContainerManagementProtocolServiceStub <: AbstractProtoServiceStub{false}
    impl::ProtoServiceStub
    ContainerManagementProtocolServiceStub(channel::ProtoRpcChannel) = new(ProtoServiceStub(_ContainerManagementProtocolService_desc, channel))
end # mutable struct ContainerManagementProtocolServiceStub

mutable struct ContainerManagementProtocolServiceBlockingStub <: AbstractProtoServiceStub{true}
    impl::ProtoServiceBlockingStub
    ContainerManagementProtocolServiceBlockingStub(channel::ProtoRpcChannel) = new(ProtoServiceBlockingStub(_ContainerManagementProtocolService_desc, channel))
end # mutable struct ContainerManagementProtocolServiceBlockingStub

startContainers(stub::ContainerManagementProtocolServiceStub, controller::ProtoRpcController, inp::hadoop.yarn.StartContainersRequestProto, done::Function) = call_method(stub.impl, _ContainerManagementProtocolService_methods[1], controller, inp, done)
startContainers(stub::ContainerManagementProtocolServiceBlockingStub, controller::ProtoRpcController, inp::hadoop.yarn.StartContainersRequestProto) = call_method(stub.impl, _ContainerManagementProtocolService_methods[1], controller, inp)

stopContainers(stub::ContainerManagementProtocolServiceStub, controller::ProtoRpcController, inp::hadoop.yarn.StopContainersRequestProto, done::Function) = call_method(stub.impl, _ContainerManagementProtocolService_methods[2], controller, inp, done)
stopContainers(stub::ContainerManagementProtocolServiceBlockingStub, controller::ProtoRpcController, inp::hadoop.yarn.StopContainersRequestProto) = call_method(stub.impl, _ContainerManagementProtocolService_methods[2], controller, inp)

getContainerStatuses(stub::ContainerManagementProtocolServiceStub, controller::ProtoRpcController, inp::hadoop.yarn.GetContainerStatusesRequestProto, done::Function) = call_method(stub.impl, _ContainerManagementProtocolService_methods[3], controller, inp, done)
getContainerStatuses(stub::ContainerManagementProtocolServiceBlockingStub, controller::ProtoRpcController, inp::hadoop.yarn.GetContainerStatusesRequestProto) = call_method(stub.impl, _ContainerManagementProtocolService_methods[3], controller, inp)

export ContainerManagementProtocolService, ContainerManagementProtocolServiceStub, ContainerManagementProtocolServiceBlockingStub, startContainers, stopContainers, getContainerStatuses

ext: jl


text:
@testset "9ketamineoriginalcropped.tif" begin
filename = joinpath(testbase, "9ketamineoriginalcropped.tif")
    # $ ScanImageTiffReader image shape 9ketamineoriginalcropped.tif
    # Shape: 119 x 114 x 300 @ u16
    tsize = ScanImageTiffReader.open(filename) do io
        size(io)
    end
    @test tsize == (119, 114, 300)
    ttype = ScanImageTiffReader.open(filename) do io
        pxtype(io)
    end
    @test ttype == UInt16

    # $ ScanImageTiffReader image bytes 9ketamineoriginalcropped.tif
    # 7.76253 MB
    dat = ScanImageTiffReader.open(filename) do io
        data(io)
    end
    @test abs(sizeof(dat)/2^20 - 7.76253) < 1e-5
    # @test ScanImageTiffReader.open(filename, length) == size(dat, 3) # TODO: check this assumption

    # $ ScanImageTiffReader descriptions --frame 0 9ketamineoriginalcropped.tif
    # "ImageJ=1.49v"
    # "images=300"
    # "slices=300"
    # "cf=0"
    # "c0=-32768.0"
    # "c1=1.0"
    # "vunit=Gray Value"
    # "loop=false"
    # "min=32556.0"
    # "max=34238.0"
    desc = ScanImageTiffReader.open(filename) do io
        split(description(io, 1), "\n")
    end
    @test desc[1] == "ImageJ=1.49v"
    @test desc[2] == "images=300"
    @test desc[3] == "slices=300"
    @test desc[4] == "cf=0"
    @test desc[5] == "c0=-32768.0"
    @test desc[6] == "c1=1.0"
    @test desc[7] == "vunit=Gray Value"
    @test desc[8] == "loop=false"
    @test desc[9] == "min=32556.0"
    @test desc[10] == "max=34238.0"

    # $ ScanImageTiffReader metadata 9ketamineoriginalcropped.tif
    #
    @test ScanImageTiffReader.open(metadata, filename) == ""
end


text:
###########################################################################
# Example: approximate posterior of power law described by two parameters #
###########################################################################

# Extra packages needed for this example
# PyPlot to plot the posterior approximation
# Random to fix the random seed
# Statistics for calculating the `uncertainty tube` when plotting results
# If missing, they must be installed explicitly!

using PyPlot, Random, Statistics

# ensure repeatability

seed = 15

rg = MersenneTwister(seed)

# ensure repeatability
seed = 101

rg = MersenneTwister(seed)

# Define model and parameters responsible for generating the data below

trueparameter = [0.5; 0.5]

decayfunction(x,p) = p[1] .+ exp.(-p[2]*x)

# Generate some data

N = 20

σ = 0.05

x = rand(rg, N)*10.0

y = decayfunction(x, trueparameter) .+ σ*randn(rg, N)

# define log-likelihood function

decaylogpdf(x, y, p) = -sum(abs2.(decayfunction(x, p).-y))/(2*σ*σ)

# Approximate posterior with Gaussian

posteriorfull, =  VI( p->decaylogpdf(x,y,p), randn(rg, 2), S = 100, iterations = 50, show_every=1)

posteriormvi,  = MVI( p->decaylogpdf(x,y,p), randn(rg, 2), S = 100, iterations = 50, show_every=1)

# Plot data

figure(1)

cla()

plot(x, y, "bo", label="observations")

xrange = collect(LinRange(minimum(x), maximum(x), 100))

plot(xrange, map(xi->decayfunction(xi, trueparameter), xrange), "k-", label="true power law")

legend()

# Plot true unnormalised posterior

figure(2)

cla()

prange = 0.3:0.001:0.9

contourf( repeat(prange,1,length(prange)),  repeat(prange',length(prange),1), map(p -> exp(decaylogpdf(x,y,p)), [[x;y] for x in prange, y in prange]), cmap=plt.cm.binary)

xlabel("p1")

ylabel("p2")

# Plot our Gaussian approximation

ApproximateVI.plot_ellipse(posteriorfull, "b", "full")

ApproximateVI.plot_ellipse(posteriormvi,  "g", "mvi")

# Add the true parameter on the plot

plot(trueparameter[1], trueparameter[2], "mo", label="true parameter")

legend()

# Add predictions from approximations

figure(1)

# plot full approximation with mean and standard deviation tube

samples = [decayfunction(xrange, rand(posteriorfull)) for _ in 1:1000]

fill_between(xrange, mean(samples) - std(samples), mean(samples) + std(samples), "b", alpha = 0.1)

plot(xrange, mean(samples), "b--", label = "mean prediction full")

# plot mvi approximation with mean and standard deviation tube

samples = [decayfunction(xrange, rand(posteriormvi)) for _ in 1:1000]

fill_between(xrange, mean(samples) - std(samples), mean(samples) + std(samples), "g", alpha = 0.1)

plot(xrange, mean(samples), "g--", label = "mean prediction mvi")

legend()

nothing


text:
"""
    nsoli(F!, x0, FS, FPS, Jvec=dirder; rtol=1.e-6, atol=1.e-12,
               maxit=20, lmaxit=-1, lsolver="gmres", eta=.1,
               fixedeta=true, Pvec=nothing, pside="right",
               armmax=10, dx = 1.e-7, armfix=false, pdata = nothing,
               printerr = true, keepsolhist = false, stagnationok=false)
)

C. T. Kelley, 2021

Julia versions of the nonlinear solvers from my SIAM books.
Herewith: nsoli

You must allocate storage for the function and the Krylov basis in advance
--> in the calling program <-- ie. in FS and FPS

Inputs:\n
- F!: function evaluation, the ! indicates that F! overwrites FS, your
    preallocated storage for the function.\n
    So FS=F!(FS,x) or FS=F!(FS,x,pdata) returns FS=F(x)

- x0: initial iterate\n

- FS: Preallocated storage for function. It is an N x 1 column vector\n
  You may store it as (n,) or (n,1), depending on what F! likes to see.

- FPS: preallocated storage for the Krylov basis. It is an N x m matrix where
       you plan to take at most m-1 GMRES iterations before a restart. \n

- Jvec: Jacobian vector product, If you leave this out the
    default is a finite difference directional derivative.\n
    So, FP=Jvec(v,FS,x) or FP=Jvec(v,FS,x,pdata) returns FP=F'(x) v. \n
    (v, FS, x) or (v, FS, x, pdata) must be the argument list,
    even if FP does not need FS.
    One reason for this is that the finite-difference derivative
    does and that is the default in the solver.

- Precision: Lemme tell ya 'bout precision. I designed this code for
    full precision functions and linear algebra in any precision you want.
    You can declare FPS as Float64 or Float32 and nsoli
    will do the right thing. Float16 support is there, but not working well.

    If the Jacobian is reasonably well conditioned, you can cut the cost
    of orthogonalization and storage (for GMRES) in half with no loss.
    There is no benefit if your linear solver is not GMRES or if
    othogonalization and storage of the Krylov vectors is only a
    small part of the cost of the computation. So if your preconditioner
    is good and you only need a few Krylovs/Newton, reduced precision won't
    help you much.

    BiCGSTAB does not benefit from reduced precsion.

----------------------

Keyword Arguments (kwargs):\n

rtol and atol: relative and absolute error tolerances\n

maxit: limit on nonlinear iterations\n

lmaxit: limit on linear iterations. If lmaxit > m-1, where FPS has
m columns, and you need more than m-1 linear iterations, then GMRES
will restart.

The default is -1 for GMRES. This means that you'll take m-1 iterations,
where size(V) = (n,m), and get no restarts. For BiCGSTAB the default is 10.

lsolver: the linear solver, default = "gmres"\n
Your choices will be "gmres" or "bicgstab". However,
gmres is the only option for now.

eta and fixed eta: eta > 0 or there's an error

The linear solver terminates when ||F'(x)s + F(x) || <= etag || F(x) ||

where

etag = eta if fixedeta=true

etag = Eisenstat-Walker as implemented in book if fixedeta=false

The default, which may change, is eta=.1, fixedeta=true

Pvec: Preconditioner-vector product. The rules are similar to Jvec
    So, Pv=Pvec(v,x) or Pv=Pvec(v,x,pdata) returns P(x) v where
    P(x) is the preconditioner. You must use x as an input even
    if your preconditioner does not depend on x

pside: apply preconditioner on pside, default = "right". I do not
      recommend "left". See Chapter 3 for the story on this.

armmax: upper bound on step size reductions in line search\n

dx: default = 1.e-7\n
difference increment in finite-difference derivatives
      h=dx*norm(x,Inf)+1.e-8

armfix: default = false\n
The default is a parabolic line search (ie false). Set to true and
the step size will be fixed at .5. Don't do this unless you are doing
experiments for research.\n

pdata:\n
precomputed data for the function/Jacobian-vector/Preconditioner-vector
products.  Things will go better if you use this rather than hide the data
in global variables within the module for your function/Jacobian

If you use pdata in any of F!, Jvec, or Pvec, you must use in in all of them.

printerr: default = true\n
I print a helpful message when the solver fails. To suppress that
message set printerr to false.

keepsolhist: default = false\n
Set this to true to get the history of the iteration in the output
tuple. This is on by default for scalar equations and off for systems.
Only turn it on if you have use for the data, which can get REALLY LARGE.

stagnationok: default = false\n
Set this to true if you want to disable the line search and either
observe divergence or stagnation. This is only useful for research
or writing a book.

Output:\n
- A named tuple (solution, functionval, history, stats, idid,
               errcode, solhist)
where

   -- solution = converged result

   -- functionval = F(solution)

   -- history = the vector of residual norms (||F(x)||) for the iteration

   -- stats = named tuple of the history of (ifun, ijac, iarm, ikfail), the
number of functions/Jacobian-vector prods/steplength reductions/linear solver
failures at each iteration. Linear solver failures DO NOT mean that the
nonlinear solver will fail. You should look at this stat if, for example,
the line search fails. Increasing the size of FPS and/or lmaxit might
solve the problem.

I do not count the function values for a finite-difference derivative
because they count toward a Jacobian-vector product.

  -- idid=true if the iteration succeeded and false if not.

  -- errcode = 0 if if the iteration succeeded

        = -1 if the initial iterate satisfies the termination criteria

        = 10 if no convergence after maxit iterations

        = 1  if the line search failed

   -- solhist:\n
      This is the entire history of the iteration if you've set
      keepsolhist=true\n

solhist is an N x K array where N is the length of x and K is the number of
iteration + 1. So, for scalar equations, it's a row vector.

------------------------

### Example for nsoli

#### Simple 2D problem.
You should get the same results as for nsol.jl because
GMRES will solve the equation for the step exactly in two iterations. Finite
difference Jacobians and analytic Jacobian-vector products for full precision
and finite difference Jacobian-vector products for single precision.

BiCGSTAB converges in 5 itertions and each nonlinear iteration costs
two Jacobian-vector products. Note that the storage for the Krylov
space in GMRES (jvs) is replace by a single vector (fpv) when BiCGSTAB
is the linear solver.

```jldoctest
julia> function f!(fv,x)
       fv[1]=x[1] + sin(x[2])
       fv[2]=cos(x[1]+x[2])
       end
f! (generic function with 1 method)

julia> function JVec(v, fv, x)
       jvec=zeros(2,);
       p=-sin(x[1]+x[2])
       jvec[1]=v[1]+cos(x[2])*v[2]
       jvec[2]=p*(v[1]+v[2])
       return jvec
       end
JVec (generic function with 1 method)

julia> x0=ones(2,); fv=zeros(2,); jv=zeros(2,2); jv32=zeros(Float32,2,2);

julia> jvs=zeros(2,3); jvs32=zeros(Float32,2,3);

julia> nout=nsol(f!,x0,fv,jv; sham=1);

julia> kout=nsoli(f!,x0,fv,jvs,JVec; fixedeta=true, eta=.1, lmaxit=2);

julia> kout32=nsoli(f!,x0,fv,jvs32; fixedeta=true, eta=.1, lmaxit=2);

julia> [nout.history kout.history kout32.history]
5×3 Array{Float64,2}:
 1.88791e+00  1.88791e+00  1.88791e+00
 2.43119e-01  2.43120e-01  2.43119e-01
 1.19231e-02  1.19231e-02  1.19231e-02
 1.03266e-05  1.03261e-05  1.03273e-05
 1.46416e-11  1.40862e-11  1.45457e-11

julia> fpv=zeros(2,);

julia> koutb=nsoli(f!,x0,fv,fpv,JVec; fixedeta=true, eta=.1, lmaxit=2,
       lsolver="bicgstab");

julia> koutb.history
6-element Vector{Float64}:
 1.88791e+00
 2.43120e-01
 1.19231e-02
 4.87500e-04
 7.54236e-06
 3.84646e-07
```

"""
function nsoli(
    F!,
    x0,
    FS,
    FPS,
    Jvec = dirder;
    rtol = 1.e-6,
    atol = 1.e-12,
    maxit = 20,
    lmaxit = -1,
    lsolver = "gmres",
    eta = 0.1,
    fixedeta = true,
    Pvec = nothing,
    pside = "right",
    armmax = 10,
    dx = 1.e-7,
    armfix = false,
    pdata = nothing,
    printerr = true,
    keepsolhist = false,
    stagnationok = false,
)
    itc = 0
    idid = true
    iline = false
    #
    #   If I'm letting the iteration stagnate and turning off the
    #   linesearch, then the line search cannot fail.
    #
    stagflag = stagnationok && (armmax == 0)
    #=
    Named tuple with the iteration data. This makes communiction
    with the linear solvers and the line search easier.
    =#
    (ItRules, x, n, solhist) = Newton_Krylov_Init(
        x0,
        dx,
        F!,
        Jvec,
        Pvec,
        pside,
        lsolver,
        eta,
        fixedeta,
        armmax,
        armfix,
        maxit,
        lmaxit,
        printerr,
        pdata,
        keepsolhist,
    )
    #    keepsolhist ? (solhist = solhistinit(n, maxit, x)) : (solhist = [])
    #
    # First Evaluation of the function. Initialize the iteration stats.
    # Fix the tolerances for convergence and define the derivative FPF
    # outside of the main loop for scoping.
    #
    FS = EvalF!(F!, FS, x, pdata)
    resnorm = norm(FS)
    tol = rtol * resnorm + atol
    FPF = []
    ItData = ItStatsK(resnorm)
    newiarm = -1
    newfun = 0
    newjac = 0
    newikfail = 0
    ke_report = false
    residratio = 1.0
    armstop = true
    etag = eta
    #
    # Preallocate a few vectors for the step, trial step, trial function
    #
    step = copy(x)
    xt = copy(x)
    FT = copy(x)
    #
    # If the initial iterate satisfies the termination criteria, tell me.
    #
    toosoon = (resnorm <= tol)
    #
    # The main loop stops on convergence, too many iterations, or a
    # line search failure after a derivative evaluation.
    #
    while resnorm > tol && itc < maxit && (armstop || stagnationok)
        #
        newfun = 0
        newjac = 0
        newikfail = 0
        #
        #
        # The GMRES solver will do the orthogonalization in lower
        # precision. I've tested Float32, but see the docstrings
        # for all the caveats. This is not the slam dunk it was
        # for Gaussian elimination on dense matrices.
        #
        step .*= 0.0
        etag = forcing(itc, residratio, etag, ItRules, tol, resnorm)
        kout = Krylov_Step!(step, x, FS, FPS, ItRules, etag)
        step .= kout.step
        #
        # For GMRES you get 1 jac-vec per iteration and there is no jac-vec
        # for the initial inner iterate of zero. For BiCGSTAB it's two
        # jac-vecs per iteration.
        #
        newjac = kout.Lstats.lits
        (lsolver == "gmres") || (newjac *= 2)
        linok = kout.Lstats.idid
        linok || (ke_report = Krylov_Error(lmaxit, ke_report); newikfail = 1)
        #
        # Compute the trial point, evaluate F and the residual norm.
        # The derivative is never old for Newton-Krylov
        #
        AOUT = armijosc(xt, x, FT, FS, step, resnorm, ItRules, false)
        #
        # update solution/function value
        #
        x .= AOUT.ax
        FS .= AOUT.afc
        #
        # If the line search fails
        # stop the iteration. Print an error message unless
        # stagnationok == true
        #
        armstop = AOUT.idid
        iline = ~armstop && ~stagflag
        #
        # Keep the books.
        #
        residm = resnorm
        resnorm = AOUT.resnorm
        residratio = resnorm / residm
        updateStats!(ItData, newfun, newjac, AOUT, newikfail)
        newiarm = AOUT.aiarm
        itc += 1
        keepsolhist && (@views solhist[:, itc+1] .= x)
        #        ~keepsolhist || (@views solhist[:, itc+1] .= x)
    end
    #    solution = x
    #    functionval = FS
    (idid, errcode) = NewtonOK(resnorm, iline, tol, toosoon, itc, ItRules)
    newtonout = CloseIteration(x, FS, ItData, idid, errcode, keepsolhist, solhist)
    return newtonout
end


text:
module RocketFilterOperatorTest

using Test
using Rocket

include("../test_helpers.jl")

@testset "operator: filter()" begin

    println("Testing: operator filter()")

    run_proxyshowcheck("Filter", filter(_ -> true))

    run_testset([
        (
            source      = from(1:5) |> filter(d -> d % 2 == 0),
            values      = @ts([ 2, 4, c ]),
            source_type = Int
        ),
        (
            source      = from(1:5) |> filter(d -> d % 2 == 1),
            values      = @ts([ 1, 3, 5, c ]),
            source_type = Int
        ),
        (
            source = completed() |> filter(d -> d % 2 == 1),
            values = @ts(c)
        ),
        (
            source = faulted("e") |> filter(d -> d % 2 == 1),
            values = @ts(e("e"))
        ),
        (
            source = never() |> filter(d -> d % 2 == 1),
            values = @ts()
        ),
        (
            source = of(1) |> filter(d -> false),
            values = @ts(c)
        ),
        (
            source = from(1:5) |> filter(d -> true),
            values = @ts([ 1:5, c ])
        )
    ])

end

end


text:
#=
  should probably just use an integer as hash: p*root of poly

  so a is an integral element where the norm is almost smooth, it means
  norm(a) = prod over factor base * p
  where p is a prime.
  This means
  <p, a> is a prime ideal of norm p hence of degree 1
  if p is no index divisor, then <p,a> = <p, b> where
  b is a linear factor of f, the defining polynomial, mod p.
  I can compute b as gcd(a, f) of course.
  =#
function special_prime_ideal(p::fmpz, a::nf_elem)
  K = parent(a)
  f = K.pol
  R = parent(f)
  Zx = PolynomialRing(FlintZZ)[1]
  Zpx = PolynomialRing(GF(UInt(p), cached=false), "\$x_p", cached=false)[1]
  g = Zpx(a)
  ff = Zpx(f)
  gcd!(g, g, ff)
  return lift(Zx, g)
end

function class_group_add_relation(clg::ClassGrpCtx{T}, a::nf_elem; orbit::Bool = true, integral::Bool = false) where T
  return class_group_add_relation(clg, a, norm(a), fmpz(1), orbit = orbit, integral = integral)
end

#deal with integral and non-integral elements differently. Computing the order
#denominator is expensive (and mostly unnecessary)
function class_group_add_relation(clg::ClassGrpCtx{T}, a::nf_elem, n::fmpq, nI::fmpz; orbit::Bool = true, integral::Bool = true) where T
  if iszero(a)
    return false
  end
  if hash(a) in clg.RS
    return false
  end

  nb = div(nbits(numerator(n)), 2)
  if haskey(clg.normStat, nb)
    clg.normStat[nb] += 1
  else
    clg.normStat[nb] = 1
  end

  O = order(clg.FB.ideals[1])
  easy = isdefining_polynomial_nice(parent(a))
  @vprint :ClassGroup 3 "trying relation of length $(Float64(length(a))) and norm $(Float64(n*nI)), effective $(Float64(n))\n"
  if integral #element is known to be integral
    fl, r = issmooth!(clg.FB.fb_int, numerator(n*nI))
  else
    fl, r = issmooth!(clg.FB.fb_int, numerator(n*nI)*denominator(a, O))
  end
  if !fl
    @vprint :ClassGroup 3 "not int-smooth\n"
#    println("not int-smooth");
    # try for large prime?
    if easy && abs(r) < clg.B2 && isprime(r) && !isindex_divisor(O, r)
      @vprint :ClassGroup 3 "gives potential large prime\n"
      i = special_prime_ideal(r, a)
      #TODO: check Galois orbit of special ideal
      if haskey(clg.largePrime, (r, i))
        lp = clg.largePrime[(r, i)]
        fl, r1 = _factor!(clg.FB, a, false, n*nI)
        fl, r2 = _factor!(clg.FB, lp[1], false, norm(lp[1]))
        b = FacElem(Dict([(a,1), (lp[1],-1)]))
        fl = class_group_add_relation(clg, b, r1 - r2)
        if fl
          clg.largePrime_success += 1
        else
          clg.largePrime_no_success += 1
        end
      else
        clg.largePrime[i] = (a, n*nI)
      end
      clg.largePrimeCnt += 1
    else
      clg.bad_rel += 1
    end
    #println(" -> fail")
    return false
  end
  fl, res = _factor!(clg.FB, a, false, n*nI, integral)
  if fl
    if res in clg.M.rel_gens || res in clg.M.bas_gens
      return false
    end
    @vprint :ClassGroup 3 "adding $res\n"
    if add_gen!(clg.M, res)
      push!(clg.R_gen, a)
    else
      push!(clg.R_rel, a)
    end
    push!(clg.RS, hash(a))
    if orbit && isdefined(clg, :aut_grp)
      n = res
      o = clg.aut_grp
      function op_smat(n::SRow, p::Nemo.Generic.perm)
        r = [(p[i], v) for (i,v) = n]
        sort!(r, lt = (a,b)->a[1]<b[1])
        return typeof(n)(r)
      end

      @v_do :ClassGroup 1 println(" adding orbit with $(length(o)) elements")
      for (b, m) in o
        nn = op_smat(n, m)
        if nn != n
          ba = b(a)
          if nn in clg.M.bas_gens || nn in clg.M.rel_gens
            break
          end
          if add_gen!(clg.M, nn, false)
            push!(clg.R_gen, ba)
            clg.rel_cnt += 1
            push!(clg.RS, hash(ba))
          else
            #push!(clg.R_rel, ba)
          end
        end
      end
    end

    clg.rel_cnt += 1
#    @assert clg.rel_cnt < 2*ncols(clg.M)
    @v_do :ClassGroup 1 println(" -> OK, rate currently ",
           clg.bad_rel/clg.rel_cnt, " this ", clg.bad_rel - clg.last,
           " rank now ", rank(clg.M), " of ", length(clg.FB.ideals))
    clg.last = clg.bad_rel
    return true
  else
    @vprint :ClassGroup 3 "not alg-smooth\n"
    clg.bad_rel += 1
    return false
  end
end

function class_group_add_relation(clg::ClassGrpCtx{SMat{fmpz}}, a::FacElem{nf_elem, AnticNumberField})
  R = SRow{fmpz}()
  for i = 1:length(clg.FB.ideals)
    p = clg.FB.ideals[i]
    v = valuation(a, p)
    if !iszero(v)
      push!(R.values, v)
      push!(R.pos, i)
    end
  end
  return class_group_add_relation(clg, a, R)
end

function class_group_add_relation(clg::ClassGrpCtx{SMat{fmpz}}, a::FacElem{nf_elem, AnticNumberField}, R::SRow{fmpz})

  if hash(a) in clg.RS
    return false
  end

  O = order(clg.FB.ideals[1])

  @vprint :ClassGroup 3 "adding $R\n"

  if add_gen!(clg.M, R)
    push!(clg.R_gen, a)
  else
    push!(clg.R_rel, a)
  end
  push!(clg.RS, hash(a))

  clg.rel_cnt += 1
#    @assert clg.rel_cnt < 2*ncols(clg.M)
  @v_do :ClassGroup 1 println(" -> OK, rate currently ",
         clg.bad_rel/clg.rel_cnt, " this ", clg.bad_rel - clg.last)
  clg.last = clg.bad_rel
  return true
end


text:
# BBI Total/Section Summary
# =========================

# Supplemental Table 7.
struct Summary
    # base coverage
    cov::UInt64
    # minimum value
    min::Float64
    # maximum value
    max::Float64
    # sum of values
    sum::Float64
    # sum of squared values
    ssq::Float64
end

const SUMMARY_SIZE = 40

function Base.read(io::IO, ::Type{Summary})
    return Summary(
        read(io, UInt64),
        read(io, Float64), read(io, Float64),
        read(io, Float64), read(io, Float64))
end

function Base.write(stream::IO, summary::Summary)
    return write(
        stream,
        summary.cov,
        summary.min,
        summary.max,
        summary.sum,
        summary.ssq)
end

# Summary per section.
# These are used to create the header, total summary, and data index.
struct SectionSummary
    # data range
    chromid::UInt32
    chromstart::UInt32
    chromend::UInt32

    # file offset to data and its size
    offset::UInt64
    datasize::UInt64
end


text:
#find out the test statistic when m0=1
import LatentGaussianMixtureModel
import Distributions, StatsBase
import Yeppp

@everywhere using LatentGaussianMixtureModel
@everywhere using Distributions, Yeppp, StatsBase, Random

#Brun calculate the statistic for one data set;
#b is the the random number seed, from 1 to 100
@everywhere function Brun(b::Integer, X::Matrix{Float64},
    groupindex::Vector; debuginfo::Bool=false)

    mu_true = [-0.979]
    wi_true =  [1.0]
    sigmas_true = [0.26]
    betas_true=[0.019498134778826972,0.007110782868926634,0.030919753747641616,-0.26957245719168493,-0.5257588591653842,-0.6314135703754931,-0.799509687993618,0.07786105411635978,0.12059632526169425,0.22504173963135948]

    #Randomly data generation based on the setting on datagen.jl
    Random.seed!(b * 100)
    n = Int(maximum(groupindex))
    m = MixtureModel(map((u, v) -> Normal(u, v), mu_true, sigmas_true), wi_true)
    gamma_true = rand(m, n)

    prob = exp(gamma_true[groupindex] .+ X*betas_true)
    prob= prob ./ (1 .+ prob)
    Y = Bool[rand(Binomial(1, prob[i])) == 1 for i in 1:length(groupindex)];

    m = LGMModel(X, Y, groupindex, 2)
    LatentGaussianMixtureModel.fit!(m)
    LatentGaussianMixtureModel.ranef!(m)
    println("Mission $b finished.")
    return [sort(m.gammaprediction)[1:10], m.p, m.μ, m.σ;]
end


text:
"""
    dm_vals(e::Array{<:Real,2}, closest::Array{<:Real,2}, max_f_vals::Array{<:Real,2}, min_f_vals::Array{<:Real,2})

Compute dm values for each feature (see reference, auxiliary function).
"""
function dm_vals(e::Array{<:Real,2}, closest::Array{<:Real,2}, max_f_vals::Array{<:Real,2}, min_f_vals::Array{<:Real,2})

    # Allocate matrix for results.
	results =  Array{Float64, 2}(undef, size(closest))

    # Compute diff values and dm values.
	diff_vals = abs.(e .- closest)./(max_f_vals .- min_f_vals .+ eps(Float64))
	@inbounds for i = 1:length(e)
		results[:,i] = Statistics.mean(diff_vals[:, 1:end .!= i], dims=2)
	end

    # Return dm value for each feature.
	return results, diff_vals
end

"""
    reliefmss(data::Array{<:Real,2}, target::Array{<:Integer,1}, m::Signed=-1,
                   k::Integer=10, dist_func::Any=(e1, e2) -> sum(abs.(e1 .- e2), dims=2);
                   f_type::String="continuous")

Compute feature weights using ReliefMSS algorithm. The f_type argument specifies whether the features are continuous or discrete
and can either have the value of "continuous" or "discrete".

---
# Reference:
- Salim Chikhi and Sadek Benhammada. ReliefMSS: A variation on a
feature ranking ReliefF algorithm. IJBIDM, 4:375–390, 2009.
"""
function reliefmss(data::Array{<:Real,2}, target::Array{<:Integer,1}, m::Signed=-1,
                   k::Integer=10, dist_func::Any=(e1, e2) -> sum(abs.(e1 .- e2), dims=2);
                   f_type::String="continuous")

    # Initialize feature weights vector.
    weights = zeros(Float64, 1, size(data, 2))

    # Compute vectors of maximum and minimum feature values.
    max_f_vals = maximum(data, dims=1)
    min_f_vals = minimum(data, dims=1)

    # Sample m examples without replacement.
    sample_idxs = StatsBase.sample(1:size(data, 1), if (m==-1) size(data,1) else m end, replace=false)
    if (m == -1) m = size(data, 1) end

    # Compute probabilities of class values in training set.
    classes_map = countmap(target)
    num_samples = length(target)
    p_classes = Array{Float64}(undef, length(keys(classes_map)), 2)
    @inbounds for (idx, p) = enumerate(classes_map)  # Compute matrix representation of probabilities.
        p_classes[idx, :] = [p[1], p[2]/num_samples]
    end

    # Compute pairwise distances between samples (vector form).
    dists = Array{Float64}(undef, Int64((size(data, 1)^2 - size(data, 1))/2))

    # Construct pairwise distances vector using vectorized distance function.
    top_ptr = 1
    @inbounds for idx = 1:size(data,1)-1
        upper_lim = top_ptr + size(data, 1) - idx - 1
        dists[top_ptr:upper_lim] = dist_func(data[idx:idx, :], data[idx+1:end, :])
        top_ptr = upper_lim + 1
    end

    # Go over sampled indices.
    @inbounds for idx = sample_idxs

        # Row and column indices for querying pairwise distance vector.
        row_idxs = repeat([idx - 1], size(data, 1))
        col_idxs = collect(0:size(data, 1)-1)

        # Get indices in distance vector (from square form indices).
        neigh_idx = square_to_vec(row_idxs[(target .== target[idx]) .& (1:length(target) .!= idx)],
                                  col_idxs[(target .== target[idx]) .& (1:length(target) .!= idx)], size(data, 1)) .+ 1
        idx_k_nearest_same = partialsortperm(dists[neigh_idx], 1:k)

        # Get k nearest hits.
        k_nearest_same = data[(target .== target[idx]) .& (1:length(target) .!= idx), :][idx_k_nearest_same, :]

        # Allocate matrix for storing the k nearest misses.
        k_nearest_other = Array{Float64}(undef, k * (length(keys(classes_map)) - 1), size(data, 2))

        # Go over class values not equal to class value of currently sampled sample.
        top_ptr = 1
        @inbounds for cl = keys(classes_map)
            if cl != target[idx]
                # If class not equal to sampled example, find indices in distance vector of examples with this class.
                neigh_idx_nxt = Int64.(square_to_vec(row_idxs[target .== cl], col_idxs[target .== cl], size(data, 1))) .+ 1
                idx_k_nearest_other_nxt = partialsortperm(dists[neigh_idx_nxt], 1:k)

                # Get k closest examples from this class and store in matrix.
                k_nearest_other_nxt = data[target .== cl, :][idx_k_nearest_other_nxt, :]
                k_nearest_other[top_ptr:top_ptr+k-1, :] = k_nearest_other_nxt
                top_ptr += k
            end
        end

        ### MARKING CONSIDERED FEATURES ###

        # Compute DM values and DIFF values for each feature of each nearest hit and nearest miss.
        dm_vals_same, diff_vals_same = dm_vals(data[idx:idx, :], k_nearest_same, max_f_vals, min_f_vals)
        dm_vals_other, diff_vals_other = dm_vals(data[idx:idx, :], k_nearest_other, max_f_vals, min_f_vals)

        # Compute masks for considered features of nearest hits and nearest misses.
        features_msk_same = diff_vals_same .> dm_vals_same
        features_msk_other = diff_vals_other .> dm_vals_other

        ###################################

        # Get probabilities of classes not equal to class of sampled example.
        p_classes_other = p_classes[p_classes[:, 1] .!= target[idx], 2]

        # Compute diff sum weights for closest examples from different classes.
        p_weights = p_classes_other./(1 .- p_classes[p_classes[:, 1] .== target[idx], 2])
        weights_mult = reshape(repeat(p_weights, inner=k), :, 1)

        ### Weights Update ###
        #
        if f_type == "continuous"
            # If features continuous.

            # Penalty term
            penalty = (abs.(data[idx:idx, :] .- k_nearest_same)./(max_f_vals .- min_f_vals .+ eps(Float64))) .- dm_vals_same
            penalty[.!features_msk_same] .= 0.0
            penalty = sum(penalty, dims=1)

            # Reward term
            reward = weights_mult .* (abs.(data[idx:idx, :] .- k_nearest_other)./(max_f_vals .- min_f_vals .+ eps(Float64))) - dm_vals_other
            reward[.!features_msk_other] .= 0.0
            reward = sum(reward, dims=1)

            # Weights update
            weights = weights .- penalty./(m*k) .+ reward./(m*k)

        elseif f_type == "discrete"
            # If features discrete.

            # Penalty term
            penalty = Int64.(data[idx:idx, :] .!= k_nearest_same) .- dm_vals_same
            penalty[.!features_msk_same] .= 0.0
            penalty = sum(penalty, dims=1)

            # Reward term
            reward = weights_mult .* Int64.(data[idx:idx, :] .!= k_nearest_other) - dm_vals_other
            reward[.!features_msk_other] .= 0.0
            reward = sum(reward, dims=1)

            # Weights update
            weights = weights .- penalty./(m*k) .+ reward./(m*k)

        else
            throw(DomainError(f_type, "f_type can only be equal to \"continuous\" or \"discrete\"."))
        end

        ######################
    end

    # Return computed feature weights.
    return vec(weights)
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct VirtualMachineScaleSetProperties <: SwaggerModel
    upgradePolicy::Any # spec type: Union{ Nothing, UpgradePolicy } # spec name: upgradePolicy
    virtualMachineProfile::Any # spec type: Union{ Nothing, VirtualMachineScaleSetVMProfile } # spec name: virtualMachineProfile
    provisioningState::Any # spec type: Union{ Nothing, String } # spec name: provisioningState
    overprovision::Any # spec type: Union{ Nothing, Bool } # spec name: overprovision
    uniqueId::Any # spec type: Union{ Nothing, String } # spec name: uniqueId
    singlePlacementGroup::Any # spec type: Union{ Nothing, Bool } # spec name: singlePlacementGroup

    function VirtualMachineScaleSetProperties(;upgradePolicy=nothing, virtualMachineProfile=nothing, provisioningState=nothing, overprovision=nothing, uniqueId=nothing, singlePlacementGroup=nothing)
        o = new()
        validate_property(VirtualMachineScaleSetProperties, Symbol("upgradePolicy"), upgradePolicy)
        setfield!(o, Symbol("upgradePolicy"), upgradePolicy)
        validate_property(VirtualMachineScaleSetProperties, Symbol("virtualMachineProfile"), virtualMachineProfile)
        setfield!(o, Symbol("virtualMachineProfile"), virtualMachineProfile)
        validate_property(VirtualMachineScaleSetProperties, Symbol("provisioningState"), provisioningState)
        setfield!(o, Symbol("provisioningState"), provisioningState)
        validate_property(VirtualMachineScaleSetProperties, Symbol("overprovision"), overprovision)
        setfield!(o, Symbol("overprovision"), overprovision)
        validate_property(VirtualMachineScaleSetProperties, Symbol("uniqueId"), uniqueId)
        setfield!(o, Symbol("uniqueId"), uniqueId)
        validate_property(VirtualMachineScaleSetProperties, Symbol("singlePlacementGroup"), singlePlacementGroup)
        setfield!(o, Symbol("singlePlacementGroup"), singlePlacementGroup)
        o
    end
end # type VirtualMachineScaleSetProperties

const _property_map_VirtualMachineScaleSetProperties = Dict{Symbol,Symbol}(Symbol("upgradePolicy")=>Symbol("upgradePolicy"), Symbol("virtualMachineProfile")=>Symbol("virtualMachineProfile"), Symbol("provisioningState")=>Symbol("provisioningState"), Symbol("overprovision")=>Symbol("overprovision"), Symbol("uniqueId")=>Symbol("uniqueId"), Symbol("singlePlacementGroup")=>Symbol("singlePlacementGroup"))
const _property_types_VirtualMachineScaleSetProperties = Dict{Symbol,String}(Symbol("upgradePolicy")=>"UpgradePolicy", Symbol("virtualMachineProfile")=>"VirtualMachineScaleSetVMProfile", Symbol("provisioningState")=>"String", Symbol("overprovision")=>"Bool", Symbol("uniqueId")=>"String", Symbol("singlePlacementGroup")=>"Bool")
Base.propertynames(::Type{ VirtualMachineScaleSetProperties }) = collect(keys(_property_map_VirtualMachineScaleSetProperties))
Swagger.property_type(::Type{ VirtualMachineScaleSetProperties }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_VirtualMachineScaleSetProperties[name]))}
Swagger.field_name(::Type{ VirtualMachineScaleSetProperties }, property_name::Symbol) =  _property_map_VirtualMachineScaleSetProperties[property_name]

function check_required(o::VirtualMachineScaleSetProperties)
    true
end

function validate_property(::Type{ VirtualMachineScaleSetProperties }, name::Symbol, val)
end


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

export
    av_dv_frame_profile,
    av_dv_codec_profile,
    av_dv_codec_profile2

function av_dv_frame_profile(sys, frame, buf_size::Integer)
    ccall((:av_dv_frame_profile, libavcodec), Ptr{AVDVProfile}, (Ptr{AVDVProfile}, Ptr{UInt8}, UInt32), sys, frame, buf_size)
end

function av_dv_codec_profile(width::Integer, height::Integer, pix_fmt::AVPixelFormat)
    ccall((:av_dv_codec_profile, libavcodec), Ptr{AVDVProfile}, (Cint, Cint, AVPixelFormat), width, height, pix_fmt)
end

function av_dv_codec_profile2(width::Integer, height::Integer, pix_fmt::AVPixelFormat, frame_rate::AVRational)
    ccall((:av_dv_codec_profile2, libavcodec), Ptr{AVDVProfile}, (Cint, Cint, AVPixelFormat, AVRational), width, height, pix_fmt, frame_rate)
end


text:
for file in ["comm.jl", "msg.jl", "execute_request.jl", "stdio.jl"]
    println(file)
    include(file)
end


text:
module Testeachrow

using Test
using DataFrames
using DataFramesMeta
using Statistics

const ≅ = isequal

df = DataFrame(a = [1, 2, missing], b = [4, 5, 6])

@testset "passmissing by row" begin

    no_missing(x::Int, y::Int) = x + y

    df = DataFrame(a = [1, 2, missing], b = [4, 5, 6])

    d = @transform df @byrow @passmissing :c = no_missing(:a, :b)
    @test d.c ≅ [5, 7, missing]

    d = @transform df @passmissing @byrow :c = no_missing(:a, :b)
    @test d.c ≅ [5, 7, missing]

    d = @transform df @byrow begin
        :c = :a + :b
        @passmissing :d = no_missing(:a, :b)
    end
    @test d.c ≅ d.d

    d = @rselect df @passmissing :c = no_missing(:a, :b)
    @test d.c ≅ [5, 7, missing]

    d = @rselect df @passmissing :c = no_missing(:a, :b)
    @test d.c ≅ [5, 7, missing]

    d = @rselect df begin
        :c = :a + :b
        @passmissing :d = no_missing(:a, :b)
    end
    @test d.c ≅ d.d

    d = @rorderby df @passmissing no_missing(:a, :b)
    @test d ≅ df

    @test_throws LoadError @eval @transform df @passmissing :c = :a + :b
end

end # module


text:
using Revise, StaticArrays, ChainRules, FiniteDifferences, ForwardDiff, Zygote
using OptiMode
using GeometryPrimitives
SHM3 = SHermitianCompact{3,Float64,6}   # static Hermitian 3×3 matrix Type, renamed for brevity

###############################################################################
#                                                                             #
#                                  Rules                                      #
#                                                                             #
###############################################################################

# @Zygote.adjoint (T::Type{<:SMatrix})(xs::Number...) = T(xs...), dv -> (nothing, dv...)
# @Zygote.adjoint (T::Type{<:SMatrix})(x::AbstractMatrix) = T(x), dv -> (nothing, dv)
# @Zygote.adjoint (T::Type{<:SVector})(xs::Number...) = T(xs...), dv -> (nothing, dv...)
# @Zygote.adjoint (T::Type{<:SVector})(x::AbstractVector) = T(x), dv -> (nothing, dv)

@Zygote.adjoint reinterpret(::Type{T}, x::T) where T = x, Δ->(nothing, Δ)
@Zygote.adjoint reinterpret(::Type{NTuple{N,T}}, xs::AbstractArray{T}) where {N,T} = reinterpret(NTuple{N,T}, xs), Δ->(nothing, reinterpret(T, Δ))
@Zygote.adjoint reinterpret(::Type{T}, xs::AbstractArray{NTuple{N,T}}) where {N,T} = reinterpret(T, xs), Δ->(nothing, reinterpret(NTuple{N,T}, Δ))
@Zygote.adjoint reinterpret(::Type{T}, x) where T = reinterpret(T, x), _ -> error("Non-trivial reinterpreting is not supported")

Zygote._zero(xs::AbstractArray{<:AbstractArray}, T=eltype(xs)) = eltype(xs)[_zero(x,eltype(T)) for x in xs]

Zygote.refresh()

###############################################################################
#                                                                             #
#                                 Parameters                                  #
#                                                                             #
###############################################################################

w               =               1.7
t_core          =               0.7
edge_gap        =               0.5             # μm
n_core          =               2.4
n_subs          =               1.4

Δx              =               6.              # μm
Δy              =               4.              # μm
Nx              =               60
Ny              =               40

###############################################################################
#                                                                             #
#                               Generate Data                                 #
#                                                                             #
###############################################################################

g   =   MaxwellGrid(Δx,Δy,Nx,Ny)
s   =   ridge_wg(w,t_core,edge_gap,n_core,n_subs,g)
ε   =   εₛ(s,g);
plot_ε(s,g)

##

function circ_wg2(w,t_core,edge_gap::Float64,n_core::Float64,n_subs::Float64,g::MaxwellGrid)::Array{GeometryPrimitives.Shape{2,4,SHM3},1}
    t_subs = (g.Δy -t_core - edge_gap )/2.
    c_subs_y = -g.Δy/2. + edge_gap/2. + t_subs/2.
    ε_core = ε_tensor(n_core)
    ε_subs = ε_tensor(n_subs)
    ax1,ax2 = GeometryPrimitives.normalize.(([1.,0.], [0.,1.]))
    b_core = GeometryPrimitives.Sphere(					# Instantiate N-D sphere, here N=2 (circle)
                    SVector(0.,t_core),			# c: center
                    w,						# r: "radii" (half span of each axis)
                    ε_core,					        # data: any type, data associated with box shape
                )
    return [b_core,b_subs]
end

function foo(p) #;g=MaxwellGrid(6.,4.,30,20))
    s   =   ridge_wg(p[1],p[2],edge_gap,n_core,n_subs,g)
    tree = KDTree(s)
    ε   =    [εₛ(s,tree,g.x[i],g.y[j],g.δx,g.δy)[k,l] for i=1:g.Nx,j=1:g.Ny,k=1:3,l=1:3]
    return vec(ε)
end

function goo(r) #;g=MaxwellGrid(6.,4.,30,20))
    s   =   [ circ_wg2(r,t_core,edge_gap,n_core,n_subs,g)[1] ]
    tree = KDTree(s)
    ε   =    [εₛ(s,tree,g.x[i],g.y[j],g.δx,g.δy)[k,l] for i=1:g.Nx,j=1:g.Ny,k=1:3,l=1:3]
    return vec(ε)
end

function unflatten(v,Nx,Ny)
    v_rs = reshape(v,(Nx,Ny,3,3))
    return [SMatrix{3,3,Float64,9}(v_rs[i,j,:,:]) for i=1:Nx,j=1:Ny]
end

##
ε = unflatten(foo((1.5,0.8)),Nx,Ny); # returns Array{Float64,4} with size(out)= (Nx,Ny,3,3);
δε_fd = FiniteDifferences.jacobian(central_fdm(3,1),foo,(1.5,0.8))[1];
δw_ε_fd = unflatten(δε_fd[:,1],Nx,Ny);
δt_ε_fd = unflatten(δε_fd[:,2],Nx,Ny);
##

plot_ε(δw_ε_fd,g)
# plot_ε(δt_ε_fd,g)
# plot_ε(ε,g)

##
using ForwardDiff
# ForwardDiff.jacobian(foo,[1.5,0.8])
ForwardDiff.derivative(goo,1.1)

##

ReverseDiff.jacobian(foo,[1.5,0.8])
##
using Zygote
Zygote.gradient(1.1) do a
    sum(goo(a))
end
# Zygote.gradient(1.1) do a
#     Zygote.forwarddiff(a) do a
#         goo(a)
#     end
# end

##

##
ε = unflatten(goo(1.1),Nx,Ny); # returns Array{Float64,4} with size(out)= (Nx,Ny,3,3);
δε_fd = FiniteDifferences.jacobian(central_fdm(3,1),goo,1.1)[1];
δr_ε_fd = unflatten(δε_fd[:,1],Nx,Ny);
##

plot_ε(δr_ε_fd,g)


text:
# Some constant parameters for the Jaynes-Cummings model

# Convenience symbol
const ⊗ = kron

# Physical quantities
const ħ = Complex128(1.05457173E-34) # in Joules*seconds

# Matrices
const σ_z      = Array{Complex128}([1 0 ; 0 -1])
const σ_plus   = Array{Complex128}([0 0 ; 1 0])
const σ_minus  = Array{Complex128}([0 1 ; 0 0])

# Wigner function space
const MAXDISP = 2.5 # maximum displacement operator value, real and imaginary


text:
"Main logging function"
function log(lgr::AbstractLogger, fqmn::AbstractString, lvl::Level.EventLevel, mkr::MARKER, msg, params...)
    if isenabled(lgr, lvl, mkr, msg, params...)
        log(lgr, fqmn, lvl, mkr, message(lgr)(msg, params...))
    end
    return
end
log(lgr::AbstractLogger, lvl::Level.EventLevel, mkr::MARKER, msg, params...) =
    log(lgr, string(current_module()), lvl, mkr, msg, params...)
log(lgr::AbstractLogger, lvl::Level.EventLevel, msg, params...) =
    log(lgr, lvl, MARKER(), msg, params...)

"Logger wrapper for `LoggerConfig`"
type Logger <: AbstractLogger
    name::AbstractString
    message::DataType
    config::LoggerConfig
    filter::FILTER  # This filter comes from Configuration
end
show(io::IO, lgr::Logger) = print(io, name(lgr), ":", level(lgr))
name(lgr::Logger) = lgr.name
message(lgr::Logger) = lgr.message
level(lgr::Logger) = level(lgr.config)
level!(lgr::Logger, lvl::Level.EventLevel) = level!(lgr.config, lvl)

"Logs a message"
log(lgr::Logger, fqmn::AbstractString, lvl::Level.EventLevel, mkr::MARKER, msg::Message) =
    log(lgr.config, name(lgr), fqmn, lvl, mkr, msg)

""" Check if message is valid for logging

Context-wide (configuration) filters and level checks are execuded.
"""
function isenabled(lgr::Logger, lvl::Level.EventLevel, mkr::MARKER, msg, params...)
    if !isnull(lgr.filter)
        r = filter(get(lgr.filter), lvl, mkr, msg)
        r != FilterResult.NEUTRAL && return r == FilterResult.ACCEPT
    end
    level(lgr.config) >= lvl
end

"Simple IO logger"
type SimpleLogger <: AbstractLogger
    name::AbstractString
    message::DataType
    level::Level.EventLevel

    showdatetime::Bool
    showname::Bool
    dateformat::Dates.DateFormat
    io::IO

    function SimpleLogger(nm::AbstractString, msg::DataType, lvl::Level.EventLevel)
        return new(nm, msg, lvl, true, true, Dates.DateFormat("yyyy-mm-ddTHH:MM:SS.sss"), STDERR)
    end
end
show(io::IO, lgr::SimpleLogger) = print(io, name(lgr), ":", lgr |> level |> string)
name(lgr::SimpleLogger) = lgr.name
message(lgr::SimpleLogger) = lgr.message
level(lgr::SimpleLogger) = lgr.level
level!(lgr::SimpleLogger, lvl::Level.EventLevel) = lgr.level = lvl

isenabled(lgr::SimpleLogger, lvl, mkr, msg, params...) = level(lgr) >= lvl

function log(lgr::SimpleLogger, fqmn::AbstractString, lvl::Level.EventLevel, mkr::MARKER, msg::Message)
    lgr.showdatetime && print(lgr.io, Dates.format(Dates.unix2datetime(time()), lgr.dateformat), " ")
    print(lgr.io, string(lvl), " ")
    lgr.showname && print(lgr.io, name(lgr), " ")
    !isnull(mkr) && print(lgr.io, "[", get(mkr), "] ")
    println(lgr.io, msg |> formatted, " ")
end