text:
# These use PascalCase to avoid clashing with source filenames.
const LOG_GROUP_PARSING = :Parsing
const LOG_GROUP_RECORDER = :Recorder
const LOG_GROUP_SERIALIZATION = :Serialization
const LOG_GROUP_SYSTEM = :System
const LOG_GROUP_SYSTEM_CHECKS = :SystemChecks
const LOG_GROUP_TIME_SERIES = :TimeSeries

# Try to keep this updated so that users can check the known groups in the REPL.
const LOG_GROUPS = (
    LOG_GROUP_PARSING,
    LOG_GROUP_RECORDER,
    LOG_GROUP_SERIALIZATION,
    LOG_GROUP_SYSTEM,
    LOG_GROUP_TIME_SERIES,
)
const SIIP_LOGGING_CONFIG_FILENAME =
    joinpath(dirname(pathof(InfrastructureSystems)), "utils", "logging_config.toml")

const LOG_LEVELS = Dict(
    "Debug" => Logging.Debug,
    "Progress" => ProgressLevel,
    "Info" => Logging.Info,
    "Warn" => Logging.Warn,
    "Error" => Logging.Error,
)

"""
Contains information describing a log event.
"""
mutable struct LogEvent
    file::String
    line::Int
    id::Symbol
    message::String
    level::Logging.LogLevel
    count::Int
    suppressed::Int
end

function LogEvent(file, line, id, message, level)
    if isnothing(file)
        file = "None"
    end

    return LogEvent(file, line, id, message, level, 1, 0)
end

struct LogEventTracker
    events::Dict{Logging.LogLevel, Dict{Symbol, LogEvent}}

    # Defining an inner constructor to prohibit creation of a default constructor that
    # takes a parameter of type Any. The outer constructor below causes an overwrite of
    # that method, which results in a warning message from Julia.
    LogEventTracker(events::Dict{Logging.LogLevel, Dict{Symbol, LogEvent}}) = new(events)
end

"""
Returns a summary of log event counts by level.
"""
function report_log_summary(tracker::LogEventTracker)
    text = "\nLog message summary:\n"
    # Order by criticality.
    for level in sort!(collect(keys(tracker.events)), rev = true)
        num_events = length(tracker.events[level])
        text *= "\n$num_events $level events:\n"
        for event in
            sort!(collect(get_log_events(tracker, level)), by = x -> x.count, rev = true)
            text *= "  count=$(event.count) at $(event.file):$(event.line)\n"
            text *= "    example message=\"$(event.message)\"\n"
            if event.suppressed > 0
                text *= "    suppressed=$(event.suppressed)\n"
            end
        end
    end

    return text
end

"""
Returns an iterable of log events for a level.
"""
function get_log_events(tracker::LogEventTracker, level::Logging.LogLevel)
    if !_is_level_valid(tracker, level)
        return []
    end

    return values(tracker.events[level])
end

"""
Increments the count of a log event.
"""
function increment_count!(tracker::LogEventTracker, event::LogEvent, suppressed::Bool)
    if _is_level_valid(tracker, event.level)
        if haskey(tracker.events[event.level], event.id)
            tracker.events[event.level][event.id].count += 1
            if suppressed
                tracker.events[event.level][event.id].suppressed += 1
            end
        else
            tracker.events[event.level][event.id] = event
        end
    end
end

function _is_level_valid(tracker::LogEventTracker, level::Logging.LogLevel)
    return level in keys(tracker.events)
end

struct LoggingConfiguration
    console::Bool
    console_stream::IO
    console_level::Base.LogLevel
    progress::Bool
    file::Bool
    filename::Union{Nothing, String}
    file_level::Base.LogLevel
    file_mode::String
    tracker::Union{Nothing, LogEventTracker}
    set_global::Bool
    group_levels::Dict{Symbol, Base.LogLevel}
end

function LoggingConfiguration(;
    console = true,
    console_stream = stderr,
    console_level = Logging.Error,
    progress = true,
    file = true,
    filename = "log.txt",
    file_level = Logging.Info,
    file_mode = "w+",
    tracker = nothing,
    set_global = true,
    group_levels = Dict{Symbol, Base.LogLevel}(),
)
    return LoggingConfiguration(
        console,
        console_stream,
        console_level,
        progress,
        file,
        filename,
        file_level,
        file_mode,
        tracker,
        set_global,
        group_levels,
    )
end

function LoggingConfiguration(config_filename)
    config = open(config_filename, "r") do io
        return TOML.parse(io)
    end

    console_stream_str = get(config, "console_stream", "stderr")
    if console_stream_str == "stderr"
        config["console_stream"] = stderr
    elseif console_stream_str == "stdout"
        config["console_stream"] = stdout
    else
        error("unsupport console_stream value: {console_stream_str")
    end

    config["console_level"] = get_logging_level(get(config, "console_level", "Info"))
    config["file_level"] = get_logging_level(get(config, "file_level", "Info"))
    config["group_levels"] =
        Dict(Symbol(k) => get_logging_level(v) for (k, v) in config["group_levels"])
    config["tracker"] = nothing
    return LoggingConfiguration(; Dict(Symbol(k) => v for (k, v) in config)...)
end

function make_logging_config_file(filename = "logging_config.toml"; force = false)
    cp(SIIP_LOGGING_CONFIG_FILENAME, filename, force = force)
    println("Created $filename")
    return
end

"""
Tracks counts of all log events by level.

# Examples
```Julia
LogEventTracker()
LogEventTracker((Logging.Info, Logging.Warn, Logging.Error))
```
"""
function LogEventTracker(levels = (Logging.Info, Logging.Warn, Logging.Error))
    return LogEventTracker(Dict(l => Dict{Symbol, LogEvent}() for l in levels))
end

"""
Creates console and file loggers per caller specification and returns a MultiLogger.

Suppress noisy events by specifying per-event values of `maxlog = X` and
`_suppression_period = Y` where X is the max number of events that can occur in Y
seconds. After the period ends, messages will no longer be suppressed. Note that
if you don't specify `_suppression_period` then `maxlog` applies for the for the
duration of your process (standard Julia logging behavior).

**Note:** Use of log message suppression and the LogEventTracker are not thread-safe.
Please contact the package developers if you need this functionality.

**Note:** If logging to a file users must call Base.close() on the returned MultiLogger to
ensure that all events get flushed.

# Arguments
- `console::Bool=true`: create console logger
- `console_stream::IOStream=stderr`: stream for console logger
- `console_level::Logging.LogLevel=Logging.Error`: level for console messages
- `progress::Bool=true`: enable progress logger
- `file::Bool=true`: create file logger
- `filename::Union{Nothing, String}=log.txt`: log file
- `file_level::Logging.LogLevel=Logging.Info`: level for file messages
- `file_mode::String=w+`: mode used when opening log file
- `tracker::Union{LogEventTracker, Nothing}=LogEventTracker()`: optionally track log events
- `set_global::Bool=true`: set the created logger as the global logger

# Example
```Julia
logger = configure_logging(filename="mylog.txt")
@info "hello world"
@info "hello world" maxlog = 5 _suppression_period = 10
```
"""
function configure_logging(;
    console = true,
    console_stream = stderr,
    console_level = Logging.Error,
    progress = true,
    file = true,
    filename = "log.txt",
    file_level = Logging.Info,
    file_mode = "w+",
    tracker = LogEventTracker(),
    set_global = true,
)
    config = LoggingConfiguration(
        console = console,
        console_stream = console_stream,
        console_level = console_level,
        progress = progress,
        file = file,
        filename = filename,
        file_level = file_level,
        file_mode = file_mode,
        tracker = tracker,
        set_global = set_global,
    )
    return configure_logging(config)
end

function configure_logging(config_filename::AbstractString)
    return configure_logging(LoggingConfiguration(config_filename))
end

function configure_logging(config::LoggingConfiguration)
    if !config.console && !config.file && !config.progress
        error("At least one of console, file, or progress must be true")
    end

    loggers = Vector{Logging.AbstractLogger}()
    if config.console
        console_logger = TerminalLogger(config.console_stream, config.console_level)
        push!(loggers, console_logger)
    end

    if config.progress
        progress_logger = ProgressLogger(config.console_stream, ProgressLevel)
        push!(loggers, progress_logger)
    end

    if config.file
        io = open(config.filename, config.file_mode)
        file_logger = FileLogger(io, config.file_level)
        push!(loggers, file_logger)
    end

    logger = MultiLogger(
        loggers,
        config.tracker,
        Dict{Symbol, Base.LogLevel}(),
        LogSuppressionTracker(),
    )
    if config.set_global
        Logging.global_logger(logger)
    end

    return logger
end

"""
Specializes the behavior of SimpleLogger by adding timestamps and process and thread IDs.
"""
struct FileLogger <: Logging.AbstractLogger
    logger::Logging.SimpleLogger
end

function FileLogger(stream::IO, level::Base.CoreLogging.LogLevel)
    return FileLogger(Logging.SimpleLogger(stream, level))
end

function Logging.handle_message(
    file_logger::FileLogger,
    level,
    message,
    _module,
    group,
    id,
    file,
    line;
    maxlog = nothing,
    kwargs...,
)
    return Logging.handle_message(
        file_logger.logger,
        level,
        "$(Dates.now()) [$(getpid()):$(Base.Threads.threadid())]: $message",
        _module,
        group,
        id,
        file,
        line;
        maxlog = maxlog,
        kwargs...,
    )
end

function Logging.shouldlog(logger::FileLogger, level, _module, group, id)
    return Logging.shouldlog(logger.logger, level, _module, group, id)
end

Logging.min_enabled_level(logger::FileLogger) = Logging.min_enabled_level(logger.logger)
Logging.catch_exceptions(logger::FileLogger) = false
Base.flush(logger::FileLogger) = flush(logger.logger.stream)
Base.close(logger::FileLogger) = close(logger.logger.stream)

"""
Opens a file logger using Logging.SimpleLogger.

# Example
```Julia
open_file_logger("log.txt", Logging.Info) do logger
    global_logger(logger)
    @info "hello world"
end
```
"""
function open_file_logger(
    func::Function,
    filename::String,
    level = Logging.Info,
    mode = "w+",
)
    stream = open(filename, mode)
    try
        logger = FileLogger(stream, level)
        func(logger)
    finally
        close(stream)
    end
end

mutable struct LogEventSuppressionStats
    "The event being tracked"
    event_id::Symbol
    "The time we started tracking an event."
    tracking_start_time::Float64
    "Number of times an event has occurred since tracking started."
    count::Int
    "Length in seconds of one tracking period"
    period::Int
    "Number of times an event has been suppressed since tracking started."
    num_suppressed::Int
    "Whether tracking is active"
    is_tracking_active::Bool
    "Whether the event is being suppressed"
    is_suppression_enabled::Bool
end

function LogEventSuppressionStats(event_id::Symbol, cur_time::Float64, period::Int)
    return LogEventSuppressionStats(event_id, cur_time, 0, period, 0, false, false)
end

function start_tracking!(stats::LogEventSuppressionStats, cur_time::Float64)
    stats.is_tracking_active = true
    stats.tracking_start_time = cur_time
    stats.count = 1
    stats.is_suppression_enabled = false
    stats.num_suppressed = 0
    return
end

function increment!(stats::LogEventSuppressionStats, maxlog::Int)
    @assert stats.is_tracking_active
    stats.count += 1
    if stats.count > maxlog
        if !stats.is_suppression_enabled
            stats.is_suppression_enabled = true
        end
        stats.num_suppressed += 1
    end
    return
end

function should_suppress!(stats::LogEventSuppressionStats, cur_time::Float64, maxlog::Int)
    diff = cur_time - stats.tracking_start_time
    num_suppressed = 0
    if !stats.is_tracking_active
        start_tracking!(stats, cur_time)
    elseif cur_time - stats.tracking_start_time < stats.period
        increment!(stats, maxlog)
        num_suppressed = stats.num_suppressed
    else
        # The suppression period has ended.
        num_suppressed = stats.num_suppressed
        start_tracking!(stats, cur_time)
    end

    return stats.is_suppression_enabled, num_suppressed
end

struct LogSuppressionTracker
    event_stats::Dict{Symbol, LogEventSuppressionStats}
end

function LogSuppressionTracker()
    return LogSuppressionTracker(Dict{Symbol, LogEventSuppressionStats}())
end

should_suppress!(tracker::LogSuppressionTracker, ::Any, ::Any, ::Nothing) = false, 0

function should_suppress!(
    tracker::LogSuppressionTracker,
    event_id::Symbol,
    maxlog::Int,
    suppression_period::Int,
)
    cur_time = time()
    stats = get!(
        tracker.event_stats,
        event_id,
        LogEventSuppressionStats(event_id, cur_time, suppression_period),
    )
    return should_suppress!(stats, cur_time, maxlog)
end

"""
Redirects log events to multiple loggers. The primary use case is to allow logging to
both a file and the console. Secondarily, it can track the counts of all log messages.

# Example
```Julia
MultiLogger([TerminalLogger(stderr), SimpleLogger(stream)], LogEventTracker())
```
"""
mutable struct MultiLogger <: Logging.AbstractLogger
    loggers::Array{Logging.AbstractLogger}
    tracker::Union{LogEventTracker, Nothing}
    group_levels::Dict{Symbol, Base.LogLevel}
    suppression_tracker::LogSuppressionTracker
end

"""
Creates a MultiLogger with no event tracking.

# Example
```Julia
MultiLogger([TerminalLogger(stderr), SimpleLogger(stream)])
```
"""
function MultiLogger(loggers::Array{T}) where {T <: Logging.AbstractLogger}
    return MultiLogger(
        loggers,
        nothing,
        Dict{Symbol, Base.LogLevel}(),
        LogSuppressionTracker(),
    )
end

function MultiLogger(
    loggers::Array{T},
    tracker::LogEventTracker,
) where {T <: Logging.AbstractLogger}
    return MultiLogger(
        loggers,
        tracker,
        Dict{Symbol, Base.LogLevel}(),
        LogSuppressionTracker(),
    )
end

function Logging.shouldlog(logger::MultiLogger, level, _module, group, id)
    return get(logger.group_levels, group, level) <= level
end

function Logging.min_enabled_level(logger::MultiLogger)
    return minimum([Logging.min_enabled_level(x) for x in logger.loggers])
end

Logging.catch_exceptions(logger::MultiLogger) = false

function Logging.handle_message(
    logger::MultiLogger,
    level::Int,
    message,
    _module,
    group,
    id,
    file,
    line;
    kwargs...,
)
    return Logging.handle_message(
        logger::MultiLogger,
        Logging.LogLevel(level),
        message,
        _module,
        group,
        id,
        file,
        line;
        kwargs...,
    )
end

function Logging.handle_message(
    logger::MultiLogger,
    level::Logging.LogLevel,
    message,
    _module,
    group,
    id,
    file,
    line;
    maxlog = nothing,
    _suppression_period = nothing,
    kwargs...,
)
    suppressed, num_suppressed =
        should_suppress!(logger.suppression_tracker, id, maxlog, _suppression_period)
    if !suppressed
        # Takeover maxlog if suppression_period is set.
        maxlog = _suppression_period === nothing ? maxlog : nothing
        for _logger in logger.loggers
            if level >= Logging.min_enabled_level(_logger)
                if Logging.shouldlog(_logger, level, _module, group, id)
                    if num_suppressed > 0
                        kwargs =
                            merge(Dict(kwargs), Dict(:num_suppressed => num_suppressed))
                    end
                    # @show typeof(_logger)
                    Logging.handle_message(
                        _logger,
                        level,
                        message,
                        _module,
                        group,
                        id,
                        file,
                        line;
                        maxlog = maxlog,
                        kwargs...,
                    )
                end
            end
        end
    end

    if !isnothing(logger.tracker)
        id = isa(id, Symbol) ? id : :empty
        event = LogEvent(file, line, id, string(message), level)
        increment_count!(logger.tracker, event, suppressed)
    end

    return
end

"""
Empty the minimum log levels stored for each group.
"""
function empty_group_levels!(logger::MultiLogger)
    empty!(logger.group_levels)
    return
end

"""
Set the minimum log level for a group.

The `group` field of a log message defaults to its file's base name (no extension) as a
symbol. It can be customized by setting `_group = :a_group_name`.

The minimum log level stored for a console or file logger supercede this setting.
"""
function set_group_level!(logger::MultiLogger, group::Symbol, level::Base.LogLevel)
    logger.group_levels[group] = level
    return
end

"""
Set the minimum log levels for multiple groups. Refer to [`set_group_level`](@ref) for more
information.
"""
function set_group_levels!(logger::MultiLogger, group_levels::Dict{Symbol, Base.LogLevel})
    merge!(logger.group_levels, group_levels)
    return
end

"""
Return the minimum logging level for a group or nothing if `group` is not stored.
"""
get_group_level(logger::MultiLogger, group::Symbol) =
    get(logger.group_levels, group, nothing)

"""
Return the minimum logging levels for groups that have been stored.
"""
get_group_levels(logger::MultiLogger) = deepcopy(logger.group_levels)

"""
Returns a summary of log event counts by level.
"""
function report_log_summary(logger::MultiLogger)
    if isnothing(logger.tracker)
        error("log event tracking is not enabled")
    end

    return report_log_summary(logger.tracker)
end

"""
Flush any file streams.
"""
function Base.flush(logger::MultiLogger)
    return _handle_log_func(logger, Base.flush)
end

"""
Ensures that any file streams are flushed and closed.
"""
function Base.close(logger::MultiLogger)
    return _handle_log_func(logger, Base.close)
end

function get_logging_level(level::String)
    if !haskey(LOG_LEVELS, level)
        error("Invalid log level $level: Supported levels: $(values(LOG_LEVELS))")
    end

    return LOG_LEVELS[level]
end

function _handle_log_func(logger::MultiLogger, func::Function)
    for _logger in logger.loggers
        if isa(_logger, Logging.SimpleLogger)
            func(_logger.stream)
        elseif isa(_logger, FileLogger)
            func(_logger)
        end
    end
end

mutable struct ProgressLogger <: Logging.AbstractLogger
    logger::TerminalLogger
end

ProgressLogger(args...) = ProgressLogger(TerminalLogger(args...))

function Logging.shouldlog(logger::ProgressLogger, level, _module, group, id)
    return level == ProgressLevel
end

Logging.min_enabled_level(x::ProgressLogger) = Logging.min_enabled_level(x.logger)
Logging.catch_exceptions(x::ProgressLogger) = Logging.catch_exceptions(x.logger)

function Logging.handle_message(x::ProgressLogger, args...; kwargs...)
    Logging.handle_message(x.logger, args...; kwargs...)
end


text:
abstract type
    HydroGen <: Generator
end

struct TechHydro <: TechnicalParams
    installedcapacity::Float64
    activepower::Float64 # [MW]
    activepowerlimits::NamedTuple{(:min, :max),Tuple{Float64,Float64}} # [MW]
    reactivepower::Union{Float64,Nothing} # [MVAr]
    reactivepowerlimits::Union{NamedTuple{(:min, :max),Tuple{Float64,Float64}},Nothing} # [MVAr]
    ramplimits::Union{NamedTuple{(:up, :down),Tuple{Float64,Float64}},Nothing} #MW/Hr
    timelimits::Union{NamedTuple{(:up, :down),Tuple{Float64,Float64}},Nothing} # Hrs
    function TechHydro(installedcapacity, activepower, activepowerlimits, reactivepower, reactivepowerlimits, ramplimits, timelimits)

        new(installedcapacity, activepower, PowerSystems.orderedlimits(activepowerlimits, "Real Power"), reactivepower, PowerSystems.orderedlimits(reactivepowerlimits, "Reactive Power"), ramplimits, timelimits)

    end
end

TechHydro(;installedcapacity = 0.0,
          activepower = 0.0,
          activepowerlimits = (min = 0.0, max = 0.0),
          reactivepower = nothing,
          reactivepowerlimits = nothing,
          ramplimits = nothing,
          timelimits = nothing
        ) = TechHydro(installedcapacity, activepower, activepowerlimits, reactivepower, reactivepowerlimits, ramplimits, timelimits)

struct EconHydro <: TechnicalParams
    curtailpenalty::Float64 # [$/MWh]
    variablecost::Union{Float64,Nothing} # [$/MWh]
end

EconHydro(; cost = 0.0, curtailcost = 0.0) = EconHydro(cost, curtailcost)

struct HydroFix <: HydroGen
    name::String
    available::Bool
    bus::Bus
    tech::TechHydro
    scalingfactor::TimeSeries.TimeArray
end

HydroFix(; name="init",
        status = false,
        bus = Bus(),
        tech = TechHydro(),
        scalingfactor = TimeSeries.TimeArray(Dates.today(),ones(1))) = HydroFix(name, status, bus, tech, scalingfactor)

struct HydroCurtailment <: HydroGen
    name::String
    available::Bool
    bus::Bus
    tech::TechHydro
    econ::Union{EconHydro,Nothing}
    scalingfactor::TimeSeries.TimeArray # [0-1]
    function HydroCurtailment(name, status, bus, tech, curtailcost::Float64, scalingfactor)
        econ = EconHydro(curtailcost, nothing)
        new(name, status, bus, tech, econ, scalingfactor)
    end
end

HydroCurtailment(; name = "init",
                status = false,
                bus= Bus(),
                tech = TechHydro(),
                curtailcost = 0.0,
                scalingfactor = TimeSeries.TimeArray(Dates.today(),ones(1))) = HydroCurtailment(name, status, bus, tech, curtailcost, scalingfactor)

struct HydroStorage <: HydroGen
    name::String
    available::Bool
    bus::Bus
    tech::TechHydro
    econ::Union{EconHydro,Nothing}
    storagecapacity::Float64 #[m^3]
    scalingfactor::TimeSeries.TimeArray
end

HydroStorage(; name = "init",
                status = false,
                bus= Bus(),
                tech = TechHydro(),
                econ = EconHydro(),
                storagecapacity = 0.0,
                scalingfactor = TimeSeries.TimeArray(Dates.today(),ones(1))) = HydroStorage(name, status, bus, tech, econ, storagecapacity, scalingfactor)

ext: jl


text:
# Horizontal generation #
# ========================== #
function generate!(stochasticprogram::StochasticProgram{N}, structure::HorizontalStructure{N}) where N
    # Generate all stages
    for stage in 1:N
        generate!(stochasticprogram, structure, stage)
    end
    return nothing
end

function generate!(stochasticprogram::TwoStageStochasticProgram, structure::HorizontalStructure{2}, stage::Integer)
    if stage == 1
        # Check generators
        has_generator(stochasticprogram, :stage_1) || error("First-stage problem not defined in stochastic program. Consider @stage 1.")
        # Prepare decisions
        structure.proxy.ext[:decisions] = structure.decisions[1]
        add_decision_bridges!(structure.proxy)
        # Generate first stage
        generator(stochasticprogram, :stage_1)(structure.proxy, stage_parameters(stochasticprogram, 1))
    else
        stage == 2 || error("Stage $stage not available in two-stage model.")
        # Check generators
        has_generator(stochasticprogram, :stage_1) || error("First-stage problem not defined in stochastic program. Consider @stage 1.")
        has_generator(stochasticprogram, :stage_2) || error("Second-stage problem not defined in stochastic program. Consider @stage 2.")
        # Sanity check on scenario probabilities
        if num_scenarios(structure, stage) > 0
            p = stage_probability(structure, stage)
            abs(p - 1.0) <= 1e-6 || @warn "Scenario probabilities do not add up to one. The probability sum is given by $p"
        end
        # Generate
        generate_horizontal!(scenarioproblems(structure, stage),
                             generator(stochasticprogram, :stage_1),
                             generator(stochasticprogram, :stage_2),
                             stage_parameters(stochasticprogram, 1),
                             stage_parameters(stochasticprogram, 2),
                             structure.decisions[stage - 1],
                             subproblem_optimizer(stochasticprogram))
    end
    return nothing
end

function generate_horizontal!(scenarioproblems::ScenarioProblems,
                              stage_one_generator::Function,
                              stage_two_generator::Function,
                              stage_one_params::Any,
                              stage_two_params::Any,
                              decisions::Decisions,
                              optimizer)
    for i in num_subproblems(scenarioproblems)+1:num_scenarios(scenarioproblems)
        push!(scenarioproblems.problems, _WS(stage_one_generator,
                                             stage_two_generator,
                                             stage_one_params,
                                             stage_two_params,
                                             scenario(scenarioproblems,i),
                                             decisions,
                                             optimizer))
    end
    return nothing
end
function generate_horizontal!(scenarioproblems::DistributedScenarioProblems,
                              stage_one_generator::Function,
                              stage_two_generator::Function,
                              stage_one_params::Any,
                              stage_two_params::Any,
                              decisions::Decisions,
                              optimizer)
    @sync begin
        for w in workers()
            @async remotecall_fetch(
                w,
                scenarioproblems[w-1],
                stage_one_generator,
                stage_two_generator,
                stage_one_params,
                stage_two_params,
                scenarioproblems.decisions[w-1],
                optimizer) do sp, gen_one, gen_two, one_params, two_params, decisions, opt
                    generate_horizontal!(fetch(sp),
                                         gen_one,
                                         gen_two,
                                         one_params,
                                         two_params,
                                         fetch(decisions),
                                         opt)
                end
        end
    end
    return nothing
end

function clear!(structure::HorizontalStructure{N}) where N
    # Clear all stages
    for stage in 1:N
        clear_stage!(structure, stage)
    end
    return nothing
end

function clear_stage!(structure::HorizontalStructure{N}, s::Integer) where N
    1 <= s <= N || error("Stage $s not in range 1 to $N.")
    clear!(scenarioproblems(structure, s))
    return nothing
end


text:
function update_weights(data, e, closest_same, closest_other, weights, m, max_f_vals, min_f_vals)
	weights = weights .- (abs.(e .- closest_same)./(max_f_vals .- min_f_vals .+ eps(Float64)))./m .+
		(abs.(e .- closest_other)./(max_f_vals .- min_f_vals .+ eps(Float64)))./m
	return weights
end

max_stars_repo_name: Architecton/diplomska-naloga-koda

max_issues_repo_name: Architecton/diplomska-naloga-koda

max_forks_repo_name: Architecton/diplomska-naloga-koda


text:
if isdefined(Base, :exponent_bits)
    using Base: exponent_bits
else
    exponent_bits(::Type{Float64}) = 11
    exponent_bits(::Type{Float32}) =  8
    exponent_bits(::Type{Float16}) =  5
end
if isdefined(Base, :significand_bits)
    using Base: significand_bits
else
    significand_bits(::Type{Float64}) = 52
    significand_bits(::Type{Float32}) = 23
    significand_bits(::Type{Float16}) = 10
end


text:
using DataVoyager
using DataValues
using Electron
using Test

Electron.prep_test_env()

@testset "DataVoyager" begin

    source = [(a = 1, b = 1), (a = 2, b = 2)]
    source2 = [(a = DataValue(1), b = DataValue{Int}()), (a = DataValue{Int}(), b = DataValue(2))]

    v = Voyager()
    @test typeof(v.w) == Electron.Window

    v = Voyager(source)
    @test typeof(v.w) == Electron.Window

    v(source)
    @test typeof(v.w) == Electron.Window

    source |> v
    @test typeof(v.w) == Electron.Window

    v = Voyager(source2)
    @test typeof(v.w) == Electron.Window

    v(source2)
    @test typeof(v.w) == Electron.Window

    source2 |> v
    @test typeof(v.w) == Electron.Window

end


text:
macro plot(plottable_table, plot_call = :(plot()), err_style = :(:yerr))
    if @capture(plot_call, fun_(as__; kws__))
        return esc(Expr(:call, :(GroupedErrors.plot_helper), plottable_table, fun, err_style, as..., kws...))
    elseif @capture(plot_call, fun_(as__))
        return esc(Expr(:call, :(GroupedErrors.plot_helper), plottable_table, fun, err_style, as...))
    else
        return esc(plottable_table)
    end
end

function plot_helper(t::ProcessedTable, f::Function, err_style, args...; kwargs...)
    s = t.table
    axis_type = t.kw[:axis_type]
    group = Tuple(columns(t.table, s) for s in listsplits(t))
    group_list = unique(Columns(group))
    plot_kwargs = [(plot_kwarg[1], reshape(plot_kwarg[2].(group_list), 1, :))
        for plot_kwarg in t.kw[:plot_kwargs]]
    if axis_type == :pointbypoint
        err_kwargs = []
        if eltype(columns(s, :x))<:Tuple
            x = getindex.(columns(s, :x), 1)
            push!(err_kwargs, (:xerr, getindex.(columns(s, :x), 2)))
        else
            x = columns(s, :x)
        end
        if eltype(columns(s, :y))<:Tuple
            y = getindex.(columns(s, :y), 1)
            push!(err_kwargs, (:yerr, getindex.(columns(s, :y), 2)))
        else
            y = columns(s, :y)
        end
        f(args..., x, y; group = group, err_kwargs...,kwargs..., plot_kwargs...)
    else
        x, y = columns(s, :x), columns(s, :y)
        err = :err in colnames(s) ? columns(s, :err) : nothing
        (axis_type == :continuous) && (err_style = :ribbon)
        kw_itr = Iterators.flatten([[(err_style, err)], kwargs, plot_kwargs])
        f(args..., x, y; group = group, kw_itr...)
    end
end

plot_helper(t, f::Function, err_style, args...; kwargs...) =
    plot_helper(ProcessedTable(t), f, err_style, args...; kwargs...)


text:
using Plots, MAT

nprocs = 4

@views function vizme1D_mpi(nprocs)
    H = []
    for ip = 1:nprocs
        file = matopen("H_$(ip-1).mat"); H_loc = read(file, "H"); close(file)
        nx_i = length(H_loc)-2
        i1   = 1 + (ip-1)*nx_i
        if (ip==1)  H = zeros(nprocs*nx_i)  end
        H[i1:i1+nx_i-1] .= H_loc[2:end-1]
    end
    fontsize = 12
    display(plot(H, legend=false, framestyle=:box, linewidth=3, xlims=(1, length(H)), ylims=(0, 1), xlabel="nx", title="diffusion 1D MPI", yaxis=font(fontsize, "Courier"), xaxis=font(fontsize, "Courier"), titlefontsize=fontsize, titlefont="Courier"))
    return
end

vizme1D_mpi(nprocs)


text:
# VEM computes the virtual element solution of a Poisson problem on a polygonal mesh
# Translation to Julia of MATLAB code from paper:
# O. Suttton, The virtual element method in 50 lines of MATLAB, Numer Algor(2017) 75:1141-1159

using jFEMTools

import Tensors
using SparseArrays
using LinearAlgebra

mesh = rectangle_mesh(RectangleCell, (2,2), Tensors.Vec{2}((0.0,0.0)), Tensors.Vec{2}((1.0,1.0)));

function compute_u(mesh)
	# forcing function
	rhs(x::Tensors.Vec{2}) = 2*π^2*sin(π*x[1])*sin(π*x[2])
	# Boundary condition
	g(x::Tensors.Vec{2}) = sin(π * x[2])*sin(π * x[1]);
	n_dofs = getnvertices(mesh)
	n_polys = 3; # Method has 1 degree of freedom per vertex
	assembler = start_assemble(n_dofs)
	F = zeros(n_dofs); # Forcing vector
	u = zeros(n_dofs); # Degrees of freedom of the virtual element solution
	linear_polynomials = ((0,0), (1,0), (0,1)); # Impose an ordering on the linear polynomials
	for el_id = 1:getncells(mesh)
		vert_ids = getverticesidx(mesh, el_id) # Global IDs of the vertices of this element
		verts = getverticescoords(mesh, el_id) # Coordinates of the vertices of this element
		n_sides = getnedges(mesh.cells[el_id]) # Start computing the geometric information
		area = cell_volume(mesh, el_id)
		centr = cell_centroid(mesh, el_id)
		diameter = cell_diameter(mesh, el_id)
		D = zeros(n_sides, n_polys); D[:, 1] .= 1;
		B = zeros(n_polys, n_sides); B[1, :] .= 1/n_sides;
		for vertex_id = 1:n_sides
			vert = verts[vertex_id]; #This vertex and its neighbours
			prev = verts[mod1(vertex_id - 1, n_sides)];
			nextv = verts[mod1(vertex_id + 1, n_sides)];
			vertex_normal = [nextv[2] - prev[2], prev[1] - nextv[1]]; # Average of edge normals
			for poly_id in 2:n_polys # Only need to loop over non-constant polynomials
				poly_degree = linear_polynomials[poly_id];
				monomial_grad = poly_degree ./ diameter; # Gradient of a linear is constant
				D[vertex_id, poly_id] = dot(vert - centr, poly_degree) / diameter;
				B[poly_id, vertex_id] = 0.5 * dot(monomial_grad, vertex_normal);
			end
		end
		projector = (B*D) \ B; # Compute the local Ritz projector to polynomials
		stabilising_term = (I - D * projector)' * (I - D * projector);
		G = B*D; G[1, :] .= 0;
		local_stiffness = projector' * G * projector + stabilising_term;
		assemble!(assembler,vert_ids,local_stiffness)
		assemble!(F, [dof for dof in vert_ids], rhs(centr)*ones(length(vert_ids)) * area / n_sides)
	end
	K = end_assemble(assembler)
	boundary_nodes = getvertexset(mesh,"boundary")
	boundary_vals = Vector{Float64}(undef, length(boundary_nodes))
	for (i,node) in enumerate(boundary_nodes)
		boundary_vals[i] = g(getvertexcoords(mesh, node))
	end

	internal_dofs = setdiff(1:n_dofs, boundary_nodes)
	F = F - K[:, collect(boundary_nodes)] * boundary_vals; # Apply the boundary condition
	u[internal_dofs] = K[collect(internal_dofs), collect(internal_dofs)] \ F[collect(internal_dofs)]; # Solve
	u[collect(boundary_nodes)] = boundary_vals; # Set the boundary values
	u
end
u = compute_u(mesh)

#plot solution
using Makie
import AbstractPlotting
using ColorSchemes
include("../src/plot_recipes.jl")
scene = Scene(resolution = (400, 200), colormap = ColorSchemes.viridis.colors)
plot!(scene, mesh, color = u)
display(scene)

max_stars_repo_name: Paulms/jFEMToools

max_issues_repo_name: Paulms/jFEMToools

max_forks_repo_name: Paulms/jFEMToools


text:
const LAYERINCREMENTS = (0.0, 0.025, 0.05, 0.1, 0.2, 0.3, 0.5, 1.0) .* m
const LAYERNAMES = ("0", "2.5", "5", "10", "20", "30", "50", "100")
const NEXTLAYER = 2.0m
const LAYERRANGE = (0.01, 1.20) .* m

" Load multiple years for an index "
load_point(basepath::String, years::AbstractVector, shade::Int, i::CartesianIndex, skip=()) =
    MicroclimPoint{shade,LAYERINCREMENTS,NEXTLAYER,LAYERRANGE}(load_microclim(basepath, years, shade, skip, i)...)

" Load grids for a single year "
load_grid(basepath::String, years::AbstractVector, shade::Int, skip=()) =
    MicroclimGrid{shade,LAYERINCREMENTS,NEXTLAYER,LAYERRANGE}(load_microclim(basepath, years, shade, skip)...)

ifhasdir(loader, basepath, file, empty) =
    if isdir(joinpath(basepath, file))
        loader()
    else
        println("$file does not exits. Skipping and returning $empty.")
        empty
    end

"""
Generic microcliate loader. An Int year and no args will return a MicroclimGrid
with full raster layers for asingle year. A range of years within 1990:2017 and a
CartesianIndex in args will return a MicroclimPoint.
"""
load_microclim(basepath, years, shade, skip, args...) = begin
    check_year(years)

    radiation = :radiation in skip ? [] : load_file(basepath, "SOLR", years, args...)

    snowdepth = :snowdepth in skip ? [] :
        load_file((n, y) -> "$(n)/$(shade)pctShade/$(n)_$(shade)pctShade_$(y).nc", basepath, "SNOWDEP", years, args...)

    airtemperature = :airtemperature in skip ? () :
        load_file((n, y) -> "Tair1cm/$(shade)pctShade/$(n)_$(shade)pctShade_$(y).nc", basepath, "TA1cm", years, args...),
        load_file(basepath, "TA120cm", years, args...)

    relhumidity = :relhumidity in skip ? () :
        load_file((n, y) -> "$(n)/$(shade)pctShade/$(n)_$(shade)pctShade_$(y).nc", basepath, "RH1cm", years, args...),
        load_file(basepath, "RH120cm", years, args...)

    windspeed = :windspeed in skip ? () :
        load_file(basepath, "V1cm", years, args...),
        load_file(basepath, "V120cm", years, args...)

    soiltemperature = :soiltemperature in skip ? () : load_folder(basepath, "soil", LAYERNAMES, years, shade, args...)
    soilwaterpotential = :soilwaterpotential in skip ? () : load_folder(basepath, "pot", LAYERNAMES, years, shade, args...)
    soilwatercontent = :soilwatercontent in skip ? () : load_folder(basepath, "moist", LAYERNAMES, years, shade, args...)

    radiation, snowdepth, airtemperature, relhumidity, windspeed,
    soiltemperature, soilwaterpotential, soilwatercontent
end

check_year(year::Integer) = year >= 1990 && year <= 2017 || error("$year is not between 1990 and 2017")
check_year(years) = check_year(years.start) && check_year(years.stop)

load_file(formatter::Function, basepath::String, name::String, years, args...) =
    load_variable(formatter, basepath, name, name, years, args...)
load_file(basepath::String, name::String, years, args...) = begin
    formatter = (n, y) -> "$(n)/$(n)_$(y).nc"
    load_variable(formatter, basepath, name, name, years, args...)
end

load_folder(basepath, name, layers, years, shade, args...) =
    ((load_variable(basepath, name, "$(name)$(l)cm", years, args...) do n, y
         "$(n)$(l)cm/$(shade)pctShade/$(n)$(l)cm_$(shade)pctShade_$(y).nc"
     end for l in layers)...,)

"load timeseries for a single index"
load_variable(formatter::Function, basepath::String, name::String, varname::String, years,  i::CartesianIndex) = begin
    println("Adding $varname for $years at $i")
    combined_data = Float64[]
    for year in years
        path = joinpath(basepath, formatter(name, year))
        year_data = [get_yeardata(path, varname)[i, :]...]
        combined_data = vcat(combined_data, year_data)
    end
    combined_data
end

"load timeseries for the whole grid"
load_variable(formatter::Function, basepath::String, name::String, varname::String, years) = begin
    println("Adding $varname for $years for complete grid")

    # Initialise for the first year, as the eltype is variable
    path1 = joinpath(basepath, formatter(name, years[1]))
    years_vector = [get_yeardata(path1, varname)]

    # Add the rest of the years, if there are any
    for i in 2:length(years)
        path = joinpath(basepath, formatter(name, years[i]))
        push!(years_vector, get_yeardata(path, varname))
    end

    # Return a vector of NCDatasets for each year
    years_vector
end

function get_yeardata(path, varname)
    println("Retrieving ", varname, " from ", path)
    Dataset(path)[varname][:,:,:]
end


text:
export alazar_exception
export alazar_no_error
const alazar_no_error = 512

exceptions = [
    "ApiNoError",                             # 512
    "ApiFailed",                              # 513
    "ApiAccessDenied",                        # 514
    "ApiDmaChannelUnavailable",               # 515
    "ApiDmaChannelInvalid",                   # 516
    "ApiDmaChannelTypeError",                 # 517
    "ApiDmaInProgress",                       # 518
    "ApiDmaDone",                             # 519
    "ApiDmaPaused",                           # 520
    "ApiDmaNotPaused",                        # 521
    "ApiDmaCommandInvalid",                   # 522
    "ApiDmaManReady",                         # 523
    "ApiDmaManNotReady",                      # 524
    "ApiDmaInvalidChannelPriority",           # 525
    "ApiDmaManCorrupted",                     # 526
    "ApiDmaInvalidElementIndex",              # 527
    "ApiDmaNoMoreElements",                   # 528
    "ApiDmaSglInvalid",                       # 529
    "ApiDmaSglQueueFull",                     # 530
    "ApiNullParam",                           # 531
    "ApiInvalidBusIndex",                     # 532
    "ApiUnsupportedFunction",                 # 533
    "ApiInvalidPciSpace",                     # 534
    "ApiInvalidIopSpace",                     # 535
    "ApiInvalidSize",                         # 536
    "ApiInvalidAddress",                      # 537
    "ApiInvalidAccessType",                   # 538
    "ApiInvalidIndex",                        # 539
    "ApiMuNotReady",                          # 540
    "ApiMuFifoEmpty",                         # 541
    "ApiMuFifoFull",                          # 542
    "ApiInvalidRegister",                     # 543
    "ApiDoorbellClearFailed",                 # 544
    "ApiInvalidUserPin",                      # 545
    "ApiInvalidUserState",                    # 546
    "ApiEepromNotPresent",                    # 547
    "ApiEepromTypeNotSupported",              # 548
    "ApiEepromBlank",                         # 549
    "ApiConfigAccessFailed",                  # 550
    "ApiInvalidDeviceInfo",                   # 551
    "ApiNoActiveDriver",                      # 552
    "ApiInsufficientResources",               # 553
    "ApiObjectAlreadyAllocated",              # 554
    "ApiAlreadyInitialized",                  # 555
    "ApiNotInitialized",                      # 556
    "ApiBadConfigRegEndianMode",              # 557
    "ApiInvalidPowerState",                   # 558
    "ApiPowerDown",                           # 559
    "ApiFlybyNotSupported",                   # 560
    "ApiNotSupportThisChannel",               # 561
    "ApiNoAction",                            # 562
    "ApiHSNotSupported",                      # 563
    "ApiVPDNotSupported",                     # 564
    "ApiVpdNotEnabled",                       # 565
    "ApiNoMoreCap",                           # 566
    "ApiInvalidOffset",                       # 567
    "ApiBadPinDirection",                     # 568
    "ApiPciTimeout",                          # 569
    "ApiDmaChannelClosed",                    # 570
    "ApiDmaChannelError",                     # 571
    "ApiInvalidHandle",                       # 572
    "ApiBufferNotReady",                      # 573
    "ApiInvalidData",                         # 574
    "ApiDoNothing",                           # 575
    "ApiDmaSglBuildFailed",                   # 576
    "ApiPMNotSupported",                      # 577
    "ApiInvalidDriverVersion",                # 578
    "ApiWaitTimeout",                         # 579
    "ApiWaitCanceled",                        # 580
    "ApiBufferTooSmall",                      # 581
    "ApiBufferOverflow",                      # 582
    "ApiInvalidBuffer",                       # 583
    "ApiInvalidRecordsPerBuffer",             # 584
    "ApiDmaPending",                          # 585
    "ApiLockAndProbePagesFailed",             # 586
    "ApiWaitAbandoned",                       # 587
    "ApiWaitFailed",                          # 588
    "ApiTransferComplete",                    # 589
    "ApiPllNotLocked",                        # 590
    "ApiNotSupportedInDualChannelMode",       # 591
    "ApiNotSupportedInQuadChannelMode",       # 592
    "ApiFileIoError",                         # 593
    "ApiInvalidClockFrequency",               # 594
    "ApiInvalidSkipTable",                    # 595
    "ApiInvalidDspModule",                    # 596
    "ApiDESOnlySupportedInSingleChannelMode", # 597
    "ApiInconsistentChannel",                 # 598
    "ApiLastError"                            # Do not add API errors below this line
]

alazar_exception(r) = exceptions[r-511]


text:
#!/usr/bin/env julia

# record start time
tstart = time()

# set root path of the ND repo
NDPATH = begin
    path = splitpath(abspath(@__DIR__))
    i = findlast(contains("NDPrototype"), path)
    abspath(path[begin:i]...)
end
BMPATH = joinpath(NDPATH, "benchmark")

# activate the benchmark folder as current environment
io = IOBuffer() # hide output
import Pkg; Pkg.activate(BMPATH; io); Pkg.develop(path=NDPATH; io);

using PkgBenchmark
using LibGit2
using Dates
using Random
using ArgParse

s = ArgParseSettings()
@add_arg_table! s begin
    "--target", "-t"
        help = "Specify branch, commit id, tag, … for target benchmark. If `directory` use the current state of the directory. If *.date file use prev. exporte raw results."
        default = "directory"
    "--baseline", "-b"
        help = "Same for baseline benchmark. If `none` or `nothing` skip the baseline benckmark and coparison step."
        default = "main"
    "--command"
        help = "Julia command for benchmarks"
        default = "julia"
    "--tcommand"
        help = "Julia command for target benchmarks. If nothing use `--command`"
    "--bcommand"
        help = "Julia command for baseline benchmarks. If nothing use `--command`"
    "--threads"
        help = "Set number of threads."
        arg_type = Int
        default = 4
    "--prefix", "-p"
        help = "Prefix for filenames, defaults to timestamp"
        default = Dates.format(now(), "yyyy-mm-dd_HHMMSS")
    "--verbose", "-v"
        help = "Print out the current benchmark."
        action = :store_true
    "--retune"
        help = "Force retuneing of parameters befor the first benchmark. (Second benchmark will use the tune file)"
        action = :store_true
    "--exportraw"
        help = "Export raw data of trials. I.e. to use benchmarks results again als baseline or target."
        action = :store_true
end
args = parse_args(s; as_symbols=true)
args[:prefix] *= '_'

@info "Run benchmarks on $(args[:target]) and compare to $(args[:baseline])"

#=
PkgBenchmark will do stuff to the repo (i.e. check out the desired branch),
therefore we copy the repo to tmp/NDPrototyp
We want to have the same benchmarks for target and base therefore we copy
NDPrototyp/benchmark to tmp/benchmark so it won't change if the repo
state changes.
=#
tmp = tempname(); mkpath(tmp)
ndpath_tmp = joinpath(tmp, "NDPrototyp")
@info "Copy repo to $ndpath_tmp"
cp(NDPATH, ndpath_tmp)
cd(ndpath_tmp) # the tmp repo is a good place to chill

# copy bm scripts to separate folder (so it won't change if the state of the nd_temp repo chagnes)
bmpath_tmp = joinpath(tmp, "benchmark")
cp(BMPATH, bmpath_tmp)
script_path = joinpath(bmpath_tmp, "benchmarks.jl")

# if directory is dirty, commit everything (otherwise we cant switch do a differen commit/branch)
isdirty = with(LibGit2.isdirty, GitRepo(ndpath_tmp))
if isdirty
    @info "Dirty directory, add everything to new commit!"
    @assert realpath(pwd()) == realpath(ndpath_tmp) "Julia is in $(pwd()) not it $ndpath_tmp"
    run(`git checkout -b $(randstring(15))`)
    run(`git add -A`)
    run(`git commit -m "tmp commit for benchmarking"`)
    # assert that repo is clean now
    with(LibGit2.isdirty, GitRepo(ndpath_tmp)) && throw(error("Repository is still dirty!"))
end

# activate the environment in the tmp copy of the scripts
# env will be used for the julia processes which do the benchmarking
Pkg.activate(bmpath_tmp)
Pkg.develop(path=ndpath_tmp)

"""
Create a runnable `Cmd` object from a string like "juila-1.5".
Due to a bug in PkgBenchmark this also adds the default sysimge
to the command. See

https://github.com/JuliaCI/PkgBenchmark.jl/issues/136
"""
function string_to_command(cmd::AbstractString)
    @assert !contains(cmd, " ") "julia commands with arguments not allowed, got $cmd"
    out = IOBuffer()
    capturecommand = `$cmd --startup-file=no -e "println(Base.julia_cmd()[3])"`
    run(pipeline(capturecommand; stdout=out))
    sysimg = String(take!(out))[1:end-1] # remove \n
    @assert contains(sysimg, r"^-J.*[dylib|so]$")  "Captured sysimg looks weird? $sysimg"
    return Cmd([cmd, sysimg])
end

"""
    benchmark(; name, id, cmd)

- name determines the filenames
- id determines the git id, i.e. branch name. If id=="directory" use currend
directory even if dirty!
- cmd if empty use default julia command
"""
function benchmark(; name, id, cmd, retune)
    # magic id string directory will benchmark dirty state
    id = id=="directory" ? nothing : id

    # PkgBenchmarks does not call Pkg.resolve after checking out a different
    # state of ndpath_tmp! Therefore we need to check out the desired commit in
    # the tmp directory and resolve the manifest of the currently activated
    # environment (-> bmpath_tmp)
    if !isnothing(id)
        @assert realpath(pwd()) == realpath(ndpath_tmp) "Julia is in $(pwd()) not it $ndpath_tmp"
        run(`git checkout $id`)
        Pkg.resolve()
    end

    # choose specific command only if given, elso the default command
    cmd = isnothing(cmd) ? args[:command] : cmd
    @info "Run $name benchmarks on $id with $cmd"

    config = BenchmarkConfig(;id,
                             env = Dict("JULIA_NUM_THREADS" => args[:threads]),
                             juliacmd = string_to_command(cmd))
    results = benchmarkpkg(ndpath_tmp, config;
                           script=script_path,
                           retune,
                           verbose=args[:verbose])

    args[:exportraw] && writeresults(joinpath(BMPATH, args[:prefix]*name*".data"), results)
    export_markdown(joinpath(BMPATH, args[:prefix]*name*".md"), results)
    return results
end

if contains(args[:target], r".data$")
    path = joinpath(BMPATH, args[:target])
    target = readresults(joinpath(BMPATH, args[:target]))
else
    target = benchmark(; name="target", id=args[:target], cmd=args[:tcommand], retune=args[:retune])
end

if args[:baseline] ∉ ["nothing", "none"]
    if contains(args[:baseline], r".data$")
        path = joinpath(BMPATH, args[:baseline])
        baseline = readresults(joinpath(BMPATH, args[:baseline]))
    elseif args[:baseline] ∉ ["nothing", "none"]
        baseline = benchmark(; name="baseline", id=args[:baseline], cmd=args[:bcommand], retune=false)
    end
    # compare the two
    judgment = judge(target, baseline)
    export_markdown(joinpath(BMPATH, args[:prefix]*"judgment.md"), judgment)
end

# copy tune file over to real repo if there is none yet or it was retuned
if !isfile(joinpath(BMPATH, "tune.json")) || args[:retune]
    println("Update tune.json in $BMPATH")
    cp(joinpath(ndpath_tmp, "benchmark", "tune.json"), joinpath(BMPATH, "tune.json"); force=args[:retune])
end

s = round(Int, time()-tstart)
m, s = s ÷ 60, s % 60
@info "Benchmarking endet after $m min $s seconds"


text:
using BVHFiles
using Test

@testset "Graph" begin
    include("graph.jl")
end

@testset "File" begin
    include("file.jl")
end

@testset "Interpolation" begin
    include("interpolation.jl")
end


text:
import XLSX

#
# examples from docstrings
#

v = XLSX.getdata("myfile.xlsx", "mysheet", "A1:B4")

f = XLSX.read("myfile.xlsx")

sheet = f["mysheet"]

sheet = XLSX.getsheet("myfile.xlsx", "mysheet")

cell = XLSX.getcell(sheet, "A1")

left = 1  # 1st column
right = 4 # 4th column
for sheetrow in XLSX.eachrow(sheet)
    for column in left:right
        cell = XLSX.getcell(sheetrow, column)

        # do something with cell
    end
end

using DataFrames, XLSX

df = DataFrame(XLSX.gettable("myfile.xlsx", "mysheet")...)

XLSX.decode_column_number("D")
XLSX.encode_column_number(4)

cn = XLSX.CellRef("AB1")
println( XLSX.row_number(cn) ) # will print 1
println( XLSX.column_number(cn) ) # will print 28
println( string(cn) ) # will print out AB1

cn = XLSX.ref"AB1"
println( XLSX.row_number(cn) ) # will print 1
println( XLSX.column_number(cn) ) # will print 28
println( string(cn) ) # will print out AB1
cr = XLSX.range"A1:C4"

XLSX.getcellrange("myfile.xlsx", "mysheet", "A1:B4")

#
# examples from README.md
#

xf = XLSX.read("myfile.xlsx")
XLSX.sheetnames(xf)
sh = xf["mysheet"]
sh["B2"] # access a cell value
sh["A2:B4"] # access a range
XLSX.getdata("myfile.xlsx", "mysheet", "A2:B4") # shorthand for all above
sh[:] # all data inside worksheet's dimension
XLSX.getdata(sh) # same as sh[:]

using DataFrames, XLSX

df = DataFrame(XLSX.gettable("myfile.xlsx", "mysheet")...)


text:
using PyCall
using FlightSims
using FSimROS
using UnPack
using DifferentialEquations

rclpy = pyimport("rclpy")
rosNode = pyimport("rclpy.node")
fsim_msg = pyimport("fsim_interfaces.msg")
std_msg = pyimport("std_msgs.msg")

struct Multicopter_ZOH_Input <: AbstractEnv
    multicopter::Multicopter
end

function FSimBase.State(env::Multicopter_ZOH_Input)
    State(env.multicopter)
end

function FSimBase.Dynamics!(env::Multicopter_ZOH_Input)
    @Loggable function dynamics!(dx, x, input, t)
        @nested_log FSimZoo._Dynamics!(env.multicopter)(dx, x, nothing, t; u=input)
    end
end

@pydef mutable struct StateNode <: rosNode.Node
    function __init__(self)
        rosNode.Node.__init__(self, "state_node")
        # publisher
        self.publisher_ = self.create_publisher(fsim_msg.PoseTwist, "state", 5)
        timer_period = 1/500
        self.i = 0
        self.t = nothing
        self.time_received = false
        function timer_callback(self)
            x = copy(self.simulator.integrator.u)
            msg = state_to_msg(self.env.multicopter, x)
            self.publisher_.publish(msg)  # publish
            self.get_logger().info("time: $(self.t), state: $(x)")
            if self.control_received && self.t != nothing
                step_until!(self.simulator, self.t)
            end
        end
        self.timer = self.create_timer(timer_period, () -> timer_callback(self))
        # env
        multicopter = LeeHexacopter()
        self.env = Multicopter_ZOH_Input(multicopter)
        state0 = State(self.env)()
        @unpack m, g = self.env.multicopter
        u0 = (m*g/6) * ones(6)  # initial control input (actually meaningless; for initialisation of simulator)
        self.simulator = Simulator(state0, Dynamics!(self.env), u0; tf=100)  # second
        step_until!(self.simulator, timer_period)  # to reduce startup latency
        reinit!(self.simulator)  # to reduce startup latency
        # subscriber
        self.control_received = false
        function listener_callback(self, msg_control)
            if !self.control_received
                self.control_received = true
            end
            input = [msg_control.u1, msg_control.u2, msg_control.u3, msg_control.u4, msg_control.u5, msg_control.u6]
            self.simulator.integrator.p = input
        end
        function listener_callback_time(self, msg)
            t = msg.data
            if !self.time_received
                self.time_received = true
                self.simulator.integrator.t = t
            end
            self.t = t
        end
        self.subscription = self.create_subscription(fsim_msg.RotorRateHexa, "control", msg -> listener_callback(self, msg), 10)
        self.subscription_time = self.create_subscription(std_msg.Float64, "time", msg -> listener_callback_time(self, msg), 10)
    end
end

function main(args=nothing)
    rclpy.init(args=args)
    println("state node started")
    state_node = StateNode()
    rclpy.spin(state_node)
    state_node.destroy_node()
    rclpy.shutdown()
end

main()


text:
# AUTO GENERATED FILE - DO NOT EDIT

export 'feffery'_antdribbon

"""
    'feffery'_antdribbon(;kwargs...)
    'feffery'_antdribbon(children::Any;kwargs...)
    'feffery'_antdribbon(children_maker::Function;kwargs...)

An AntdRibbon component.

Keyword arguments:
- `children` (a list of or a singular dash component, string or number; optional)
- `id` (String; optional)
- `className` (String; optional)
- `color` (String; optional)
- `loading_state` (optional): . loading_state has the following type: lists containing elements 'is_loading', 'prop_name', 'component_name'.
Those elements have the following types:
  - `is_loading` (Bool; optional): Determines if the component is loading or not
  - `prop_name` (String; optional): Holds which property is loading
  - `component_name` (String; optional): Holds the name of the component that is loading
- `placement` (a value equal to: 'start', 'end'; optional)
- `style` (Dict; optional)
- `text` (String; optional)
"""
function 'feffery'_antdribbon(; kwargs...)
        available_props = Symbol[:children, :id, :className, :color, :loading_state, :placement, :style, :text]
        wild_props = Symbol[]
        return Component("'feffery'_antdribbon", "AntdRibbon", "feffery_antd_components", available_props, wild_props; kwargs...)
end

'feffery'_antdribbon(children::Any; kwargs...) = 'feffery'_antdribbon(;kwargs..., children = children)
'feffery'_antdribbon(children_maker::Function; kwargs...) = 'feffery'_antdribbon(children_maker(); kwargs...)

max_stars_repo_name: ReilYoung/feffery-antd-components

max_issues_repo_name: ReilYoung/feffery-antd-components

max_forks_repo_name: ReilYoung/feffery-antd-components


text:
#Signal sources
#-------------------------------------------------------------------------------

#==
===============================================================================#

struct Sine
	A::Float64
	f::Float64
end
Sine(;A=1, f=1) = Sine(A, f)

#Last


text:
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.6.2
#     language: julia
#     name: julia-1.6
# ---

# %% [markdown]
# * https://discourse.julialang.org/t/for-monte-carlo-simulation-with-same-code-same-algorithm-how-fast-is-julia-compared-with-fortran/67021/41
# * https://discourse.julialang.org/t/for-monte-carlo-simulation-with-same-code-same-algorithm-how-fast-is-julia-compared-with-fortran/67021/43
# * https://github.com/genkuroki/public/blob/main/0018/%40defkwargs.ipynb
# * https://github.com/genkuroki/public/blob/main/0018/%40defunpack.ipynb
# * https://github.com/genkuroki/public/blob/main/0018/How%20to%20define%20unpacking%20macros.ipynb
# * https://github.com/genkuroki/public/blob/main/0018/How%20to%20define%20unpacking%20macros%20Part%202.ipynb

# %%
using ConcreteStructs

@concrete struct Foo a; b; c; d; e; f; g; h; i; j; k; l; m; n; o; p; q; r; s; t; u; v; w; x; y; z end

names = fieldnames(Foo)
kwargs = Expr(:parameters, (Expr(:kw, name, v) for (v, name) in enumerate(names))...)
@eval Foo($kwargs) = Foo($(names...))

foo = Foo(a = "meow", m = 'm', z = 99.99)

# %%
"""
`@unpackall_Foo(obj)` unpacks all fields of the object `obj` of type `Foo`.
"""
macro unpackall_Foo(obj)
    names = fieldnames(Foo)
    Expr(:(=),
        Expr(:tuple, names...),
        Expr(:tuple, (:($obj.$name) for name in names)...)
    ) |> esc
end

let
    @unpackall_Foo foo
    @show a b c d e f g h i j k l m n o p q r s t u v w x y z
end;

# %%
using Parameters

let
    @unpack a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q, r, s, t, u, v, w, x, y, z = foo
    @show a b c d e f g h i j k l m n o p q r s t u v w x y z
end;

# %%
methods(Foo)

# %%
?@unpackall_Foo

# %%
?@unpack

# %%

max_stars_repo_path: 0018/How to define unpacking macros Part 2.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_name: genkuroki/public

max_forks_repo_name: genkuroki/public


text:
img = SquareImage(rows=4, cols=4, l=2, calibration_value=1, background=0)
@testset "Testing Radon transform" begin
    @test img ≈ [
        [0 0 0 0];
        [0 1 1 0];
        [0 1 1 0];
        [0 0 0 0]
    ]
    @testset "Testing diag algorithm" begin
        sinog = radon(img, Radon())
        @test sinog ≈ [
            [0.0       0.0       0.0       0.0       0.0];
            [0.178452  0.145332  0.164798  0.164798  0.145332];
            [0.785413  0.763186  0.745526  0.745525  0.763186];
            [0.785413  0.763186  0.745526  0.745525  0.763186];
            [0.178452  0.145332  0.164798  0.164798  0.145332];
            [0.0       0.0       0.0       0.0       0.0]
        ]
    end
    @testset "Testing square algorithm" begin
        sinog = radon(img, RadonSquare())
        @test sinog ≈ [
            [0.0       0.05596    0.044909  0.064763   0.0647629  0.044909  0.05596];
            [0.888889  0.764706   0.846043  0.804828   0.804828   0.846043  0.764706];
            [0.888889  0.764706   0.846043  0.804828   0.804828   0.846043  0.764706];
            [0.0       0.0559601  0.044909  0.0647629  0.0647629  0.044909  0.05596]
        ]
    end
end


text:
module TestDuplicates
    using Base.Test
    using DataFrames

    df = DataFrame(a = [1, 2, 3, 3, 4])
    udf = DataFrame(a = [1, 2, 3, 4])
    @test isequal(nonunique(df), [false, false, false, true, false])
    @test isequal(udf, unique(df))
    unique!(df)
    @test isequal(df, udf)

    pdf = DataFrame(a = NullableCategoricalArray(Nullable{String}["a", "a", Nullable(),
                                             Nullable(), "b", Nullable(), "a", Nullable()]),
                    b = NullableCategoricalArray(Nullable{String}["a", "b", Nullable(),
                                                              Nullable(), "b", "a", "a", "a"]))
    updf = DataFrame(a = NullableCategoricalArray(Nullable{String}["a", "a", Nullable(), "b", Nullable()]),
                     b = NullableCategoricalArray(Nullable{String}["a", "b", Nullable(), "b", "a"]))
    @test isequal(nonunique(pdf), [false, false, false, true, false, false, true, true])
    @test isequal(nonunique(updf), falses(5) )
    @test isequal(updf, unique(pdf))
    unique!(pdf)
    @test isequal(pdf, updf)
end


text:
using Kinetic
using KitBase.ProgressMeter, KitBase.OffsetArrays

cd(@__DIR__)
D = read_dict("briowu_2d.txt")
for key in keys(D)
    s = Symbol(key)
    @eval $s = $(D[key])
end

begin
    γ = heat_capacity_ratio(inK, 3)
    set = Setup(matter, case, space, flux, collision, nSpecies, interpOrder, limiter, boundary, cfl, maxTime)
    pSpace = PSpace1D(x0, x1, nx, nxg)

    ue0 = umin * sqrt(mi / me)
    ue1 = umax * sqrt(mi / me)
    ve0 = vmin * sqrt(mi / me)
    ve1 = vmax * sqrt(mi / me)
    kne = knudsen * (me / mi)

    vSpace = MVSpace2D(umin, umax, ue0, ue1, nu, vmin, vmax, ve0, ve1, nv, vMeshType, nug, nvg)
    plasma = Plasma1D([knudsen,kne], mach, prandtl, inK, γ, mi, ni, me, ne, lD, rL, sol, echi, bnu)

    begin
        # upstream
        primL = zeros(5, 2)
        primL[1, 1] = 1.0 * mi
        primL[2, 1] = 0.0
        primL[3, 1] = 0.0
        primL[4, 1] = 0.0
        primL[5, 1] = mi / 1.0
        primL[1, 2] = 1.0 * me
        primL[2, 2] = 0.0
        primL[3, 2] = 0.0
        primL[4, 2] = 0.0
        primL[5, 2] = me / 1.0

        wL = mixture_prim_conserve(primL, γ)
        h0L = mixture_maxwellian(vSpace.u, vSpace.v, primL)

        h1L = similar(h0L)
        h2L = similar(h0L)
        for j in axes(h0L, 3)
            h1L[:, :, j] .= primL[4, j] .* h0L[:, :, j]
            h2L[:, :, j] .= (primL[4, j]^2 + 1.0 / (2.0 * primL[end, j])) .* h0L[:, :, j]
        end

        EL = zeros(3)
        BL = zeros(3)
        BL[1] = 0.75
        BL[2] = 1.0

        # downstream
        primR = zeros(5, 2)
        primR[1, 1] = 0.125 * mi
        primR[2, 1] = 0.0
        primR[3, 1] = 0.0
        primR[4, 1] = 0.0
        primR[5, 1] = mi * 1.25
        primR[1, 2] = 0.125 * me
        primR[2, 2] = 0.0
        primR[3, 2] = 0.0
        primR[4, 2] = 0.0
        primR[5, 2] = me * 1.25

        wR = mixture_prim_conserve(primR, γ)
        h0R = mixture_maxwellian(vSpace.u, vSpace.v, primR)

        h1R = similar(h0R)
        h2R = similar(h0R)
        for j in axes(h0R, 3)
            h1R[:, :, j] .= primR[4, j] .* h0R[:, :, j]
            h2R[:, :, j] .= (primR[4, j]^2 + 1.0 / (2.0 * primR[end, j])) .* h0R[:, :, j]
        end

        ER = zeros(3)
        BR = zeros(3)
        BR[1] = 0.75
        BR[2] = -1.0

        lorenzL = zeros(3, 2)
        lorenzR = zeros(3, 2)
        bcL = zeros(5, 2)
        bcR = zeros(5, 2)
    end

    ib = IB3F(
            wL,
            primL,
            h0L,
            h1L,
            h2L,
            bcL,
            EL,
            BL,
            lorenzL,
            wR,
            primR,
            h0R,
            h1R,
            h2R,
            bcR,
            ER,
            BR,
            lorenzR,
        )

    outputFolder = pwd()

    ks = SolverSet(set, pSpace, vSpace, plasma, ib, outputFolder)
    KS = ks

    ctr = OffsetArray{ControlVolume1D3F}(undef, axes(KS.pSpace.x, 1))
    face = Array{Interface1D3F}(undef, KS.pSpace.nx + 1)
end

begin
    idx0 = (eachindex(pSpace.x) |> collect)[1]
    idx1 = (eachindex(pSpace.x) |> collect)[end]

    for i in eachindex(ctr)
        if i <= KS.pSpace.nx ÷ 2
            ctr[i] = ControlVolume1D3F( KS.pSpace.x[i], KS.pSpace.dx[i], KS.ib.wL, KS.ib.primL,
            KS.ib.h0L, KS.ib.h1L, KS.ib.h2L, KS.ib.EL, KS.ib.BL, KS.ib.lorenzL )
        else
            ctr[i] = ControlVolume1D3F( KS.pSpace.x[i], KS.pSpace.dx[i], KS.ib.wR, KS.ib.primR,
            KS.ib.h0R, KS.ib.h1R, KS.ib.h2R, KS.ib.ER, KS.ib.BR, KS.ib.lorenzR )
        end
    end

    face = Array{Interface1D3F}(undef, KS.pSpace.nx+1)
    for i=1:KS.pSpace.nx+1
        face[i] = Interface1D3F(KS.ib.wL, KS.ib.h0L, KS.ib.EL)
    end
end

begin
    simTime = 0.
    dt = Kinetic.timestep(KS, ctr, simTime)
    nt = Int(floor(ks.set.maxTime / dt))+1
    res = zeros(5, 2)
end

@showprogress for iter in 1:nt
    #dt = Kinetic.timestep(KS, ctr, simTime)
    #Kinetic.reconstruct!(ks, ctr)

    Kinetic.evolve!(ks, ctr, face, dt; mode=:kcu, isPlasma=:true)
    # it's equivalent to the following process
    #=
    @inbounds Threads.@threads for i = 1:KS.pSpace.nx+1
        flux_kfvs!(
            face[i].fw,
            face[i].fh0,
            face[i].fh1,
            face[i].fh2,
            ctr[i-1].h0,
            ctr[i-1].h1,
            ctr[i-1].h2,
            ctr[i].h0,
            ctr[i].h1,
            ctr[i].h2,
            KS.vSpace.u,
            KS.vSpace.v,
            KS.vSpace.weights,
            dt,
            1.0,
        )

        flux_em!(
            face[i].femL,
            face[i].femR,
            ctr[i-2].E,
            ctr[i-2].B,
            ctr[i-1].E,
            ctr[i-1].B,
            ctr[i].E,
            ctr[i].B,
            ctr[i+1].E,
            ctr[i+1].B,
            ctr[i-1].ϕ,
            ctr[i].ϕ,
            ctr[i-1].ψ,
            ctr[i].ψ,
            ctr[i-1].dx,
            ctr[i].dx,
            KS.gas.A1p,
            KS.gas.A1n,
            KS.gas.D1,
            KS.gas.sol,
            KS.gas.χ,
            KS.gas.ν,
            dt,
        )
    end
    =#

    Kinetic.update!(ks, ctr, face, dt, res, isMHD=true)
end

soluiton = zeros(ks.pSpace.nx, 10, 2)
for i in 1:ks.pSpace.nx
    soluiton[i, 1, 1] = ctr[i].prim[1,1]
    soluiton[i, 1, 2] = ctr[i].prim[1,2] / ks.gas.me
    soluiton[i, 2:4, 1] .= ctr[i].prim[2:4,1]
    soluiton[i, 2:4, 2] .= ctr[i].prim[2:4,2]
    soluiton[i, 5, 1] = 1. / ctr[i].prim[5,1]
    soluiton[i, 5, 2] = ks.gas.me / ctr[i].prim[5,2]

    soluiton[i, 6, 1] = ctr[i].B[2]
    soluiton[i, 6, 2] = ctr[i].E[1]
end
using Plots
plot(ks.pSpace.x[1:ks.pSpace.nx], soluiton[:,1,1:2])
plot(ks.pSpace.x[1:ks.pSpace.nx], soluiton[:,6,1:2])

using JLD2
@save "refsol.jld2" ks ctr


text:
#===================================================================================================
  Generic kernel function
===================================================================================================#

call{T}(κ::Kernel{T}, x::T, y::T) = kernel(κ, x, y)
call{T}(κ::Kernel{T}, x::AbstractArray{T}, y::AbstractArray{T}) = kernel(κ, x, y)

kernel{T}(κ::StandardKernel{T}, x::T, y::T) = pairwise(κ, x, y)
kernel{T}(κ::StandardKernel{T}, x::AbstractArray{T}, y::AbstractArray{T}) = pairwise(κ, x, y)

kernel{T}(κ::KernelComposition{T}, x::T, y::T) = phi(κ.phi, pairwise(κ.kappa, x, y))
function kernel{T}(κ::KernelComposition{T}, x::AbstractArray{T}, y::AbstractArray{T})
    phi(κ.phi, pairwise(κ.kappa, x, y))
end

#==========================================================================
  Generic kernelmatrix Functions
==========================================================================#

function kernelmatrix{T}(
        σ::Union{Type{Val{:row}},Type{Val{:col}}},
        κ::Kernel{T},
        X::AbstractMatrix{T},
        symmetrize::Bool = true
    )
    kernelmatrix!(σ, init_pairwisematrix(σ, X), κ, X, symmetrize)
end

function kernelmatrix{T}(
        κ::Kernel{T},
        X::AbstractMatrix{T},
        symmetrize::Bool = true
    )
    kernelmatrix(Val{:row}, κ, X, symmetrize)
end

function kernelmatrix{T}(
        σ::Union{Type{Val{:row}},Type{Val{:col}}},
        κ::Kernel{T},
        X::AbstractMatrix{T},
        Y::AbstractMatrix{T}
    )
    kernelmatrix!(σ, init_pairwisematrix(σ, X, Y), κ, X, Y)
end

function kernelmatrix{T}(
        κ::Kernel{T},
        X::AbstractMatrix{T},
        Y::AbstractMatrix{T}
    )
    kernelmatrix(Val{:row}, κ, X, Y)
end

#== Standard Kernel ==#

function kernelmatrix!{T}(
        σ::Union{Type{Val{:row}},Type{Val{:col}}},
        K::Matrix{T},
        κ::StandardKernel{T},
        X::AbstractMatrix{T},
        symmetrize::Bool = true
    )
    pairwisematrix!(σ, K, κ, X, symmetrize)
end

function kernelmatrix!{T}(
        σ::Union{Type{Val{:row}},Type{Val{:col}}},
        K::Matrix{T},
        κ::StandardKernel{T},
        X::AbstractMatrix{T},
        Y::AbstractMatrix{T}
    )
    pairwisematrix!(σ, K, κ, X, Y)
end

#== Composition Kernel ==#

function phi_rectangular!{T}(ϕ::CompositionClass{T}, K::AbstractMatrix{T})
    for i in eachindex(K)
        @inbounds K[i] = phi(ϕ, K[i])
    end
    K
end

function phi_symmetric!{T}(ϕ::CompositionClass{T}, K::AbstractMatrix{T}, symmetrize::Bool = true)
    if !((n = size(K,1)) == size(K,2))
        throw(DimensionMismatch("Kernel matrix must be square."))
    end
    for j = 1:n, i = (1:j)
        @inbounds K[i,j] = phi(ϕ, K[i,j])
    end
    symmetrize ? LinAlg.copytri!(K, 'U') : K
end

function kernelmatrix!{T}(
        σ::Union{Type{Val{:row}},Type{Val{:col}}},
        K::Matrix{T},
        κ::KernelComposition{T},
        X::AbstractMatrix{T},
        symmetrize::Bool = true
    )
    pairwisematrix!(σ, K, κ.kappa, X, false)
    phi_symmetric!(κ.phi, K, symmetrize)
end

function kernelmatrix!{T}(
        σ::Union{Type{Val{:row}},Type{Val{:col}}},
        K::Matrix{T},
        κ::KernelComposition{T},
        X::AbstractMatrix{T},
        Y::AbstractMatrix{T}
    )
    pairwisematrix!(σ, K, κ.kappa, X, Y)
    phi_rectangular!(κ.phi, K)
end

#== Kernel Operation ==#

for (kernel_object, scalar_op, identity, scalar) in (
        (:KernelProduct, :*, :1, :a),
        (:KernelSum,     :+, :0, :c)
    )
    @eval begin
        function kernelmatrix!{T}(
                σ::Union{Type{Val{:row}},Type{Val{:col}}},
                K::Matrix{T},
                κ::$kernel_object{T},
                X::AbstractMatrix{T}
            )
            kernelmatrix!(σ, K, κ.kappa1, X)
            broadcast!($scalar_op, K, kernelmatrix(σ, similar(K), κ.kappa2, X))
            κ.$scalar == $identity ? K : broadcast!($scalar_op, K, κ.$scalar)
        end

        function kernelmatrix!{T}(
                σ::Union{Type{Val{:row}},Type{Val{:col}}},
                K::Matrix{T},
                κ::$kernel_object{T},
                X::AbstractMatrix{T},
                Y::AbstractMatrix{T}
            )
            kernelmatrix!(σ, K, κ.kappa1, X, Y)
            broadcast!($scalar_op, K, kernelmatrix(σ, similar(K), κ.kappa2, X, Y))
            κ.$scalar == $identity ? K : broadcast!($scalar_op, K, κ.$scalar)
        end
    end
end

#==========================================================================
  Kernel Matrix Transformation
==========================================================================

# Centralize a kernel matrix K
function centerkernelmatrix!{T<:AbstractFloat}(K::Matrix{T})
    (n = size(K, 1)) == size(K, 2) || throw(DimensionMismatch("Kernel matrix must be square"))
    μ_row = zeros(T,n)
    μ = zero(T)
    @inbounds for j = 1:n
        @simd for i = 1:n
            μ_row[j] += K[i,j]
        end
        μ += μ_row[j]
        μ_row[j] /= n
    end
    μ /= n^2
    @inbounds for j = 1:n
        @simd for i = 1:n
            K[i,j] += μ - μ_row[i] - μ_row[j]
        end
    end
    K
end
centerkernelmatrix{T<:AbstractFloat}(K::Matrix{T}) = centerkernelmatrix!(copy(K))

==========================================================================
  Kernel Operation Matrix Functions
==========================================================================

for (kernel_object, matrix_op, array_op, identity, scalar) in (
        (:KernelProduct, :matrix_prod!, :scale!,     :1, :a),
        (:KernelSum,     :matrix_sum!,  :translate!, :0, :c)
    )
    @eval begin

        function kernelmatrix!{T<:AbstractFloat}(K::Matrix{T}, κ::$kernel_object{T}, X::Matrix{T},
                                                 is_trans::Bool, store_upper::Bool,
                                                 symmetrize::Bool)
            nk = length(κ.k)
            kernelmatrix!(K, κ.k[1], X, is_trans, store_upper, false)
            if nk > 1
                K_factor = similar(K)
                for i = 2:nk
                    ($matrix_op)(K, kernelmatrix!(K_factor, κ.k[i], X, is_trans, store_upper, false), store_upper, false)
                end
            end
            K = symmetrize ? (store_upper ? syml!(K) : symu!(K)) : K
            κ.$scalar == $identity ? K : ($array_op)(κ.$scalar, K)
        end

        function kernelmatrix!{T<:AbstractFloat}(K::Matrix{T}, κ::$kernel_object{T}, X::Matrix{T},
                                                 Y::Matrix{T}, is_trans::Bool)
            nk = length(κ.k)
            kernelmatrix!(K, κ.k[1], X, Y, is_trans)
            if nk > 1
                K_factor = similar(K)
                for i = 2:nk
                    ($matrix_op)(K, kernelmatrix!(K_factor, κ.k[i], X, Y, is_trans))
                end
            end
            κ.$scalar == $identity ? K : ($array_op)(κ.$scalar, K)
        end

    end
end

=#


text:
using Pkg

if !in("Distributions", keys(Pkg.installed()))
	Pkg.add("Distributions")
end

if !in("Convex", keys(Pkg.installed()))
	Pkg.add("Convex")
end

if !in("ECOS", keys(Pkg.installed()))
	Pkg.add("ECOS")
end

if !in("DataFrames", keys(Pkg.installed()))
	Pkg.add("DataFrames")
end

if !in("CSV", keys(Pkg.installed()))
	Pkg.add("CSV")
end

max_stars_repo_name: Airspace-Encounter-Models/em-model-uam

max_issues_repo_name: Airspace-Encounter-Models/em-model-uam

max_forks_repo_name: Airspace-Encounter-Models/em-model-uam


text:
# GR.inline("pdf")
# upscale = 2 #upscaling in resolution
# fntsm = font("serif", pointsize=round(8.0*upscale))
# fntlg = font("serif", pointsize=round(12.0*upscale))
# default(titlefont=fntlg, guidefont=fntlg, tickfont=fntsm, legendfont=fntsm)
# default(size=(800*upscale,600*upscale)) #Plot canvas size
#
# gr()

function plot_results_dict_line(
      results_dict_allplots,
      subnet_arr,
      x_axis,
      output_folder,
      vert_line,
      horiz_line,
      xlimits,
      ylimits,
      output_plot_label,
      xvals=[],
      points_dict=Dict();
      plot_infeasible_boundaries=true,
      series_labels=[],
      color_palette=:tab10
   )

   # GR.inline("pdf")

   upscale = 2 #upscaling in resolution
   fntsm = font("serif", pointsize=round(8.0*upscale))
   fntlg = font("serif", pointsize=round(12.0*upscale))
   default(titlefont=fntlg, guidefont=fntlg, tickfont=fntsm, legendfont=fntsm)
   default(size=(800*upscale,600*upscale)) #Plot canvas size

   # gr()

   # Given each specific key for the x axis values, collect the values into x_arr and the indices into indices array
   # println("axis: $axis")

   if isa(results_dict_allplots, Dict)
      results_dict_allplots = [results_dict_allplots]
   end
   if length(series_labels) == 0
      series_labels = [[] for i in results_dict_allplots]
   end

   function make_x_array(results_dict)
      ax_key = x_axis

      if !(ax_key[1] in keys(results_dict))
         if ax_key[1] == "sn"
            x_label = join(ax_key, " ")
            # x_label = ax_key
            # println("x_label: $x_label")
            x_arr = Array{Float64,1}()
            indices_arr = Array{Int64,1}()
            for (param,res) in results_dict[x_label]
               push!(x_arr, res)
               push!(indices_arr, param)
            end
            # Sort values by index of results dictionary (index is the parameter which was varied)
            perm = sortperm(indices_arr)
            x_sorted = x_arr[perm]
            indices_sorted = indices_arr[perm]
         else # custom x axis
            x_label = ax_key[1]
            x_arr = ax_key[2]
            indices_arr = Array{Int64,1}()
            for (param,res) in results_dict["subnet"]
               push!(indices_arr, param)
            end
            # Sort values by index of results dictionary (index is the parameter which was varied)
            perm = sortperm(indices_arr)
            indices_sorted = indices_arr[perm]
            x_sorted = x_arr[indices_sorted]
         end
      else
         (res_label_x,subnet_x) = ax_key

         x_arr = Array{Float64,1}()
         x_label = ""
         indices_arr = Array{Int64,1}()
         for (param,res) in results_dict[res_label_x]
            if length(res) == 1
               x_result = (res isa Array) ? res[1] : res
               x_label = res_label_x
            elseif subnet_x == 0
               x_result = sum(res)
               x_label = "total $res_label_x"
            else
               x_result = res[findfirst(subnet_arr.==subnet_x)]
               x_label = "subnet $subnet_x $res_label_x"
            end
            push!(x_arr, x_result)
            param_num = (typeof(param) == String) ? parse(Int64, param) : param
            push!(indices_arr, param_num)
         end
         # Sort values by index of results dictionary (index is the parameter which was varied)
         perm = sortperm(indices_arr)
         x_sorted = x_arr[perm]
         indices_sorted = indices_arr[perm]
      end

      # Now do the same for the points dictionary
      pts_x_arr = Array{Float64,1}()
      if length(points_dict) > 0
         ax_key = x_axis
         if !(ax_key[1] in keys(points_dict))
            if ax_key[1] == "sn"
               x_label = join(ax_key, " ")
               # x_label = ax_key
               # println("x_label: $x_label")
               pts_indices_arr = Array{Int64,1}()
               for (param,res) in points_dict[x_label]
                  push!(pts_x_arr, res)
                  push!(pts_indices_arr, param)
               end
            else # custom x axis
               x_label = ax_key[1]
               pts_x_arr = ax_key[2]
               pts_indices_arr = Array{Int64,1}()
               if "subnet" in keys(points_dict)
                  for (param,res) in points_dict["subnet"]
                     push!(pts_indices_arr, param)
                  end
               end
            end
         else
            (res_label_x,subnet_x) = ax_key
            pts_x_arr = Array{Float64,1}()
            x_label = ""
            pts_indices_arr = Array{Int64,1}()
            for (param,res) in points_dict[res_label_x]
               if length(res) == 1
                  x_result = (res isa Array) ? res[1] : res
                  x_label = res_label_x
               elseif subnet_x == 0
                  x_result = sum(res)
                  x_label = "total $res_label_x"
               else
                  x_result = res[findfirst(subnet_arr.==subnet_x)]
                  x_label = "subnet $subnet_x $res_label_x"
               end
               push!(pts_x_arr, x_result)
               push!(pts_indices_arr, param)
            end
         end
      end
      return x_sorted, indices_sorted, pts_x_arr, perm, x_label
   end

   #    if (label == "status") || (label == "subnet")
   #       continue
   #    end
   #    for subnet_i in 0:n_subnets
   #       # println("subnet_i: $subnet_i")
   #       br_results[1] = Array{Float64,1}()
   #       br_indices[1] = Array{String,1}()
   #       base_res = -9999.0
   #       if subnet_i == 0
   #          # Don't plot total for frequency
   #          if label == "frequency (Hz)"
   #             continue
   #          end
   #          plotlabel = "total $label"
   #          plotdir = "total"
   #       else
   #          plotlabel = "subnet $(subnet_array[subnet_i]) $label"
   #          plotdir =  "subnet$(subnet_array[subnet_i])"
   #       end
   #       for (folder,res) in results
   #          if results_dict["status"][folder] in [LOCALLY_SOLVED, ALMOST_LOCALLY_SOLVED]
   #             if subnet_i == 0
   #                subnet_result = sum(res)
   #             else
   #                try
   #                   subnet_result = res[subnet_i]
   #                catch
   #                   subnet_result = 0
   #                end
   #             end
   #          else
   #             subnet_result = -9999.0
   #             println("Subnet $subnet_i, network $folder: result is -9999.0")
   #          end
   #          if occursin("br", folder)
   #             push!(br_results[1], subnet_result)
   #             push!(br_indices[1], folder[3:end])
   #          # Remove base value to plot separately
   #          elseif occursin("base", folder)
   #             base_res = subnet_result
   #          end
   #       end

   # Collect all values in the dictionary and plot them against the x values
   # Some dictionary entries contain values from multiple subnetworks, while others are specific to a single bus, gen, or branch
   function plot_series!(plots, results_dict, x_sorted, indices_sorted, pts_x_arr, perm, x_label, series_label)
      output_plot = nothing

      for (label,results) in results_dict
         res_vals = collect(values(results))
         if (label == "status") || (label == "subnet") || (length(res_vals) < 1) || (any(res_vals.==nothing))
            continue
         end
         if length(points_dict) > 0
            pts_results = points_dict[label]
         end
         pts_y_arr = Array{Float64,1}()
         pts_labels = Array{Any,1}()
         # println("results: $results")
         n_subnets = length(res_vals[1])
         y_arr = Array{Float64,1}()
         status_pairs = sort(collect(results_dict["status"]), by = x -> x[1])
         statuses = map(x -> x[2], status_pairs)
         solved = indexin(statuses, [LOCALLY_SOLVED, ALMOST_LOCALLY_SOLVED, "LOCALLY_SOLVED", "ALMOST_LOCALLY_SOLVED"]) .!= nothing
         if plot_infeasible_boundaries
            # println("solved: $solved")
            infeas_boundary_l = (solved[2:end] - solved[1:end-1]) .< 0
            append!(infeas_boundary_l, [0])
            infeas_boundary_r = zeros(Bool, 1)
            append!(infeas_boundary_r, (solved[2:end] - solved[1:end-1]) .> 0)
            # println("infeas_boundary_l: $(infeas_boundary_l)")
            # println("infeas_boundary_r: $(infeas_boundary_r)")
            infeas_boundary_idx = x_sorted[infeas_boundary_l .| infeas_boundary_r]
         else
            infeas_boundary_idx = []
         end
         # println("infeas_boundary_idx: $infeas_boundary_idx")
         if n_subnets == 1
            plt = plots[label][1]
            plotlabel = label
            dirs = split(label," ")
            if dirs[1] == "sn"
               plotdir = join(dirs[1:end-1],"/")
            else
               plotdir = ""
            end
            for (param,res) in results
               y_result = (res isa Array) ? res[1] : res
               push!(y_arr, y_result)
            end
            # Sort values by index of results dictionary (index is the parameter which was varied)
            y_sorted = y_arr[perm]
            y_sorted[.!solved] .= NaN
            # println("y_sorted: $y_sorted")
            if length(points_dict) > 0
               for (param,res) in pts_results
                  pts_y_result = (res isa Array) ? res[1] : res
                  push!(pts_y_arr, pts_y_result)
                  pts_y_label = (param isa Array) ? param[1] : param
                  push!(pts_labels, string(pts_y_label))
               end
            end
            # println("$label: single subnet: $y_sorted")
            plotdir = "vs "*x_label*"/"*plotdir
            p = plot_line!(plt,x_sorted, y_sorted, x_label, plotlabel, output_folder, plotdir, vert_line, horiz_line, xlimits, ylimits, xvals, infeas_boundary_idx, pts_x_arr=pts_x_arr, pts_y_arr=pts_y_arr, pts_labels=pts_labels, series_label=series_label)
            if (x_label == output_plot_label[1]) && (plotlabel == output_plot_label[2])
               output_plot = p
               # println("= = = = = = = =\noutput_plot: $output_plot\n= = = = = = = =")
            end
         else
            for subnet_i in 0:n_subnets
               # println("subnet_i: $subnet_i")
               plt = plots[label][subnet_i]
               y_arr = Array{Float64,1}()
               pts_y_arr = Array{Float64,1}()
               pts_labels = Array{String,1}()
               # indices = Array{Int64,1}()
               if subnet_i == 0
                  # Don't plot total for frequency
                  if label == "frequency (Hz)"
                     continue
                  end
                  plotlabel = "total $label"
                  plotdir = "total"
               else
                  plotlabel = "subnet $(subnet_arr[subnet_i]) $label"
                  plotdir =  "subnet$(subnet_arr[subnet_i])"
               end
               for (folder,res) in results
                  if subnet_i == 0
                     subnet_result = sum(res)
                  else
                     subnet_result = res[subnet_i]
                  end
                  push!(y_arr, subnet_result)
               end
               # Sort values by index of results dictionary (index is the parameter which was varied)
               y_sorted = y_arr[perm]
               y_sorted[.!solved] .= NaN
               # println("y_sorted: $y_sorted")
               if length(points_dict) > 0
                  for (param,res) in pts_results
                     if subnet_i == 0
                        subnet_result = sum(res)
                     else
                        subnet_result = res[subnet_i]
                     end
                     push!(pts_y_arr, subnet_result)
                     push!(pts_labels, param)
                  end
               end

               plotdir = "vs "*x_label*"/"*plotdir
               p = plot_line!(plt, x_sorted, y_sorted, x_label, plotlabel, output_folder, plotdir, vert_line, horiz_line, xlimits, ylimits, xvals, infeas_boundary_idx, pts_x_arr=pts_x_arr, pts_y_arr=pts_y_arr, pts_labels=pts_labels, series_label=series_label)
               if (x_label == output_plot_label[1]) && (plotlabel == output_plot_label[2])
                  output_plot = p
                  # println("= = = = = = = =\noutput_plot: $output_plot\n= = = = = = = =")
               end
            end
         end
      end
      return output_plot
   end

   x_sorted_allplots = []
   indices_sorted_allplots = []
   pts_x_arr_allplots = []
   perm_allplots = []
   x_label_allplots = []
   for results_dict in results_dict_allplots
      (x_sorted, indices_sorted, pts_x_arr, perm, x_label) = make_x_array(results_dict)
      push!(x_sorted_allplots,x_sorted)
      push!(indices_sorted_allplots,indices_sorted)
      push!(pts_x_arr_allplots,pts_x_arr)
      push!(perm_allplots,perm)
      push!(x_label_allplots,x_label)
   end

   plots_dict = Dict{String,Any}()
   for (dict_idx,results_dict) in enumerate(results_dict_allplots)
      for (label,results) in results_dict
         # println("label: $label")
         if (label in keys(plots_dict))
            continue
         end
         if (label == "status") || (label == "subnet")
            continue
         end
         if any(values(results) .== nothing)
            continue
         end
         plots_dict[label] = Dict{Int64,Any}()
         n_subnets = length(collect(values(results))[1])
         for subnet_i in 0:n_subnets
            if !(subnet_i in keys(plots_dict[label]))
               plots_dict[label][subnet_i] = Plots.plot(palette=color_palette)
            end
         end
      end
   end

   output_plot = Plots.Plot()
   for (i,results_dict) in enumerate(results_dict_allplots)
      output_plot = plot_series!(
      plots_dict,
      results_dict,
      x_sorted_allplots[i],
      indices_sorted_allplots[i],
      pts_x_arr_allplots[i],
      perm_allplots[i],
      x_label_allplots[i],
      series_labels[i]
      )
   end

   # outprint = output_plot != nothing ? output_plot : "nothing"
   # println("\n***IN plot_results_dict_line, output_plot = $outprint")
   return output_plot
end

function plot_line(
      x_sorted,
      y_sorted,
      x_label,
      y_label,
      output_folder,
      plotdir,
      vert_line=([],""),
      horiz_line=([],""),
      xlimits=[],
      ylimits=[],
      xvals=[],
      infeasible_boundaries=[];
      pts_x_arr,
      pts_y_arr,
      pts_labels
   )
   # GR.inline("pdf")

   upscale = 2 #upscaling in resolution
   fntsm = font("serif", pointsize=round(8.0*upscale))
   fntlg = font("serif", pointsize=round(12.0*upscale))
   default(titlefont=fntlg, guidefont=fntlg, tickfont=fntsm, legendfont=fntsm)
   default(size=(800*upscale,600*upscale)) #Plot canvas size

   # gr()
   plt = Plots.plot()
   plot_line!(plt,x_sorted, y_sorted, x_label, y_label, output_folder, plotdir, vert_line, horiz_line, xlimits, ylimits, xvals, infeasible_boundaries; pts_x_arr=pts_x_arr, pts_y_arr=pts_y_arr, pts_labels=pts_labels)

   return plt
end

function plot_line!(
      plt::Plots.Plot,
      x_sorted,
      y_sorted,
      x_label,
      y_label,
      output_folder,
      plotdir,
      vert_line=([],""),
      horiz_line=([],""),
      xlimits=[],
      ylimits=[],
      xvals=[],
      infeasible_boundaries=[];
      pts_x_arr,
      pts_y_arr,
      pts_labels,
      series_label=""
   )
   # GR.inline("pdf")

   upscale = 2 #upscaling in resolution
   fntsm = font("serif", pointsize=round(8.0*upscale))
   fntlg = font("serif", pointsize=round(12.0*upscale))
   default(titlefont=fntlg, guidefont=fntlg, tickfont=fntsm, legendfont=fntsm)
   default(size=(800*upscale,600*upscale)) #Plot canvas size

   # gr()
   # println("x_sorted: $x_sorted")
   # println("y_sorted: $y_sorted")
   if length(series_label) == 0
      plot!(plt, x_sorted, y_sorted,
         legend=false,
         left_margin=(10*upscale)*mm,
         right_margin=(10*upscale)*mm, top_margin=(10*upscale)*mm, bottom_margin=(10*upscale)*mm,
         linewidth=2*upscale,
         # marker = (:dot, 1.5*upscale, 0.6, :green, stroke(0))
      )
   else
      plot!(plt, x_sorted, y_sorted,
         label=series_label,
         legend=true,
         left_margin=(10*upscale)*mm,
         right_margin=(10*upscale)*mm, top_margin=(10*upscale)*mm, bottom_margin=(10*upscale)*mm,
         linewidth=2*upscale,
         # marker = (:dot, 1.5*upscale, 0.6, :green, stroke(0))
      )
   end

   if length(xvals) > 0
      # println("xvals: $xvals")
      # println("x_sorted: $x_sorted")
      # println("y_sorted: $y_sorted")
      plot!(
         plt, xvals, y_sorted[indexin(xvals, x_sorted)],
         line=:scatter,
         marker=(:dot, 2*upscale, 0.8, :red, stroke(0))
         )
      # println("plotted $xvals, $(y_sorted[indexin(xvals, x_sorted)])")
   end
   if length(infeasible_boundaries) > 0
      for v in infeasible_boundaries
         # println("v: $v")
         vline!([v], line=(0.5*upscale, :dash, :red), label="")
      end
   end
   if (length(pts_x_arr) > 0) && (length(pts_y_arr) > 0)
      plot!(
         plt, pts_x_arr, pts_y_arr,
         # series_annotations=Plots.series_annotations(string.(pts_labels), Plots.font("Sans", 12)),
         # annotation=string.(pts_labels),
         annotations = [(pts_x_arr[i], pts_y_arr[i], text(pts_labels[i], font("Sans", 12), :left, :bottom)) for i in 1:length(pts_y_arr)],
         # annotationguide = :left,
         line=:scatter,
         marker=(:dot, 2*upscale, 0.8, :blue, stroke(0))
         )
      # annotate!(p, pts_x_arr, pts_y_arr, text.(pts_labels, :left), font("Sans", 12))
   end
   if length(vert_line[1]) > 0
      for v in vert_line[1]
         # println("v: $v")
         vline!(plt, [v], line=(2*upscale, :dash), label=vert_line[2])
      end
   end
   if length(horiz_line[1]) > 0
      for h in horiz_line[1]
         # println("h: $h")
         hline!([h], line=(2*upscale, :dash), label=horiz_line[2])
         miny = min(min(y_sorted...), min(horiz_line[1]...), min(pts_y_arr...,Inf))
         maxy = max(max(y_sorted...), max(horiz_line[1]...), max(pts_y_arr...,-Inf))
      end
   else
      miny = min(min(y_sorted...), min(pts_y_arr...,Inf))
      maxy = max(max(y_sorted...), max(pts_y_arr...,-Inf))
   end

   # Set x limits
   if length(xlimits) == 2
      xlims!(plt, min(xlimits...), max(xlimits...))
   end

   # Set y limits
   if length(ylimits) == 2
      ylims!(plt, min(ylimits...), max(ylimits...))
   else
      curr_ylims = ylims(plt)
      # println("Current y limits: $curr_ylims")
      diff = maxy - miny
      if diff > 0
         ylims!(plt, min(miny-0.1*diff, curr_ylims[1]), max(maxy+0.1*diff,curr_ylims[2]))
      else
         ylims!(plt, min(0.9*miny, curr_ylims[1]), max(1.1*maxy,curr_ylims[2]))
      end
   end

   # println("New y limits: $(ylims(plt))")
   xlabel!(plt, x_label)
   ylabel!(plt, y_label)
   # savefig(p,"$(folder)/$(y_label) vs $(x_label).pdf")
   if !isdir("$(output_folder)/plots/$plotdir")
      mkpath("$(output_folder)/plots/$plotdir")
   end
   # println("saving plot at $(output_folder)/plots/$plotdir/$(y_label) vs $(x_label).pdf")
   # display(plt)
   savefig(plt,"$(output_folder)/plots/$plotdir/$(y_label) vs $(x_label).pdf")
end


text:
"""
    getbranch(p, q)

Get the branch that connects node `p` to node `q`.
"""
function getbranch(p::AbstractNode, q::AbstractNode)::Branch
    for link ∈ p.links
        if link.to ≡ q
            return link.branch
        end
    end

    p ≡ q && throw(InvalidTopology("`p` and `q` are the same node"))
    throw(InvalidTopology("`p` is not linked to `q`"))
end # function getbranch

"""
    neighbours(p)

Get all the nodes connected to node `p`.
"""
neighbours(p::AbstractNode) = map(x -> x.to, p.links)

"""
    n_neighbour(p)

Get the number of neighbours of node `p`.
"""
n_neighbour(p::AbstractNode) = length(p.links)

"""
    istip(p::AbstractNode)

Return true if node `p` is a tip, i.e. it has exactly one neighbour.
"""
istip(p::AbstractNode) = length(p.links) == 1

"""
    areneighbours(p, q)

Return true if node `p` is a neighbour of node `q`.
"""
function areneighbours(p::AbstractNode, q::AbstractNode)
    p2q, q2p = false, false
    for link ∈ p.links
        if link.to ≡ q
            p2q = true
            break
        end
    end
    for link ∈ q.links
        if link.to ≡ p
            q2p = true
            break
        end
    end
    if p2q ⊻ q2p
        p2q && throw(InvalidTopology("the link from `q` to `p` is missing."))
        q2p && throw(InvalidTopology("the link from `p` to `q` is missing."))
    end

    return p2q
end # function areneighbours


text:
package ln

import (
	"math"
	"sort"
)

struct Tree
    Box::Box
end

type Tree struct {
	Box  Box
	Root *Node
}

func NewTree(shapes []Shape) *Tree {
	box := BoxForShapes(shapes)
	node := NewNode(shapes)
	node.Split(0)
	return &Tree{box, node}
}

func (tree *Tree) Intersect(r Ray) Hit {
	tmin, tmax := tree.Box.Intersect(r)
	if tmax < tmin || tmax <= 0 {
		return NoHit
	}
	return tree.Root.Intersect(r, tmin, tmax)
}

type Node struct {
	Axis   Axis
	Point  float64
	Shapes []Shape
	Left   *Node
	Right  *Node
}

func NewNode(shapes []Shape) *Node {
	return &Node{AxisNone, 0, shapes, nil, nil}
}

func (node *Node) Intersect(r Ray, tmin, tmax float64) Hit {
	var tsplit float64
	var leftFirst bool
	switch node.Axis {
	case AxisNone:
		return node.IntersectShapes(r)
	case AxisX:
		tsplit = (node.Point - r.Origin.X) / r.Direction.X
		leftFirst = (r.Origin.X < node.Point) || (r.Origin.X == node.Point && r.Direction.X <= 0)
	case AxisY:
		tsplit = (node.Point - r.Origin.Y) / r.Direction.Y
		leftFirst = (r.Origin.Y < node.Point) || (r.Origin.Y == node.Point && r.Direction.Y <= 0)
	case AxisZ:
		tsplit = (node.Point - r.Origin.Z) / r.Direction.Z
		leftFirst = (r.Origin.Z < node.Point) || (r.Origin.Z == node.Point && r.Direction.Z <= 0)
	}
	var first, second *Node
	if leftFirst {
		first = node.Left
		second = node.Right
	} else {
		first = node.Right
		second = node.Left
	}
	if tsplit > tmax || tsplit <= 0 {
		return first.Intersect(r, tmin, tmax)
	} else if tsplit < tmin {
		return second.Intersect(r, tmin, tmax)
	} else {
		h1 := first.Intersect(r, tmin, tsplit)
		if h1.T <= tsplit {
			return h1
		}
		h2 := second.Intersect(r, tsplit, math.Min(tmax, h1.T))
		if h1.T <= h2.T {
			return h1
		} else {
			return h2
		}
	}
}

func (node *Node) IntersectShapes(r Ray) Hit {
	hit := NoHit
	for _, shape := range node.Shapes {
		h := shape.Intersect(r)
		if h.T < hit.T {
			hit = h
		}
	}
	return hit
}

func (node *Node) PartitionScore(axis Axis, point float64) int {
	left, right := 0, 0
	for _, shape := range node.Shapes {
		box := shape.BoundingBox()
		l, r := box.Partition(axis, point)
		if l {
			left++
		}
		if r {
			right++
		}
	}
	if left >= right {
		return left
	} else {
		return right
	}
}

func (node *Node) Partition(size int, axis Axis, point float64) (left, right []Shape) {
	left = make([]Shape, 0, size)
	right = make([]Shape, 0, size)
	for _, shape := range node.Shapes {
		box := shape.BoundingBox()
		l, r := box.Partition(axis, point)
		if l {
			left = append(left, shape)
		}
		if r {
			right = append(right, shape)
		}
	}
	return
}

func (node *Node) Split(depth int) {
	if len(node.Shapes) < 8 {
		return
	}
	xs := make([]float64, 0, len(node.Shapes)*2)
	ys := make([]float64, 0, len(node.Shapes)*2)
	zs := make([]float64, 0, len(node.Shapes)*2)
	for _, shape := range node.Shapes {
		box := shape.BoundingBox()
		xs = append(xs, box.Min.X)
		xs = append(xs, box.Max.X)
		ys = append(ys, box.Min.Y)
		ys = append(ys, box.Max.Y)
		zs = append(zs, box.Min.Z)
		zs = append(zs, box.Max.Z)
	}
	sort.Float64s(xs)
	sort.Float64s(ys)
	sort.Float64s(zs)
	mx, my, mz := Median(xs), Median(ys), Median(zs)
	best := int(float64(len(node.Shapes)) * 0.85)
	bestAxis := AxisNone
	bestPoint := 0.0
	sx := node.PartitionScore(AxisX, mx)
	if sx < best {
		best = sx
		bestAxis = AxisX
		bestPoint = mx
	}
	sy := node.PartitionScore(AxisY, my)
	if sy < best {
		best = sy
		bestAxis = AxisY
		bestPoint = my
	}
	sz := node.PartitionScore(AxisZ, mz)
	if sz < best {
		best = sz
		bestAxis = AxisZ
		bestPoint = mz
	}
	if bestAxis == AxisNone {
		return
	}
	l, r := node.Partition(best, bestAxis, bestPoint)
	node.Axis = bestAxis
	node.Point = bestPoint
	node.Left = NewNode(l)
	node.Right = NewNode(r)
	node.Left.Split(depth + 1)
	node.Right.Split(depth + 1)
	node.Shapes = nil // only needed at leaf nodes
}